Macrofauna as a major driver of bentho-pelagic exchange in the southern North Sea

The contribution of sediments to nutrient cycling of the coastal North Sea is strongly controlled by the intensity of ﬂ uxes across the sediment water interface. Pore-water advection is one major exchange mechanism that is well described by models, as it is determined by physical parameters. In contrast, biotransport (i

Coastal sediments are important sites for organic matter turnover, amounting to about 10-20% of primary production (Heip et al. 1995) and providing important habitats and ecosystem functions like nutrient regeneration.These can, however, be compromised by human activities (Costanza et al. 1997).The North Sea (NW Europe) is a shallow coastal sea with strong human impact, as it is bordered by highly industrialized countries.Human activities are still expected to increase and interact in a yet unknown way with climate change aspects, like temperature rise (Emeis et al. 2015).
Various human activities such as eutrophication, sand removal, bottom trawling by demersal fisheries, and installation of offshore wind farms affect the seafloor, which may change intensities and ratios of biogeochemical processes (Tillin et al. 2006;Coates et al. 2014;Solan et al. 2016).However, our understanding of benthic biogeochemical processes in the North Sea is still insufficient to quantify present effects and to predict possible future changes.In particular, there is a lack of studies combining both spatial and temporal dynamics, and the influence of benthic fauna.Earlier studies (Lohse et al. 1993) contrasted winter and summer situations, but did not cover the entire seasonal cycle.Oehler et al. (2015a,b) covered the entire seasonal cycle but focused on a small area (around Heligoland, German Bight) that is strongly influenced by riverine runoff.More recently, Bratek et al. (2020) studied the spatial differences in benthic organic matter and nitrogen turnover in the German Bight across different hydrodynamic and sedimentary setting, but did not resolve the seasonal dynamics.
Apart from the supply of organic matter as an important driver of biogeochemical processes, sediment type and macrofauna have a strong effect on benthic biogeochemical processes and bentho-pelagic coupling.Sediment grain size is highly correlated with permeability, organic matter content, bottom shear stress, and fauna composition (Kröncke et al. 2004;Meyer et al. 2018;Neumann et al. 2019).Ahmerkamp et al. (2017) demonstrated that benthic oxygen consumption in advective systems is ultimately a function of sediment grain size and flow velocity.They found that coarse sediment with high permeability and more intense advective flushing had higher oxygen consumption rates than finer sediment.In permeable systems, the influence of macrofauna organisms on nutrient fluxes is also very limited (Huettel and Webster 2001).In less permeable sediment, macrofauna however strongly influences benthic biogeochemical processes through activities that affect the sediment matrix (Kristensen et al. 2012).Macrofaunal sediment reworking activities (e.g., burrowing and foraging of animals) mix the sediment, and bioirrigation activities actively move water between burrows and the water column.Bioturbation and bioirrigation thus enhance biogeochemical processes and oxygen consumption (Blackburn 1988;Braeckman et al. 2014).
Benthic communities are often very diverse and comprise various species, each with specific abundances, life cycles, behaviors, body dimensions, and food sources.This diversity implies that each species has a specific contribution to the intensity of bioturbation (Queirós et al. 2013;Meyer et al. 2019).Accordingly, it is difficult to predict the magnitude of biogeochemical processes in systems with high biodiversity and high abundance of macrofauna.Descriptive nonquantitative indices such as the community bioturbation potential (Solan et al. 2004) and the community bioirrigation potential (Renz et al. 2018;Wrede et al. 2018), aim to simplify complex community effects on exchange processes.Due to their nonquantitative nature, these biotransport proxies do not give absolute values for sediment reworking, bioirrigation, or biogeochemical cycling.Instead, both indices aim to estimate the joint potential of a macrofauna community to increase transport rates, based on the individual functional traits of the species within that community.Recent studies demonstrated that the bioturbation potential correlated with total organic carbon (TOC) content, chlorophyll concentration, oxygenation depth, sediment oxygen consumption, ammonium efflux, and denitrification (Solan et al. 2004, Braeckman et al. 2014;Gogina et al. 2020).The bioirrigation potential correlated with bioirrigation rate, nutrient flux across the sediment water interface, and bioirrigation depth (Renz et al. 2018;Wrede et al. 2018;de Borger et al. 2020).It is thus justified to apply these biotransport indices to estimate macrofaunal impact on biogeochemical processes.
The present study was carried out in the framework of the German project "North Sea -Observation and Assessment of Habitats (NOAH), which set out to quantify the influence of human activity on seafloor integrity, including biological structure and functioning of benthic communities and food webs.NOAH aimed to combine field observations of solute exchange across the sediment-water interface with numerical and statistical modeling to evaluate utilization, status, performance and provision of ecosystem services of seafloor habitats, and ecosystems in the German Bight.
The specific goals of this study were: 1. To establish the spatial and seasonal variability of benthic fluxes of oxygen and nutrients that are sustained by biotransport; 2. To identify the major drivers of the observed variability; and 3. To estimate the contribution of fluxes sustained by macrofauna to the overall bentho-pelagic exchange.
Our approach to address these aims was to sample all distinct, major sediment provinces in the German Bight of the southern North Sea repeatedly in all seasons to examine the benthic processes and environmental conditions.By employing a whole-core incubation method close to in situ conditions, but excluding the effects of waves and currents, we measured benthic fluxes that were driven predominantly by faunal activity.Our observations complement measurements of Ahmerkamp et al. (2017), who sampled the same sites as in our study, but focused on benthic fluxes that were driven by pore-water advection evoked by bottom water flow over rippled beds.Their study showed that benthic oxygen fluxes due to pore-water advection were controlled by volumetric respiration rate, grain size, and current velocity.We will combine results from the present study with the results of Ahmerkamp et al. (2017) to assess the significance of faunal activity for the overall bentho-pelagic exchange.

Study area and sampling
The study sites are located in the German Bight of the southern North Sea in water depths < 50 m (Fig. 1).NOAH sites (A, B, E, I, CCP-G) are characterized by coarse and permeable sand in shallow water along the coast and on the Dogger Bank (Fig. 1).The macrofauna community in these areas is characterized by Gonadiella spp., Spisula spp., Bathyporeia spp., and Tellina fabula (Meyer et al. 2018).The NOAH sites (C, D, F, G, H, CCP-J) in the deeper parts are characterized by much finer sediment with a high silt content.At these deeper stations, the near-bottom turbulence is low due to low tidal currents and low wave impact.Characteristic macrofauna species in these areas are Nucula nitidosa and Amphiura filiformis (Meyer et al. 2018).NOAH site C is located in the regional depocenter "Heligoland Mud Area," and is a unique habitat with muddy sediment, high sedimentation rate, and the highest TOC content of all study sites (Fig. 1, light blue hue).
The sediment characteristics and the corresponding provinces of all sites are summarized in Table 1.
These study sites were sampled during several cruises of RVs Heincke (HE), Alkor (AL), and Maria S. Merian (MSM) in the period 2013-2019 during various seasons.However, not all sites were accessible during each cruise (Table 2) due to time constraints and bad weather conditions.The surface sediment was sampled with a Multicorer (Oktopus Kiel) equipped with transparent acrylic (Poly(methyl methacrylate)) core liners (10 cm inner diameter, 60 cm length), and additional samples  Ex situ whole core batch incubation Ordinarily, 3-4 intact sediment cores per station were selected for ex situ incubation.Cores with visible perturbations such as cracks, voids, or injured animals were rejected.The sediment cores were typically 15-30 cm in length.Figure 2 illustrates the configuration of the incubation setup.The incubation was enclosed with a gas-tight lid that was adjusted to a resulting supernatant height of 15 cm (approximately 1 L volume).The water column was constantly stirred by horizontal propellers (3 blades, 20 mm diameter) that were magnetically driven by external motors.The stirring intensity was adjusted to the highest intensity that did not result in sediment resuspension, to ensure vigorous mixing of the supernatant, and to prevent already suspended particles from settling.Incubations were carried out onboard the research vessels in temperature-controlled labs that were set at in situ seafloor temperatures.Primary production was inhibited by wrapping the cores in aluminum foil.Incubations were terminated when the oxygen saturation dropped below approximately 80%, which was normally after 12-24 h.Depending on the oxygen consumption rate, water samples were drawn from the core supernatant in 2-4 h intervals (typically 10 time steps per incubation) through the use of syringes connected with PVC-tubing.Water samples were then filtered through 0.45 μm syringe filters and stored frozen until analysis in the land-based laboratory.After incubation, samples of uppermost 1 cm sediment layer were retrieved for chlorophyll measurement.The remaining sediment was screened through 0.5 mm sieves to collect the macrofauna.

Oxygen measurement
The oxygen concentration of the core supernatant was measured continuously with fiber-optical optodes (PreSens, NTH-PSt1).A two-point calibration was applied with airsaturated water (100%) and oxygen-free water (0%) that was prepared by addition of sodium sulfite (Na 2 SO 3 ).

Nutrient analysis
Sampled water from incubations was analyzed for ammonium, nitrite, nitrate, phosphate, and silicate with a segmented flow analyzer (Seal Analytical AA3), employing methods based on Grasshoff et al. (1983) and Kèrouel and Aminot (1997).

Flux calculation
The measured oxygen and nutrient concentrations in the core supernatant were corrected for the dilution with replacement water.Net fluxes were calculated as concentration change over time by linear regression of multiple time steps (usually 10).Water samples for nutrient measurements were drawn from the supernatant.The ammonium that is released during remineralization is partially oxidized to NO x either in the upper sediment layer or in the supernatant, thus we cannot differentiate between nitrification within the sediment and water column nitrification.Therefore, we combined the concentrations of ammonium, nitrite, and nitrate into dissolved inorganic nitrogen (DIN) and present only DIN fluxes.

Macrofauna sampling and biotransport indices
The macrofauna was sampled directly from incubated cores during the cruises HE-447, 470, 471, 511 and AL-520, while during cruises HE-395, HE-412, HE-432, and MSM-50, macrofauna was sampled only from additional (not incubated) cores.The sediment content of the cores was sieved through 0.5 mm mesh size, fixated, and stored in 4% buffered formaldehyde seawater solution.Taxa were determined down to species level, counted, and weighed; biomass is given as wet weight.The community bioturbation potential of the benthic community was calculated according to Queirós et al. (2013, eq. 1), where BP c = community bioturbation potential, B i = biomass of species i, A i = abundance, M i = mobility score, and R i = reworking score.Likewise, the community bioirrigation potential was calculated according to Wrede et al. (2018, eq. 2), where IP c = community bioirrigation potential, B i = biomass of species i, A i = abundance, BT i = burrow type, FT i = feeding type, and ID i = injection pocket depth.Both indices were calculated from species-specific biomass and abundance, from specific trait scores to account for differences in life form and behavior, and from the depth of the injection pocket.The values for the individual traits were taken from data compilations (e.g., Queiro´s et al. 2013;Wrede et al. 2018), or were deduced from the literature on macrofauna life cycles and behavior in cases when no published traits scores were available.The method of trait score assignment is described in detail in Queiro ´s et al. ( 2013) and Wrede et al. (2018).Additionally, we further calculated the core term (ct) of bioturbation and bioirrigation potential, but without the traits according to Eq. 3,

Sediment analysis
Sediment samples from incubated cores were stored frozen and in darkness until further processing.The sediment was freeze-dried under light exclusion, and chlorophyll was extracted with 90% acetone at 5 C overnight.Concentrations of chlorophyll a were measured photometrically (Hach Lange DR-6000) applying the method of Lorenzen (1967).

Data treatment, analysis, and availability
Large differences exist among the sites for bottom water chlorophyll, sediment chlorophyll, and community bioturbation potential.To characterize the sampling sites (Fig. 3), we used the relative seasonality by normalizing the with June values as reference (June = 1).
The carbon stock to oxygen flux ratio (n C : J O2 ) was calculated to relate the instantaneous oxygen flux (J O2 in μmol O 2 m −2 d −1 ) to the depth-integrated amount of available organic carbon (n C in μmol C m −2 ).The chlorophyllassociated carbon stock was calculated from depth-integrated chlorophyll concentration (in μg chl a m −2 ), a carbon conversion factor of 50 (Jakobsen and Markager 2016), and a molar weight of 12 g mol −1 to estimate the chlorophyll-associated organic carbon stock (in μmol C m −2 ).The pelagic carbon stock was calculated from depth profiles measured with the CTD-mounted fluorometer.The chlorophyll concentrations of the lower half of the profiles were integrated to consider the stratification during summer.The benthic carbon stock was calculated from direct chlorophyll measurements of sediment samples that were integrated over the uppermost 5 cm layer, as the majority of macrofauna were found herein.By assuming that all oxygen is consumed for the oxidation of organic carbon to CO 2 , the ratio of estimated carbon stock to measured oxygen flux (in μmol O 2 m −2 d −1 ) directly gives the hypothetical period (in days) that the carbon stock can fuel the oxygen flux.
The data of this study were compared to results from in situ measurements of oxygen fluxes by Ahmerkamp et al. (2017).For this comparison, we used results from RV Heincke cruises HE432 and HE447, because sites were sampled simultaneously.We averaged results from sites A, CCP-G for the coarse sand province, results from sites B, I for medium sand province, results from sites D, F, H for fine sand province, and results from site C for the silt province.
Statistical analyses were performed with the Matlab software (version R2019b).We employed the non-parametric Spearman's rank correlation coefficient to assess the independence of parameters, as the Spearman's correlation coefficient is robust against nonlinearity and outliers.Additionally, the variance inflation factors for the identification of collinear parameters were calculated from the diagonal elements of the inverse of the correlation matrix (Belsley et al. 1980).The threshold for collinearity is a variance inflation factor > 5.A partial least squares analysis (PLS) was employed to separate and quantify the specific contributions of individual site parameter (water depth, salinity, temperature, bottom water chlorophyll, sediment chlorophyll, sediment permeability, bioturbation potential, bioirrigation potential, core term), using the site characteristics as independent parameters, and benthic fluxes as dependent parameters, respectively.For PLS analysis, measurements from parallel incubations were averaged, and all parameters were then z-transformed (Freedman et al. 2007) to standardize the data set.The values of permeability were log-transformed prior to standardization to account for the wide range of magnitude.The principal components were then rotated orthogonally employing the varimax method.The PLS scores of the principal components were further employed for a clustering analysis.
A two-way ANOVA was performed with Matlab (version R2019b) to test site characteristics and measured fluxes for significant seasonal and spatial differences.For a complete factorial analysis, only measurements from cruises HE-395/432/447/471 and AL-520 were included.From these cruises, data from stations C, D, E, F, H, and I were used.The excluded cruises and stations did not provide a complete factorial data set.
For the examination of the chlorophyll concentration of coastal and open North Sea sites, we performed the Student's t-test with the Matlab software (version R2019b).
The complete data of fluxes and sediment characteristics are available via the Pangaea data repository.

Estimation of the relative temperature effect
The relative change of respiration rates (F T ) due to a change in temperature (T) with respect to the reference temperature (T Ref ) was calculated with Eq. 4.
where F T = relative temperature change, T = temperature, T Ref = reference temperature, and Q 10 = temperature coefficient.
The equation describes that a respiration rate changes by the factor Q 10 for every 10 C temperature change.Here, we assumed a temperature coefficient of 2, which is a typical value for many benthic systems (Provoost et al. 2013;Wrede et al. 2018).

Results
In the following section, we focus on general trends and present results to highlight the interplay of benthic fluxes and independent site characteristics such as sediment type, benthic community, and temperature.The synoptic Fig. 3 provides an overview, and the basic data of all relevant parameters are available as Supporting Information.

Seasonal and spatial variability of site characteristics
Each of the examined site parameters (temperature, chlorophyll, and bioturbation potential) had an individual seasonal pattern (Fig. 3A-D).In the German Bight, the observed bottom water temperatures had a pronounced seasonal cycle with a minimum of approximately 5 C in early spring (March) and a maximum of approximately 20 C in late summer (September).During any season, the spatial differences in bottom water temperatures were in the range of 5 C (Fig. 3A).
The bottom water chlorophyll concentration had a clear maximum during the spring bloom in May (Fig. 3B), which was more pronounced at the coastal stations (A-E) than at the distal stations (F-I).The highest bottom water chlorophyll was measured at site NOAH-C during a spring bloom (41.6 μg L −1 ), which is one order of magnitude above the average concentration at that time (approximately 6 μg L −1 ).After the peak in Table 3. Spearman's rank correlation coefficients of site characteristics water depth, salinity, temperature, permeability, sediment chlorophyll (Sediment chl.), bottom water chlorophyll (Bot.water chl.), bioturbation potential (bioturbation), bioirrigation potential (bioturbation), and core term.Significant correlations (α = 0.01) are indicated with an asterisk.May, bottom water chlorophyll remained at a low plateau during summer, and decreased back to minimum values between November and February (Fig. 3B).For the sediment chlorophyll data, we found a significant seasonality (ANOVA, df = 2, F = 11.3,p < 0.001) and spatial variability (ANOVA, df = 5, F = 78.5, p < 0.001).However, the seasonal pattern of sediment chlorophyll was less clear than the pattern of bottom water chlorophyll (Fig. 3C).Yet the monthly median values, which are more robust to the present combination of low relative variability and high spatial variability, were higher by a factor of 2 in summer (May-September) compared to winter (January-March) (ANOVA, df = 1, F = 38.1,p < 0.001).The highest concentrations were measured in fine-grained sediment of which the NOAH-C site had the highest concentrations of up to 13.2 μg chlorophyll per gram of dry sediment.Concentrations in coarser, sandy sediment (A, B, E, I, CCP-G) were in the range of 0.1-4.1 μg g −1 .
Both macrofauna biomass and abundance had substantial seasonal variability with a minimum during winter and spring (January-May) followed by a continuous increase during summer to a peak in August and September (see Supporting Information).Consequently, all three biotransport indices that were calculated from abundance and biomass reflected the same general seasonality: low values during winter and spring (January-May), and a rapid increase during summer (June-September).This is also evident in Table 3 in the high correlation of all three indices (Spearman's ρ = 0.74-0.95).As the seasonal variability was more pronounced in fine-grained sediment than in coarse sediment, we normalized the values of bioturbation potential in Fig. 3D with June values as references (June = 1).
In contrast to the above site characteristics, salinity and sediment permeability had no substantial seasonal pattern.However, we observed substantial spatial differences in sediment permeability: Sites in deep water with fine sediment had lower permeability (e.g., sites F, H) than sites with coarse sediment along the coast.Coarser sediments generally had a permeability > 10 −12 m 2 and was as high as 10 −10 m 2 at sites A and E (see Supporting Information).The lowest permeability of 10 −15 m 2 was measured in sediment from the silt province at site C. Sampling stations closer to the coast (A-E, CCP-G) had low salinities in the range of 31.5-33.5PSU, while the more distant stations (F-I) had slightly higher salinities in the range of 34-35 PSU (see Supporting Information).

Seasonal and spatial variability of benthic fluxes
Oxygen and nutrients fluxes had a similar seasonal pattern (Fig. 3E-H).Oxygen fluxes were significantly correlated with DIN fluxes (Spearman's ρ = 0.87, p < 0.001), silicate (ρ = 0.82, p < 0.001), and phosphate (ρ = 0.60, p < 0.001).The variation of flux magnitudes between NOAH sites was remarkably low as compared to the substantial seasonal variability.Fluxes were low during winter (January-March) with oxygen fluxes typically of order −6 mmol m −2 d −1 .Fluxes of silicate and DIN were very low (< 0.5 mmol Si m −2 d −1 , < 0.2 mmol N m −2 d −1 ).Starting in late spring (May), fluxes increased substantially, peaked in late summer (August, September), while the peak of the median flux magnitudes was in August and declined to September.A Student's t-test suggests that the coastal sites A-D, CCPG had significantly higher fluxes than the open German Bight sites E-I, CCPJ (t-test, p < 0.05, N = 9).Maximum fluxes observed in late summer were −70 mmol O 2 m −2 d −1 , 10 mmol Si m −2 d −1 , and 6 mmol DIN m −2 d −1 , respectively.Phosphate fluxes remained generally low with typically < 0.2 mmol P m −2 d −1 as a consequence of oxygen being above 80% saturation in the core supernatant, and also displayed a seasonality pattern similar to the patterns of oxygen and silicate (Fig. 3H).In the course of autumn, benthic fluxes decreased rapidly and remained low during fall and winter (Fig. 3E-G).
In order to compare our oxygen fluxes that are predominantly sustained by faunal activity with the results on oxygen fluxes driven by pore water advection by Ahmerkamp et al. (2017), we calculated average oxygen fluxes for each sediment province.Based on our data from the joint Heincke cruises HE-432 and HE-447, the average fauna-sustained oxygen fluxes were The weighted average (by province area) of the fauna-sustained oxygen flux was −12.3 ± 2.0 mmol O 2 m −2 d −1 .The sediment provinces with very fine sand and very coarse sand, that together represent less than 11% of the study area (Fig. 1), have not been sampled during these two cruises.
We calculated C stock to oxygen flux ratios based on benthic and pelagic chlorophyll concentrations and benthic oxygen fluxes.These ratios indicate the hypothetical periods during which the chlorophyll-associated carbon stock can sustain the observed oxygen fluxes.The benthic C stock to flux ratios had a substantial variation among sediment types most of the time.Especially site C had a stock to flux ratio of > 100 d, while the other sites generally had benthic C stock to flux ratios in the range of 20-60 d.During the peak oxygen consumption in summer, the benthic C stock to oxygen flux ratios dropped to 5-10 d (Fig. 4A).By contrast, the pelagic C stock to oxygen flux ratios had substantially less spatial variation, and the pelagic C stock to flux ratios were generally in the range of 5-10 d.During the spring plankton bloom in April and May, pelagic C stock to oxygen flux ratios increased to 30-70 d, and exceeded 100 d at sites A and C (Fig. 4B).

Partial least squares analysis
We employed a PLS analysis to separate and quantify the specific effect of each site parameter on the observed benthic fluxes.However, pair-wise calculated Spearman rank correlation coefficients show that not all site characteristics were perfectly independent (Table 3).For examples, we found moderate correlations for salinity and water depth (Spearman's ρ = 0.63, p < 0.01), and permeability and sediment chlorophyll content (Spearman's ρ = −0.61,p < 0.01).The community potentials for bioturbation and bioirrigation, and the trait-less core term were strongly correlated (Spearman's ρ = 0.74-0.95).The calculated variance inflation factors for the whole ensemble were < 4 for all parameters except bioturbation potential and core term, which had variance inflation factors > 20.Since factor loadings of the bioturbation potential (component 1: 4.4, component 2: 2.7) were similar to those of the core term, and were substantially higher than those of the bioirrigation potential (component 1: 0.9, component 2: 1.1), we subsequently only show results for the bioirrigation potential.
The PLS component 1 is mainly loaded by temperature and bioturbation potential and explains 63% of variability of oxygen fluxes, 48% of silicate, 36% of phosphate, and 39% of DIN fluxes, respectively.The PLS component 2 additionally explains approximately 6% of variability in observed fluxes.It is mainly positively loaded by salinity and water depth, and negatively loaded by the sediment chlorophyll content.The loadings of the first two components are visualized in Fig. 5.The loadings for the flux parameters oxygen, silicate, phosphate, and DIN are very similar and are comprised by gray circles in Fig. 5. Component 3 is loaded by sediment permeability, sediment chlorophyll content, and the bottom water chlorophyll concentration, and explains only approximately 2% of the observed flux variability.
In summary, we found the highest degree of similarity between the seasonal patterns of benthic fluxes, and the biotransport indices (bioturbation potential, bioirrigation potential, core term).This is also reflected by the dominant loading of the biotransport indices onto the first PLS component (Fig. 5).The second dominant parameter that loads onto the first PLS component is bottom water temperature, as it peaked simultaneously with the biotransport indices.However, temperature decreased gradually over the course of fall and winter, while the biotransport apparently decreased faster: By the beginning of winter (January), the biotransport indices and benthic fluxes were already close to the minimum values, while the water temperature continued to decrease until March (Fig. 3).
For further disentangling the effects of temperature and biotransport, we estimated the relative changes in temperature and bioturbation potential, and compared these with the  relative change in observed oxygen fluxes.By using Eq. 4 and the median temperature from March (5.5 C) as the reference, we estimated an increase of respiration rates by the factor of 2.4 in September (18 C).By further assuming a linear relationship of biotransport proxy and actual transport rates (Wrede et al. 2018) and using the median of normalized bioturbation potential from March (0.6), we estimated a multiplication of the biotransport rates by 4.1 in September.Applying the combined factor of 6.5 (2.4 + 4.1) to the median oxygen fluxes of −5 mmol m −2 d −1 as measured in March (Fig. 3E), we estimated a theoretical oxygen flux of approximately −33 mmol m −2 d −1 for September (compare with Fig. 3E).
The Scores of PLS components 1 and 2 for the sampled sites were utilized for a cluster analysis.Based on our data, stations assigned to the provinces "coarse sand," "medium sand," and "fine sand" are indistinguishable.By contrast, samples from the "silt" province with very fine sediment (Fig. 6) are significantly different and are grouped in a separate cluster.

Comparison of oxygen consumption by biotransport and advective pumping
In parallel to our incubations, Ahmerkamp et al. (2017) employed an in situ technique during two joint cruises (HE-432, HE-447) for the measurement of oxygen fluxes that are associated with pore water advection.As faunal activity and pore-water advection represent the two major drivers of benthic exchange in the North Sea, we assume that the sum of both processes approximates the total benthic oxygen flux.During June and September, total oxygen consumption rates in coarse sand amounted to 33.6 ± 12.0 mmol O 2 m −2 d −1 (Fig. 7 ) with biotransport contributing approximately 32%.The contribution of biotransport increased toward finer sediments, while the flux sustained by pore-water advection decreased.The total oxygen consumption was 28.8 ± 2.4 mmol O 2 m −2 d −1 in medium sand and 16.8 ± 4.8 mmol O 2 m −2 d −1 in fine sand, to which biotransport contributed 44% and 86%, respectively (Fig. 7).Sediment from the silt province has not been sampled by Ahmerkamp et al. (2017), but since the permeability of silt sediment is very low, we assume that the advectionsustained oxygen flux was negligible.

Discussion
Identifying the drivers of benthic fluxes is essential for understanding the seasonal dynamics of oxygen and nutrients exchange across the seafloor, which includes the evaluation of macrofaunal effects on the one hand, and the effects of porewater advection evoked by flow over rippled beds on the other hand.Therefore, we sampled all important sediment provinces of the German Bight repeatedly during all seasons from 2013 to 2019, and compiled a unique data set of physical,  biogeochemical, and macrofaunal parameters.This facilitated the analysis of the seasonal variability of benthic fluxes connected to benthic fauna data, and further identified the major driving factors.
Our PLS analysis suggests that bioturbation potential and temperature had the highest explanatory power, and explained much of the substantial seasonal variability in observed benthic fluxes of oxygen and nutrients.In contrast, the site characteristics water depth, sediment permeability, and sediment chlorophyll content had considerably lower explanatory power.Their combined effects account for the relatively low spatial variability of observed fluxes.
Temperature and the bioturbation potential loaded approximately equally on PLS component 1 (Fig. 5), but this ratio does not necessarily reflect the absolute contribution of the two parameters to the observed variability in fluxes, as the values were standardized prior to PLS analysis.Temperature is a significant driver of reaction rates in general, including the metabolism of benthic organisms (Provoost et al. 2013;Wrede et al. 2018).For further disentangling the effects of temperature and biotransport, we estimated the hypothetical increase of the benthic oxygen consumption from March to September due to the relative change in temperature and bioturbation potential.The theoretical estimate is close to the 3 rd quartile of actually observed oxygen fluxes (Fig. 3E).With this agreement of extrapolated and observed oxygen consumption rates, we conclude that temperature accounts for approximately 1/3 of the variability explained by PLS component 1, while the biotransport proxy accounts for approximately 2/3.As for the absolute variability in oxygen fluxes of which 63% are explained by PLS component 1, approximately 40% can be attributed to biotransport, and 20% can be attributed to temperature.

Biotransport proxies
The biotransport proxies emerged as a dominant parameter to explain seasonal variability.We calculated the following indexes (1) community bioturbation potential, (2) community bioirrigation potential, and (3) the biomass-abundance product (core term).
We initially assumed that a detailed, trait-based implementation of biotransport by utilizing the potentials for bioturbation and bioirrigation would substantially improve estimates of benthic fluxes in benthic models.Therefore, we analyzed the predictive power of both trait-based proxies in comparison with the biomass-abundance product (core term) that ignores any species-specific traits.All three biotransport proxies were significantly correlated (Table 3).Although the bioirrigation potential had a variance inflation factor < 4, which would have made it plausible to include it in the PLS model, the factor loadings were much lower than the loadings of the bioturbation potential.Additionally, the trait-based bioturbation potential and the trait-less core term were highly correlated (Table 3) and had variance inflation factors > 20, which implies that both proxies are equally well suited for the prediction of benthic fluxes.This conclusion is counterintuitive, as the macrofauna communities of the sediment provinces were substantially different with respect to species composition and feeding strategies, such as surface feeding, suspension feeding, or sediment regeneration (Meyer et al. 2018(Meyer et al. , 2019;;Meyer and Kröncke 2019).The high degree of correlation among the different biotransport proxies is remarkable.We suggest that the specific traits of individual species in sufficiently diverse communities complement each other with the consequence that abundance and biomass alone may be sufficient to predict the bioturbation/bioirrigation potential with minor contributions of the specific traits.However, an in-depth discussion of this result is beyond the scope of this study.
A potential bias might have been introduced into our data by using a Multicorer, which is not a typical sampling device for macrofauna surveys.However, using Multicores enabled us to correlate the nutrient and oxygen fluxes directly to the fauna within the cores.Sediment cores retrieved by the Multicorer were typically 15-30 cm in length.This is where the majority (80-90%) of the fauna has been previously observed in the German Bight (Flach and Heip 1996;Oehler et al. 2015a,b).In agreement with the results, we found the highest abundance and biomass in the upper 5 cm sediment layer, and both were decreasing steeply downward.From 10 to 25 cm depth, abundance and biomass were close to zero (Dauwe et al. 1998;Oehler et al. 2015a,b;Zhang et al. 2019).Higher abundances of deep-burrowing species such as Callianassa subterranea, burrowing up to 90 cm (Witbaard and Duineveld 1989;Dauwe et al. 1998;Rowden et al. 1998), or Upogebia deltaura, generally burrowing up to 30 cm (Tunberg 1986;Dauwe et al. 1998), are mainly restricted to areas with muddy sand in the central southern North Sea (Witbaard and Duineveld 1989).We conclude that the calculated biotransport indexes are applicable to the study area.

Synthesis: Benthic remineralization dynamics
After establishing the seasonal variability of benthic fluxes, identifying driving parameters, and analyzing the significance and suitability of the applied biotransport proxies, we combine these results to evaluate the seasonal dynamics of benthic processes.For the interpretation of our results, we introduced the carbon to oxygen flux ratio.The ratio of the depthintegrated organic carbon in water column or sediment, and the benthic oxygen flux, directly indicates how long the instantaneous benthic oxygen consumption can be sustained by the available carbon stock.
From the repeated surveys in a 7 yr period, we have synthesized an average year.Beginning in January, the benthic C stock to oxygen flux ratios were high (Fig. 4A) and exceeded 100 d in fine sediment (sites C, F, I), while the pelagic C stock to oxygen flux ratios were substantially lower and were in the range of 3-5 d.This suggests that the low remineralization rates were not carbon-limited and that the benthic organic matter was probably the major food source.The benthic carbon stock was gradually consumed over the winter (January-March) as indicated by the decrease in sediment chlorophyll concentrations (Fig. 3C).During the spring plankton bloom in April and May, water column chlorophyll attained the annual maximum, and the sediment chlorophyll concentrations increased concomitantly.Our samples from May probably represent a late stage of the spring bloom, as the peak of spring bloom in the southern North Sea is commonly in April (Sharples et al. 2006).However, biomass and abundance of the macrofauna were still low, and the resulting faunasustained fluxes were likewise low (Fig. 3E-H).The corresponding pelagic C stock to oxygen flux ratio peaked at typical values in the range of 30-70 d, and exceeded 100 d at sites A and C. The abundant plankton biomass available from the spring bloom sinking to the sediment is exceeding the benthos's capacity to consume it.Moreover, many larval stages of benthic organisms are part of the plankton bloom, and have to further develop and settle on the seafloor to become part of the benthic community later.Therefore, an increase in bioturbation activity due to secondary production was not detected prior to July.Because of the delayed benthic response, chlorophyll and chlorophyll-associated carbon continued to accumulate in the sediment until August, which concurs with previous studies (Boon and Duineveld 1998).After the spring bloom, phytoplankton biomass decreased but productivity remained on a high level (Joint & Pomroy 1993).As pelagic nutrients were depleted by the spring bloom, the continuing productivity relied on an efficient regeneration of nutrients (Bratek et al. 2020).The phytoplankton biomass remained at relatively low levels as reflected by the decrease in bottom water chlorophyll concentrations.During summer, the benthic remineralization rates continued to increase, while the chlorophyll-associated carbon stock of the sediment declined after a peak in August.The benthic macrofauna community continued to grow until the late summer, and in pace with the macrofauna growth, benthic remineralization rates increased until approximately August and September.Due to the decreasing resupply from the water column and the increasing consumption within the sediment, the average pelagic carbon stock to oxygen flux ratios gradually decreased to just 5-15 d.This suggests that the benthic community might have been carbon-limited during its peak turnover rate, especially in coarse sediment.The benthic carbon to oxygen flux ratios were particularly low at coarse sediment sites A (5 d) and CCP-G (12 d).Additionally, microbial respiration sustained by pore-water advection was in the same magnitude as the fauna-sustained fluxes at these sites (Ahmerkamp et al. 2017), which further reduced the carbon to oxygen flux ratio to just a few days.This highlights that the benthic community depends on a rapid resupply of organic matter from the pelagic primary production during late summer.Yet, not all of the pelagic carbon stock was available for the macrobenthos.On average, about 90% of the organic matter turnover takes place in the water column of the German Bight (Heip et al. 1995).
The fall season (October-December) is covered by a relatively low number of observations due to frequently bad weather conditions and related difficulties to carry out surveys.In early fall, water temperatures were still high, but observations indicate a rapid decline in biomass and abundance that approached the minimum values in January (Fig. 3D).Likewise, macrofauna-sustained fluxes were already at minimum levels at late fall/early winter (Fig. 3E-H).This observation is contrary to results of Saulnier et al. (2019), who presented a multidecadal time series of benthic biomass seasonality from the southern North Sea, with a gradual population decrease over the course of 5 months.However, benthos abundance data from the Dutch Wadden Sea and German Bight illustrates that the macrozoobenthos can collapse to approximately 1/3 of the peak abundance within 1-2 months in individual years (Beukema 1974;Reiss and Kröncke 2005).On the ecosystem scale, observations from the nearby Wadden Sea during September show a fast transition from a system that takes up nutrients to a nutrient-releasing system (Loebl and van Beusekom 2008).Loebl and van Beusekom (2008) suggested that the collapse of the productivity during September was due to the combined effect of increased wind events and decreasing insolation.Taken together, our limited number of observations from the fall season suggest that the macrozoobenthos and hence macrofauna-sustained fluxes rapidly decline during every fall in German Bight sediments, but more observations are needed to corroborate this.The scarcity of German Bight benthos data from fall and winter, in combination with potentially rapid transition processes, call for dedicated sampling campaigns that specifically target the fall and winter seasons.

Significance of benthic fluxes sustained by physical and biological transport
In parallel to our incubations, Ahmerkamp et al. ( 2017) measured the oxygen fluxes evoked by pore-water advection.The comparison of the results of Ahmerkamp et al. (2017) with ours enables us to evaluate the relative contribution of the benthic macrofauna to the overall exchange of oxygen between bottom water and sediment.Figure 7 illustrates the gradual transition from the dominance of biotransport in impermeable sediment to the dominance of pore-water advection in coarse and permeable sediment.Interestingly, the total oxygen fluxes were quite similar among the different sediment types, even though the corresponding sediment characteristics varied considerably between impermeable and permeable sediment.Based on the combined results of our study and Ahmerkamp et al. (2017), we estimate that the weighted average (by province area) of the total oxygen consumption for the whole study area in summer, was approximately 23 ± 5 mmol O 2 m −2 d −1 , to which faunal transport contributed 54%.
The carbon stock to oxygen flux ratio introduced above offers an additional perspective.We estimated the amount of chlorophyll-associated carbon in the water column, and integrated the pelagic carbon in the lower half of the water column to account for the stable thermocline during summer.Pelagic C stock to oxygen flux ratios representing carbon decay times of 3-5 d for macrofauna-sustained fluxes do not necessarily indicate that the sediment is a major sink for the pelagic primary production.In fact, approximately 90% of the pelagic primary production is consumed within the pelagic community (Heip et al. 1995).Nonetheless, the benthic community is not only consuming oxygen to degrade organic matter but is also recycling the associated nutrients back into the water column.When assuming a steady state, this implies that the benthic nutrient fluxes hypothetically could provide nutrients for 1/3 (33%) to 1/5 (20%) of the daily subthermocline primary production in August, which is supported by the observation of deep chlorophyll maxima in the North Sea during summer stratification (Fernand et al. 2013).These approximations agree with estimates by Heip et al. (1995), for the whole water column (10%), and is within the range of recent estimates of Bratek et al. (2020).They concluded that benthic nitrogen recycling contributed 7-59% to the subthermocline primary production in the German Bight in late summer.

Conclusions
The major sediment provinces of the North Sea based on grain size classification were repeatedly sampled over 7 yr from 2013 onward at nine locations to analyze the seasonal and spatial variability of fauna-sustained fluxes of oxygen and nutrients between sediment and bottom water.From this unique data set, we conclude: 1.The intensity of fauna-sustained bentho-pelagic exchange is inversely related to the intensity of fluxes sustained by the advective transport evoked by flow over a rippled bed.Biotransport is most intense in fine-grained sediment where pore-water advection is virtually absent.Weighted by the respective area of each sediment province, we estimate an average total benthic oxygen consumption of 23 ± 5 mmol O 2 m −2 d −1 during summer, to which biotransport contributed approximately 54%.The parallel recycling of nutrients from the sediment back into the water column was hypothetically sufficient to fuel 20-33% of the daily subthermocline primary production in August.2. Abundance and biomass of the macrozoobenthos, expressed as trait-based bioturbation potential or trait-less core term, in combination with bottom water temperature, were the dominant drivers of fauna-sustained fluxes that have clear seasonality.The spring plankton bloom provided the benthos with organic matter.Starting from low levels of abundance and biomass, the macrozoobenthos community started to grow in response to the peak of primary production.At this stage, the majority of the biomass flux toward the sediment was temporally accumulating in the sediment until benthic turnover exceeded the biomass supply from the water column.Macrobenthos abundance and biomass reached a maximum in late summer, and started to decline when pelagic primary production and temporally stored benthic biomass did no longer sustain the benthic community.3. The trait-based community bioturbation potential and the trait-less core term (ct) were significantly correlated, and were equally well suited for the prediction of benthic fluxes evoked by macrofauna.Along with temperature, biotransport indexes emerged as the parameters with the strongest explanation power.Their combined effect accounted for up to 63% of the observed variability in benthic fluxes with temperature accounting for approximately 1/3 of the explained variability and biotransport indexes for approximately 2/3 of fluxes, respectively.4. Site characteristics such as grain size, permeability, and bottom water chlorophyll concentration have only minor contributions to the observed variance in fauna-sustained benthic fluxes of oxygen and nutrients.Biotransport partially compensates for lower advective exchange in finegrained sediments (Fig. 7) and thereby masks the effects of grain size, permeability, and current velocity on total oxygen and nutrient fluxes.

Fig. 1 .
Fig. 1.Map of the study area in the German Bight of the southern North Sea.Colors indicate distinct sediment provinces; red circles indicate sampled sites (NOAH A-I, CCP-G/J).The solid, black outline marks the German Exclusive Economic Zone.Sediment classification by Walter Puls, Helmholtz Center Geesthacht (unpublished).

Fig. 4 .
Fig. 4. Carbon stock to oxygen flux ratios of (A) chlorophyll-associated benthic carbon and (B) chlorophyll-associated pelagic carbon.The color bars indicate the spring equinox (blue), fall equinox (red), and the spring plankton bloom (green).Boxes indicate 2 nd and 3 rd quartiles, median is indicated by a red line in the boxes, and extreme values are indicated by crosses.Note that y-axes are log-scaled.

Fig. 5 .
Fig. 5. Biplot of the first two principal components of partial least squares regression.The gray circles comprise the loadings of individual parameters for the fluxes of oxygen, silicate, phosphate, and DIN.Colors and numbers in brackets of axis labels indicate the vectors and explained variance for oxygen (blue), silicate (red), phosphate (black), and DIN (orange).

Fig. 6 .
Fig. 6.PLS scores for the oxygen model of principal components 1 and 2 for individual sites to illustrate the clustering of the NOAH sediment provinces.N = 37.The PLS component 1 explains 63% of variance in oxygen fluxes; PLS component 2 explains 6% of variance in oxygen fluxes.

Fig. 7 .
Fig. 7. Averages of benthic oxygen fluxes into the sediment at four sediment provinces that are sustained by fauna (gray) and flow (white) during cruises HE-432 and HE-447 (N = 3 sites per province).Data on faunasustained oxygen consumption rates from this study, flow-sustained oxygen consumption data from Ahmerkamp et al. (2017).

Table 1 .
Summary of sediment characteristics: median of grain size distribution (D 50 ), intrinsic permeability, and TOC content of the surface sediment.The assigned sediment province refers to Fig. 1.

Table 2 .
Summary of sampling campaigns of the vessels Alkor (AL), Heincke (HE), and Maria S. Merian (MSM) with sampled sites and number of incubated sediment cores used for calculations of monthly averages.