Mixing and Phytoplankton Growth in an Upwelling System

Previous studies focused on understanding the role of physical drivers on phytoplankton bloom formation mainly used indirect estimates of turbulent mixing. Here we use weekly observations of microstructure turbulence, dissolved inorganic nutrients, chlorophyll a concentration and primary production carried out in the Ría de Vigo (NW Iberian upwelling system) between March 2017 and May 2018 to investigate the relationship between turbulent mixing and phytoplankton growth at different temporal scales. In order to interpret our results, we used the theoretical framework described by the Critical Turbulent Hypothesis (CTH). According to this conceptual model if turbulence is low enough, the depth of the layer where mixing is active can be shallower than the mixed-layer depth, and phytoplankton may receive enough light to bloom. Our results showed that the coupling between turbulent mixing and phytoplankton growth in this system occurs at seasonal, but also at shorter time scales. In agreement with the CTH, higher phytoplankton growth rates were observed when mixing was low during spring-summer transitional and upwelling periods, whereas low values were described during periods of high mixing (fall-winter transitional and downwelling). However, low mixing conditions were not enough to ensure phytoplankton growth, as low phytoplankton growth was also found under these circumstances. Wavelet spectral analysis revealed that turbulent mixing and phytoplankton growth were also related at shorter time scales. The higher coherence between both variables was found in spring-summer at the ~16–30 d period and in fall-winter at the ~16–90 d period. These results suggest that mixing could act as a control factor on phytoplankton growth over the seasonal cycle, and could be also involved in the formation of occasional short-lived phytoplankton blooms.

Previous studies focused on understanding the role of physical drivers on phytoplankton bloom formation mainly used indirect estimates of turbulent mixing. Here we use weekly observations of microstructure turbulence, dissolved inorganic nutrients, chlorophyll a concentration and primary production carried out in the Ría de Vigo (NW Iberian upwelling system) between March 2017 and May 2018 to investigate the relationship between turbulent mixing and phytoplankton growth at different temporal scales. In order to interpret our results, we used the theoretical framework described by the Critical Turbulent Hypothesis (CTH). According to this conceptual model if turbulence is low enough, the depth of the layer where mixing is active can be shallower than the mixed-layer depth, and phytoplankton may receive enough light to bloom. Our results showed that the coupling between turbulent mixing and phytoplankton growth in this system occurs at seasonal, but also at shorter time scales. In agreement with the CTH, higher phytoplankton growth rates were observed when mixing was low during spring-summer transitional and upwelling periods, whereas low values were described during periods of high mixing (fall-winter transitional and downwelling). However, low mixing conditions were not enough to ensure phytoplankton growth, as low phytoplankton growth was also found under these circumstances. Wavelet spectral analysis revealed that turbulent mixing and phytoplankton growth were also related at shorter time scales. The higher coherence between both variables was found in spring-summer at the ∼16-30 d period and in fall-winter at the ∼16-90 d period. These results suggest that mixing could act as a control factor on phytoplankton growth over the seasonal cycle, and could be also involved in the formation of occasional short-lived phytoplankton blooms.

INTRODUCTION
Marine phytoplankton is responsible for about half of the primary production in the biosphere (Field et al., 1998), and therefore plays a key role in the cycling of matter and energy on Earth. The two main resources limiting phytoplankton growth, light, and nutrients availability, are strongly dependent on turbulent mixing conditions in the water column. By controlling the vertical displacement of cells, mixing determines their exposure to solar photosynthetic active radiation. In addition, turbulent mixing is one of the mechanisms responsible for the supply of inorganic nutrients from deep waters to the surface, where they can be taken up by phytoplankton. For this reason, turbulent mixing has been frequently invoked in the formulation of theoretical models, either to explain the dominance of different phytoplankton functional groups (Margalef, 1978) or the behavior of individual phytoplankton cells (Sverdrup, 1953).
The Critical Depth Hypothesis (CDH) formulated by Sverdrup (1953) to explain the onset of the North Atlantic spring bloom, pioneered the development of conceptual models relating phytoplankton growth, mixing conditions, and light availability. The CDH concluded that deep mixed-layers during winter keep phytoplankton in an unfavorable light environment and therefore limit their production. When solar heating thins the non-stratified mixed-layer, phytoplankton could have the potential to bloom because growth outweighs loses due to the increase of solar exposure. He defined the "critical depth" as the lower limit of the water column in which depth-integrated production of organic matter equals its oxidation by respiratory processes. Since its formulation, several studies have attempted to verify the CDH in the field with controversial results. Some of them found evidence supporting the CDH (Semina, 1960;Menzel and Ryther, 1961;Nelson and Smith, 1991;Obata et al., 1996;Siegel et al., 2002;Chiswell, 2011;Brody et al., 2013;Chiswell et al., 2013Chiswell et al., , 2015Lozier, 2014, 2015;Wihsgott et al., 2019;Hopkins et al., 2021), whereas others rejected it based on the observation that phytoplankton growth rate was positive during deep winter mixing (Acuña et al., 2010;Behrenfeld, 2010;Boss and Behrenfeld, 2010;Behrenfeld et al., 2013;Behrenfeld and Boss, 2014;Arteaga et al., 2020). A recent study revealed that although phytoplankton starts growing in early winter at weak rates, a proper bloom initiates only in spring when atmospheric cooling subsides and the mixedlayer rapidly shoals (Mignot et al., 2018). The observation that phytoplankton blooms sometimes occur in the apparent absence of water column stratification (Townsend et al., 1992;Ellertsen, 1993;Backhaus et al., 1999;Körtzinger et al., 2008) led to propose alternative bloom formation processes (Franks, 2015). The Critical Turbulence Hypothesis (CTH, Huisman et al., 1999a,b) proposed that if turbulence is low enough, phytoplankton in the well-lit surface layer could bloom independently of the thickness of the mixed-layer. Sverdrup (1953) was aware of the distinction between a mixed-layer, defined as a subsurface layer of relatively uniform temperature or density, and the active mixed-layer, defined by the intensity of turbulent diffusivity (K). However, by explicitly assuming that it was large enough to be ignored, he avoided the need to include K in his model. Franks (2015) emphasized that using the theoretical background of the CDH requires observations of microstructure turbulence, rather than mixed-layer depth estimates derived from thermohaline properties. He also noted that it is crucial to determine, not only the intensity of turbulence, but also its vertical structure and temporal variability. In a recent study Hopkins et al. (2021) employed 2 weeks of sub-hourly observations of turbulent kinetic energy dissipation rate to demonstrate the critical role that the strength and structure of turbulent mixing play in governing the development of spring phytoplankton blooms in the Celtic Sea. According to these authors their results could be applicable to any region where wind-driven mixing can modify nutrient and light availability.
Eastern boundary upwelling systems (EBUS) are complex regions where the interaction between hydrographic conditions and phytoplankton growth occurs within a broad range of temporal scales (Pitcher et al., 2010). Several studies have investigated the role of the major physical processes that may control biological productivity in these regions (Messié and Chavez, 2015). Patti et al. (2008) suggested that several driving factors, as nutrients concentration, light availability, shelf extension, and surface turbulence estimated from wind speed, must be considered when investigating the phytoplankton biomass distribution. By using Finite Size Lyapunov Exponents and satellite data, Rossi et al. (2009) described a global negative correlation between surface horizontal mixing and chlorophyll. Fearon et al. (2020) used a 1D modeling approach to predict vertical mixing from wind speed in St Helena Bay (Bengala upwelling system), and emphasized the role of land-sea breeze in the development of phytoplankton blooms. As far as we know direct observations of microstructure turbulence have been never used to investigate the role of mixing in phytoplankton bloom formation in EBUS.
The Ría de Vigo is a long narrow embayment located in the northern boundary of the Canary Current upwelling ecosystem. Intense and intermittent upwelling events occur mainly in spring and summer (Fraga, 1981). In these seasons, the prevailing northerly winds cause offshore Ekman transport of surface water, and the rise of cold nutrient-enriched subsurface waters, which stimulate phytoplankton growth and support a highly productive pelagic ecosystem (Fréon et al., 2009). During winter, southerly winds driving downwelling conditions are dominant. The annual cycle of phytoplankton biomass corresponds to a typical temperate shelf sea, with the development of spring and autumn blooms (Álvarez-Salgado et al., 1996;Nogueira et al., 1997;Moncoiffe et al., 2000;Cermeño et al., 2006). On the other hand, short-term variability of the upwelling regime and runoff water pulses induce changes on phytoplankton biomass and composition over shorter time-scales (Nogueira et al., 2000;Nogueira and Figueiras, 2005). The main goal of this study is to investigate the relationship between vertical turbulent mixing and phytoplankton growth at the different temporal scales involved in the coupling between physical and biological processes in this coastal upwelling system.

Sampling Site
In the framework of the REMEDIOS project (RolE of Mixing on phytoplankton bloom initiation, maintEnance, and DIssipatiOn in the Galician ríaS, http://proyectoremedios.com/inicio), 52 samplings were carried out on board R/V Kraken at station EF located in the inner part of the Ría de Vigo (8.778 o W 42.235    2018, approximately once a week. During each sampling, hydrographic, and microstructure turbulence profiles, as well as samples for the determination of inorganic nutrients, chlorophyll a, and primary production were collected. A continuous recording of current velocity profiles during the sampling period was acquired with an upward-looking bottom-moored Acoustic Doppler Current Profiler (ADCP). The weekly samplings were clustered into three groups by virtue of the predominance of upwelling (U), downwelling (D), or the transition between both conditions (T). This classification was based on the upwelling index (UI), the horizontal currents, and hydrographic information. The upwelling index was calculated as the offshore Ekman transport (Bakun, 1973) in the direction perpendicular to the shoreline (Gomez-Gesteira et al., 2006):  (Morel, 1988).

Hydrography and Turbulence
Hydrographic and turbulent data were collected with a microstructure turbulence profiler MSS90 (Prandke and Stips, 1998). The profiler was equipped with two microstructure shear sensors (type PNS06), a microstructure temperature sensor (FP07), a high-precision CTD (Conductivity-Temperature-Depth) probe, a fluorescence sensor, and an accelerometer. On each sampling day, 10 profiles were conducted and then averaged. The profiler was balanced to have negative buoyancy in the water column and a sinking velocity in the range 0.4-0.7 m s −1 . Dissipation rates of turbulent kinetic energy (ǫ) were computed in 512 data point segments, with 50% overlap, from the vertical shear (∂ z u) variance using the Taylor (1935) equation assuming isotropic turbulence: where ν is the kinematic viscosity of seawater and denotes the ensemble average. The shear variance was computed by integrating the shear power spectrum. The lower integration limit was set to 2 cpm. The upper cut-off wavenumber for the integration of the shear spectrum was set as the Kolmogoroff number k c = (2π) −1 ǫ 1 4 ν − 3 4 . An iterative procedure was applied to determine k c . The maximum upper cut-off was not allowed to exceed 30 cpm to avoid the noisy part of the spectrum. Assuming a universal form of the shear spectrum, ǫ was corrected for the loss of variance below and above the integration limits, using the polynomial functions reported by Prandke et al. (2000). ǫ values were then averaged in 1 m bins. Peaks due to particle collisions were removed by comparing the dissipation computed simultaneously from the two shear sensors. The turbulent diffusivity K was estimated from the Osborn (1980) formula: where N 2 is the squared buoyancy frequency and γ is the mixing efficiency, here considered to be 0.2 (Oakey, 1982). Although a growing body of evidence suggests that γ is not always constant, we follow the recommendation of Gregg et al. (2018) to use γ = 0.2 for microstructure studies until observations, laboratory experiments, and numerical simulations converge on a more accurate formulation. Furthermore, the Ría de Vigo is a system characterized by being marginally unstable to sheardriven turbulence (Fernández-Castro et al., 2018), which justifies the choice of this γ value (Smyth, 2020). The mixed-layer depth (MLD) was determined as the depth where an increase of 0.125 kg m −3 was observed with respect to the surface values (∼ 2 m). The turbulent layer depth (TLD) was computed using one-dimensional lagrangian simulations forced with the observed diffusivity profiles (Ross and Sharples, 2004). The new particle depth (z n+1 ) was computed in based of the present depth (z n ) after: z n+1 = z n + t∂ z K(z n ) + R 2 tK z n + 1 2 t∂ z K(z n ) r where R ∈ [−1, 1] is a random number of variance r = 1 3 and t the time-step chosen by ensuring that it was much lower than the minimum value of |(∂ zz K) −1 |. In these simulations, 1,000 particles were released at the surface and let evolve in the diffusivity field for 24 h. The TLD was defined as the lower limit of the depth range containing 95% of 1,000 particles at the end of the simulation. The light availability was computed as the dailyaverage surface PAR in the turbulent layer computed following the expression reported in (Vallina and Simó, 2007): where κ is the light attenuation coefficient determined from PAR profiles obtained with a Licor PAR sensor on each sampling day.

Current Velocities
The RD Instruments Acoustic Doppler Current Profiler (ADCP) was bottom-moored in the central part of the Ría de Vigo close to the EF station (8.761 o W 42.241 o N, ∼ 45 m depth, Figure 1). Profiles of current velocity were acquired continuously from 3 March 2017 to 29 May 2018. A gap in data collection occurred from 30 May 2017 to 8 June 2017 due to maintenance operations. The measurements were made in 89 layers with a bin thickness of 0.5 m and the first bin located at 2 m above the ADCP transducers. Every 10 min, 85 samples of 3D currents were acquired at 2 Hz, averaged in ensemble and recorded.
Resulting zonal and meridional velocity vector components were projected into the axis along the main channel (35 o north of east) of the Ría de Vigo. Hence, the along-Ría component was positive into the Ría. Finally, these time-series were smoothed with the A 2 24 A 25 operator (Godin, 1972), which applies three consecutive moving averages with window size of 24, 24, and 25 h, with the aim of filtering out the tidal and supertidal frequencies. The residual current velocity along the main axis of the Ría at 32 m depth (u 32 ) was selected to characterize the deep circulation in the Ría. Positive (negative) values indicate the deep water inflow (outflow) into the Ría and surface water outflow (inflow), corresponding to upwelling and downwelling conditions, respectively.

Inorganic Nutrients Concentration
Samples for the determination of dissolved inorganic nutrients (nitrate, nitrite, ammonium, phosphate, and silica) were collected from 8 to 9 depths, in the water column by using 5 l Niskin bottles. Samples were frozen at −20 o C until further determination at the laboratory, following the methods described by Hansen and Koroleff (2007).

Nutrient Supply
Nutrient input into the surface layer of the Rías occurs mainly by coastal upwelling, while continental runoff and precipitation  and turbulent diffusion (Moreira-Coello et al., 2017) represent minor inputs. An estimate of the supply of nitrate due to upwelling pulses was computed as the nitrate vertical advection induced by convergence of upwelled waters inside the Ría: where l and A are the length of the mouth and surface area of the Ría (∼ 10 km and ∼ 174 km 2 , respectively), and NO − 3 (40) is the nitrate concentration at the deepest sampling depth around 40 m depth. UI is the upwelling index averaged over the 3 days prior each sampling (Gilcoto et al., 2017). A null value was assigned to those samplings under downwelling conditions.

Chlorophyll a Concentration and Primary Production Rates
Chlorophyll a concentration (Chl a) and primary production rates were determined from samples collected at surface (∼0.5-1 m) and 10 m depth. For the determination of total Chl a, 250 ml of seawater were filtered through polycarbonate filters of 0.2 µm, which were later frozen at −20 o C until analysis. The fluorescence (Fluo) emitted by the Chl a was measured from pigments extracted in 90% acetone at 4 o C overnight using a Turner designs Trilogy fluorometer previously calibrated with pure Chl a. The measured Chl a concentrations were used to calibrate the MSS fluorometer by using a linear regression (Chl a = 1.47 × Fluo − 0.16, R 2 = 0.89, n = 22) for Chl a concentration ranging between 0.26 and 5.83 mg m −3 .
Primary production rates were determined by running incubations with the radioisotope 14 C as described in Cermeño et al. (2016). Briefly, four 72 ml acid-washed polystyrene bottles (three light and one dark bottle) were filled with seawater from each depth. Each bottle was inoculated with ∼ 5 µCi of NaH 14 CO 3 and then incubated for 2 h starting at noon. An incubator equipped with a set of blue and neutral density plastic filters was used to simulate irradiance conditions at the original sampling depths. Temperature conditions during the incubation period were kept similar to those observed at surface (±1.5 o C) by employing a closed refrigerated water system. Immediately after incubation, samples were filtered through 0.2 µm polycarbonate filters under low-vacuum pressure. Non-assimilated radioactive inorganic carbon retained in the filters was removed by exposing them to concentrated HCl fumes overnight. Radioactivity signal on each sample was determined on a 1409-012 Wallac scintillation counter, which used an internal standard for quenching correction.
Growth rates (µ) were computed as the ratio primary production to phytoplankton biomass. Carbon phytoplankton biomass was estimated assuming the C:Chl a ratios for three hydrographic conditions in the Ría de Vigo: stratification (55 ± 13), winter mixing (54 ± 27) and upwelling (35 ± 18) at the same station (Cermeño et al., 2005). In our case 35 ± 18 was used for upwelling and 55 ± 30 for downwelling and transitional periods. Daily growth rates were computed from hourly values considering the number of light hours on the sampling day.

Wavelet Analysis
Time-series of phytoplankton growth rates and turbulent layer depth were analyzed with the wavelet transform, a technique in which the evolution of the times-series variance is analyzed at different frequencies. The contribution of variance at different times and periods was determined by the wavelet power spectrum (WPS), which was represented in the time-period plane, the periodogram (Cazelles et al., 2008). The finite length of timeseries provokes discontinuities at their borders delimited by the cone of influence, in which the edge effects are negligible (Torrence and Compo, 1998). The relation of two time-series can be identified in time and period in terms of coherence (Ŵ 2 , Liu, 1994), from 0 to 1 due to a smoothing in both time and period following the recommendation of Torrence and Webster (1999). Peaks in the spectrum can be assumed to be a true feature with a certain confidence by using a test against an assumed background level (Torrence and Compo, 1998). In this case, the background was simulated with 1000 replicas of surrogate times-series by using a hidden Markov model (HMM), in which the shortterm temporal correlation of the original time-series is preserved (Cazelles and Stone, 2003). The null hypothesis tested is that the original and surrogate time-series had same value distribution and short-term autocorrelation (Cazelles et al., 2014) with 95% confidence level.
To apply the wavelet transform, the time-series have to be spaced in a regular time interval (δt). Hence, phytoplankton growth rates and turbulent layer depth were linearly interpolated onto an equi-spatially time-grid with δt ≈ 8 d and N = 52 points. The comparison between times-series was ensured by using the standardization of the variables, so they had nil mean and unit variance. An specific code was built on Python in order to conduct the wavelet analysis.

Environmental and Hydrographic Conditions
The temporal variability of the smoothed (fortnight moving average) upwelling index and the smoothed current velocity sampled at 32 m depth allowed us to establish upwelling, downwelling, or transitional periods (Figure 2). Based on the bidirectional circulation of the Ría de Vigo, positive values of UI and deep current velocity correspond with the input of deep water, and negative values with the output of deep water. From sampling 1 to 11 (9 Mar 2017-18 May 2018) a spring transitional period (T 1 ) was sampled characterized by the alternation of relatively short upwelling and downwelling events, which resulted in low averaged upwelling index of UI = 0.3 ± 0.1 m 2 s −1 (mean ± standard uncertainty of the mean, Table 1). Samplings 12-26 (25 May 2017 to 14 Sep 2017), when the averaged upwelling index was positive (UI = 1.7 ± 0.1 m 2 s −1 ), were characterized by springsummer upwelling conditions (U 1 ). Samplings 27-36 (21 Sep 2017-21 Nov 2017) recorded a fall transitional period (T 2 ), with alternation of upwelling and downwelling events and null averaged upwelling index (UI = 0.0 ± 0.1 m 2 s −1 ). Between late fall and early spring (samplings 34-47, 30 Nov 2017-27 Mar 2018, downwelling conditions were dominant (D) and averaged upwelling index was negative (UI = −1.1 ± 0.1 m 2 s −1 ). Finally, in spring 2018 (samplings 48-52, 12 Apr 2018-10 May 2018) the beginning of the spring upwelling season (U 2 ) was sampled, characterized by a mean upwelling index (UI = 1.9 ± 0.4 m 2 s −1 ) similar to the previous upwelling period (U 1 ).
The turbulent layer depth (TLD), an indication of the depth reached by currently active mixing (Brainerd and Gregg, 1995) computed by simulations forced with K profiles (see section 2), was significantly correlated with the mixed-layer depth (TLD = 0.48 × MLD + 5.4, R 2 = 0.66, p < 0.001, n = 52), and exhibited similar seasonal variability (Figure 3). As a result of the enhanced K during samplings 36 and 38, maximum TLD values were computed during the fall transitional T 2 (13 ± 2 m) and winter downwelling D (12 ± 3 m) periods, whereas TLD was shallower during the spring-summer upwelling U 1 (7.1 ± 0.4 m). Most observations (n = 40) corresponded to shallow (< 9 m) MLD coinciding with slightly deeper (3-4 m) TLD. This could be the result of the limited number of microstructure turbulence data near the surface used for the calculation of both TLD and MLD. For samplings showing relatively deep mixed-layers (MLD > 9 m, n = 12), TLD was 8 ± 2 m shallower than MLD in average. In fact, in four sampling days, TLD was at least 12 m shallower than MLD (samplings 3, 31, 36, and 37), which suggests the occurrence of mixed-layers that were not actively mixing at the time of sampling.
The temporal variability of Chl a concentration, primary production rates and phytoplankton growth rates determined at the surface and 10 m depth from water samples collected from Niskin bottles are shown in Figure 4. All these variables exhibited a distinct seasonal variability, which allowed fitting the data to a sinusoidal curve to highlight the seasonal cycle. The amplitude of the seasonal variability for Chl a and primary production was larger at 10 m. In agreement with the chlorophyll-derived fluorescence data from the microstructure profiler, maximum values of extracted Chl a and primary production were measured during spring-summer upwelling, whereas minimum values were obtained mainly during winter downwelling. On average, Chl a and primary production were lower at the surface than at 10 m depth during the spring transition T 1 and spring-summer upwelling U 1 , whereas the opposite pattern was found during the fall transition T 2 , the winter downwelling and the spring upwelling U 2 .
As the result of the observed variability, phytoplankton growth rates were in general higher at the surface than at 10 meters depth. On average, surface rates were higher during spring-summer upwelling U 1 (1.4 ±0.2 d −1 ), downwelling (1.1 ± 0.2 d −1 ) and spring upwelling U 2 (1.1 ± 0.1 d −1 ), whereas lower rates were computed during transitional periods T 1 (0.8 ± 0.2 d −1 ) and T 2 (0.9 ± 0.1 d −1 ). At 10 m, growth rates were higher during the spring transitional period T 1 (1.4 ± 0.2 d −1 ) and spring-summer upwelling U 1 (1.0 ± 0.1 d −1 ), whereas FIGURE 4 | Time-series of primary production rates (PP, red), chlorophyll a (Chl a, green) and growth rates (µ, blue) at (A) surface (z ∼ 0m) and (B) 10 m depth. Data for each sampling are represented by circles connected by a dashed line, and the seasonal mode (computed by fitting to a sinusoidal curve), with solid line. The ticks and numbers on the top axis indicate the sampling number and the colors, the periods for upwelling (U 1 and U 2 , in red), downwelling (D, in blue), and transition (T 1 and T 2 , in white). lower values were computed for the transitional T 2 (0.72 ± 0.08 d −1 ) and downwelling period (0.7 ± 0.2 d −1 ). Since higher phytoplankton growth rates were estimated at the surface layer and a similar seasonal pattern was observed at the two sampling depths, we focus on surface rates to investigate the relationship between phytoplankton growth and mixing.

Relationship Between Turbulent Mixing and Phytoplankton Growth
We first investigated the relationship between turbulent mixing and surface phytoplankton growth rate over seasonal scales by plotting the variability of the TLD and phytoplankton growth in Figure 5. The annual mode, computed for both variables by fitting the observations to a sinusoidal curve, showed that in general maximum growth rates coincided with shallow turbulent layers, whereas minimum values corresponded to deep turbulent layers. Similarly, high phytoplankton growth was coincident with periods of high values of the light and nitrate availability indices, and vice versa.
In addition to the seasonal mode described above, Figure 5 evidenced the existence of short-term variability both in TLD and phytoplankton growth. In fact, 44 and 60% of the variance of TLD and surface growth rates, respectively, were distributed in periods shorter than 128 d. The wavelet analysis, used to investigate this shorter-term variability (Figure 6), showed that TLD variance was larger during the fall transition T 2 and winter downwelling through all the periods with maxima at 34 and 91 days ( Figure 6A). The variance of surface growth rates displayed two main components extending mostly during the springsummer upwelling U 1 centered at 26 and 56 days ( Figure 6B). Growth rates at 10 m (Supplementary Figure 2A) exhibited variance maxima during the spring transition T 1 (centered at 25 days), and during the fall transition T 2 and winter downwelling (109 days period).
The relationship between the variance in TLD and surface phytoplankton growth rates was investigated via the coherence (Ŵ 2 ) computed using the wavelet analysis ( Figure 6C). Also, the phase difference (δϕ) between both signals was calculated to characterize the delay between the two time-series. The variables are in phase (δϕ = 0 o ) when their maxima coincide in time, whereas they are in phase opposition (δϕ = 180 o ) when the maximum of one variable coincides with the minimum of the other, and vice versa. TLD and surface growth rates showed significant coherence during the end of the spring transitional T 1 and the beginning of the spring-summer upwelling U 1 (∼16-30 d period) and from the end of the fall transitional period T 2 to the middle of the winter downwelling (∼16-90 d periods).
The highest coherence (Ŵ 2 > 0.8) occurred for periods ∼16-20 d, for which the averaged phase difference was 145 ± 4 o . The coherence between TLD and growth rates at 10 m depth (Supplementary Figure 2B) was significant during the end of the spring transitional T 1 and the beginning of the spring-summer upwelling U 1 (∼16-56 d period) and from the fall transitional period T 2 to the spring upwelling U 2 with maximum values at ∼ 16 and ∼ 35 d.

DISCUSSION
As far as we know this study represents the first attempt to investigate the relationship between mixing and phytoplankton growth, by using coincident observations of microstructure turbulence and phytoplankton growth, in a coastal upwelling system. In general, higher turbulent mixing was observed during fall-winter transitional and downwelling periods, whereas weaker mixing corresponded to spring-summer transition and upwelling. However, this pattern was mainly driven by the enhanced turbulent diffusivity quantified during two samplings (36 and 38 sampled on 21 November 2017 and 5 December 2017, respectively). This increase in diffusivity was mainly due to a decrease in the buoyancy frequency, as the consequence of the weak vertical gradient observed in both temperature and salinity. According to the analysis of hydrographic and current velocity data, these two samplings were under the influence of fall upwelling conditions. Although this system is characterized by predominant upwelling conditions during the spring and summer (Fraga, 1981), previous studies described the existence of an important short-term variability in wind conditions, which generates a succession of upwellingdownwelling episodes in periods of 15 days (Álvarez-Salgado et al., 1993;Nogueira et al., 1997;Figueiras et al., 2002). Recently, using direct observations of currents, the mean duration of these upwelling-downwelling episodes has been estimated in ∼ 3 days (Gilcoto et al., 2017). The vertical homogenization of the thermohaline conditions, due to the advection of salty oceanic waters inside the Rías, during fall upwelling events erodes the frequent haline stratification driven by increasing precipitation and runoff (Rosón et al., 2008), and enhances mixing.
The seasonality of turbulent mixing in the upper ocean, both in coastal and open-ocean regions, remains poorly assessed due to the limited number of microstructure observations spanning entire annual cycles, such that the limited present knowledge derived mostly from indirect methods (Whalen et al., 2012;Warner et al., 2016;Inoue et al., 2017;Evans et al., 2018;Thakur et al., 2019;Cherian et al., 2020). Previous studies carried out in this region provided information about the variability of turbulent diffusivity on short-time scales during specific seasons (Barton et al., 2016;Villamaña et al., 2017;Fernández-Castro et al., 2018;Broullón et al., 2020) or along seasonal scales with coarser temporal resolution (Cermeño et al., 2016;Moreira-Coello et al., 2017).
Our results show that surface growth rates, which ranged from 0.11 to 2.42 d −1 (0.008-0.226 h −1 ), were in general higher during spring-summer transitional and upwelling periods, and minimum during fall-winter transitional and downwelling. These values are consistent with previous estimates determined for this system, by using a similar methodology, but with coarser temporal resolution. An early study carried out in the Ría de Arousa during an upwelling event (Hanson et al., 1986), reported depth-averaged phytoplankton growth rates ranging from < 0.1 to 1.89 d −1 (from < 0.004 to 0.0788 h −1 ). Cermeño et al. (2016) estimated slightly lower phytoplankton growth rates at the same station in the Ría de Vigo (∼0.004-0.114 h −1 ). Our data, derived from short (2 h) incubations experiments, are also consistent with apparent phytoplankton growth rates estimated from the dilution technique which ranged from 0.31 to 1.75 d −1 (Teixeira and Figueiras, 2009). However, maybe due to aliasing of short-term variability in their coarser observations, previous studies in the region did not report a clear seasonal cycle in phytoplankton growth rates. Seasonal trends in phytoplankton growth, consistent with our findings, have been described in several studies carried out in contrasting regions, mainly by using the dilution technique Gutiérrez-Rodríguez et al., 2011;Lawrence and Menden-Deuer, 2012;Anderson and Harvey, 2019).
The temporal coverage and resolution of our study allowed interpreting our results within the conceptual framework of the Critical Turbulence Hypothesis (CTH, Huisman et al., 1999a,b, Figure 7A). According to the CTH, positive net phytoplankton growth will occur when (1) the water column depth is shallower than the critical depth, which is consistent with the Critical Depth Hypothesis (CDH, Sverdrup, 1953), or (2) when turbulent diffusivity in a mixed water column is lower than a critical value, which causes the timescale for phytoplankton growth to be less than the cell's residence time in the portion of the water column where net growth can be achieved. In this case, the depth of the turbulent layer would be shallower than the mixed layer depth. Huisman et al. (1999a) described these two distinct and independent mechanisms by using simulations with variable turbulent diffusivity and water column depth. Similar to the CDH, they assumed that phytoplankton growth was only limited by light, and also constant values in time and depth for turbulence diffusivity and the phytoplankton loss term. In order to interpret our results under this theoretical framework, we performed some modifications on the original diagram ( Figure 7A). First, we used the mixed-layer depth instead of the water column depth, to indicate the depth of the water column where density was relatively vertically homogeneous. Second, we used the turbulent layer depth, instead of turbulent diffusivity, to consider both the magnitude and the vertical structure of turbulent mixing. Finally, we used surface growth rate instead of depth-integrated phytoplankton growth rates to represent phytoplankton growth in the turbulent layer ( Figure 7B). Consequently, phytoplankton growth was always positive and its variability must be interpreted in relative terms instead of positive or negative growth, as in the original diagram.
In agreement with the first scenario proposed by the CTH, our results showed that higher phytoplankton growth rates were found for periods of shallow turbulent and mixed-layers during spring-summer transitional and upwelling periods, whereas low values were observed during periods of deep turbulent and mixed-layers (fall-winter transitional and downwelling periods). Unfortunately, our dataset was not suitable to assess the second scenario described by Huisman et al. (1999a), as most observations corresponded to shallow mixed-and turbulent layers. Only two samplings (31 and 37) showed the mixed-layer to be much deeper (∼ 20 m) than the turbulent layer, suggesting that the mixed-layer was not actively mixing in <1 day. The surface phytoplankton growth rate for sampling 31, 2 days after the passage of Hurricane Ophelia, was relatively low (0.73 d −1 ), whereas this rate was slightly higher during sampling 37 (1.0 d −1 ).
Due to the difficulty to quantify turbulence in the field, previous studies using the CTH theoretical framework mainly used hydrographic and atmospheric data to infer mixing conditions in the upper layer, and its relationship with phytoplankton growth in different regions as the North Atlantic (Chiswell, 2011;Chiswell et al., 2013;Mignot et al., 2018), the North West European Shelf (Wihsgott et al., 2019), the southwest Pacific (Chiswell, 2011;Chiswell et al., 2013), the Western Mediterranean Sea (Bernardello et al., 2012;Kessouri et al., 2018), the Chukchi Sea (Lowry et al., 2018), and the Labrador Sea (Marchese et al., 2019). The decrease in turbulence that could explain the TLD being shallower than the MLD (second scenario proposed for the CTH) has been associated with a reduction in air-sea fluxes that generally precedes mixed-layer shoaling at the end of winter (Taylor and Ferrari, 2011) and lowered wind stress (Chiswell, 2011;Chiswell et al., 2013). Moreover, Brody and Lozier (2014) proposed that the mixing-length scale of the largest energy-containing eddies in the upper ocean is a better predictor for bloom initiation than the decrease in mixed-layer depth, the onset of positive heat fluxes, or the decrease in wind strength. They also found that the shift from buoyancy-driven to wind-driven mixing in late winter creates the decrease in mixing length scale, and thus the conditions necessary for blooms to begin (Brody and Lozier, 2015).
Only a limited number of studies have previously used microstructure observations to investigate the relationship between turbulence and phytoplankton activity and composition. Among those, most studies were focused on analyzing the role of mixing and nutrient supply as drivers of phytoplankton community structure (Sharples et al., 2007(Sharples et al., , 2009Machado et al., 2014;Villamaña et al., 2017Villamaña et al., , 2019. As far as we know only two studies have used microturbulence observations to investigate phytoplankton growth under the CTH framework. Huisman et al. (2004) extended the CTH including competition theory, to predict how changes in turbulent mixing affect competition for light between buoyant and sinking phytoplankton species. Consistent with the model prediction, results from a lake experiment showed that changes in turbulent mixing caused a dramatic shift in phytoplankton species composition. In coherence with the CTH, Hopkins et al. (2021) used a relatively short (2 weeks) period of high-frequency (sub-hourly) observations collected from gliders in the Northwest European Shelf to conclude that turbulent mixing, that mediate light availability, was the key process governing phytoplankton growth in spring.
Our study covers a full annual cycle which allows investigating the coupling between turbulent mixing and phytoplankton growth at different temporal scales. Our results show that shallower mixed-and turbulent layers conditions were not enough to ensure phytoplankton growth, as low phytoplankton growth was also found under these circumstances. It is important to note that, instead of following a single phytoplankton bloom over time, our weekly approach involves sampling different blooms at different stages of development. Previous studies in this region have described that the spring-summer and fall phytoplankton growth seasons are the result of a succession of blooms triggered by the short-term variability in the upwelling regime and runoff water pulses (Nogueira and Figueiras, 2005). This is consistent with observations carried out in the North Atlantic where fluctuations in the winter mixed-layer depth triggered occasional short-live blooms (Mignot et al., 2018). Using a wavelet transform analysis we showed that the coupling between turbulent mixing and phytoplankton growth in this system occurs not only over the seasonal cycle, but also at shorter temporal scales. These results suggest that mixing could act as a control factor on phytoplankton growth over the seasonal cycle, but also being involved in the formation of occasional short-lived phytoplankton blooms.

CONCLUSIONS
By analyzing a unique dataset combining weekly observations of microstructure turbulence and phytoplankton growth collected in the NW Iberian coastal upwelling we show that coupling between turbulent mixing and phytoplankton growth occurs over the seasonal cycle, but also at shorter temporal scales. In agreement with the Critical Turbulence Hypothesis (Huisman et al., 1999a,b), our results indicate that higher values of phytoplankton growth occurred during spring-summer transitional and upwelling periods, when the depth of the turbulent layer was shallower, whereas the opposite pattern was observed during fall-winter transitional and downwelling periods. However, shallower turbulent layers were not enough to stimulate phytoplankton growth, as low growth rates were observed during the same conditions. Despite our unprecedented sampling effort, we are aware that a significant part of the variability in this system was still not captured, which deterred us from unveiling the mechanisms driving the coupling between mixing and phytoplankton growth rates, specially at short time-scales. Furthermore, longer observations will be needed to discern whether the seasonal patterns described here are characteristic of the study system or determined by the specific investigated sampling period.
The Critical Turbulence Hypothesis assumes that phytoplankton growth is exclusively limited by light and that the phytoplankton loss rate is independent of depth and constant through time. However, mixing conditions in the water column influence light but also nutrient availability for phytoplankton cells. Our results indicate that, in general, higher phytoplankton growth occurs when light but also nitrate availability was higher. Although our study was performed in a nutrient-rich coastal upwelling ecosystem, previous studies in the region have shown that phytoplankton is limited by nitrogen during the summer, or at least it positively responds to its addition (Martínez-García et al., 2010). Moreover, although phytoplankton blooms have traditionally been attributed to changes in "bottom-up" environmental factors controlling phytoplankton division rates, such as light and nutrients (e.g., Hunter-Cevera et al., 2016;Wihsgott et al., 2019), changes in phytoplankton biomass are the result of the interplay between phytoplankton division rates and the sum of all loss terms.
Recent studies using satellite, in situ observations and modeling approaches highlighted the relevance of variability in loss rates (Behrenfeld, 2010;Boss and Behrenfeld, 2010;Behrenfeld et al., 2013). The Disturbance Recovery Hypothesis proposed by Behrenfeld and Boss (2014) argues that, contrary to the CDH, bloom initiation takes place in the winter because deep mixing dilutes phytoplankton cell density, thus reducing the encounter rate between predator and prey. It is important to note that our phytoplankton growth estimates, as derived from short carbon uptake experiments, are close to gross rates, and therefore do not reflect the balance between cell division and loss rates. Moreover, the time step we used for the wavelet analysis (8 days) is much longer than the timescale for phytoplankton growth (∼1 day). Future studies combing higher temporal resolution and longer observations, together with modeling approaches will be required in order to discern the mechanisms responsible for the coupling between mixing and phytoplankton growth observed in this system, and to evaluate the role of "bottom-up" and "top-down" control factors.
Despite their small surface extension (1% of the global ocean) eastern boundary upwelling systems account for 5% to the global marine primary production (Carr, 2001) and 20% of global fish catch . One unresolved question is these regions is the regulation of primary production. Substantial variability in the ratio of nitrate supply to primary production has been described , indicating that regulators other than nitrate supply must be at play. In agreement with previous studies (Hardman-Mountford et al., 2003;Rossi et al., 2009;Fearon et al., 2020), our results indicate that turbulent mixing may control phytoplankton growth at the different time-scales involved in the physical-biological coupling in these systems.
Anthropogenic global change perturbations have the potential to affect phytoplankton dynamics, including the timing and magnitude of blooms (e.g., Hunter-Cevera et al., 2016;Jena and Narayana Pillai, 2020). Therefore, identifying the mechanisms responsible for the formation of phytoplankton blooms, including turbulent mixing (Károly et al., 2020), is key to understand the role of the ocean in the present and future climate.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
AC, BM-C, and BF-C designed the overall study. PC, EF, AF-L, MG, MP-L, and BM-C carried out the observations and performed experiments in the lab. AC analyzed the dataset. AC, BM-C, and BF-C prepared the manuscript with contributions from all the co-authors. All authors contributed to the article and approved the submitted version.

FUNDING
This research was supported by grant CTM2016-75451-C2-1-R from the Spanish Ministry of Economy and Competitiveness to BM-C. AC was supported by a predoctoral fellowship for the doctoral formation (BES-2017-080935) from the Spanish Ministry of Economy and Competitiveness.