Benthic oxygen fluxes in a coastal upwelling system (Ria de Vigo, NW Iberia) measured by aquatic eddy covariance

Organic carbon mineralization and nutrient cycling in benthic environments are critically important for their biogeochemical functioning, but are poorly understood in coastal up welling systems. The main objective of this study was to determine benthic oxygen fluxes in a muddy sediment in the Ria de Vigo (NW Iberian coastal upwelling), by applying the aquatic eddy covariance (AEC) technique during 2 campaigns in different seasons (June and October 2017). The main drivers of benthic fluxes were studied and compared among days in each season and between seasons. The 2 campaigns were characterized by an upwelling-relaxation period followed by a downwelling event, the last of which was due to the extratropical cyclone Ophelia in October. The mean (±SD) seasonal benthic oxygen fluxes were not significantly different for the 2 campaigns despite differences in hydrodynamic and biogeochemical conditions (June: −20.9 ± 7.1 mmol m−2 d−1 vs. October: −26.5 ± 3.1 mmol m−2 d−1). Benthic fluxes were controlled by different drivers depending on the season. June was characterized by sinking labile organic material, which enhanced benthic fluxes in the downwelling event, whereas October had a significantly higher bottom velocity that stimulated the benthic fluxes. Finally, a comparison with a large benthic chamber (0.50 m2) was made during October. Despite methodological differences between AEC and chamber measurements, concurrent fluxes agreed within an acceptable margin (AEC:benthic chamber ratio = 0.78 ± 0.13; mean ± SD). Bottle incubations of water sampled from the chamber interior indicated that mineralization could explain this difference. These results show the importance of using non-invasive techniques such as AEC to resolve benthic flux dynamics.


INTRODUCTION
Coastal regions are active biogeochemical areas which play an important role in the global carbon cycle. Although coastal regions occupy only 7% of the total oceanic surface, between 20 and 50% of the total oceanic primary production occurs in these areas (Wollast 1998), and they support approximately 50% of the total benthic carbon mineralization (Middelburg et al. 1997). Primary production generates a large amount of labile organic matter in the coastal margins, which, to a large extent, is deposited on the sediment surface. This particulate matter can subsequently be subjected to resuspension, advection, or biogeochemical transformations, which promote the benthic mineralization of organic carbon and benthic regeneration of nutrients (Wollast 1998). Conversely, these recycled nutrients support the large primary production of coastal areas (Grenz et al. 2010). Oxygen is the most energetically favourable oxidant to react in the direct mineralization of organic carbon and also to re-oxidize reduced products of anaerobic mineralization. Thus, benthic oxygen uptake is the most common proxy used to quantify the total mineralization rate of organic carbon on the seafloor, because it is relatively easy to measure compared to dissolved inorganic carbon fluxes (Froelich et al. 1979, Glud 2008. Among coastal margins, coastal upwelling regions constitute the most biologically productive areas due to the upwelling of cold and nutrient-rich subsurface waters. Additionally, these regions are highly dynamic and are characterized by the occurrence of upwelling/downwelling events, which potentially influence benthic remineralization. Consequently, in these regions, it is particularly important to measure benthic oxygen fluxes over different time scales to integrate the most important biogeochemical processes (Berelson et al. 2013, Reimers et al. 2016. In the context of global climate change, global warming and coastal eutrophication are altering benthic dynamics of coastal zones. The synergistic effect of these stressors increases benthic organic matter mineralization, which intensifies benthic fluxes, and consequently, decreases oxygen concentrations. The number of hypoxic coastal areas have increased 10fold since the first half of the 20 th century and have triggered a severe loss of biodiversity (Diaz & Rosenberg 2008). This increase supports the idea that the distribution and exchange of oxygen between the water column and the sediment are critical parameters to determine the biological status of the coastal environment (Glud 2008).
The non-invasive aquatic eddy covariance (AEC) technique (Berg et al. 2003) to determine benthic oxygen fluxes has several advantages over traditional methods. This technique calculates fluxes from measured turbulent fluctuations of the vertical velocity and oxygen concentration (Berg et al. 2003). AEC has a greater temporal resolution than traditional techniques like benthic chambers and integrates benthic fluxes over a larger area (10−100 m 2 ) . It is a nonintrusive technique, as the benthic fluxes can be determined without disturbing natural hydrodynamic conditions and without resuspending the benthic surface. Moreover, AEC can measure benthic fluxes in permeable sediments, whereas traditional benthic chambers cannot mimic the natural porewater exchange produced by horizontal pressure gradients well (Cook et al. 2007, Berg & Huettel 2008. In fact, most of the previous coastal studies relying on the AEC technique have been carried out in permeable sediments (Huettel et al. 2014and references therein, Attard et al. 2015, Reimers et al. 2016. Only a small number of AEC measurements have been carried out in muddy coastal sediments (Berg et al. 2003, Glud et al. 2010, Holtappels et al. 2013, Attard et al. 2016.
Our study was conducted in the Ria de Vigo (Fig. 1), located in the NW Iberian coastal upwelling system. This area is characterized by high primary production during the upwelling season when nutrient-rich subsurface waters reach the sea surface and trigger substantial phytoplankton blooms (up to 5 mg chlorophyll a m −3 , Figueiras et al. 2002). This enables the Rias Baixas, including the Ria de Vigo (Fig. 1), to support the highest mussel production in Europe (Figueiras et al. 2002). In spite of the importance of benthic mineralization as a source of nutrients to sustain primary production in coastal systems, there are few previous investigations of benthic fluxes in the Rias Baixas (Alonso-Pérez & Castro 2014, Alonso-Pérez et al. 2015. These studies estimated that benthic mineralization is responsible for 40% of the total inorganic nitrogen that supports primary production in the Ria de Vigo. These studies were carried out in muddy sediments, which are abundant in the ria, playing a key role in organic carbon mineralization and in benthic nutrient regeneration (Vilas et al. 2005, Alonso-Pérez & Castro 2014. Based on this background, the main objectives of this study were to (1) determine the magnitude of benthic mineralization in a muddy area of the Ria de Vigo by implementing the AEC technique and (2) establish the principal forcing factors that control the benthic oxygen flux dynamics during the study periods. In order to do this, AEC benthic fluxes were measured during 2 campaigns, June 2017 (upwelling season) and October 2017 (transition from upwelling to downwelling season), with different hydrodynamic and biogeochemical conditions.

Study area
The Ria de Vigo is 1 of the 4 Rias Baixas (temperate coastal embayments, Evans & Prego 2003) located in the NW Iberian coastal upwelling system (Fig. 1). It is characterized by a V-shaped transversal profile with a NE−SW central channel. It has a length of 21 km from Cíes Islands to the Rande Strait. The inner part and central axis of the ria consist of predominately fine-grained sediment, compared to coarser-grained sediment that dominates the outer part of the ria (Vilas et al. 2005).
Hydrography and residual circulation in the ria are largely modulated by the seasonal pattern of shelf winds caused by the meridional migration of the Azores High. Responding to shelf winds, the Ria de Vigo presents a classical estuarine circulation. Northerly winds (prevalent from March/April to September/October) produce an offshore Ekman transport of surface waters and the upwelling of cold and nutrient-rich subsurface waters. Conversely, southerly winds generate an onshore Ekman transport triggering the accumulation of oceanic water onto the coast and the downwelling of surface waters (Fraga 1981, Nogueira et al. 1997, Figueiras et al. 2002. How ever, this seasonal cycle only explains about 10% of the variability in the wind regime off northwest Spain, whereas more than 70% of the variability arises from periods < 30 d (Álvarez-Salgado et al. 2003). In fact, the upwelling/downwelling seasons in the region are described as a succession of upwelling/ downwelling events lasting between 3 and 10 d (Nogueira et al. 1997, Gilcoto et al. 2017. Shelf winds are the main factor controlling the subtidal circulation of the ria, although other important factors can modify this circulation, such as river discharge . Our sampling site was located in the main channel in the inner part of the ria (Fig. 1). This area is characterized by unvegetated muddy sediments and is mainly exposed to NE−SW currents (Villacieros-Robineau et al. 2013). The presence of molluscs, ophiuroids, and pennatulacean anthozoans has been documented in the main channel of the Ria de Vigo (Aneiros et al. 2015).

Water sampling
Data collection occurred during 2 different seasons: summer (12−30 June) and autumn (2−20 October) of 2017. Our sampling site ( Fig. 1) was located in the inner ria (42°16.186' N, 8°41.985' W) with an average depth of 21.8 m. During each campaign, the site was sampled 3 times a week over a 3 wk period. Data collection was carried out during research campaigns on board the RV 'Mytilus'.
Vertical profiles of temperature, salinity, oxygen, and fluorescence were obtained with a CTD SBE 911 plus (Sea-Bird) coupled to a bottle rosette system with 12 Niskin bottles of 10 l capacity to obtain seawater samples for oxygen, chlorophyll a (chl a), and particulate organic carbon (POC). Three or 4 water column depths were sampled depending on the location of the subsurface chlorophyll maximum depth determined through CTD fluorescence. Dissolved oxygen was determined by Winkler potentiometric titration. The precision of the method was ± 0.32 μM.
To measure chl a, 100 ml of seawater were filtered using GF/F filters with a pore size of 0.7 μm. After filtration, samples were frozen (−20°C) until pigment extraction in 90% acetone over 24 h in the dark at 4°C. Chl a was analysed through fluorometric determination using a Turner Designs fluorometer (Yentsch & Menzel 1963). The precision of the method was ± 0.07 μg l −1 . For POC analysis, 0.5 l seawater samples were filtered on pre-weighted and pre-combusted (4 h, 450°C) GF/F filters of 0.7 μm pore size. Measurements of POC were determined with a Perkin Elmer 2400 CNH analyser. The precision was ± 0.3 μM.

Residual and total current
Total current was obtained using an upward looking 400 kHz Nortek Aquadopp Profiler moored at 25.8 m depth at the study site (42°15.870' N, 8°42. 132' W, Fig. 1), between 8 June (20:00 h UTC) and 13 July (09:45 h) 2017 for the June campaign and between 25 September (20:00 h) and 30 October (11:40 h) 2017 for the October campaign. Twelve cells (bin size 200 cm) were established through the water column. The residual current was extracted using a low-pass symmetric digital filter (PL33) with a cutoff period of 33 h (Flagg et al. 1976, Beardsley et al. 1985 to remove tidal frequencies. This analysis was carried out in MATLAB.

Aquatic eddy covariance
Benthic oxygen fluxes were measured with the AEC technique (Berg et al. 2003). Instrumentation for the AEC measurements was mounted on a light frame (Berg & Huettel 2008) and consists of a Nortek Vector acoustic Doppler velocimeter (ADV) and a dual sensor (RINKO EC, JFE Advantech), which al lows high-frequency (64 Hz) recording of temperature and oxygen concentrations at the same volume of water ). The AEC system was de ployed at our sampling site (42°16.217' N, 8°41.849'W, Fig. 1) at an average depth of 20.1 m. It was deployed by divers and care was taken to direct the ADV x-axis into the main current direction. A total of 18 deployments, lasting between 7 and 9 h to preferentially cover either tidal inflow or outflow, were carried out during the June and October campaigns (Table 1). The ADV measuring volume was ap proximately 11 cm above the bottom. A stable dual oxygen−temperature optode (miniDOT, PME) was attached to the frame for independent oxygen concentration measurements every 5 min.

Benthic chamber
Benthic oxygen fluxes were also measured in a deployed benthic chamber to complement the AEC results. Benthic chamber deployments were carried out at our field site (42°16.226' N, 8°41.860' W, Fig. 1) at an average depth of 21.1 m. The chamber was carefully deployed for the same days as the AEC deployments using a mooring line. The deployment operation was checked by divers. Details of the benthic chamber measurements are described by Ferrón et al. (2008). In short, the benthic chamber has a large volume of 140 l and covers a large total sediment area of 0.50 m 2 (one of the largest benthic chambers  (Table 1). Benthic chamber fluxes were calculated from changes in oxygen concentrations during the incubation period by an empirical linear fitting. Uncertainties of fluxes were related to the fit of the data to a linear function and the propagation of random errors (Ferrón et al. 2008). On 2 occasions during the October campaign, 6 replicate samples of seawater from the benthic chamber interior were manually collected in 125 ml Winkler flasks by divers just before the beginning of the benthic chamber incubation period. Three replicates were fixed on board to measure oxygen concentration at the starting point of the incubation, and the other 3 samples were incubated in situ during the benthic chamber deployment time. Dissolved oxygen was measured by a Winkler potentiometric titration as described in Section 2.2 for the water column. Oxygen consumption in the water column inside the benthic chamber was calculated as the oxygen concentration difference at the initial and final time of the incubation.

Sediment samples
Surface sediment samples were collected using a Van Veen grab sampler on 3 different days (8 June, 13 July, and 30 November). A total of 3 replicate samples were obtained for each day. Granulometric analysis was carried out using an LS 13 320 Beckman Coulter system on non-homogenized sub-samples of roughly 1 g after removal of organic matter and carbonates by adding 10% H 2 O 2 and 25% HCl, respectively. Organic carbon content analysis was performed on triplicates of 20 mg sub-samples after eliminating carbonates by repeatedly adding 10 μl of 25% HCl (Zúñiga et al. 2016). Organic carbon was analysed using a Perkin Elmer 2400 CNH analyser as previously explained for POC water samples.

Upwelling index
Shelf winds were obtained from the Seawatch buoy located in Cape Silleiro (42°7.2' N, 9°25.8' W, code 2248) anchored at a depth of 600 m. These hourly data were provided by Puertos del Estado (www.puertos.es). Wind data measured by the Cape Silleiro buoy at a 3 m height were converted to 10 m height winds by applying a logarithmic profile with a roughness length of 0.002 m . Shelf wind data were filtered using a low-pass filter. The upwelling index (UI) was calculated from the shelf winds following Bakun (1973): (1) where ρ air is the average density of the air (1.22 kg m −3 at 15°C), C D is the empirical drag coefficient (1.3 × 10 −3 according to Hidy 1972), f is the Coriolis factor at 42° latitude, and ρ sw is the density of the sea water (1025 kg m −3 ). |V | is the module and V y is the north component of the wind velocity at a height of 10 m. Positive values of UI correspond to upwelling conditions and negative values indicate downwelling.

Aquatic eddy covariance
Fluxes of oxygen were extracted from the raw eddy covariance data following the multi-step AEC process described below and shown in Fig. 2. First, oxygen concentration was calibrated against the stable independent dual oxygen−temperature sensor data (Fig. 2b). All 64 Hz data were then reduced to 8 Hz data, which reduces noise while still providing sufficient resolution to contain the full frequency spectrum carrying the flux signal (Berg et al. 2009). Benthic oxygen fluxes were then computed using the software package Eddy Flux version 3.2 (P. Berg unpubl.).
Benthic oxygen fluxes were calculated using highresolution data of vertical velocity (v z ) and oxygen concentration. The mean benthic oxygen flux Due to the response time of the fast-responding oxy gen sensor combined with its position down- stream from the ADV's measuring volume, a time shift correction was applied . Following the standard practice for eddy covariance data processing (Fan et al. 1990, Lorrai et al. 2010, this was done by repeating the outlined flux extraction procedure, while shifting the 8 Hz oxygen concentration data back in time, 1/8 s at a time, until the numerically largest flux was identified (Eddy Flux 3.2, P. Berg unpubl.).
A data quality check was performed to remove erroneous measurements produced by collisions between floating particles and the RINKO EC sensor. Cumulative fluxes (Fig. 2c) were calculated by a running sum of instantaneous vertical fluxes (v ' z C ') over time within each 15 min interval. In each burst, the cumulative flux was verified for any anomalous fluctuation that was not due to natural turbulence ( Fig. 2c; Berg et al. 2013). Data recorded when the  Berg et al. 2009) were also removed. The percentage of the fluxes kept after applying the data quality check was 38% (for June and October). Finally, the selected 15 min fluxes (Fig. 2d) were averaged over each deployment to calculate mean daily fluxes. Negative fluxes indicate oxygen uptake by the sediment and positive fluxes indicate oxygen release.

Statistical analyses
A Wilcoxon-matched pair test and a Mann-Whitney U-test were performed to identify differences between AEC fluxes and benthic chamber fluxes, to determine seasonal variation in the benthic fluxes, and to analyse which factors controlled them. For all statistical tests, the statistical software Statistica 6.0 (StarSoft) was used with a confidence level of α = 0.05. All mean values are reported ±1 SD.

Hydrodynamic conditions
During the summer campaign in June 2017, 2 hydrodynamic periods were identified. The first period (12−25 June) was characterized by an intense positive upwelling pulse (12−16 June), followed by a wind relaxation (Fig. 3a). In the second period (26− 30 June), there was a downwelling event with negative UI values.
Water column temperature responded to these upwelling and downwelling events (Fig. 3b). During the upwelling event, the water column was weakly stratified, with temperatures of 17°C at the sea surface and colder than 15°C at the bottom. However, during the wind relaxation, the water column was more stratified and the thermocline strengthened. This situation was maintained until the downwelling event when there was a thermal homogenization and warming of the entire water column.
Total current velocity had values between −30 and 30 cm s −1 with a strong tidal signal (Fig. 3c). Residual current ranged between −10 and 10 cm s −1 (Fig. 3d), and it responded to the observed upwelling and downwelling events. During the upwelling event and subsequent wind relaxation, the estuarine circulation was positive with an inflow through the bottom layer and outflow at the surface layer. After that, dur-ing the downwelling (26−29 June), residual current changed to a negative circulation pattern with inflow through the upper layer and outflow through the bottom layer.
Similar to the summer campaign, 2 contrasting hydrodynamics periods were discerned for the autumn campaign. During the first period (2−14 October), there was an upwelling event with positive UI values (Fig. 4a) followed by a wind relaxation. The second period (15−20 October) was characterized by a strong downwelling forced by extratropical cyclone Ophelia during 15−16 October, which caused very negative UI.
The temperature profile showed changes in the water column related to these 2 different periods (Fig. 4b). During the upwelling-relaxation period (2−14 October), the water column was slightly stratified, with values of 17°C at the sea surface and 14°C at the bottom, followed by notable thermal homogeneity. In the downwelling period, there was a warming of the entire water column.
Residual currents followed the wind regime, with a positive estuarine circulation during the upwellingrelaxation period (Fig. 4d). These conditions persisted until the circulation reversed due to downwelling forced by extratropical cyclone Ophelia.

Biogeochemical conditions
The biogeochemical variables responded to the hydrodynamic and thermohaline conditions previously explained. During the first summer period, maximum oxygen concentrations were reached at the sea surface with values up to 240 μM (Fig. 3e) and minima at the bottom (12−25 June). There were high surface concentrations of POC and a subsurface maximum located at approximately 10 m depth with values up to 40 μM (Fig. 3f) during the relaxation period. Chl a (not shown) had a similar pattern with sea surface maximum of 15 mg m −3 showing a high significant positive correlation with POC (R 2 = 0.86, p < 0.05), which indicates that this particulate or ganic matter was highly labile. Hence, the first summer period was characterized by the sinking of this labile organic matter. During the subsequent downwelling event, there was a ventilation of the water column increasing the oxygen concentration in the entire water column, with concentrations higher than 230 μM. Associated with this homogenized water column, a decrease in POC concentrations was observed (Fig. 3e,f). In the first period of the autumn campaign, there was a maximum oxygen concentration at the sea surface and minima beginning at 10 m (Fig. 4e). POC showed the same pattern with a sea surface maximum (Fig. 4f). These conditions persisted until wind relaxation, when there was an intense mixing (6− 13 October) that disrupted the previous surface maxima of oxygen and POC. During the subsequent downwelling period (15−20 October), oxygen concentrations increased, reaching values higher than 200 μM (Fig. 4e). There was also a slight increase in POC in the last days of this period (18−20 October, Fig. 4f), but these POC concentrations were consistently lower than during the summer campaign.

Aquatic eddy covariance fluxes
The AEC system was deployed at a site above sandy mud sediment (clay 20%, silt 59%, sand 21%) according to Folk's classification (Folk 1954), with an organic carbon content of 3.1 ± 0.2% by weight. Due to these sediment characteristics, this site is subject to the formation of a fluff layer over the seabed. Daily mean benthic fluxes were obtained from high resolution 15 min interval flux data (Fig. 2). All of the fluxes determined in this study were negative, indicating oxygen uptake by the sediment. During the summer campaign, bottom velocity ranged from 1.1 to 3.5 cm s −1 (Fig. 5a). In the first period of the campaign (12− 25 June), bottom velocity reached a minimum value on 19 June and increased during the following downwelling event (26−30 June). Bottom oxygen concentrations increased progressively during the summer campaign (Fig. 5b). Temperature showed a similar pattern with a minimum value on 19 June and maximum value on 28 June (Fig. 5b). Benthic fluxes were low during the summer upwelling-relaxation period with a minimum of −15.1 mmol m −2 d −1 on 12 June (Fig. 5c). In the subsequent downwelling period, fluxes were higher, reaching a maximum of −35.0 mmol m −2 d −1 on 26 June. Mean benthic oxygen flux was −16.4 ± 1.2 mmol m −2 d −1 during the first period and −28.6 ± 5.7 mmol m −2 d −1 during the second period, During the autumn campaign, bottom velocity ranged from 0.7 to 6.0 cm s −1 with a maximum on 6 October (Fig. 6a). Bottom oxygen concentrations increased progressively during the campaign from values of 131 μM on 4 October to 200 μM on the last day (20 October). Temperatures showed a similar pattern to oxygen concentrations, with lower temperatures during the upwelling-relaxation period and higher temperatures during the downwelling event. Benthic fluxes were not computed for 11, 13, and 16 October, as bottom velocities were very low for these deployments (<1 cm s −1 ). The lowest benthic flux (−23.6 mmol m −2 d −1 ) was measured on 9 October. Mean benthic oxygen fluxes were −25.7 ± 1.9 and −27.7 ± 5.2 mmol m −2 d −1 for the first and second autumn period, respectively. There were no large changes in benthic fluxes, with a low monthly coefficient of variation of 12% for October (Fig. 6c).
Daily coefficients of variation (7−47%) were of the same order of magnitude compared to monthly variation (34%) on most of the days for June. In contrast, in October daily variability (13−37%) was much higher than monthly variability (12%). In fact, only on 18 October was daily variability of the same order of magnitude as monthly variability.

Benthic chamber fluxes
In situ benthic chamber fluxes showed similar patterns to those obtained by AEC for the October campaign. The initial temperatures and oxygen concentrations inside the benthic chamber (Fig. 7a) were similar to daily mean values measured by the AEC system. Similarly, oxygen concentrations and temperatures were higher for the second downwelling period than during the previous upwelling-relaxation period, reaching the maxima on 20 October.
The benthic oxygen fluxes measured by in situ benthic chamber incubations produced higher fluxes than the AEC system, with a minimum of −29.1 ± 1.2 mmol m −2 d −1 on 9 October (Fig. 7b). Benthic chamber fluxes had greater variability than those measured by AEC mean daily fluxes, with a coefficient of variation of 30%. The higher variability was mainly due to the large fluxes found on 13 and 20 October. For comparison, mean benthic oxygen fluxes using the same deployment days as AEC were −31.4 ± 3.0 mmol m −2 d −1 during the upwelling-relaxation event and −39.8 ± 5.3 mmol m −2 d −1 during the downwelling event. Seawater collected from the benthic chamber interior and incubated in situ in the Winkler flasks during the deployment time had an average oxygen uptake of 10.2 ± 1.8 mmol m −2 d −1 .

Benthic fluxes
The mean benthic oxygen flux of −23.1 ± 6.4 mmol m −2 d −1 obtained by AEC from 13 deployments over the muddy sediments of Ria de Vigo (Table 1) is comparable to those (−32.9 mmol m −2 d −1 ) derived from global-scale sediment−water flux data for coastal regions (Boynton et al. 2018). In addition, the mean in our study is close to the average oxygen demand (−34 ± 10 mmol m −2 d −1 ) obtained previously by Alonso-Pérez & Castro (2014) using benthic chambers during 4-season campaigns in the Ria de Vigo. To our knowledge, the AEC technique has only been applied in a handful of muddy coastal marine sediments (Berg et al. 2003, Glud et al. 2010, Holtappels et al. 2013, Attard et al. 2016. Our mean value for the Ria de Vigo is within the range of values reported for fine-grained muddy sediments in the shallow bay of Aarhus (Berg et al. 2003) and the glacial fjord Loch Etive (Holtappels et al. 2013). In fact, all of these values followed the simple power relationship between benthic oxygen uptake and water depth described by Glud (2008) based on a global database of benthic oxygen fluxes (Fig. 8).
A comparison between results from AEC and the benthic chamber was carried out in October 2017 using the same deployment days. The mean benthic flux was −26.5 ± 3.1 mmol m −2 d −1 for AEC and −34.8 ± 5.7 mmol m −2 d −1 for the benthic chamber ( Fig. 9). Although these fluxes are significantly different (p = 0.04, Wilcoxon matched pair test), they agree within an acceptable margin given the fundamental differences between the 2 approaches and possible heterogeneity of the sediment surface at our site. Specifically, the mean daily ratio between AEC and benthic chamber fluxes was 0.78 ± 0.13, which is comparable to other ratios ranging from 1.4 to 0.75 found in previous studies in muddy sediments (Holtappels et al. 2013 and references compiled by Attard et al. 2015). A likely explanation of the difference is that mineralization processes occurred in the bottom water en closed in the chamber. Indeed, bottle incubations of water sampled from the chamber interior after it was deployed indi cated an average oxygen consumption of 10.2 ± 1.8 mmol m −2 d −1 . This result emphasizes the importance of deploying the benthic chamber with the utmost care, especially when the sediment is a fluffy muddy sediment as was the case at our study site.

Main forcing factors of the benthic fluxes
The intensive hydrographic and biogeochemical sampling during the June and October 2017 campaigns allowed us to investigate the controls of changes in benthic oxygen fluxes during 2 different seasonal periods for the NW Iberian upwelling system. June is part of a substantial upwelling season, when northerly winds predominate, causing the uplift of subsurface nutrient-rich waters and triggering intense phytoplankton blooms inside the Rias Baixas. In contrast, October usually corresponds to the transition period between the upwelling and downwelling seasons (Nogueira et al. 1997, Figueiras et al. 2002. In June 2017, there was a strong upwelling event that produced a proliferation of large diatoms (M. Froján et al. unpubl. data). The stratified water column of the following days favoured this phytoplankton growth, with chl a levels as high as 15 mg m −3 . During the subsequent downwelling, this labile organic matter (senescent phytoplankton cells) began to sink, and conveyed by the water circulation it was deposited on the sediment surface. In October 2017, there was an upwelling event with a slightly stratified water column followed by a wind relaxation during which there was a thermal homogenization. Under these conditions, phytoplankton proliferation was much lower than during the phytoplankton bloom of June 2017. Subsequently, the occurrence of extratropical cyclone Ophelia triggered an intense downwelling, which caused the inflow of warm and well-oxygenated oceanic surface waters into the entire water column.
Changes in hydrodynamics and biogeochemical conditions of the water column affect benthic oxygen fluxes. Studies have shown the influence of oxygen, temperature, and organic matter supply on the magnitude of benthic oxygen fluxes (Cowan et al. 1996, Lomas et al. 2002, Hopkinson & Smith 2005, Murrell & Lehrter 2011, Alonso-Pérez & Castro 2014. Additionally, the way in which hydrodynamics control benthic oxygen fluxes has been analysed in several papers , Berg & Huettel 2008. Based on the distinct water column conditions measured during the June and October 2017 sampling periods, we expected different magnitudes of benthic oxygen fluxes for these 2 periods. Daily benthic oxygen fluxes in June were characterized by high variability, distinguished by significantly smaller benthic oxygen fluxes for the upwellingrelaxation period (average of −16.4 ± 1.2 mmol m −2 d −1 ; Fig. 10a) compared to the mean value of −28.6 ± 5.7 mmol m −2 d −1 during the downwelling event at the end of June (26−30 June, Mann-Whitney U-test, p = 0.03, Fig. 10a). The large benthic fluxes at the end of June are associated with significantly higher temperatures (Mann-Whitney U-test, p = 0.03, Fig. 10b), and relatively higher, though not significant, oxygen concentrations (Mann-Whitney U-test, p = 0.1, Fig. 10c) and bottom velocity (Mann-Whitney U-test, p = 0.5, Fig. 10d). Interestingly, it was during this summer downwelling event that the largest benthic oxygen fluxes for the entire study were measured, larger than the mean benthic oxygen flux for the autumn downwelling event due to extratropical cyclone Ophelia. Even though the stronger bottom velocity in the October campaign (4.3 ± 1.2 cm s −1 ) intensified benthic oxygen exchange, the highest mean flux during the summer downwelling had relatively lower bottom velocities (2.3 ± 0.9 cm s −1 ) and similar temperature and oxygen concentrations. organic matter triggered the high benthic oxygen fluxes for the summer downwelling. Consequently, during the summer campaign, the principal forcing factor that modulated the benthic oxygen fluxes was the supply of organic matter. In October, no significant differences occurred in benthic oxygen fluxes between the upwelling-relaxation period (−25.7 ± 1.2 mmol m −2 d −1 ) and the strong downwelling after 14 October (−27.7 ± 5.2 mmol m −2 d −1 , Mann-Whitney U-test, p = 0.6, Fig. 10a). Additionally, no significant differences were observed in any of the studied variables (temperature, oxygen concentration , bottom velocity, and POC) between both periods. Comparison of the mean benthic fluxes for the June (−20.9 ± 7.1 mmol m −2 d −1 ) and October (−26.5 ± 3.1 mmol m −2 d −1 ) campaigns resulted in similar values, even though available POC concentrations were much higher in June (Fig. 10e), and bottom temperature and oxygen concentrations were similar (Fig. 10b,c). Thus, we suggest that the much higher bottom velocity (Fig. 10d) during the October campaign is the key factor enhancing benthic oxygen fluxes compared to the June campaign. Intensification of bottom shear stress produced by strong currents can resuspend surface sediments, mixing the organic matter and reduced compounds into the overlying water and favouring their oxidation (Almroth et al. 2009, Niemistö et al. 2018, Camillini et al. 2021. In fact, intense resuspension events in the Ria de Vigo are more likely to happen in autumn (Villacieros-Robineau et al. 2013). Additionally, an increase in bottom velocity reduces the thickness of the diffusive boundary layer, producing an increase in diffusive gradients and deeper oxygen penetration into the sediment ).
Using a global sediment−water flux data set for coastal zones, in which only flux measurements obtained through benthic chambers and sediment cores were considered, Boynton et al. (2018) concluded that benthic oxygen fluxes were modulated by temperature, depth, and the organic matter supply. In contrast to other techniques, AEC is not an invasive technique that interferes with the natural environmental conditions. Thus, through use of this technique, the important role hydrodynamic conditions can play on benthic fluxes , Attard et al. 2015, Reimers et al. 2016) has been documented. In this first study applying AEC in the Ria de Vigo, we document that benthic fluxes are modulated by different key factors, most notably, hydrodynamics (October campaign) and the supply of labile organic matter fueling benthic oxygen uptake (June campaign).