Control of a phytoplankton bloom by wind-driven vertical mixing and light availability

The balance of physical and biological processes governing phytoplankton growth rates and the accumulation of biomass is widely debated in the literature, notably during the winter – spring transition. Here we show, in a temperate shelf sea that variability in the depth of the actively mixing surface layer is the leading order control. During a 2-week period preceding the peak of the spring bloom we observe two distinct regimes; ﬁ rst, growth within the euphotic zone during the day and re-distribution of new biomass to the seasonal pycnocline at night by convective mixing; then, more rapid biomass accumulation trapped within a shallower, wind-driven actively mixing layer that was decoupled from the pycnocline below. Our observations of the bloom in the Celtic Sea, Northwest European Shelf, were made using ocean gliders and include measurements of the dissipation of turbulent kinetic energy. A 1-D phytoplankton growth model driven by our measurements of dissipation and incident irradiance replicates the observed bloom and reinforces the conclusion that physical processes that mediate light availability were key. Day-to-day variability in cloud cover and the ability of phytoplankton to acclimate to their light environment were also important factors in determining growth rates, and the timing of the biomass peak. Our results emphasize the need for accurate turbulent mixing parameterizations in coupled hydrodynamic-ecosystem models. Our ﬁ ndings are

Marine phytoplankton are responsible for $ 50% of primary production on Earth and form the base of the global ocean food web (Field et al. 1998). Seasonal phytoplankton blooms make a significant contribution to oceanic primary production, air-sea CO 2 fluxes and carbon sequestration, as well as being key to the life cycles and trophic interactions of marine organisms (Lutz et al. 2007;Koeller et al. 2009;Signorini et al. 2012). Shelf seas contribute between 10% and 30% of global marine primary production, a disproportionately large contribution relative to their size (Mackenzie et al. 2005), they are a net sink for atmospheric carbon dioxide (Laruelle et al. 2018) and are responsible for up to 50% of the organic carbon supplied to the deep ocean (Jahnke 2010).
A bloom occurs when phytoplankton growth rates exceed losses such that a sustained period of growth leads to a (net) accumulation of biomass. Controls on growth rates include light availability, nutrient supply, and temperature, while losses may occur through mortality (e.g., grazing, viral lysis), advection, sinking or detrainment out of the euphotic zone. Hence, a range of biological and physical controls govern the balance between phytoplankton production and loss (see Lindemann and John 2014 for a review). Several hypotheses concerning the mechanisms for bloom initiation have emerged, with three main themes: critical depth, critical turbulence, and dilution re-coupling. While we do not implicitly test these hypotheses, we briefly review them here in order to provide the reader with an overview of the relevant physical and biological drivers behind increases in phytoplankton growth rates and the accumulation of biomass. Many of these drivers are relevant throughout the year, not just for spring bloom initiation. Sverdrup (1953) introduced the concept of a critical depth, based upon the assumption that phytoplankton growth was proportional to light, such that it decayed exponentially with depth, and that losses were constant throughout the water column. He defined a critical depth that needed to be deeper than the surface mixed layer for vertically integrated phytoplankton growth to outweigh losses. This evolved into the widely accepted view that the spring bloom is triggered once a seasonal surface mixed layer is established above the critical depth (e.g., Bishop et al. 1986;Siegel et al.. 2002). Conventional temperature or density thresholds for the depth of the surface mixed layer (de Boyer Montegut et al. 2004) however, do not always accurately represent the depth to which cells are actively mixed (Brainerd and Gregg 1995), which violates Sverdrup's assumption of a thoroughly well mixed layer.
The critical turbulence hypothesis proposes that a surface bloom can start in an arbitrarily deep surface mixed layer if the turbulence is weak enough in the well-lit surface waters for phytoplankton to receive sufficient light before being mixed down beneath the critical depth, or when an intensely mixed surface layer shoals, exposing phytoplankton to favorable light conditions (Huisman et al. 1999;Waniek 2003). A bloom can therefore potentially start in the winter following the shutdown of deep convective mixing and before onset of the mixed layer in spring (Townsend et al. 1994;Taylor and Ferrari 2011;Ferrari et al. 2015). Decreases in turbulence are typically attributed to net surface warming and reduced wind stress (Chiswell 2011;Taylor and Ferrari 2011;Chiswell et al. 2013). Later, Brody and Lozier (2014) proposed that a decrease in the dominant mixing length scale from deep winter convection to a shallower wind mixing regime above the compensation depth (where net phytoplankton growth is zero) is a reliable condition for the initiation of the spring bloom. Such conditions can occur when surface heat fluxes are still negative (i.e., convective), the seasonal thermocline is deep or when stratification is still weak. In subtropical regions, where nutrient availability rather than light is the limiting factor on phytoplankton growth, mixed layer deepening during the winter provides an essential injection of nutrients that may stimulate high enough growth rates to support an increase in biomass throughout the entire, deepening mixed layer (Zarubin et al. 2017).
In addition to a thoroughly well-mixed surface layer, Sverdrup also assumed a constant loss rate with depth (encompassing respiration, grazing, viral lysis, and sinking). The dilution re-coupling hypothesis (and later the disturbance-recovery hypothesis) argues that grazing pressure is the leading order control on bloom timing (Behrenfeld 2010;Behrenfeld et al. 2013), and allows a bloom to start before stratification is established. Behrenfeld (2010) argues that within a deep winter mixed layer, replete with nutrients, phytoplankton cell density and predator-prey encounter rates are reduced, stimulating an increase in phytoplankton biomass prior to the onset of stratification. As the surface mixed layer develops in early spring, phytoplankton and the grazer populations become concentrated into a thin surface layer, which increases encounter rates and eventually limits bloom biomass.
Existing studies have therefore identified key mechanisms controlling bloom initiation, and these processes are also important for understanding the subsequent development and dynamics of a bloom. However, many studies are based on open ocean data and/or coarse spatial ($ 10s km) and temporal ($ week) scale climatologies and satellite imagery (Henson et al. 2009;Behrenfeld 2010;Brody and Lozier 2014). These coarse time-and basin-scale approaches are generally unable to resolve short-lived (sub-daily) changes in wind stress, wave climate, incidental irradiance or diurnal cycles of heating and cooling that all impose a significant adjustment to stratification, turbulent mixing and light availability during a bloom. Furthermore, they do not consider the role of tidal mixing, an important forcing in shelf sea environments.
Like the open ocean, in shelf seas sub-daily changes in meteorological forcing impose significant adjustments to the length of the growing season and have implications that cascade to other functions of the ecosystem (Powley et al. 2020). However, the relatively shallow bathymetry of shelf seas and their proximity to land means that they are biologically and dynamically different to the open ocean. The penetration of light can be strongly attenuated by particulate and dissolved matter that is re-suspended from the seabed, or originates from rivers (Babin et al. 2003). Tidal currents drive strong vertical mixing, such that the seasonal pycnocline is established once the rate of stabilizing heat input is able to overcome the de-stabilizing effect of wind and tidally-generated turbulent mixing (Simpson and Bowers 1981). The seasonal pycnocline is therefore always established from the sea surface downward, above the critical depth. Further, the maximum winter mixing length scale is constrained by the bathymetry, which limits the likelihood of phytoplankton and their grazers becoming sufficiently decoupled to allow seasonal cycles of dilutionrecoupling as may occur in the open ocean.
In this paper we present high-resolution observations from moorings and gliders collected in the central Celtic Sea, Northwest European Shelf, combined with analysis of a 1-D physics and phytoplankton growth model. Seasonal stratification at our study site started to build on the 26th March 2015 coincident with the beginning of sustained positive heat flux (Ruiz-Castillo et al. 2019a). Surface chlorophyll a (Chl a) concentrations also started to increase above background winter levels on this day. The observations presented here start 10 d later, providing detailed information on the evolution of the bloom over the 2-weeks leading up to its peak. These 2-weeks accounted for 90% of the observed increase in surface Chl a concentration between bloom initiation and the biomass maximum. We demonstrate the critical role that the strength and structure of turbulent mixing from winds, nighttime convection and tidally generated shear played in governing the development of the bloom, via controlling exposure of phytoplankton cells to light. Uniquely, we have coincident measurements of all the key physical variables that underpin phytoplankton growth rates: temperature, irradiance, stratification, mixed layer depth, vertical turbulent structure, and active mixing layer depths. Unlike previous studies that have relied on parameterizations and length scales to estimate the depth of active mixing, we employ sub-hourly profiles of the observed turbulent kinetic energy dissipation rate. The conclusions we draw are widely applicable to any shelf or open ocean region where wind-driven mixing has the ability to modify nutrient and light availability, especially across subpolar shelves in the northern hemisphere where light rather than nutrients is typically the limiting factor on phytoplankton growth.

Study site
Our observations are from a central location in the Celtic Sea ( Fig. 1a) on the Northwest European Shelf. The Central Celtic Sea study site was 120 km from the shelf edge and 200 km from the nearest coast or tidal mixing front. The total water depth was 145 m. Dynamics in the Central Celtic Sea are typical of a seasonally stratifying temperate shelf with vertical ocean-atmosphere heat fluxes principally controlling the development and breakdown of seasonal stratification (Simpson and Bowers 1981;Wihsgott et al. 2019).

Gliders
Two gliders were deployed and maintained a position in close proximity to a mooring array in the Central Celtic Sea for 21 d between the 4th and 25th April 2015 (Fig. 1a). The gliders remained within 10 km of the moorings and of each other throughout the deployment, which is within one tidal excursion and so considered quasi-stationary for the purposes of this study (Fig. 1b). Each glider was equipped with the following sensors:

A Teledyne Slocum glider carrying a SeaBird CTD and a
Rockland Scientific Microrider ocean microstructure package, providing coincident measurements of the dissipation rate of turbulent kinetic energy (ϵ). A total of 1642 profiles were completed, averaging one profile every 16 min. Data were processed following Palmer et al. (2015) to obtain coincident profiles of temperature, salinity and ϵ from near-surface (typically with 5 m) to 110 m depth at approximately 1 m vertical resolution. 2. An iRobot Seaglider carrying a SeaBird CT sensor, Paine pressure sensor and Wetlabs ECO-puck for Chl a fluorescence (470/695 nm) completed 1547 profiles, providing information on vertical water column density structure and Chl a fluorescence on average every 19 min. Thermal lag corrections were applied following the methods of Garau et al. (2011). Glider temperature, salinity and fluorescence were calibrated against CTD casts from a coincident ship-based survey (see CTD profiles and nutrient sampling section) and each profile interpolated onto a 1 m resolution grid. Corrections were applied to minimize the impact of fluorescence quenching (see Data S1).

Moorings and surface forcing
Wind speed and direction, air temperature, relative humidity, air pressure, sea surface temperature, salinity, Chl a fluorescence and photosynthetically active radiation (PAR, 400-700 nm) were measured from a mooring at the Central Celtic Sea study site. Full water column currents were provided by an upward looking 150 kHz Flowquest acoustic current profiler mounted in a bed frame.
Wind stress, τ = C d ρ a W 2 , was calculated where W is the wind speed (m s −1 ) at 10 m above the sea surface, ρ a the air density (1.22 kg m −3 ), and C d the drag coefficient. The bottom tidal stress, τ b = k b ρ 0û j jû , was calculated whereû is the depth mean tidal current velocity (m s −1 ), k b a drag coefficient (0.0025) and ρ 0 (kg m −3 ) the near bed density.
Neglecting any changes in salinity, the net surface heat flux (Q net ; W m −2 ), was obtained by subtraction of the longwave back radiation, sensible and latent heat fluxes from the net shortwave radiation entering the ocean. The sensible and latent heat fluxes were calculated following Fairall et al. (1996). The net longwave and net shortwave radiation were taken four-times a day from the National Center for Environmental Prediction and the National Center for Atmospheric Research reanalysis product (Kalnay et al. 1996). Time series were extracted from the grid point covering the Central Celtic Sea mooring location. The heat flux is defined to be positive when the ocean is gaining heat from the atmosphere.
CTD profiles, Chl a, and nutrient sampling Between 3rd and 28th April, 42 full depth CTD casts were taken within 10 km of the Central Celtic Sea mooring array (Fig. 1b). A SeaBird 911plus CTD recorded temperature, conductivity, and pressure. Derived salinity was calibrated against in situ samples analyzed on a Guildline 8400B autosal. Chl a fluorescence was recorded by a Chelsea Technology Group Aquatracka MKIII and calibrated against > 200 samples of extracted Chl a, including 24 deep (> 60 m) samples from beneath the seasonal pycnocline on the shelf (see Fig. S1). A Biospherical QCP Cosine PAR sensor measured down-welling surface irradiance. Water samples for the determination of Chl a and nutrient concentrations (including nitrate + nitrite) were collected on each CTD cast. Chl a samples were collected from six depths and fluorescence was measured on a Turner Design Trilogy fluorometer using a nonacidification module and calibrated against a solid standard and a pure Chl a standard (Sigma-Aldrich, UK) (see details in Mayers et al. 2019). Nutrient samples were analyzed onboard using a Bran and Luebbe segmented flow colorimetric auto-analyzer using techniques described in Woodward and Rees (2001).

Seasonal pycnocline and mixing layer definitions
We define the base of the active mixing layer (z mix ), as the depth at which ϵ approached background values, chosen here as ϵ = 1 × 10 -8.5 W kg −1 . This fixed threshold is midway between the two values used by Sutherland et al. (2014) and indicative of a region that, from 9th April onwards, consistently separated upper and mid-water turbulent processes. Between the 5th and 9th April this fixed dissipation threshold was less reliable, a recognized problem when the water column is experiencing diurnal stratification, nighttime convection and wind-driven mixing (see e.g., Sutherland et al. 2014).
During this spring transition period we were unable to reliably identify the seasonal pycnocline by assessing changes in density from a near surface reference. The commonly used threshold of 0.125 kg m −3 (de Boyer Montegut et al. 2004;Henson et al. 2009;Chiswell 2011) was too high and only correctly identified the seasonal pycnocline on the 21st April, clearly after the observed change in vertical structure and after the main peak in phytoplankton biomass. We therefore define the base of the seasonal pycnocline (z pyc ) as the depth at which the density decreased by 0.015 kg m −3 below the bottom density. This is comparable to the value of 0.02 kg m −3 used by Hickman et al. (2012) in nearby waters and is aligned with the deepest local maxima in buoyancy frequency.
The buoyancy frequency, N 2 = − g ρ 0 dρ dz , was calculated as a measure of the strength of stratification and the depths of multiple mixed layers. Here, g = 9.81 m s −2 is acceleration due to gravity, dρ dz is the vertical potential density gradient and ρ 0 the average density. We also calculate the potential energy anomaly, Φ = 1 , between the sea surface and depth of the pycnocline (z pyc ), whereρ is the mean potential density of the water column. Φ is a measure of the work required to bring about complete vertical mixing. It is zero for a fully mixed water column and increases as stratification strengthens.

Light attenuation and euphotic zone depth
The vertical light attenuation coefficient of photosynthetically active radiation (PAR), K d , was estimated from 42 PAR profiles taken within 10 km of the Central Celtic Sea site (Fig. 1b). K d was estimated from the slope of a linear fit through log(PAR) vs. depth within the surface mixed layer. The depth of the surface mixed layer was defined according to a density threshold of Δσ θ = 0.025 kg m −3 . Where there was strong vertical variability in Chl a and PAR with depth and/or high near surface Chl a, the fit was performed on the nearsurface portion of the data only, avoiding the top 5-10 m where light data were unreliable. A significant linear relationship (r 2 = 0.78, p < 0.01, n = 35) was established between log (K d ) and log(chl sml ) where chl sml is the mean Chl a concentration within the surface mixed layer. Based on this relationship, for each glider profile we calculated the depth of the euphotic zone (z eu ), defined as the 1% irradiance level, from a measure of the mean Chl a concentration within the surface mixed layer.

The phytoplankton growth model
To help diagnose how and why phytoplankton biomass accumulation and carbon production rates were changing during the glider deployment we setup a coupled 1-D Lagrangian random walk and phytoplankton growth model. The model moves particles around a 1-D water column. The trajectories of the particles are determined by the vertical eddy diffusivity values that we derived from the observed dissipation profiles. The phytoplankton cells "contained" within each of the particles photosynthesize and respire according to the availability of light and nutrients. The model uses fixed grazing and cell mortality rates, a simplicity that allows the physical controls of turbulence and light to be isolated. In this section we describe the setup and forcing of the model and justify our parameter choices.

Lagrangian random walk
A Lagrangian random walk model, appropriate for a spatially nonuniform turbulent mixing environment with sharp vertical gradients in diffusivity (Visser 1997;Ross and Sharples 2004), was used to track 10,000 model particles through a vertical eddy diffusivity field estimated from the dissipation values measured by the glider. Each new particle position z n + 1 was calculated from its current position z n based on the following random walk: where K z (m 2 s −1 ) is the vertical turbulent eddy diffusivity and K 0 z = dKz = dz (m s −1 ) its derivative. Δt is the model time step (1 s) and R a random process of zero mean and variance v.
The vertical eddy diffusivities, K z z ð Þ = κ ε z ð Þ N 2 z ð Þ , were calculated assuming a mixing efficiency of κ = 0.2. Dissipation (ϵ) and buoyancy frequency (N 2 ) were taken from the glider profiles. The Slocum glider did not always collect reliable microstructure data within the top 5-10 m of the water column. In these instances, to achieve a realistic logarithmic decay of K z near the surface, near-surface dissipation was estimated using ε z ð Þ = w 3 Ã κz , where w * is the friction velocity derived from 10 m wind speeds (W 10 ) using the relation w 2 Ã = ρ a ρ 0 C d W 2 10 , and the von Karman constant, κ = 0.41. Wind stress was the dominant control on the mixing length scale, so this approach was considered acceptable.
Individual K z profiles were averaged using a 4-h running mean window, interpolated onto a 5-min grid and then filtered using a 4 m × 20 min Wiener noise removal filter. To avoid artificial particle accumulation, a piecewise cubic smoothing spine with 17 breakpoints was fitted to these averaged and filtered profiles of log 10 (K z ) to ensure that both the first and second derivatives of K z were smooth and continuous. An example is shown in Fig. 2a. To avoid particles accumulating at the boundaries the first derivative of K z was forced to zero at the surface and seabed by applying the spline to the original profile and its extensions, such that the profile was symmetric about the boundaries (Ross and Sharples 2004). Averaging the number of particles in each depth bin across the whole 24-d simulation, we find that the maximum deviation from the original uniform particle distribution is < 1.5% within 5 m of the surface and bottom, and < 5% within the euphotic zone (31 m).

Phytoplankton growth, light, and nutrient dynamics
The light exposure of phytoplankton cells is affected on hourly to daily time-scales by: (1) day-night, (2) cloud cover, and (3) vertical mixing. Each day a phytoplankton cell may experience both extreme light and fully dark conditions. It will continually try to acclimate (adapt) to its ambient light conditions to achieve optimal growth rates (Marra 1978;MacIntyre et al. 2000). A cells photosynthetic performance is therefore dependent upon its previous light history and its physiological state. We adapted the phytoplankton-growth model used by Ross and Sharples (2008) which was based on the production-acclimation model of Denman and Marra (1986) to simulate the instantaneous production of a cell being moved through a vertical light gradient.
Typically, photosynthesis-irradiance (P-I) curves are used to determine the average rate of photosynthetic carbon production at varying levels of light intensity (Jassby and Platt 1976;MacIntyre et al. 2002). The light intensity at which the photosynthesis of a cell starts approaching its maximum value is called the saturation onset. At higher light intensities, the rate of carbon fixation levels out, or even decreases, as the high light damages the phytoplankton photosystem. This is known as photo-inhibition and starts once light levels exceed the inhibition onset. Photo-acclimation of a cell to continually varying light levels throughout the day however means that the instantaneous rate of carbon production cannot be described by a single P-I curve (Denman and Marra 1986).
Within the model, we describe a cell that has adapted to achieve optimal growth rates under high light levels as fully "light acclimated." In contrast, a cell adapted to achieve maximum growth rates under extreme light limitation (very low light availability) is said to be "dark acclimated." Changes in the photosynthetic carbon production rate of a cell that is fully dark acclimated (P d ) or fully light acclimated (P l ) is described by two different curves: where P d=l m is the maximum dark or light-acclimated production, I d/l is the saturation onset for the dark or light-acclimated curve, and I is the light intensity at the particle depth. The instantaneous production, P, of a cell continually adapting to varying light levels, always falls between these two limiting curves, and depends on the acclimation status of the cell, i.e., its degree of adaption from the dark to the light acclimated curve. As a cell moves through light gradients its cellular acclimation status, Y, changes according to: where γ is the acclimation timescale (1 h) and Δt the model time step (1 s). Y reflects the light history of the cell that is always moving towards the acclimation status appropriate for the ambient irradiance level, described by the exponential function X(I) (see Denman and Marra (1986) for full details). The instantaneous light-dependent production (photosynthetic rate, P; d −1 ), is then calculated from: When the cell is fully adapted to a high light environment (Y = 1), the production rate is P l . If the cell has been in darkness for a sustained period the cellular acclimation status will be zero (Y = 0) and the production rate is P d . This photoacclimation model reproduces the observed hysteresis in photosynthetic rates over a day: high rates in the morning, a midday depression and values recovering in the afternoon but lower than at the same irradiances experienced in the morning (Marra 1978).
The incident irradiance (I 0 ) decays exponentially with depth and was taken from the PAR sensor on the mooring in the Central Celtic Sea. The model absorption coefficient is the sum of a background light absorption (0.09 m −1 ) and the contribution from self-shading, a depth dependent variable determined by the carbon biomass of cells in the water column above and their cell-specific absorption coefficient (see Ross 2004). Accurate replication of the light field is therefore dependent upon growth within the model, which is itself a function of light availability.
The model takes into account one key limiting nutrient, appropriate for the Celtic Sea (Holligan et al. 1984), dissolved inorganic nitrogen (DIN). The model cells take up DIN, at concentration N (mg m −3 ), store it up to a maximum nutrientcarbon ratio, and consume it during photosynthesis. The uptake of nutrients, U, follows: where Q is the cellular nutrient-to-carbon ratio which needs to remain above a subsidence quota (Q min ) and below a maximum storage quota (Q max ). U m is the maximum uptake rate and κ N the half-saturation constant. The cells take up as much nutrient as possible based on the nutrient concentration available in the water column and their cellular nutrient-to-carbon ratio. Water column nutrients are modeled on an Eulerian grid using a standard diffusion equation. A fixed N concentration of 92 mg N m −3 is maintained at the seabed to account for nutrients released from the sediments (see Ross (2004) and Ross and Sharples (2008) for full details).
The cellular increase in carbon is calculated according to where r is the cellular respiration rate (d −1 ), G is the grazing rate (d −1 ), C is the cellular carbon (mg C) and P is the instantaneous photosynthetic production rate (d −1 ). The model contains a lower carbon threshold (C starve ) at which cells begin to die and an upper level (C divide ) at which the cells divide (so the number of cells per Lagrangian particle doubles). Losses due to grazing and cell mortality act by reducing the number of cells within a particle. Lr is a carbon dependent mortality rate and is nonzero when C < C starve . All model parameters are given in Table 1. Grazing (G), respiration (r), and mortality rates (Lr) are all held constant. This is a defensible simplification for this 2-week period and allows us to isolate where there are shifts in the balance between physical and biological controls. It is not appropriate however, over longer time-scales where losses must vary in time in order to recreate annual cycles in phytoplankton biomass that have realistic minima and maxima (Behrenfeld and Boss 2018). To test this simplification, we ran sensitivity experiments with a range of fixed grazing rates. We also allowed for a linear six-fold increase as the bloom developed, equivalent to the increase in micro-zooplankton grazing pressure reported by Mayers et al. (2019) during the lead up to the biomass maximum for the same bloom. In all cases our conclusions remain the same. The importance of variable loss terms in the context of model-observation mismatch is discussed later.

Parameter choices and model initialization
The values assigned to each model parameter (see Table 1) were guided by similar model setups for temperate shelf sea environments including the Celtic Sea (Ross and Sharples 2008;Sharples 2008;Marsh et al. 2015). Maximum dark and light acclimated production (P d=l m ), the dark and light saturation onset (I d/l ) and the inhibition onset (I b ) were chosen based on the results of 14 C-uptake photosynthesis-irradiance experiments in the Celtic Sea (Hickman et al. 2012) and previous successful implementations of the model (Ross and Sharples 2008). The maximum production rate, P d m = 1.6 d −1 , lies comfortably within the range of carbon specific photosynthesis rates (0.27-3.83 d −1 ) reported in a review by MacIntyre et al. (2002).
Small flagellates (2-20 μm) as opposed to larger diatoms (> 20 μm) dominated the 2015 spring bloom in the Celtic Sea in terms of both Chl a and primary production (Poulton & Hickman, unpubl.). Although there was some succession within the flagellate groups present (cryptophytes and other nano-eukaryotic taxa; G. Tarran pers. comm.) as the bloom evolved, nano-eukaryotes remained the dominant group in terms of carbon biomass throughout April (Poulton & Tarran, unpubl.). Hence, we include only one type of cell in this basic model and assume that differences in the photo-physiology between the different nano-plankton taxa observed are considerably less than between nano-plankton and large diatoms.
The model was seeded with 10,000 model particles evenly distributed across a 110 m deep model water column. The initial DIN profile was based on fitting a hyperbolic tangent to discrete nitrate + nitrite samples taken at the Central Celtic Sea site on the 3rd and 4th April (Fig. 2b). Maximum (97 mg N m −3 ) and minimum (84 mg N m −3 ) concentrations in bottom and surface waters were separated by a 12 m thick nutricline centered at 60 m depth. Each model particle initially represented 800,000 "real" cells. Based on the observed biomass at the Central Celtic Sea during spring 2015 (Poulton et al. 2019), cells in the surface and bottom waters were given initial carbon contents of 1.5 × 10 −6 and 1.01 × 10 −6 mg C, respectively, and a vertical structure that mirrored the observed structure in Chl a fluorescence (Fig. 2c). These values are comparable to the mean phytoplankton carbon for autotrophic dinoflagellates (1.02 × 10 −6 mg C cell −1 ) reported in Poulton et al. (2019). The cellular nutrient-to-carbon ratio of each particle was initially set at the maximum (i.e., nutrient replete). The initial acclimation status (Y) was set to zero (fully dark acclimated) for cells deeper than 20 m and increased to 0.7 (approaching light acclimated) towards the surface.

Definition of a phytoplankton bloom
There is no consensus in the literature on the quantitative definition of a bloom, which leads to differing conclusions on when, where and why a bloom occurs (Brody et al. 2013). Here, following Behrenfeld and Boss (2018), we consider the specific rate of change of biomass and calculate it within the actively mixing surface layer. Increases and decreases in phytoplankton biomass result from the imbalance between carbon production (i.e., the growth rate, μ; d −1 ) and loss (l; d −1 ) rates. This imbalance can be expressed by the specific rate of change of biomass, or accumulation rate (a; d −1 ): where a is normalized to C phyto (mg C m −3 ), the depth average biomass between the sea surface and z mix . Carbon production (μ) is defined as the difference between a cells photosynthetic (P; d −1 ) and respiration rates (r; d −1 ). When light and nutrient availability support photosynthetic production that exceeds cell respiration (P > r) carbon is produced, and μ is positive. For phytoplankton biomass to increase (a > 0), and a phytoplankton population to be considered "blooming," μ must exceed l.
Correctly determining the balance between μ and l requires an appreciation of how the accumulation rate in a depth varying actively mixing layer, within which phytoplankton are uniformly distributed, is not always accurately represented by the rate of change of C phyto . In a deepening mixing layer for example, dilution of a phytoplankton population may be equal to or greater than the excess in μ over l. In this situation, the carbon biomass (in mg C m −3 ), may remain constant or even decrease, despite production rates exceeding losses. However, the depth-integrated carbon (ΣC phyto ) increases during the deepening and it is the specific rate of change of ΣC phyto that should be used to identify the balance between μ and l.
Mortality rate of cells if C < C starve 0. Quantified as the difference between photosynthesis (P) and respiration (r) and normalized by the total carbon biomass. † Quantifiable from changes in C phyto and ΣC phyto normalized by C phyto .
The inverse of this situation can also occur: as the actively mixing layer shoals, an increase in light availability may stimulate an increase in μ over l and an increase in biomass (a > 0). However, the depth-integrated biomass (ΣC phyto ) of the population blooming in the much shallower mixing layer is not necessarily greater than ΣC phyto in the previously deeper layer.
We therefore quantify the specific accumulation rate within the active mixing layer (a zmix ; d −1 ) from the rate of change of depth averaged C phyto (normalized by C phyto ) when the active mixing layer is constant or shoaling, and from normalized changes in ΣC phyto when z mix is deepening. Specifically, These definitions help us to understand the conditions under which phytoplankton biomass above z mix changes. We have calculated C phyto and ΣC phyto for the model (in mg C m −3 and mg C m −2 ) and Chl phyto and ΣChl phyto (in mg Chl m −3 and mg Chl m −2 ) for the observations.
From the model we also calculate the specific carbon production rate within the active mixing layer (μ zmix ; d −1 ). This is defined as the difference between photosynthetic carbon production and respiration rates (P − r), normalized by the total carbon biomass. The average between the sea surface and z mix is then taken. Table 1 summarizes our definition of key variables.

Physical observations Wind stress and tidal forcing
Over the deployment period, wind speeds varied between 2 and 13 m s −1 equivalent to a 0.4 N m −2 range in surface wind stress (Fig. 3a). Between the 4th and 9th April the wind stress remained relatively stable, fluctuating around 0.1 N m −2 , with subsequent peaks occurring on the 11th, 12th, 16th, 18th, 19th, 21st, and 22nd April separating otherwise low stress levels. Episodes of weak wind stress (< 0.05 N m −2 ), equivalent to wind speeds < 5 m s −1 , lasting more than 24 h, occurred between 9th and 10th April and between 13 and 15 th April. There were also shorter episodes of low wind stress on the 17th, 20th, and 23rd April.
Spring tides occurred on the 6th and 20th April (Fig. 3a) when tidal bed stress peaked at 0.75 and 1 N m −2 , respectively. There was a marked asymmetry between the two spring periods with the second spring tide being the stronger of the two. Only one neap tide was observed, occurring on the 13th April however, the final stages of the study observed relatively low tidal flows, close to a second neap period. On sub-daily time scales, the bottom tidal stress was dominated by quarterdiurnal (flood-ebb) variability.

Incident radiation
The maximum irradiance each day ranged between 430 and 1315 μE m −2 s −1 , reflecting the highly variable nature of cloud cover and fog across the Celtic Sea (Fig. 3b). During the afternoon of the 18th April, the incoming light dropped to 200 μE m −2 s −1 . Maximum midday irradiance was close to or below 500 μE m −2 s −1 on eight individual days, including a period of four consecutive days between the 14th and 18th April. High levels of radiation (> 1000 μE m −2 s −1 ) between the 5th and 13th April were interrupted by a cloudy day (< 500 μE m −2 s −1 ) on the 10th April. Similarly, a cloudy day on the 23rd April interrupted four consecutive clear sky days and was followed by strong incident radiation on the 24th and 25th April. The average depth of the euphotic zone over the glider deployment was 31 m (range of 21-46 m).

Heat flux and stratification
There was a diurnal cycle of the ocean gaining heat during the day and then loosing heat to the atmosphere during the night, indicated by alternating positive and negative surface heat flux respectively (Fig. 3c). The average daily surface heat flux between the 4th and 25th April was 89 W m −2 , representing a net heat flux into the ocean throughout the deployment period. This resulted in a gradual increase in near surface temperature from 10 C on 5th April to over 11 C on the 25th (Fig. 4a). The surface (1-5 m average) to bottom (105-110 m average) difference in potential density subsequently increased from 0.05 to 0.20 kg m −3 during this time. For wider seasonal context, at the peak of stratification in August 2014 the top-bottom density difference at the Central Celtic Sea was > 2.00 kg m −3 (Wihsgott et al. 2019).
Throughout the deployment period, increases in nearsurface temperature were observed to coincide with short (hourly to 1-2 d) episodes of reduced wind stress (notably on the 9th-10th and 13th-14th April) that allowed shallow, warm surface layers to form over the upper 10-20 m (Fig. 4a). These shallow, thin layers were subsequently eroded by the following input of increased wind stress (or by nighttime convection), and the heat retained within the near surface then re-distributed down to the base of the seasonal pycnocline. This cycle of surface warming followed by mixing of the upper ocean, increasing and then lowering the potential energy anomaly between the sea surface and z pyc , respectively (Fig. 4c), resulted in an increasing strength of stratification over the whole water column and at the main pycnocline, where the buoyancy frequency (N 2 ) gradually increased throughout April, from < 1 × 10 −5 to > 2 × 10 −4 s −2 (Fig. 4b).
The periodic establishment and breakdown of warm layers created a hydrographically multi-layered water column that is reflected in the buoyancy frequency (Fig. 4b). For example, between the 13th and 15th April there were 4 hydrographically well mixed layers; a 5-m thick warm (> 10.8 C) near surface layer, cold (< 10 C) bottom water beneath the seasonal pycnocline and two more weakly defined 10-20 m thick layers in-between.
Between the 4th and 8th April, when the wind stress was relatively low and stable, the potential energy anomaly (Φ) between the sea surface and z pyc became weakly positive during the day, reflecting diurnal warming and the establishment of weak near surface stratification (Fig. 4c). Φ then dropped to zero in the early hours of the following morning, consistent with nighttime convection fully mixing the upper layer of the water column.

Dissipation of turbulent kinetic energy
Throughout the deployment, the depth of the active surface mixing layer (z mix ) varied between 5 and 56 m (Fig. 4e) and was observed to deepen every night relative to the preceding day. For example, in the early hours of the morning on 6th April z mix reached a depth of 56 m, then shoaled during the day to 18 m depth as a warm surface layer developed. Between the 5th and 9th of April, the potential energy anomaly (Fig. 4c) conclusively shows that the water column between the surface and z pyc became fully mixed each night. However this is not always mirrored by z mix reaching the pycnocline (Fig. 4e), a consequence of the fixed dissipation threshold not always reliably identifying the full depth of active mixing taking place within the surface waters during this period.
From the 9th April onwards the most pronounced deepening (shoaling) events were associated with increases (decreases) in wind stress. The linear correlation coefficient between wind stress and z mix from the 9th April onwards was 0.60 (p < 0.01, n = 1188). Prior to the 9th April, the correlation coefficient was 0.2 (p < 0.01, n = 359). Furthermore, the turbulent Langmuir number, La = (v * /U s ) 1/2 (Belcher et al. 2012), a measure of the relative influences of directly wind-driven shear (with friction velocity v * ) and wave-driven stokes drift (U s ) was > 0.3 for 55% of the deployment (> 0.2 for 98% of the time). This further suggests that, as a first approximation, wind stress was the leading order control on the depth and structure of near surface dissipation.
The measurements of ε reveal that the mixed layers, identifiable in Fig. 4b by regions of N 2 = 0, were not all regions of active mixing. For example, the multiple mixed layers between the 13th and 15th April had coincident levels of turbulence close to or at background levels. During this period, ε < 1 × 10 −9 W kg −1 was observed throughout the region immediately beneath the surface constrained active mixing layer and was only observed to increase towards the base of the seasonal pycnocline, z pyc . The previously identified period of reduced wind stress, starting on the 13th April, left behind remnant mixed layers, which were characterized by only weak levels of vertical mixing.

Hopkins et al. Turbulent control of phytoplankton bloom
Dissipation beneath z pyc was dominated by tidally driven quarter-to semi-diurnal and fortnightly spring-neap variability (Fig. 4e) attributable to tidally induced bed shear stress. During spring tides, the turbulent bottom boundary layer extended to the base of the pycnocline. In contrast, during less active tidal mixing an area of low-level turbulence was maintained between the extent of tidal mixing and z pyc indicating that the pycnocline depth takes some considerable time to increase relative to the spring-neap transition.

Observed phytoplankton bloom
The glider deployment between 4th and 25th April coincided with (1) maximum near surface (1 m) concentrations of Chl a recorded at the Central Celtic Sea site in 2015 (Fig. 1c), (2) the 2015 peak in Chl a concentration derived from satellite remote sensing (Fig. 1d) and (3) the annual maximum in depth integrated Chl a biomass calculated in 2015 from CTD casts (Wihsgott et al. 2019). This independently confirms that the gliders were in the water for the annual biomass maximum in the Central Celtic Sea.
The vertical structure of Chl a fluorescence was tightly coupled to that of the hydrographic structure with the depth of the maximum vertical gradient in Chl a closely correlated with the depth of z mix , especially from 9th April onwards (Fig. 5a,b). Phytoplankton were therefore well mixed throughout the active mixing layer.
Between the 5th and 8th April, the Chl a concentration above the pycnocline remained low (Fig. 5a). Depth averaged Chl a values increased gradually from 1.5 to 2.5 mg m −3 (Fig. 5c). Over these 4 d, Chl a concentration increased during the day in the warming, sunlit near surface waters (within the euphotic zone). During the early hours of the morning (between 04:00 and 06:00 on 5th-8th April) the depth average Chl a concentration within the active mixing layer (above z mix ) became the same as between the sea surface and the pycnocline (Fig. 5c inset), coincident with the potential energy anomaly dropping to zero (Fig. 4c). Therefore, during the night new biomass was being redistributed down to the seasonal pycnocline by convective mixing, homogenizing the Chl a concentration. From the 9th April onwards the average Chl a biomass within the active mixing layer and above the seasonal pycnocline diverged, reflecting the decoupling between the active mixing layer and seasonal pycnocline, and a bias towards a higher relative biomass accumulation above z mix (Fig. 5a,c). Figure 6a shows the depth-integrated Chl a between the surface and seabed (145 m). The glider did not observe the bottom 20 m but, due to its well-mixed nature, we assume that the Chl a concentration in this near-bed layer was homogenized and so the same as was observed at 125 m. There was a maximum in the full depth integrated biomass of 154 mg Chl m −2 on 17 th April and a smaller, second peak (136 mg Chl m −2 ) on the 21st. On these two dates, 82% and 78% of the total water column biomass was contained within the active mixing layer (Fig. 6b). Following this period, there was a steady decline in the total water column biomass until the end of the observational period on the 25th April.
There were four periods, each coinciding with or immediately following a reduction in wind stress and shoaling of z mix , where the depth average Chl a within the active mixing layer increased (dashed black boxes in Fig. 5c). Maxima in Chl a during these times occurred on the 10th, 14th, 17th, and 20th April. The highest observed depth average (11 mg Chl m −3 ) and depth integrated (160 mg Chl m −2 ) biomass peaks were on the 17th April (Fig. 5c,d). Figure 5e shows the specific biomass accumulation within the active mixing layer, calculated from the observed Chl (c) Depth averaged Chl a concentration (C phyto ) from the glider within the active mixing layer (red) and above the seasonal pycnocline (blue). Top left inset is a zoom between 4th and midnight on 9th April. Dashed black boxes mark periods coinciding with or immediately following a reduction in wind stress and shoaling of z mix . (d) Depth integrated (ΣC phyto ) Chl a concentration from the glider (left) and carbon from the model (right) between the surface and z mix . Dashed back line is ΣC phyto smoothed with a running 14 h averaging window to remove variability associated with semidiurnal tidal advection. (e) Observed specific biomass accumulation rate (a zmix ; day −1 ) within the active mixing layer. a concentration. Until the 18th April it is almost exclusively positive, indicating that phytoplankton growth rates were always greater than losses during this time (μ > l). Local maxima in biomass accumulation rates occurred on the 7th and 10th April. An extended period of positive accumulation took place between the 13th and 18th April, with a local maximum on the 15th. On the 18th April, the specific accumulation was negative suggesting that biomass losses within the active mixing layer were greater than phytoplankton growth rates. This was followed by a final episode of positive accumulation between the 19th and 21st, and then a switch to primarily negative rates from the 21st April onwards. Between 5th and 24th April phytoplankton were blooming (accumulating in biomass, a zmix > 0) for 75% of the time.

Modeled phytoplankton bloom
Before presenting the modeled specific accumulation and carbon production rates we first assess the model performance against independent measurements of DIN concentration, PAR and Chl a profiles, and the variability of observed biomass within the active mixing layer.

Model validation
Between the 4th and 15th April nutrient samples taken from CTD casts at the Central Celtic Sea study site show that the surface water nitrate + nitrite (DIN) concentrations decreased from 84 mg N m −3 (6 μM) to 17 mg N m −3 (1.2 μM). This is well reproduced by the model which predicts that DIN within the active mixing layer was fully depleted by the end of the 17 th April (Fig. 7; red squares). In situ samples from the 20th and 21st show that the DIN concentration above z mix increased again, reaching 30 mg N m −3 (2.1 μM). This pattern is also reproduced by the model which shows a sudden increase in DIN on the 18th April (Fig. 7). There was a gradual decline in bottom water DIN over the deployment, which is well simulated by the model.
The top panel in Fig. 8 compares individual profiles of modeled carbon biomass against vertical profiles of Chl a from the glider. Each of these times coincides with a deployment of the CTD (times indicated in Fig. 5a by pink squares). The depth and gradient of the observed chlorocline matches the vertical structure of biomass in the model well. Given the tight coupling between turbulent structure and biomass ( Fig. 5a) we compare the modeled depth integrated carbon between the surface and z mix with the equivalent depth integrated Chl a from the glider (Fig. 5d). These time series compare remarkably well, particularly up until the final biomass peak on the 21st April. From the 21st onwards the depthintegrated Chl a started to reduce, whereas carbon biomass in the model did not. From the 18th April onwards there is semidiurnal variability in the observations that is not captured by the model. This is probably the result of tidal advection moving patchy growth back-and-forth. Since the 1D model does not include an advective term, it was unable to simulate this signal.
The carbon-chlorophyll ratio within the euphotic zone during the deployment, as determined from discrete fluorometric Chl a measurements and phytoplankton cell counts (Poulton et al. (2019), flow cytometry and microscope counts converted to carbon using literature cell carbon conversions), ranged from 14 to 37 (g C : g Chl a), with a cruise average of 27 ± 8 (mean ± SD) (Poulton et al. 2019). The equivalent ratio within the euphotic zone implied by the model-glider comparison is 39 ± 11 (mean ± SD), which is in general agreement with the cruise estimates and well within the range of values reported in the literature (MacIntyre et al. 2002).
The bottom panels in Fig. 8 compare modeled profiles of PAR (based on the models own self-shading algorithm) against the nearest (in time) profile of PAR directly measured during a CTD cast (within 10 km of the Central Celtic Sea site). This confirms that the model attenuated light realistically throughout the whole timeseries.
Light availability, biomass accumulation and production Figure 9a shows the incidental PAR at the sea surface (I 0 , as in Fig. 3b) and the average light dose experienced by cells within the active mixing layer (I zmix ; μE m −2 s −1 ). The maximum depth average light available to cells above z mix at midday varied by a factor of 13, from 50 μE m −2 s −1 on the 18th to 650 μE m −2 s −1 on the 13th April. In contrast, the maximum midday irradiance at the sea surface varied by a factor of three. Between the 5th and 23rd April, light levels exceeded the inhibition onset on 10 individual days for between 2 and 9 h. On 9 d, the maximum midday light availability ranged between the saturation and inhibition onset (90-200 μE m −2 s −1 ). The 18th April was the only day where the average light within the active mixing layer did not reach saturation.
Day-to-day changes in the average light levels within the active mixing layer did not always correlate with increases and decreases in the surface irradiance. This was driven by changes   Fig. 6a). (Bottom row) Model PAR (black dashed) and PAR profile from the nearest CTD cast (black solid). Black dot is PAR at the sea surface measured at the mooring array.
in the depth of the active mixing layer. For example, on the 9th April, despite the incoming radiation at the surface being the same as on the 8th April, phytoplankton received twice the amount of light. Similarly, light levels at the surface on the 21st and 22nd April were comparable but deepening of z mix on the 22nd reduced the light availability.
In Fig. 9b the modeled and observed specific biomass accumulation rates within the active mixing layer are compared. There are five main modeled peaks on the 7th, 10th, 15th, 18th, and 21st April. These compare well to the observations, particularly on the 7th and 10th April. The daily average (midnight to midnight) specific carbon production rate within the active mixing layer is plotted in Fig. 9c. From this depth average perspective, maximum carbon production rates (> 0.25 d −1 ) occurred on the 14th April. On the 18th April, there was no net carbon production. Other local maxima in production occurred on the 6th, 10th, 17th, and 20th April. Figure 9d shows the depth of the active mixing layer shaded by the (model) specific biomass accumulation rate. The peaks in biomass accumulation on the 10th and 15th April occurred after the active mixing layer had been stable at around 10 m depth for the previous 24 h. The peaks on the 7th, 18th, and 21st April coincided with z mix deepening. In Fig. 9e the depth of the active mixing layer is shaded according to the instantaneous depth averaged specific carbon production rate. As z mix shoaled, typically during the day, carbon production increased. During the night the carbon production is negative (respiration > photosynthesis). Note that  (black) and observed (dashed gray) specific biomass accumulation rate within the active mixing layer (a zmix ; d −1 ). (c) Daily average carbon production rates (μ; d −1 ) within the active mixing layer (z mix ) and between z mix and the seasonal pycnocline layer (z mix to z pyc ). (d) Depth of the active mixing layer shaded according to the modeled specific biomass accumulation rate within it. (e) Depth of the active mixing layer shaded according to the instantaneous carbon production rate within it. Note that the color scale saturates at zero, so periods of negative production, typically during the night, are all dark blue. the color scale saturates at zero in order to highlight the periods of growth (positive production rates). The maximum carbon production rates were on the 10th and 14th April (Fig. 9c) when the active mixing layer was at its shallowest.
There was a maximum in phytoplankton biomass on the afternoon of 17th April (Fig. 5). This coincided with the surface water becoming fully depleted of nutrients (Fig. 7). During the 18th, one of the cloudiest days on record (Fig. 9a), the wind stress increased (Fig. 3a) to 0.4 N m −2 (> 13 m s −1 wind speed) and the active mixing layer deepened from 10 to 44 m. This resulted in the average light levels above z mix being lower than the saturation onset (Figs. 9a, 10a). Nutrient depletion and reduced light availability resulted in the lowest carbon production rates of the month (Fig. 9c) and signified the end of the spring bloom. The wind-driven deepening of the mixing layer on the 18th injected new nutrients into the surface waters from depth (Fig. 7), which triggered another period of positive biomass accumulation, and resulted in a smaller, secondary biomass maximum between the 20th and 21st April.

Nighttime convection and wind-mixing regimes
Our glider observations and model results have revealed important aspects of the drivers, timing, magnitude and rates of phytoplankton growth on a temperate continental shelf during the spring. We identify two distinct regimes that followed the initial onset of stratification and initiation of the bloom (Fig. 10). During the first, a diurnal cycle of near surface growth followed by nighttime convective mixing dominated. During the second, the largest and most rapid changes in phytoplankton biomass took place within the wind-driven active mixing layer.
From previous work, we know that positive heat input initiated seasonal stratification on 26 th March but, vertical differences in salinity (not temperature) controlled 50%-100% of water column stability between the 28th March and 7th April (Ruiz-Castillo et al. 2019a). Our gliders were deployed during this time (on 4th April) and we initially observe (and model) modest increases in phytoplankton biomass within the nutrient replete euphotic zone during the day, followed by redistribution of this biomass to the depth of the pycnocline by convective overturning at night (Fig. 10a).
The regime changed on the 9th April, coincident with temperature regaining full control of water column stability (Ruiz-Castillo et al. 2019a) and wind stress as opposed to convection becoming the dominant driver of the mixing length scale (Wihsgott et al. 2019). The increased importance of wind-driven mixing is consistent with the increase in correlation between z mix and wind stress that we find here, from 0.2 during the 5 d before the 9th April to 0.6 for the remainder of the deployment. On the 9th April, reduced wind stress allowed the warm near surface cap established during the day to be maintained throughout the night, decoupling the active mixing layer from the seasonal pycnocline by more than 30 m. Henceforth periods of low wind stress resulted in shoaling of the active mixing layer and stimulated rapid increases in biomass that were retained within it (Fig. 10b). The model confirms that maximum carbon production and biomass accumulation rates occurred after periods of reduced wind stress and shoaling of z mix , which drove changes in light availability to the phytoplankton cells. plankton growth within the sunlit surface waters during the day. At night convective mixing evenly redistributes new biomass to the seasonal pycnocline (z pyc ). Nutrient availability does not limit phytoplankton growth. (b) Wind-driven mixing. The active mixing layer (z mix ) and the seasonal pycnocline (z pyc ) are decoupled. Wind stress is the dominant control on the depth of active mixing. Phytoplankton growth rates and biomass accumulation above z mix increase when the wind stress is weak and the active mixing layer shoals. Initially, when nutrients are replete, light availability mediated by changes in the depth of active mixing dominates phytoplankton dynamics (10th-17th April, 2a). Once nutrients have been depleted within the active mixing layer, both light and nutrient availability are important controls on growth (18th-25th April, 2b).
The dynamics that we observe during the wind-dominated regime are well aligned with Brody and Lozier (2014) who predict phytoplankton growth during the winter-spring transition once a shallow wind-mixing regime can be established, and also with Chiswell (2011) who shows the spring bloom forming in shallow, weakly stratified near-surface layers, above the seasonal pycnocline.
In summary, we observed two regimes, where (1) convective overturning initially dominated changes in the depth of active mixing and re-distributed new growth within the euphotic zone towards the base of the seasonal pycnocline each night (Fig. 10a), to (2) a wind-mixing regime where phytoplankton and new biomass were trapped within the wind-driven layer, only being mixed to the pycnocline during periods of high wind stress (Fig. 10b). The wind-mixing regime could further be separated into (a) a nutrient replete condition where light availability and photoacclimation dominated the phytoplankton dynamics, and (b) when nutrients began to deplete and nutrient resupply from beneath the pycnocline also became important.
Although we identify two distinct periods of time in our data set where the dominant drivers of vertical mixing were different, the evolution of the spring bloom in temperate seas progresses differently each year owing to interannual and spatial variability in, for example, wind and tidal forcing, insolation and water column stability. Not every year will exactly mirror the progression between the regimes described here in 2015, but in situations where light rather than nutrients is the limiting factor for phytoplankton growth, convective and/or wind-driven mixing are likely to be of leading order importance during the spring bloom period.
We now consider further details behind these key findings, specifically; the importance of distinguishing between the depth of active mixing and the pycnocline, the role of variable cloud clover and photo acclimation, the balance between physical and biological controls on phytoplankton growth, the impact of wind and tidally driven advection, and the role of the barotropic tide in entraining biomass into bottom water, a process unique to shelf sea environments.

Active mixing and pycnocline depth
The base of the active mixing layer and the seasonal pycnocline was decoupled during most of our study period, separated by a region of low turbulent kinetic energy dissipation and vertical eddy diffusivities (Fig. 4e). We did not see any significant biomass increases within this weak mixing layer (Fig. 5a), even though z pyc was close to or shallower than estimates of Sverdrup's critical depth ($ 50 m; Wihsgott et al. 2019) and surface waters were not depleted in nutrients (Fig. 7). Although daily depth average light between the sea surface and z pyc varied from 12 to 60 μE m −2 s −1 , phytoplankton cells were not being thoroughly mixed over the whole light gradient. Those trapped between z mix and z pyc received an average daily light dose of only 0.15-17 μE m −2 s −1 , much less than those retained within the active mixing layer (Fig. 11a). The daily carbon production rate of cells between z mix and z pyc was less than zero (except on the 9th and 13th) (Fig. 9c). Therefore, even though there were sufficient nutrients to support growth, cells remained light-limited beneath z mix . Light-limitation was only fully relieved within the nearsurface mixing layers that kept cells shallower in the euphotic zone and above z eu (Fig. 5a). On the 9th and 13-14th April the average biomass between z mix and z pyc was at its highest (2.5 mg Chl m −3 ). At these times of minimum wind stress, z mix was very shallow and sat 20 m above z eu (Fig. 5a). Growth was therefore possible in a 10-20 m layer below z mix , which is reflected in the positive daily production rates between z mix and z pyc on these dates (Fig. 9c).
In contrast to the canonical two-layer view of a stratifying shelf sea, we propose that during the winter-spring transition and build-up to the spring biomass peak it is more appropriate to consider a three-layer system comprising: a wind-driven actively mixing surface layer where there are sufficient nutrients and optimal light conditions to support rapid phytoplankton growth and biomass accumulation; bottom water, beneath the seasonal pycnocline, where there is no light to support net growth; and, in between, a layer where nutrients are available but where low vertical eddy diffusivities hold phytoplankton cells at low light levels that are typically insufficient to allow for any sustained net growth and biomass accumulation.

Cloud cover and photoacclimation
In addition to the close coupling between z mix , carbon production rates and changes in phytoplankton biomass accumulation (Fig. 9d,e), we observed the maximum midday irradiance to vary by a factor of three and the average light available to cells within the surface mixing layer to vary by a factor of 13. Furthermore, based on a significant log(K d -Chl) relationship established using PAR and Chl a fluorescence profiles at the Central Celtic Sea, we know that K d varied by a factor of three due to changes in in situ biomass and the effect of (chlorophyll) self-shading.
Our results suggest that the exact timing and magnitude of increases in production and accumulation rates following shoaling of z mix was moderated by day-to-day variability in cloud cover. For example, rather than observing significant increases in production rates and biomass accumulation on the 9th and 13th, when z mix shoaled and incident radiation was high, maximum daily carbon production actually occurred on the 10th and 14th April, two of the cloudiest days (Fig. 9a,c). In the early hours of the morning on the 9th April, one of the sunniest days observed (> 1000 μE m −2 s −1 at midday), z mix shoaled from 29 to 8 m depth. Rather than increase, the daily average carbon production on the 9 th remained constant at 0.15 d −1 (Fig. 9c). Production increased to over 0.23 d −1 the next day (10th April). The active mixing layer remained above 15 m until late afternoon on the 10th, but there was an increase in cloud cover that reduced the maximum midday irradiance to < 500 μE m −2 s −1 and therefore also the daily light dose received by cells in the active mixing layer (Figs. 9a, 11a). A similar sequence of events took place between the 13th and 14th April. This implies that variability in cloud cover can make a significant difference to the timing and magnitude of phytoplankton blooms by moderating the daily light dose.
To more clearly separate the effect of variable cloud cover from changes in z mix , we ran the model with two different fixed levels of maximum midday irradiance, representative of the brightest (1330 μE m −2 s −1 ) and cloudiest days (416 μE m −2 s −1 ). All other variables, including the vertical eddy diffusivities driving the particle trajectories, were kept the same. The average (and total) light dose experienced by cells within the active mixing layer over each day, for fixed values of midday irradiance, is shown in Fig. 11a. In these runs the day-to-day variability in light exposure was driven only by changes in the depth of z mix . Figure 11b shows modeled daily carbon production rates for the observed light field (as in Fig. 9c) and for the synthetic light time series. Figure 11c shows the specific biomass accumulation rate anomalies (observed minus synthetic). The acclimation status of cells with the active mixing layer is shown in Fig. 11d.
Carbon production on the 9th April when the maximum observed daily irradiance reached 1105 μE m −2 s −1 (Fig. 9a) did not increase as the mixing layer shoaled because (1) in the high light levels phytoplankton cell acclimation status was high (Y = 0.8 at midday) (Fig. 11d), reducing the maximum possible instantaneous production rates and (2) the average light experienced by the cells was close to the inhibition threshold of 200 μE m −2 s −1 (Fig. 11a), indicating that they may have become photo-inhibited. If the 9th April had instead been a cloudy day, then the daily production rate would have increased to 0.3 d −1 (yellow line in Fig. 11b), and the biomass accumulation rate would have been greater (Fig. 11c). Conversely, if there had been slightly more light then the average light exposure would have exceeded the inhibition threshold (Fig. 11a), cells would have been photoinhibited and carbon production rates would have decreased (red line in Fig. 11b), in spite of the shoaling event.
On the 13th and 14th April the depth of the active mixing layer remained above 10 m. The maximum midday irradiance however decreased from 1300 μE m −2 s −1 on the 13th to 600 μE m −2 s −1 on the 14th (Fig. 9a). This reduced the average daily light availability above z mix from close to the inhibition threshold (200 μE m −2 s −1 ), to 100 μE m −2 s −1 (blue line in Fig. 11a) which is much closer to the light levels required for cells to achieve their maximum photosynthetic performance (the saturation onset of 90 μE m −2 s −1 ). If the 13th had been cloudy, then carbon production rates and biomass accumulation would have increased instead (yellow lines in Fig. 11b,c). Conversely, if light levels had remained high on the 14th then carbon production would not have increased and biomass accumulation would have reduced.
Variability in production rates associated with day-to-day variability in cloud cover can affect the timing and magnitude of the biomass peak. The observed and modeled (using the observed light levels) biomass peak within the active mixing layer occurred during the afternoon of the 17th April (11 and 189 mg C m −3 , respectively). If each prior day had been sunny, then the biomass peak would still have occurred on the 17th, but it would have been 26% higher. In contrast, persistently cloudy days would have resulted in an equivalent peak in biomass (189 mg C m −3 ) but it would have occurred 3 d earlier, on the 14th April.
The production-acclimation model was specifically chosen for its ability to allow a phytoplankton cells instantaneous photosynthetic performance to continually vary between two limiting P-I curves, representative of fully dark-and lightacclimated states. An acclimation time-scale of 1 h (γ = 1 h) was chosen. It is reasonable to ask whether this representation of short-term acclimation was important. We address this question by running the model with a range of acclimation time scales: 1 s, 1, 3, 6, 9, 12, 24 h and 60 d. For γ = 1 s, the cell acclimates immediately at each model time step to its ambient light environment. For γ = 60 d, acclimation is effectively turned off, and the cells retain the acclimation status that they are initially prescribed (Y = 0 below 20 m, Y = 0.7 at the surface). When the acclimation time scale was 6 h or more, the model bloom occurred 3 d earlier than in our observations, and the surface water nutrient concentration decreased below 0.5 mg N m −3 more than 3 d too soon. For time scales of 3 h or less the model biomass peak occurred within a few hours of what we observed. We could not therefore have reproduced the timing of the 2015 spring bloom biomass peak in the Celtic Sea without an appropriate production-acclimation model and a suitable choice of photoacclimation time scale. Overall, we have shown that the effect of light on bloom dynamics is the result of complex feedbacks between active mixing depth, cloud cover, short-term photoacclimation, and in situ biomass.

Physical vs. biological control
While we find physical processes to be the most important control on phytoplankton growth during the spring bloom, biological processes, particularly losses through predator-prey relationships, are also relevant when considering biomass accumulation rates and termination of phytoplankton blooms. Since the model was run with fixed respiration, grazing, and mortality rates (Table 1), differences between our observations and model are indicators of when biological rather than physical controls were more important in driving changes in biomass.
The modeled and observed specific biomass accumulation rates (Fig. 9b) are different signs late in the day on the 18th April and then again from late afternoon on the 21st until the 23rd April. On both occasions the observed specific biomass accumulation rate is negative, indicating that loss rates exceeded phytoplankton growth at these times, whereas the model showed weakly positive rates. Divergence between the model and observations is also seen in Fig. 5d; the modeled depth integrated biomass within the active mixing layer between the 21st and 23rd April remains high whereas the observed biomass steadily decreases. These periods coincide with observed changes in the abundance, biomass, net growth and grazing rates of zooplankton at the Central Celtic Sea site during spring 2015 (Giering et al. 2019;Mayers et al. 2019), suggesting that at specific times grazing pressure exerted topdown control on the bloom. Between midday on the 21st and 22nd April a 20% reduction in depth integrated biomass within the active mixing layer was observed. A grazing rate of 0.2 d −1 , an order of magnitude larger than the fixed rate imposed in the model, would have been required at this time to achieve the equivalent 20% loss in model biomass. An order of magnitude increase in micro-zooplankton grazing rate throughout April was reported by Mayers et al. (2019) for the same bloom. Viral mortality, although not measured, may also have become more important following the peak of the bloom. Variable grazing and mortality rates are therefore reasonable explanations of the observation-model mismatches.
Despite the fixed rates of respiration, grazing and mortality and the imposed coupling between biomass and cell division, such convincing replication of the bloom leads us to conclude that variable grazing rates are not a first order control during the build-up of the spring bloom in shallow, temperate shelf seas. However, these loss terms become increasingly important once the peak in biomass has been reached and grazing pressure becomes significant.

Wind driven and tidal advection
The data presented here provides only a 1-D Eulerian perspective of phytoplankton growth. Therefore, it is reasonable to ask whether advective processes, moving patchy phytoplankton growth around, could contribute to the changes in Chl a biomass that we observed.
Over timescales of days to months the wind sets up shelfscale circulations and is able to strain horizontal temperature and salinity gradients across the Celtic Sea (Ruiz-Castillo et al. 2019a, 2019b. During early spring, prevailing westerly winds frequently drive a surface layer of relatively cold and freshwater from the inner shelf southwards, and a compensatory northward return flow of warmer and more saline water at depth. This process plays an important role in delivering nutrients with an oceanic origin to the inner shelf (Ruiz-Castillo et al. 2019b) and in determining the onset of permanent seasonal stratification (Ruiz-Castillo et al. 2019a). The O (0.1 C) changes in bottom water temperature seen in Fig. 4a are consistent with wind driven advection. A wind driven surface current of 0.01 m s −1 , sustained over 1 month, could advect a parcel of water up to 30 km. Based on 4 km resolution ocean color satellite images, within a 30 km range of the Central Celtic Sea, the maximum horizontal Chl a gradient was 0.07 mg Chl m −3 km −1 . A 0.01 m s −1 current could therefore have introduced no more than a 2.5 × 10 −3 mg Chl m −3 h −1 change in Chl a at a fixed sampling location. This is two orders of magnitude less than the mean rate of change of surface Chl a recorded by the glider (0.15 mg Chl m −3 h −1 ) and three orders of magnitude less than the maximum rate of increase (1.2 mg Chl m −3 h −1 ). Hence, wind driven surface circulation could not introduce significant variability into our timeseries of Chl a biomass.
The Celtic Sea has a strong semi-diurnal tide (Fig. 3a). This is reflected in the bottom water temperature (Fig. 4a) where a 5-10 km tidal excursion drove 0.025 C amplitude temperature oscillations. A 0.35 m s −1 tidal current advecting a 0.07 mg Chl m −3 km −1 gradient in surface Chl a past a fixed point could introduce increases and decrease in Chl a concentration of 0.1 mg Chl m −3 h −1 . In fact, we see this throughout the observational record. This known variability however is superimposed on top of much greater changes in surface Chl a that cannot be attributed to advection. On the 10th, 14th, and 17th April increases in surface Chl a concentration of 0.4, 1.2, and 1.2 mg Chl m −3 h −1 , respectively, were observed. Thus, during the periods of reduced wind stress and shoaling of z mix , when we observe (and model) increases in growth rate and the accumulation of biomass, tidal advection is not the leading order driver of the signal.
The 1D model used in this study cannot account for advection. This is particularly noticeable after the biomass peak on the 17th April when there is a strong semi-diurnal signal in the observed depth integrated Chl a that is not replicated by the model (Fig. 5d). The strong spring tidal currents in combination with high levels of (presumably spatially patchy) biomass during this period would accentuate the mismatch. Importantly though, the model still captures the overall decrease in biomass within the active mixing layer between the 18th and 20th April, followed by an increase on the 21st April.

Entrainment of phytoplankton into bottom waters
Although the biomass specific accumulation rates between the 5th and 9th April were positive within the active mixing layer, including a local maximum on the 7th April, there was no local biomass maximum above either z mix or z pyc (Fig. 5c). Instead, the depth-integrated biomass above the pycnocline decreased from 88 mg Chl m −2 on 5th April to 78 mg Chl m −2 on the 8th April (not shown). So, if phytoplankton were blooming (actively growing) above z mix , where did the biomass accumulate? Between the 5th and 8th April, there was a 3 mg Chl m −2 increase in the total depth integrated biomass (from 118 to 121 mg Chl m −2 ), additional evidence that phytoplankton growth was occurring (Fig. 6a). Over the same period, the depth integrated Chl a below the pycnocline increased from 30 to 43 mg Chl m −2 , which is a 13 mg Chl m −2 increase and accounts for the 10 mg Chl m −2 lost from above the pycnocline, as well as the 3 mg Chl m −2 of extra (Chl a) biomass produced. Figure 6b shows the percentage of the total water column integrated Chl a that was contained within the active mixing layer, above the pycnocline and in bottom waters. On the 5th April, 75% of the total water column biomass was above the seasonal pycnocline while the remaining 25% was contained within the bottom layer. On the 8th April, only 64% of the total biomass was above z pyc , while the remaining 36% was in the bottom water. The loss of biomass from above the pycnocline and the coincident increase in biomass below coincided with spring tides (Fig. 3a), a weak early season pycnocline (Fig. 4b) and tidally driven dissipation in the bottom boundary layer that extended to the base of the pycnocline (Fig. 4e). Figure 12a shows biomass being stripped from the base of the pycnocline and drawn into bottom water between the 5th and 9th April. Each erosion event happened just after peaks in the depth-integrated dissipation beneath the pycnocline (vertical dashed lines in Fig. 12). Since there was no light available to support new growth beneath the pycnocline, either within the bottom water or at the seabed, the source of high Chl a concentration (bottom layer depth average of 0.3-0.5 mg m −3 , see Fig. S2) must have been the surface water above.
Although dissipation during the second spring tide on 20th April was greater (Fig. 12b), comparatively little surface biomass was eroded. At this time, the strength of the pycnocline was an order of magnitude stronger (Fig. 4b), and most of the phytoplankton biomass was contained within the active mixing layer, 10 m or more above z pyc (Figs. 5a, 6b). The decoupling of z mix and z pyc therefore prevented biomass within the active mixing layer from being eroded and entrained into the bottom water. Although very little biomass was stripped from the sunlit surface water at this time, the average Chl a concentration within the bottom layer remained relatively high (0.2 mg m −3 , see Fig. S2). This is because phytoplankton cells stripped from the surface layer between 5th and 9th April continued to fluoresce as they were homogenously mixed across the bottom layer by strong tidal currents (Figs. 3a,4e).
Tidally driven shear stress propagating upwards from the seabed and eroding biomass from subsurface phytoplankton populations has been observed before in shelf systems (Sharples et al. 2001), but never for such a prolonged period. It is a physical erosion process unique to shallow, tidally active shelves where the turbulent bottom boundary layer is able to extend to the base of the pycnocline. It is difficult to quantify what impact this event may have had on the timing and peak magnitude of the bloom, although given that seasonal stratification was established on 26th March, 4 days after maximum spring tides it is probable that the erosion event we observed here was the first of the season.

Conclusions
We have presented observations and a model that show how the growth of phytoplankton and accumulation of biomass during the spring bloom period in a temperate shelf sea is tightly coupled to changes in the depth of the surface active mixing layer and day-to-day variability in cloud cover. During a 2-week period leading up to the blooms' biomass peak we observed two different mixing regimes. First, a diurnal regime of growth within the sunlit surface waters, followed by nighttime convection evenly redistributing new biomass to the base of the seasonal pycnocline (Fig. 10a). Then, a transition to a shallower wind-mixing regime where more rapid phytoplankton growth and biomass accumulation took place within the active mixing layer. Periods of reduced wind stress and shoaling of the active mixing layer stimulated peaks in growth rate and accumulation (Fig. 10b). Our results demonstrate that physical processes that mediate light availability are the leading order control on early spring phytoplankton growth. We also find that the ability of phytoplankton cells to acclimate to their light environment on hourly time scales to be an important control on the timing and magnitude of phytoplankton biomass peaks during the spring bloom. Lastly, we observed tidally-driven turbulence eroding phytoplankton biomass from the base of the seasonal pycnocline, entraining it into the bottom water, a process unique to shelf sea environments. These new insights into the significance of sub-daily changes in meteorological forcing, irradiance and turbulent mixing to bloom dynamics provide a detailed process level of understanding that has not been possible from coarser temporal and spatial resolution data sets. While it is inevitable that the spring bloom across temperate shelves evolves differently each year, it is likely that nighttime convection and wind-driven mixing will be of leading order importance.
Our key results are relevant across both continental shelves and the open ocean. Globally, day-to-day variability in local cloud cover is an order of magnitude greater than inter-annual variability (Stubenrauch et al. 2013). The sensitivity of phytoplankton growth rates and biomass accumulation to day-to-day changes in incidence irradiance is therefore highly relevant and important on global ocean scales, notably across cloudy polar and subpolar regions, and those in the path of storm tracks. Further, the passage of weather systems and storms drives hourly to daily changes in the depth of the active mixing layer across whole ocean basins. Under nutrient replete conditions (i.e., pre-bloom) the tight coupling between the depth of wind-driven mixing, light availability and biomass accumulation that we report here is highly important.
Our findings present a number of challenges for coupled hydrodynamic-ecosystem models if we are to make accurate predictions of the timing and magnitude of phytoplankton blooms. Firstly, mixing parameterizations need to accurately resolve changes in the depth of surface mixing, on sub-daily time scales. This is particularly important when the water column is experiencing diurnal stratification, nighttime convection and episodes of wind-driven mixing. Secondly, ensuring that cloud cover and light attenuation are realistic and that phytoplankton have the ability to photo-acclimate on hourly timescales should be a priority.
As a consequence of anthropogenic global warming, many of the physical drivers of phytoplankton bloom dynamics (e.g., heat fluxes, stratification, cloud cover, precipitation, and wind stress) are changing (IPCC 2013). The canonical view is that increasing stratification will reduce nutrient supply and thus reduce net primary production (Steinacher et al. 2010). Projected regional trends in net primary production across the NW European Shelf however are both positive and negative (Holt et al. 2016). Our results suggest that on a nutrient replete shelf, which fully mixes during the winter, changes in wind stress and local weather conditions (cloud cover) may have more influence over future changes in the timing, magnitude and duration of the spring phytoplankton bloom than the strength of the seasonal pycnocline.