Net community oxygen production derived from Seaglider deployments at the Porcupine Abyssal Plain site (PAP; northeast Atlantic) in 2012-13

As part of the OSMOSIS project, a fleet of gliders surveyed the Porcupine Abyssal Plain site (Northeast Atlantic) from September 2012 to September 2013. Salinity, temperature, dissolved oxygen concentration and chlorophyll fluorescence were measured in the top 1000 m of the water column. Net community production ( N ) over an annual cycle using an oxygen-budget approach was compared to variations of several parameters (wind speed, mixing layer depth relative to euphotic depth, temperature, density, net heat flux) showing that the main theories (Critical Depth Hypothesis, Critical Turbulence Hypothesis, Heat-flux Hypothesis) can explain the switch between net heterotrophy to net autotrophy in different times of the year, The dynamics leading to an increase in productivity were related to shifts in regimes, such as the possible differences in nutrient concentration. The oxygen concentration profiles used for this study constitute a unique dataset spanning the entire productive season resulting in a data series longer than in previous studies. Net autotrophy was found at the site with a net production of (6.4±1.9) mol m -2 in oxygen equivalents (or (4.3±1.3) mol m -2 in carbon equivalents). The period exhibiting a deep chlorophyll maximum between 10 m and 40 m of depth contributed (1.5±0.5) mol m -2 in oxygen equivalent to the total N . These results are greater than most previously published estimates.


Introduction
Marine net biological production (N) is the balance between oxygen (O 2 ) production by phytoplankton during photosynthesis and O 2 consumption during respiration by the entire marine community. The seas around the world harbour almost half of the global plant production (Field et al, 1998;Williams, 1998), moving carbon and oxygen within and across compartments and reservoirs. By causing supersaturation or undersaturation of surface waters, biota is able to drive fluxes between the ocean and the atmosphere. This makes the ocean a carbon sink or a source depending on biological activity, which is important for estimating the global carbon budget and understanding how CO 2 influences climate as a greenhouse gas (Falkowksi, 1998).
Measurement of N over the entire annual cycle are important to understand the metabolic balance of the open ocean (i.e., the sign and magnitude of N), which is the focus of a long-running debate (del Giorgio et al., 1997;Duarte and Agustí,1998;Williams 1998;Duarte et al., 1999;Williams and Bowers, 1999;del Giorgio and Duarte, 2002;Karl et al., 2003;Hansell et al., 2009;Ducklow and Doney, 2013;Williams et al., 2013;Duarte et al, 2013). Uncertainty about N derives from the use of different methods for the calculation of N and its components. Several biases are known to affect in vitro measurements and their comparability with the real ocean (Williams et al., 1998;Kaiser et al., 2005). There are also challenges in separating the influence of biological and physical processes on in situ measurements (Hamme and Emerson, 2006;Emerson et al., 2008).
It is not even clear what some of the methods used are actually measuring (Regaudie-de-Gioux et al., 2014).
The analysis of N variations within the annual cycle, and their comparison with variations in other parameters., are instead useful to understand what factors are limiting or stimulating production and to investigate the validity of different mechanisms proposed so far. In regimes of nutrient limitation, gradual deepening of the mixed layer into nutrient-richer waters has been recognised as a plausible explanation for autumn blooms (Marra et al., 1990;Findlay et al., 2006). More recently, productivity peaks have been related to pulses of nutrients created by the interaction between wind and surface currents (Rumyantseva et al., 2015). The discussion about what triggers autotrophy when nutrients are not limited is more complicated. The Sverdrup Hypothesis (Critical Depth Hypothesis or CDH, Sverdrup, 1953) sees light as a driving factor: the plankton community is productive when the mixed layer is shallower than the critical depth, which is the depth above which total production exceeds total respiration. Since 1953, a long discussion has flourished to confirm or refute the Sverdrup CDH and new hypotheses have been proposed based on its weak points such as the assumption of phytoplankton behaving as a passive tracer. New hypotheses focus on the influence that turbulence has on the ability of the phytoplankton to access light. According to the Critical Turbulence Hypothesis (CTH, Huismann, 1999), high turbulence displaces the plankton at a faster rate than its growth rate. When turbulence decreases below a critical value, plankton grows faster than it is displaced and this leads to blooms (i.e., accumulation of oxygen and chlorophyll at the surface). Taylor and Ferrari (2011) linked the turbulence to the net heat flux, suggesting that the inversion from negative heat flux (water cooling) to positive heat flux (water warming) and the consequent shut down of convective mixing is a reliable parameter to predict the start of the bloom on an interannual timescale (Heat Flux Hypothesis or HFH). Enriquez and Taylor (2015) proposed another model linking the variations in turbulence induced by wind stress and water cooling (negative net heat flux leading to convective mixing) to the depth of the mixing layer. They predicted that when the mixing layer shoals, the phytoplankton respond with an increased growth rate (and then increased production). Behrenfeld (2010) suggested the Recoupling-Dilution Hypothesis, according to which phytoplankton has a positive growth rate when the mixed layer is deepening because of lower predation pressure.
The first goal of the present study is to estimate the magnitude of N in the productive layer of the water column through the analysis of variations of the oxygen inventory over time, based on oxygen concentration measured in situ by underwater gliders. Thanks to the high frequency of glider measurements and considering the length of this time series, the present study tries to overcome limitations in calculating N due to low spatial or temporal resolution. The present study surveyed an area located in the North Atlantic, in the proximity of the frequently sampled Porcupine Abyssal Plain (PAP) Sustained Observatory. This allows us to compare N estimates with results of previous studies that focused on the same area (e.g., Körtzinger et al., 2008a;Frigstad et al., 2015) along with basin-wide estimates. Furthermore, the availability of a suite of different parameters provided insights into the mechanisms that trigger increases in production. The second aim of the paper is to compare variations in N with other parameters to understand how observation fit with the different theories suggested to explain the increase in productivity.

Data acquisition
The OSMOSIS project included five cruises and three glider missions performed at PAP ( Figure 1) between September 2012 and September 2013. The cruises enabled the deployment and recovery of the gliders and the collection of a suite of in situ parameters to be used for calibration of the glider data. The three glider missions were carried out with short overlapping periods when the research vessels visited the area to swap gliders.
Details of the glider campaigns, quality control, calibration and analysis of the physical oceanographic context of the year-long time series are provided by Damerell et al. (2016). During each mission two gliders operated at the same time moving along two separate butterfly-(or hourglass-) shaped transects oriented perpendicular to each other (one glider moving north-south and the other east-west) centred around 48.7º N and 16.2º W with 15 km long edges ( Figure 1). Because of a malfunctioning oxygen sensor and deviations from the designated transect, only one glider per mission was used to create the year-long dataseries. Data from glider SG566 were considered between September 2012 and January 2013; SG502 between January 2013 and April 2013 and SG566 a second time between April 2013 andSeptember 2013. The glider data are held at the British Oceanographic Data Centre and can be accessed at https://doi.org/10/cqc6.

CTD calibration
The calibration of the data followed several steps. Ship-CTD oxygen concentration, c(O 2 ), profiles were calibrated against water samples analysed by Winkler titration. Glider c(O 2 ) profiles were first adjusted to account for the response time of the optodes and subsequently calibrated against the Winkler-calibrated ship-CTD c(O 2 ) profiles.
CTD casts were performed just before the deployment of the gliders and soon after their recovery. Oxygen concentration was measured by a Clark-type electrode (Seabird SBE43) attached to the rosette frame and water was sampled by the means of Niskin bottles attached to the same rosette. At the end of the cast, Winkler samples were collected from selected Niskin bottles and their c(O 2 ) was measured by Winkler titration following WOCE protocols (Culberson 1991;Dickson,1996). For each cast, CTD c(O 2 ) was calibrated by linear regression against Winkler-derived c(O 2 ).

Response time correction
The gliders recorded c(O 2 ) by means of optodes (Aanderaa Data Instruments; Tengberg et al., 2003), which measure the lifetime of the red light emitted by excited porphyrins in a sensing foil as a temperature-compensated phase difference (ϕ TC ), which depends on the c(O 2 ). The diffusion time of the gas through the optode foil is characterised by a certain response time (τ). When plotting ϕ TC against pressure, the response time shifts ϕ TC profiles in the same direction the glider is moving in. Gradients will appear deeper than their actual depth when gliders descend and shallower when gliders ascend through the water column. A best-fit response time (τ) was determined from each pair of consecutive ascents and descents over the top 300 m. All response times of ascent-descent pairs in the same glider mission were fitted to a normal distribution and its central value was used as τ to shift all the profiles of the same glider mission. The best-fit response times were applied to shift the profiles backwards in time, which resulted in a shift of the descent profiles upwards and the ascent profiles downwards. Variations and uncertainty in τ may be caused by variations in temperature and short-term changes in the vertical c(O 2 ) profile during subsequent glider dives.

Glider calibration and despiking
The τ-corrected ϕ TC profiles were calibrated using the Winkler-calibrated ship-CTD c(O 2 ) profiles, using the closest CTD cast (less than 4 km and less than 3 h away). ϕ TC profiles from several glider dives and several CTD profiles were used in the calibration as long as the glider and the ship were within the limits of proximity defined above. Calibrated phase differences ϕ cal (ship-CTD) were calculated from ship-CTD c(O 2 ) as follows: 1) calculate water vapour pressure p vap using potential temperature and practical salinity; 2) calculate oxygen saturation concentration c sat (O 2 ) using the parameterisation of Garcia and Gordon (1992) with the solubility coefficients of Benson and Krause (1984); 3) calculate air saturation s(O 2 ) equivalent to an atmospheric pressure of 1013.25 hPa where s(O 2 ) = c(O 2 ) / c sat (O 2 ); 4) calculate partial pressure of oxygen, Δp(O 2 ), using s(O 2 ) and p vap 5) calculate ϕ cal by inverting the manufacturer-provided sensing-foil specific polynomial that uses temperature, Δp(O 2 ) and 21 coefficients (A 0 to A 13 and B 0 to B 6 ) to derive ϕ cal ϕ cal (ship-CTD) profiles and glider ϕ TC profiles were then binned according to potential density and compared.
The offset and slope of their linear regression was used as C 0 and C 1 in Eq. 1: Since CTD casts were performed during the cruises at the deployment and recovery of the gliders, there was a linear calibration equation obtained at the beginning and at the end of each mission. In case the two calibration equations were not the same (indicating drift of the sensor over time), a time-varying C 0 and C 1 was calculated for each dive interpolating over time between the values at the beginning and end of each mission. Each dive had therefore a unique calibration equation that transformed ϕ TC into ϕ cal .
ϕ cal (glider) profiles were then transformed into calibrated glider c(O 2 ) with the same five steps used for the back-calculation of the CTD ϕ cal in reverse order.
Spikes in the profiles were automatically flagged when a data point matched any of the following criteria: 1) unrealistic c(O 2 ), i.e values <0 µmol kg -1 or >1000 µmol kg -1 ; 2) significant increase in the standard deviation of a c(O 2 ) profile due to a single point; 3) single points with anomalous c(O 2 ) within water masses with constant concentrations; 4) c(O 2 ) values apparently above the surface (due to pressure sensor inaccuracies).
Visual despiking of profiles was carried out to remove anomalous c(O 2 ) at the surface due to light hitting the foil, waves that expose the sensor to air or problems in τ correction.
At the end of the process, 527 points were flagged as spikes in SG566 (0.14 % of the total), 837 in SG502 At the end of the calibration, entire profiles were also flagged as anomalous and not considered in further calculations. In particular, the profiles recorded after the 11 th August 2013 were disregarded because biofouling affected the optode of the SG566. Biofouling was identified by anomalous readings in c(O 2 ) throughout the whole water column and its presence was confirmed by visual inspection when the glider was recovered. More details about biofouling during this survey are provided in Appendix A -Biofouling, including the explanation of why the 11 th August was chosen as cut-off date.

Mixed layer calculation
The calculation of z mix was performed for each of the 4035 c(O 2 ) profiles of the OSMOSIS time series using three consecutive glider missions (one glider per mission). c(O 2 ) at 5 m depth (calculated by interpolation) was

Production calculation
Marine biological production (N) at the top of the water column was calculated analysing the changes in the oxygen inventory per unit area (I). Only the c(O 2 ) profiles from glider descents were used for productivity calculations because some of the ascents were affected by spikes caused by sunlight hitting the optode foil near the surface (see Appendix A). The calculations focused on the top 60 m (z lim ) of the water column, which was equivalent to the mean euphotic depth (z eup ) of (60±15) m during the study. z eup was defined as the depth at which PAR falls to 1% of the level measured at the surface. PAR was measured in situ by sensor installed on the seagliders.
N was calculated as where F as is the air-sea O 2 flux (positive for O 2 outgassing), E is entrainment and N is net community production. ΔI was calculated as the c(O 2 ) inventory change above z lim between consecutive profiles. Profile inventory was calculated by integrating c(O 2 ) over depth.
F as was calculated using the correction for bubble injection (Δ) formulated by Woolf and Thorpe (1991): where k(O 2 ) is gas transfer velocity for oxygen; c(O 2 ) is the dissolved oxygen concentration in the mixed layer; and c sat (O 2 ) is the oxygen saturation concentration calculated according to the Benson and Krause (1984) fit of Garcia and Gordon (1992) using the atmospheric pressure derived by interpolation of ERA-Interim reanalysis data (http://www.ecmwf.int/en/research/climate-reanalysis/era-interim, resolution of 6 hours and 0.125º in latitude and longitude). c(O 2 ) and c sat (O 2 ) used to calculate F as were derived from the mean value for the top 10 m or as the mean above z mix when z mix < 10 m. The gas transfer velocity at a Schmidt number Sc = 600 was parameterised following Nightingale et al. (2000) using the daily averaged wind speed at 10 m derived from ERA-Interim reanalysis with the same resolution as the atmospheric pressure. F as used to calculate N between two profiles was the mean of the F as measured for these two profiles. It is worth noting that Song et al. (2015) compared ERA-Interim data with measurements from eight buoys, showing a good agreement between reanalysis and in situ measurements (regression coefficients for wind speed was above 0.7 and for direction was greater than 0.79). However, ERA-Interim overestimated wind data at the buoy stations, with the max difference of 1.8 m/s; and 13% ERA-Interim wind data below 6 m/s were flagged as not good. There is therefore the possibility that F as used here might be slightly overestimated.
Entrainment was considered as the change of c(O 2 ) due to deep water mixed all the way to the surface in the mixed layer when z mix deepens. Entrainment could be positive or negative, corresponding to an increase or a decrease of the oxygen inventory. E between any two profiles at times t 1 and t 2 was considered only when z mix deepened and when z mix at t 2 was deeper than z lim . Otherwise, E was equal to zero, e.g. when z mix deepened, but remained above z lim , a redistribution of the oxygen was assumed without any O 2 flux occurring through z lim and therefore no variation in the oxygen inventory above z lim (I(z lim )). In order to calculate z mix deepening and shoaling, z mix values were smoothed using a moving filter with 5 points span.
This was an Eulerian rather than Lagrangian study and therefore part of ΔI between adjacent profiles was the signal of geographical heterogeneity (patchiness) and horizontal advection. Advection has been considered negligible in previous studies Nicholson et al., 2015) due to the rapid effect of air-sea O 2 flux inequilibrating the concentration at the surface. However, Alkire et al. (2014) and Hull et al. (2016) showed that advection can significantly affect N estimates over time scales of days/months and for spatial scales less than 50 km. In this study, in order to consider the effect of advection, the individual N values measured between consecutive profiles were averaged over 7 days. Looking at the time-series, 7 days was in fact a period longer than the glider took to traverse the distance between large N values (sign of glider entering a O 2 -rich advected water mass) and large negative N values (gliders getting out of the advected water mass). Averaging over 7 days was considered therefore an effective way to cancel out the positive and negative contribution of any water mass advected in the area to N calculation. This is also in line with Alkire et al. (2014), which in a similar glider experiment showed that the time scale of advection processes was around 4 days. In the present study, running averages of N for overlapping 7 day-long bins were assumed to be a valid estimate of biological activity (sensitivity to averaging period is tested in Appendix B).
The period of 7 days is also the approximate time any glider took to complete its butterfly-or hourglass-shaped transects and the averaging therefore gave estimates of N for the entire surveyed area, disregarding its internal geographical heterogeneity.
The variation in c(O2) were also analysed in relation to the thermal exchange between atmosphere and the ocean. In order to so this, the timeseries of net heat flux (H) was obtained from ERA-Interim reanalysis. The resolution of 6 hours and 0.125º in latitude and longitude.

Annual time series
The distribution of oxygen measured by the gliders at PAP between September 2012 and September 2013 is plotted against time and depth in Figure 2, along with oxygen saturation (s(O 2 ), ratio of c(O 2 ) over c sat (O 2 )) and Apparent Oxygen Utilization (AOU, difference c sat (O 2 ) -c(O 2 )). The other physical parameters measured concurrently with c(O 2 ) are shown in Figure 3. The vertical distributions suggest the presence of three layers in the water column. These layers have also been described by Damerell et al. (2016) who analysed salinity and temperature measured concurrently with the oxygen concentrations analysed in the present study. The top layer was roughly 150 m deep. This layer included the ocean surface boundary layer and had a seasonal cycle in the temperature due to solar insolation. Salinity was more variable, it did not follow any seasonal cycle and varied at all time scales, probably due to horizontal advection, local air-sea interaction and vertical mixing.
The intermediate layer, between 150 m and 700 m, was characterised by a significant intra-seasonal variability in temperature and salinity, also strongly intercorrelated. This variability was mostly linked to gyre-scale and mesoscale dynamics rather than the surface forcing. Bottom layer was between 700 and 1000 m and had high variability at all timescales in temperature and salinity, strongly influenced by the Mediterranean Outflow Water (MOW) that appeared at these depths. This paper focuses on the top layer because the euphotic depth was always shallower than 100 m at any time and, therefore, the plankton blooms were restricted to in this layer.   In order to calculate N, the physical factors affecting c(O 2 ) variations (F as and E) were calculated. F as depends on the difference between c(O 2 ) and c sat (O 2 ) following Eq. 5 ( Figure 5a). There was supersaturation (c(O 2 ) > c sat (O 2 )) in September at the beginning of the timeseries and after May. There was a period of undersaturation (c(O 2 ) < c sat (O 2 )) lasting from November until March and then a period of quasi equilibrium from March to May. F as varied between -193 mmol m -2 d -1 and 155 mmol m -2 d -1 , with a mean value of (-13±53) mmol m -2 d -1 . F as showed a strong seasonality, with a short period of outgassing at the beginning of the timeseries, followed by a long period of ingassing from the end of September to the beginning of March (Figure 5c). This was in turn followed by two months of quasi equilibrium with weaker influx and, from the end of May, F as switched sign and started a period of outgassing that lasted until the end of the mission. Over nearly an annual cycle, this region of the North Atlantic is shown here to be a sink of oxygen rather than a source, with 4.8 mol m -2 of O 2 absorbed by the ocean during the surveyed period. This is driven by pulses of strong influx due to high wind that induce high bubble influx (Δ, Figure 5b), but also by the late-occurring supersaturation. The data after 11 th August 2013 were disregarded because of biofouling, making the time series one month shorter than an annual cycle (see Appendix A). This missing month was probably a productive period and therefore its inclusion would have likely increased the magnitude of the annual outgassing if taken into account. bubbles supersaturation parameterisation (Δ) according to Woolf and Thorpe (1991)

; (c) air-sea oxygen flux. Positive values indicate outgassing of oxygen in the atmosphere and negative values indicate influx of oxygen in the water column.
The other element in the calculation of N was entrainment, E. When z mix did not deepen below z lim , E was considered to be zero ( Figure 6). When z mix deepened, but remained above z lim , it was assumed that a redistribution of the I(z lim ) occurred without any O 2 flux occurring through z lim . Also, when z mix shoaled, the change in I(z lim ) was assumed not to be related to any mixing with deeper water masses below z lim and, therefore, no E was assumed to occur.
Fluctuation in z mix linked to geographical variability and to the sensitivity of the threshold used for z mix computation would affect E because deepening events are not compensated by the shoaling events in the calculation. In order to mitigate the effect of z mix variability on N calculation, z mix values were smoothed using a moving average filter over 5 datapoints (black line in Figure 6a). The values of N were calculated as the variation in c(O 2 ) (ΔI) between consecutive profiles not explainable by  Four periods (Table 1)

Autumn bloom
During the autumn period, net autotrophy alternated with net heterotrophy (Figure 9a). During the first switch from net respiration to net production on the 27 th September N reached a magnitude of (16±12) mmol m -2 d -1 in oxygen equivalent. This peak in productivity lasted until 2 nd October, and at the same time chlorophyll

Heterotrophic period
The period between 21 st November 2012 and 9th February 2013 was characterised by a long initial period of heterotrophy ( Figure 10a; N < 0 for 62 % of the time). Net consumption in the area was calculated at -0.3 mol m -2 . Despite the mean N in this period was (-3±34) mmol m -2 d -1 , the community seemed to go through a train of short peaks in production. These peaks ( Figure 11a) coincided with sharp changes in potential density (visible in Figure 11c), and this horizontal heterogeneity suggests that the glider was probably crossing a mesoscale feature, as discussed by Thompson et al. (2016). These are therefore likely false increases in N, since not related to biological activity. A valid correction for these values was not considered in this study and they were considered in the annual calculation of N. However, they represent an overestimation of biological N. Between 30 th December 2012 and 9 th January 2013 there was a gradual transition towards a shallower z mix, which eventually became shallower than z eup (Figure 10b). This coincided with high N exceeding 65 mmol m -2 d -1 . In the same way, a decrease in wind speed and the shoaling of z mix to the depth of z eup coincided with the peak between 28th January and 6th February.

Spring
The spring season is here defined as the period between 10 th February 2013 and 19 th June 2013, 130 days during which the area was autotrophic (N > 0) with cumulative O 2 production of 4.5 mol m -2 and a mean N of (34±44) mmol m-2 d -1 (Figure 11). In this period, however, there was high variability (N standard deviation = 44 mmol m -2 d -1 ) because production alternated with net respiration from the beginning of May (Figure 13).
Six 7-day bins of N estimates were above 100 mmol m -2 d -1 with a maximum of 149 mmol m -2 d -1 .
The period showed a series of N fluctuations. In February, despite minimal variations in c(Chl a), N is high.
High levels of U 10 are linked to periodic deepening of z mix , that is however usually above z eup , which indicate an increase in the amount of light experienced by the cells and, therefore, a boost in productivity. This also coincide with a period of H oscillating between negative and positive values, after a period of strong negative values. Another obvious feature is a short chlorophyll bloom happening in the upper 50 m of the water column between 24 th and 28 th February 2013, here called 'End-February Event' (EFE, Figure 12a). This event happened when z mix was very shallow (20-25 m), wind decreased and net heat flux (H, Figure 12e) became temporarily positive (heat from the atmosphere to the ocean). and z mix deepened. There was a slight increase of potential density and slight decrease of temperature showing that water from below the z eup , probably enriched in nutrients, was mixed to the surface. Wind then decreased and N increased for a brief time, followed by an increase in c(O 2 ) below z mix rather than above, which could be evidence of low nutrient concentrations at the surface.

Summer bloom and deep chlorophyll maximum
From 20 th June onwards, N was relatively high and above zero and z mix was always shallower than z eup ( Figure   14). This summer period as a whole had a mean N of (47± 36) mmol m -2 d -1 and produced 2.5 mol m -2 . The main features in this period are a surface bloom between 26 th June and 4 th July and the development of a deep chlorophyll maximum (DCM). The surface bloom was very productive with a mean N of (71±24) mmol m -2 d -1 , reaching 110 mmol m -2 d -1 . However, since it lasted only 8 days, it produced only 0.6 mol m -2 . There was also an increase in c(Chl a). The bloom ended when wind increased again and z mix deepened.
From 8 th July wind decreased (Figure 14c) and the water column transitioned to a regime of low turbulence and strong stratification. A DCM developed and the production increased significantly. Both the subsurface oxygen-and chlorophyll-rich feature were above the z eup of 60 m. During the time in which there was a DCM, the system remained productive until 4 th August, when the productivity decreased along with an increase in wind speed, potentially leading to increased turbulence in the water. z mix deepened within the DCM, eroding it and mixing it with surface waters. The decrease of N at the end of the DCM period occurred at the same time as a decrease in the c(Chl a). During the presence of DCM (30 days, 8 th July to 8 th August) 1.5 mol m -2 were produced with a mean N of (48±32) mmol m -2 d -1 .

Annual cycle of N
This study calculated the productivity of the plankton community at the Porcupine Abyssal Plain for an annual cycle based on variations in c(O 2 ). Data were acquired between surface and 1000 m depth, but the analysis focused on the production in the euphotic layer, which was always within the upper ~ 100 m of the water column. In winter, colder temperatures increased c sat (O 2 ) and triggered an influx of O 2 from the atmosphere that, considering the absence of substantial biological activity and the rapid gas exchange due to strong winds, was expected to equilibrate to saturation (Broecker and Peng, 1992;Woolf and Thorpe, 1991;Chester, 2000;Ito et al., 2004). However, the water stayed undersaturated during this period showing that the air-sea O 2 flux was not sufficient to saturate the water. These results confirm previous observations (e.g. Körtzinger et al., 2001;Russell and Dickson, 2003;Körtzinger et al., 2004;Keeling et al., 2010, Duteil et al., 2013 and model output (Ito et al., 2004), which also reported undersaturation in surface waters in several oceans. the community of the deep chlorophyll maximum seems able more efficient, being able to produce more O2 and supersaturate the water at lower c(Chl a) with respect to earlier periods of the year. Changes in the phytoplankton community, such as the succession of dominant species over time, are linked to variation in parameters in the water column; the decrease of silicates after the uptake during diatom blooms, is one of the phenomena that can drive these successions (Egge and Aksnes, 1992;Martin-ezequel et al, 2000) making diatoms bloom before autotrophic dinoflagellates (Margalef, 1978;Leterme et al., 2005;McQuatters-Gollop et al., 2007, Barton et al., 2013. Furthermore, dimensional classes within the same group show a succession related to the ability of smaller species to uptake nutrients more efficiently in oligotrophic environments (Barton et al., 2013). This might also be related to the cyclic changes in nitrate concentration shown in the area by Hartman et al. (2010). The changes in the length of day light has also been linked to changes in bacterioplankton composition, which in turn has been linked to changes in phytoplankton (Gilbert et al., 2012).
Considering that each taxon produces a different amount of oxygen per mole of chlorophyll a, a change in the dominant taxa when productivity is moved in the deep chlorophyll maximum could explain the higher amount of oxygen per unit chlorophyll. It must also be considered that the same species have variations in its amount of chlorophyll a, for example because of photoacclimation (Sakshaug et al. 1997;Goericke and Montoya 1998;Henriksen et al. 2002), which means that variations in the amount of chlorophyll a per cell do not change linearly with the amount of O 2 produced. This shift in species and/or in the physiology of the cells is influenced by many environmental factors such as light intensity, nutrient availability or the regime of turbulence (e.g., Huisman et al., 2004;Veldhuis et Kraay, 2004;Brunet et al., 2008;Dimier et al., 2009, Barton et al., 2013.
Considering that the formation of the summer deep chlorophyll maximum suggests a substantial attenuation of mixing forces and that the end of the main bloom can be related to low nutrients, it is reasonable to induce that environmental changes drove a shift in the phytoplankton community, which led to higher saturation per unit chlorophyll in the DCM. It shall be also considered that the reduced turbulence might have given time to phytoplankton to adapt to different levels of light and that photosensitivity might therefore have played a game in reducing productivity at the top of the water column.
The productive period (spring and summer together) spanned from 9 th February to the start of biofouling on 11 th August and had a time-integrated oxygen production of (7.1±2.1) mol m -2 with a mean N of (39±41) mmol m -2 d -1 . The seasonal production was converted to C equivalents using the photosynthetic quotient (PQ) of 1.5 as in Alkire et al., (2014) and resultant N C values are listed in Table 3.
Considering only the productive period, our study region produced N C = (4.8±1.4) mol m -2 . This value is lower than the (6.4±1.1) mol m -2 estimated by Körtzinger et al. (2008a)  . This result fits with previous works that show net production when incubation-free methods are used (Letscher and Moore, 2017) and support the position that this part of the ocean has positive net production. Despite the major focus of discussion about the trophic state of the ocean has been focused mostly on subtropical oligotrophic gyres (e.g., Williams et al., 2013;Duarte et al, 2013), our results should be taken in account in global budgeting to estimate the carbon cycle and to consider whether the ocean is a net sink or source of carbon.
Considering the whole time series, the PAP site was autotrophic between September 2012 and August 2013, with annually integrated net community O 2 production of (6.4 ±1.9) mol m -2 a -1 ((4.3±1.3) mol m -2 in C equivalents). This value was computed without the last month of the year, which was disregarded due to biofouling on the optode. However, the shape of the biofouled profiles showed a DCM above 60 m (data not shown). Biofouling and its progressive growth also show a productive phytoplankton community. The disregarded period can therefore be considered productive and the cumulative N of 6.4 mol m -2 is likely an underestimation of the real production in the area over the full year.
The annual production values are higher than previous annual N C estimates of Quay et al. (2012) who estimated 2.8 mol m -2 in the subpolar North Atlantic Ocean or by Neuer et al. (2007) who estimated Nc=3.3 mol m -2 as a mean between 1996 and 2000 in a more southerly area. The annual production estimated in the present study is instead similar to the 5.5 mol m -2 estimated by Ostle et al. (2015) for 2012 in their region 2 (where PAP site is located). This area was found in their study to be the most productive sector in the basin. This similarity, however, hide seasonal differences since estimates from Ostle et al. (2015) are lower during the productive period and higher during the winter.
The differences among studies are probably due to factors such as interannual variability in the area and differences in methods used for the calculations. N in the present study was an estimate of the production in the euphotic layer and, therefore, studies analysing variation at greater depths than z eup are expected to be lower because of the respiration occurring deeper. For example, some of the studies compared here (e.g., Frigstad et al., 2015;Ostle et al., 2015) analyse the changes above z mix rather than above z eup , while others (i.e. Körtzinger

Bloom initiation dynamics
Measuring N based on oxygen variations (direct by-product of photosynthesis) shows that heterotrophy and autotrophy alternate throughout the whole year. Different processes seem to trigger autotrophic peaks at different times of year. During autumn (nutrient limitation), N peaks have been related to pulses of nutrients created by the interaction between wind and surface currents (see Rumyantseva et al., 2015). However, the trigger for later N peaks seems to follow instead the gradual deepening of the mixed layer into nutrient-richer waters, a dynamic already suggested in previous papers (Marra et al., 1990;Findlay et al., 2006). This process also explains how peaks of N can develop at the end of the spring, when nutrient can be depleted as well after a big bloom.
When nutrient limitation could be excluded, N increased only when the mixed layer was shoaling while there was net heterotrophy during the winter, when the mixed layer was deepening. Our results therefore disagree with the Recoupling-Dilution Hypothesis. Instead, the present study presents evidence supporting the validity of the mechanism proposed by Enriquez and Taylor (2015). When nutrients are not limiting, the peaks in N are associated with decreasing wind speed and positive net heat flux, which in turn are linked to a shoaling mixing layer. Our results imply that the plankton community needs low turbulence conditions in order to bloom. The magnitude of the blooms also seems to be related to the relative depth of the mixing layer to the euphotic depth, rather than the critical depth used in Enriquez and Taylor (2015). The main blooms developed when z mix shoaled near or just above z eup . The increase of production when the mixing layer shoaled did not always correspond to significant increases of chlorophyll a concentration at the surface (Rumyanteva et al., 2019) but rather subsurface (e.g., the peak starting on 3 rd March).
Most of the peaks of N are associated with positive net surface heat flux i.e. ocean warming. In particular, the start of the main bloom during spring coincided with the switch between a period of mean negative net heat flux and a period of mean positive net heat flux at the beginning of April consistent with the HFH theory of Taylor and Ferrari (2011). This suggests that the time of this switch in the sign of net surface heat flux could be used as a proxy to analyse interannual variability in the starting time of the main bloom. However, N increased after a delay due to the presence of a storm, showing the need to take into consideration the turbulence induced by the wind stress in order to have more accurate bloom timing estimates, as hypothesised by Chiswell (2011) and Brody et al. (2013).
This study also highlights the presence of peaks in productivity when chlorophyll concentration showed no variations, which have to be considered along with the chlorophyll fluorescence-defined blooms in order to analyse correctly the triggering factors that increase production. It is also important to use high temporal resolution in situ data instead of climatologies to better appreciate the high variability of the system. The use of the mixing layer depth instead of the mixed layer depth is important to analyse variations in turbulence that affect the plankton and its metabolic activity.

Autumn period
The presence of increased productivity during the autumn is well known for this part of the ocean and is usually referred as the "autumn bloom" (Colebrook, 1982;Longhurst, 1995;Dandonneau et al., 2004;Lévy et al., 2005;Neuer et al., 2007;Martinez et al., 2011). The term 'bloom' however suggests a prolonged period of stable productivity, which was not observed in this data series. In the present study in fact, a series of autotrophic peaks happened in this season at the end of storms. Production enhancement after storms has already been seen in previous studies (Babin et al., 2004;Son et al., 2006;Wu et al., 2008, Rumyantseva et al., 2015. This supports the notion that autumn blooms are sustained by nutrient pulses through the pycnocline due to shear spiking (Rippeth et al., 2005, Rippeth et al., 2009Williams et al., 2013, Rumyantseva et al., 2015 generated by rapid change in wind stress (Pollard, 1980). This suggests that pulses of nutrients from below stimulate biological production in shallow and nutrient depleted mixed layers. In these post-storm blooms, wind has to decrease before N could peak.
In contrast, the last peak of this season (30 th October to 6 th November) was linked to the gradual deepening of z mix and the introduction to the surface layer of nutrients from below. The decrease of c(Chl a) at the end of this productivity peak (figure 3d) marks the passage to a less productive regime, with N not increasing even when z mix shoaled. This peak seemed therefore to follow the dynamics described by Marra et al. (1990) and Findlay et al. (2006) according to which the nutrient input fuelling the autumn bloom is caused by the gradual deepening of z mix . These two different dynamics of N peaks were discussed by Dutkiewicz et al. (2001), who showed that increasing wind speed can enhance N by bringing nutrients towards the surface as well as decrease N moving phytoplankton cells deeper, where they consume more than they produce (for example during storms).
The peaks of productivity during autumn were of lower magnitude than in spring and the total N was not significantly different from 0. This is in line with the conclusions of Martinez et al. (2011) who showed an asymmetry in the magnitude of the blooms in different seasons. According to Martinez et al. (2011) , there was a shift from the 1980s, when autumn blooms had a magnitude comparable with the spring blooms, to the present day, when autumn blooms are smaller than spring blooms, as we find here. Martinez et al. (2011) linked this change to the delayed deepening of z mix at the end of the summer that now happens later in the year than in the past.
Lateral advection, presence of mesoscale events, change in zooplankton community or even the effect of wind and storms are other possible causes for the smaller magnitude of autumn blooms proposed by Martinez et al. (2011). The present study supports the conclusions of Martinez et al. (2011) of non-symmetric blooms between seasons, and uses in situ measurements to support their hypothesis, which was based on satellite data.

Heterotrophic period
Heterotrophic periods have been already recorded in the North Atlantic (see literature in Duarte et al., 2013), however their magnitude and impact on the annual metabolic balance are debated (Duarte et al., 2013).
Multiannual In this study, pulses of positive N during the heterotrophic period were linked to the glider crossing a mesoscale feature. The averaging process was probably not able to fully eliminate the signal of this geographical heterogeneity in N because the feature stayed in the area longer than one week.
The feature crossed by the glider at the end of November -beginning of December 2013 ( Figure 11) had higher c(O 2 ) and part of this might be actually due to production. However, the density of the water was lower and an increase in c(O 2 ) was explainable by the solubility effect (higher c sat (O 2 )). This peak was therefore probably overestimating N.
Other peaks occurred when z mix stopped deepening and shoaled above z eup . The potential reduction in turbulence in these cases seems to be linked to higher productivity since the N peaks were interrupted when the wind speed increased again.
The consumption estimated in the heterotrophic period (0.3 mol m -2 ) was one order of magnitude lower than the production estimates in the rest of the year. The present study therefore shows that the presence of potentially protracted periods of net heterotrophy in this part of the North Atlantic have only a moderate impact on the production on an annual scale.

Spring
The PAP site is located in the North Atlantic between the subpolar and subtropical gyres, where, according to Longhurst (1998), blooms are expected in May. The timing of this bloom and its intensity have high interannual and geographical variability (Ueyama and Monger, 2005;Henson et al., 2006;Henson et al., 2009;Kahru et al., 2011;Zhai et al., 2013; and this explains why, despite being one of the most studied systems in oceanography, the dynamics of the North Atlantic spring bloom have not been fully understood as yet. The pattern of several parameters (wind speed, z mix in relation to z eup , temperature, density, net heat flux) was compared with variations of N, with several peaks in productivity observed before the main bloom between April and May. Following this comparison, possible explanations for the variations of production over time were suggested: the water near-surface was considered nutrient repleted at the beginning of spring, while nutrient-limitation was assumed to happen later on in the season,considering the seasonal pattern showed by Hartman et al. (2010) in this area. Variations in nutrients availability to phytoplankton could therefore be the cause of the oscillations between N > 0 and N < 0 in the second part of the spring, with phytoplankton becoming more productive when nutrients were supplied. However, the absence of direct measurements of nutrient concentrations in this study makes it difficult to confirm these speculations and extrapolate them to infer more general dynamics. At the end of this period there were rapid transitions between accumulation of oxygen at the surface and below z mix . These were probably related to geographical patchiness and demonstrate the heterogeneity of biological production at this time of the year.

Summer and deep chlorophyll maximum
Changes in nutrient concentrations may have caused the variations of N seen during the summer. Particularly interesting in this period is the DCM that lasted for over 30 days in the area thanks to a well-stratified water column with a very shallow z mix above 10 m. The presence of this feature suggests nutrient limitation in the upper water column, as shown in previous studies (Klausmeier and Litchman, 2001;Klausmeier et al., 2007, Denaro et al., 2013. When the DCM was present, N integrated stayed high, accounting for 38 % of the cumulative N estimated throughout the whole study. The formation of the DCM is usually related to increases in biomass (Beckmann and Hense, 2007) and/or to adaptation in the chlorophyll content of the cells (Fennel and Boss, 2003). This feature is a challenge for N calculations based on remote measurements or on the sampling of the plankton community for in vitro incubation. The ocean colour measured by satellite-borne sensors can be biased if the DCM is shallower than ~ 45 m depth study (Stramska and Stramski, 2005), as found in the present, de facto decoupling fluorescence readings from the real value at the surface. Annual N estimates obtained with the method used here should therefore be of higher accuracy and reliability than the ones based on remotely sensed ocean colour.
The demise of the DCM is probably related to nutrient limitation. N decreased when z mix started to deepen at the end of July; however, z mix was still above z eup and so the reduced productivity was not related to the limitation of light. Instead, wind speed increased and the ensuing vertical turbulence may have exposed the plankton to the nutrient depleted water above, lowering the production. Evidence of this is the decrease of c(Chl a) happening at the same time between 20 and 40 m. An alternative explanation could be the reduction of the photosynthetic performances in the cells due to changes in photosensitivity. The low turbulence could in fact narrow the difference between phytoplankton adaptation time and water mixing time, resulting in changes in the pigment physiology of the cells (Claustre et al., 1994).
From the end of June, F as was coupled to N values. The entrainment in this period was negligible, thanks to the strong stratification that allowed the formation of the DCM. This F as can then be considered biologically induced, as found by Kaiser et al. (2005) for systems with negligible vertical and horizontal mixing.

Conclusions
Net community production (N) above the mean euphotic depth near the PAP site from September 2012 to August 2013 has been calculated by analysing the variations in depth integrated oxygen concentration over time. The area is autotrophic, with a mean N value of 19 mmol m -2 d -1 and a total production of 6.5 mol m -2 and an estimated annual production of 7 mol m -2 . The analysis of the annual cycle of net community production shows the presence of four periods with different regimes: the autumn period, a heterotrophic period and two productive periods (spring and summer) separated by the depletion of nutrients after the spring bloom. During the summer a very productive deep chlorophyll maximum developed which was responsible for a significant portion of the annual production. The values calculated fit the range of published estimates of net community production in the North Atlantic basin and in the same area. The variations within this range are attributed in part to the differences among the methods used for the calculations and also to interannual variability.
Variations in production are associated with factors such as wind speed, net heat flux and mixing layer depth.
The theories proposed in the last decades for the explanation of the blooms (Critical Depth Hypothesis, Critical Turbulence Hypothesis, Heat-flux Hypothesis) are consistent with each other in explaining different mechanisms for how the system passes from net heterotrophy to net autotrophy when favourable conditions are matched.

Acknowledgments
This work would not have been possible without the piloting skills of several colleagues at the University of East Anglia, National Oceanography Centre (Southampton) and California Institute of Technology, plus the contribution of many people not represented in the author list. The deployment and recovery operations of gliders has also been possible thanks to the scientists, technicians, officers and crews of the RRS Discovery (cruise D381), RV Celtic Explorer (cruise CE13001), and RRS James Cook (cruises JC085, JC087 and JC090). The study is part of the project OSMOSIS, funded by NERC grants NE/I020083/1 and NE/I019905/1 and supported by NSF award OCE 1155676. UB was supported by Cefas funding during the writing of this paper. GMD and KJH were supported during the writing of this paper by the descent of each glider dive. Figure A.1 shows the concentration at 11 m as measured during descents and ascents. After 11 th August, c(O 2 ,11 m) in the ascents is higher than in the descents. The magnitude of this difference increased over time, especially in the first metres of the water column down to the deep chlorophyll maximum (data not shown). However, during the night c(O 2 ) values measured during the ascents and descents matched again ( Figure A.1c).
Sunlight seems therefore to be a possible factor causing this difference. This was possibly related to the different angle that the optode had with respect to the incident light according to the direction of the glider ( Figure A.2). The foil was virtually parallel to the surface in the ascents and more angled with respect to the incident light during the descents. This means that the probe was hit directly by the light when the glider went towards the surface, whereas it received less light when it went towards the deep. This was not enough to explain why the two phases of the dives are different in the last month of measurements, because otherwise this phenomenon would have been visible throughout the whole time series. There must have been therefore a new factor that, interacting with the foil and with the light, caused the difference between ascents and descents in this part of the year. The increasing mismatch between phases ( Figure A.1 a-c) also showed that this new factor had a growing influence on the sensor over time.
In the last month of the dataset there was also an increase in c(O 2 ) in the otherwise overall stable minimum c(O 2 ), c min (O 2 ) ( Figure 14b). Being distant from the surface and from the euphotic depth z eup , this deep water mass was expected to be stable because it was not exposed to the big perturbations due to air-sea exchange and biological productivity. After 11 th August there was a fast and un-interrupted increase of c min (O 2 ) that reached 226 μmol kg -1 . Considering that this sharp increase in c(O 2 ) at depth began at the same time as the discrepancy between ascents and descents (on 11 th August), these events were considered to be caused by the same factor. The descents seem to be less affected, while ascents show obviously unrealistic values during the day ( Figure A1a). However, descents still show an increasing pattern over time, showing that data cannot be used despite the direction of the glider movement.  Biofouling of the foil was the most likely factor behind the phenomena just described. It probably developed on top of the optode foil after the beginning of the productive period, when chlorophyll a concentration at the top of the water column was higher than in the rest of the year (beginning of July 2013, Figure 2h). This is usually a proxy for the presence of high phytoplankton biomass, which makes it plausible that phytoplankton started to grow into a biofilm on the foil. The algae, producing more O 2 when exposed to direct and stronger light (during ascents), would have caused the difference between profiles in different phases. O 2 produced by the biofilm would have given high c(O 2 ) readings not reflecting the actual c(O 2 ) in the water column.
Furthermore, the amount of gas released by the biofilm would have been proportional to its biomass -the growth of the biofilm would explain why there was an increase in the difference between phases, of c(O 2 , surface) and of c min (O 2 ). At the recovery of the glider, all the sensors were covered by a green biofilm (Stephen Woodward, personal communication). The data collected after 11 th August are therefore considered not valid for the scope of this study. As a lesson learnt, datasets should be checked, especially when missions last for several months in productive areas; discrepancies between ascents and descents appearing during the day and disappearing during the night, and the increase of values in deep water masses that are usually stable should be signs to look for to spot the possible presence of biofouling and question the validity of the data.

Appendix B -Method Sensitivity
In order to test the sensitivity of the method and determine the uncertainties associated with the N estimates discussed above, we assessed the influence of different parameters and choices made.
If glider ascents are used in the calculations, the mean N is 2 % greater than calculated with descents ( Figure   A.3). However, considering that the optode was influenced in a different way during ascents and descents, this calculation could have been potentially influenced by the initial growth of biofouling at the end of the dataseries. Since the mean and standard deviation of z eup were 60 m and 15 m respectively, mean and total cumulative N were recalculated using 45 m and 75 m as z lim (Table A.1Error! Reference source not found.). These shallower and deeper z lim were considered respectively an underestimation and overestimation of N in the euphotic zone. The shallower z lim measure in fact only the very productive top layer of the water disregarding the respiration happening deeper down. A deeper z lim based on the deepest z eup , instead, accounts for all the respiration happening in the deeper parts when the euphotic layer is thinner. The mean difference between N determined using z lim = 45, 60 and 75 m was used as a measure of the uncertainty associated with N eup (±6.3 mmol m -2 d -1 , ±2.1 mol m -2 , ±30%). Another test was performed to assess the sensitivity of the method to the 7-day length of the averaging bins. N was recalculated binning over 1 day, 3 days, 5 days, 7 days, 9 days, 11 days, 13 days and 15 days (Table A.2).
The maximum change with respect to the values averaged over 7-day bins was obtained using 15-day bins, which increased mean N by 6.5 %. This value was one order of magnitude smaller than the uncertainty related to changes in z lim so was ignored in the error budget. Therefore, the uncertainty of ±30 % estimated from the choice of z lim was used, as the uncertainty introduced by different bin lengths was considered to be negligible in comparison.

Supplementary material -Calibration of the gliders
The first two missions were calibrated using CTD casts from the cruises that visited the site to deploy and recover the gliders. These casts were performed when the ship was close to the profiling gliders (figure 2.13ab). SG599 was instead calibrated using CTD casts from three different cruises. Apart the deployment and recovery cruises, another one (JC087) visited the site in the middle of the mission (locations in figure 2.13c).
Seven different calibrations were then performed during the study, whose details are listed in Table Sup