Modelling biogeochemical processes in sediments from the north western Adriatic Sea : response to enhanced POC fluxes

This work presents the result of a study carried out in the north-western Adriatic Sea, by combining two different types of biogeochemical models with field sampling efforts. A longline mussel farm was taken as a local source of perturbation to the natural POC downward flux. This flux was first quantified by means of a pelagic model of POC deposition coupled to sediment traps data, and its effects on sediment bioirrigation capacity and OM degradation pathways were investigated by constraining an early diagenesis model, linked to new data in sediment porewaters. The measurements 15 were performed at stations located inside and outside the area affected by mussel farm deposition. Model-predicted POC fluxes showed marked spatial and temporal variability, which were mostly associated with the dynamics of the farming cycle. Sediment traps data at the two sampled stations (in and out of the mussel farm) showed average POC background flux of 20.0-24.2 mmol C m -2 d -1 . The difference of OC fluxes between the two stations was in agreement with model results, ranging between 3.3 and 14.2 mmol C m -2 d -1 , and primarily associated with mussel physiological conditions. Although 20 restricted, these changes in POC fluxes induced visible effects on sediment biogeochemistry. Observed oxygen microprofiles presented a 50 % decrease in oxygen penetration depth (from 2.3 to 1.4 mm), accompanied by an increase in the O2 influx at the station below the mussel farm (19-31 versus 10-12 mmol O2 m -2 d -1 ) characterized by higher POC flux. DIC and NH4 + concentrations had similar behavior, with a more evident effect of bioirrigation underneath the farm. This was confirmed through constraining the early diagenesis model, which calibration leads to an estimation of enhanced and shallower 25 bioirrigation underneath the farm: bioirrigation rates of 40 y -1 and irrigation depth of 15 cm were estimated inside the shellfish deposition footprint versus 20 y -1 and 20 cm outside. These findings were confirmed by independent data on macrofauna composition collected at the study site. Early diagenesis model results indicated a larger organic matter mineralization below the mussel farm (11.1 versus 18.7 mmol m -2 d -1 ), characterized by similar proportions between oxic and anoxic degradation pathways at the two stations, with an increase in the absolute values of oxygen consumed by OM 30 degradation and reduced substances re-oxidation underneath the mussel farm.


Introduction
Disturbance gradients in benthic habitats of marine soft sediment can cause shifts in species behaviours and interactions, thus affecting the biodiversity-ecosystem function relationship (Snelgrove et al., 2014;Villnäs et al., 2013).Disturbance gradients can also affect spatial habitat heterogeneity, determining changes in essential ecosystem traits such as food web functioning Shellfish farming is regarded as an extractive aquaculture activity (Barrington et al., 2009).However, the production of faeces and pseudofaeces leads to a net transfer of organic matter from the water column to the surface sediment (Tenore and Dunstan, 1973;Cranford et al., 2007).This process is expected to locally affect sediment biogeochemistry, benthic-pelagic coupling, and benthic community functioning.A range of studies peformed in the last 30 years reported on farm-induced changes in sedimentation rates (Callier et al., 2006), sulfate reduction (Dahlbäck and Gunnarson, 1981), NH 4 + and PO 4 2-5 regeneration (Hatcher et al., 1994;Nizzoli et al., 2005), porewater nutrient concentration gradients (Mesnage et al., 2007), benthic community structure (Stenon-Dozey et al., 1999;Mirto et al., 2000;Christensen et al., 2003).Biogeochemical and contaminant cycling associated to the presence of shellfish farming in a Mediterranean lagoon was the focus of a special issue (Estuaine Coastal and Shelf Science 72(3); Rabouille et al., 2007).Only recently, the quantitative understanding of these processes have received attention in relation to the assessment of "positive modifications" induced by shellfish farm 10 deposition on benthic habitats (McKindsey et al., 2011).The influence of local hydrodynamics on the fluxes of organic matter deposited by mussel farms was the focus of two modelling studies, by Hartstein and Stevens (2005) and Weise et al. (2009).Based on these works it was possible to have a clearer mechanistic understanding of the relationship between the degree of deposition and the different farming conditions (in terms of local hydrodynamics and farm characteristics -depth; geometry).To our knowledge, less attention was received by the pelagic-benthic coupling associated to the different phases 15 of the rearing cycle, and, ultimately linked to the physiology of the farmed mussel -faeces and pseudofaeces production rates are dependent on water temperature and seston quantity and quality (Tenore and Dunstan, 1973), parameters which could present remarkable variations at the annual time scale, in particular at those sites characterized by a relatively fast growth out cycle -e.g. the Mediterranean Sea.
In this work, a longline mussel farm located in the north western Adriatic Sea was regarded as a local source of perturbation 20 of natural organic matter downward fluxes.Average yearly increase in POC flux induced by the farm was first quantified by applying a biogeochemical model of POC production and deposition (mussel faeces and pseudofaeces), coupled to sediment traps data, and subsequently used to force an early-diagenesis model simulation, which was calibrated by collecting new data in sediment porewaters.The objectives of this study are a better understanding of seasonal and inter-annual variability in POC deposition fluxes from the mussel farm, and the quantification of benthic recycling of organic matter under contrasted 25 forcings linked to mussel farm, ecological mechanisms associated (bioirrigation), and relative importance of oxic versus anoxic processes.

Materials and methods
The study was articulated in four steps: 1) a pelagic biogeochemical model was used, with the aim of estimating the typical flux of POC originating from the 30 mussel farm throughout the year; 2) two sediment trap deployments were carried out at the beginning and at the end of mussels farming cycles, with the aim of corroborating model predictions; 3) sediment cores were collected at two stations, a pristine one, and an impacted one, which were localized on the basis of model application.Measurements included O 2 micro-profiling, porosity and micro-porosity, pore waters NH 4 + , SO 4 2-

35
, and Dissolved Inorganic Carbon (DIC); 4) two early diagenesis models (EDM) were set-up and constrained by the observed field data in the sampled cores at the two stations: bioirrigation parameters and ratio among degradation pathways were estimated on the basis of model application.

Study site and mussel farm description
The study was performed at a Mediterranean mussel (Mytilus galloprovincialis) farm located approximately 1.5 nautical miles off-shore the city of Jesolo (Italy), in the North-Western Adriatic Sea (see Figure 1).The area is characterized by a flat bathymetry ranging between 13 and 15 m (Colla, 2017).The farm is affected by the freshwater plume of the Sile River, which outlet is located at approximately 1.5 nm of distance from the south-western edge of the farm (see Figure 1).The 5 mean annual Sile river discharge in 2008-2009 was 10.9 m 3 s -1 , ranging between 5.2 m 3 s -1 in March 2008, and 14.7 m 3 s -1 in December 2009 (ARPAV, 2010).A grain size analysis was recently carried out in the area: sediments were classified as silty-sandy, with percentages of sand and mud respectively within the ranges 11.9- 25.5 % and 88.1-74.5 % (Colla, 2017).
Offshore mussel farming is extensively practiced in the northern Adriatic Sea, and farmed mussels account for approximately 2/3 of the national Italian production (MiPAAF, 2014).The studied mussel farm covers an area of about 2 10 km 2 and has been operating since 1990 with an average annual production ranging between 600 and 800 tons year -1 (farm manager pers.comm.).Mussels are grown on ropes approximately 4 m long, which are suspended on cables, and placed at depths comprised between 2 and 4 m.Lines are positioned parallel to the coast, along the principal current direction, at a distance of 40 m among each other.Length of each line is approximately 2 km.Mussels are normally harvested within July-September, after a rearing cycle lasting a single year, during which they are re-socked 2-3 times (farm manager pers.15 comm.).

Model theory: mussel farm deposition model and Early diagenesis model
For the purposes of this work, the model described by Brigolin et al. (2014) was modified to simulate the C, N, P biogeochemical fluxes across shellfish farms.For an introductive description of the model, the reader should refer to the original publication, here we will focus on those changes required to adapt the model to the simulation of mussel farms 20 deposition.The model (see Figure 2) combines two generic modules, respectively accounting for: i) individual growth and dynamics of the farmed population; ii) organic particle tracking and deposition.
Mussel growth and population dynamics were estimated by means of the individual-based approach by Brigolin et al. (2009).This individual model is capable to simulate physiological processes and their response to key environmental forcings, i.e. suspended particulate matter quality/quantity and water temperature.The individual growth model requires in 25 input time series of daily values of sea water temperature, and concentrations of chlorophyll-a, particulate organic carbon (POC), and total suspended solids (TSS).The model allows accounting for the daily energy intake, weight gain, and faeces and pseudofaeces production rate, these latter two representing the input for the deposition module.The individual model was up-scaled to the population level by means of a set of Monte Carlo simulations, which were used for estimating the size structure of the population (the virtual population was made up of 5000 individuals).In accordance with Bacher and 30 Gangnery (2006) such differences were accounted for by assigning to each specimen a different maximum clearance rate, reflecting variability in individual phenotypes as well as differences in the localization of specimens within the farm.The model allows to specify farm geometry based on a set of basic parameters.An interval is selected in order to account for the depth at which ropes are positioned, z fmin -z fmax , and initial position of mussel biodeposits (faeces and pseudofaeces) is assigned randomly within this interval.Longlines -cables on which mussel ropes containing mussels are attached -are 35 disposed parallel each to the other, at a fixed distance, S x , and with a defined orientation with respect to the North, D n .
Mussels are considered to have an homogeneous density, B i , within the same i-th longline.A fixed number of longlines, n L , of length l L , are assumed to be productive within a farming cycle.The bathymetry of the farmed area can either be specified through an external file, or assumed to be flat.Faeces and pseudofaeces deposition from mussel lines is simulated by means of a Lagrangian technique, consistent with the advection-diffusion equation (for details see Jusup et al., 2007) 4) in Brigolin et al. (2014).These fluxes are equally partitioned among C, N and P content to each generic particle.Particles are released homogeneously within the 24 hours, with a total of 5000 particles launched every day from each mussel line.The output of this module are daily maps of downward fluxes of organic C, N and P reaching the sea bed (in g m -2 d -1 ).The complete set of parameters used in the deposition model, values and their references, are reported in 5 Table A1 (Appendix).Settling velocities for mussel faeces and pseudofaeces were set to 1.0 ± 0.1 cm s -1 and 0.1 ± 0.01 cm s -1 (see Weise et al., 2009;Chamberlain, 2002).The settling velocity of each particle was randomly selected from a Gaussian distribution.
The early diagenesis model (EDM) was developed by means of the Biogeochemical Reaction Network Simulator -BRNS - (Regnier et al., 2002), through the Knowledge-Based Reactive Transport Model application (Aguilera et al., 2005).In the 10 present application we used a simplified version of the EDM identified for the northern Adriatic Sea by Brigolin et al. (2011).The model solves the diagenetic equations describing mass conservation for solids and dissolved species in a vertical sediment column -see Eqs. ( 1) and (2) in Table A2 (Berner, 1980;Boudreau, 1997).The advection term includes burial, compaction, and bioirrigation; the diffusion term includes molecular and ionic diffusion, as well as bioturbation.Organic matter oxidation is described by means of a multi-G model (Westrich and Berner, 1984) with three types of organic matter: 15 labile (OM1), semi-refractory (OM2), and mussel biodeposits (OM3), following the approach by van Cappellen and Wang (1996).Oxic and anoxic pathways of organic matter oxidation, as well as secondary redox reactions, are included (see Table A2).As said, the reaction network is simpler with respect to the one published in Brigolin et al. (2011), not including reactions involving Fe and Mn, which processes were not the focus of the present work.A fixed concentration was imposed at the upper boundary for all solutes, while a fixed flux was used for solids.The model was coded in Fortran.The ordinary 20 differential equations were integrated numerically by means of a 4 th order Runge-Kutta scheme (Press et al., 1987).The Lagrangian equation for the deposition model was solved following Jusup et al. (2007).The set of partial differential equations in the EDM was numerically solved by means of an implicit method; details on operator splitting techniquesequential non-iterative approach, and definition of the function residuals and the Jacobian matrix are provided in Aguilera et al. (2005) and Regnier et al. (2002).Model runs were performed on SCSCF (www.dais.unive.it/scscf),a multiprocessor 25 cluster system owned by Ca' Foscari University of Venice running under GNU/Linux.

Model application: simulations set up and EDM calibration
Based on the rearing cycle characteristics, described in section 2.1, in model simulations shellfish are stocked in September, and harvested after one year.The farm, made up of n L =25 longlines of lenght l L =2000 m each, occupies a total area of 2 x 1 km 2 .Longlines orientation with respect to the North, D n , was of 60°.A biomass density, B i = 15 ind.m -2 was considered at 30 all the longlines of the farm.A depth range, z fmin -z fmax , comprised between 1 m and 7 m was selected as representative for ropes at the farm considered in this study.
As explained at the beginning of this section, the deposition model application was aimed at constraining the typical variability of mussel deposition at the farm site.In order to simulate the average flux of POC deposited by the farm, we carried out 10 different runs, considering each one a rearing cycle under forcings for a different year within the 2002-2011 35 time frame.Time series of monthly Sea Surface Temperature (SST) and concentration of chlorophyll-a were extracted from the EMIS (http://emis.jrc.ec.europa.eu/)data base from July 2002 to December 2012, (longmin 12.5;longmax 12.6;latmin 45.4;latmax 45.5) by means of the R package EMISR v0.1 (R version 3.0.3).Chlorophyll-a and SST data were derived from the sensor Modis (Moderate Resolution Imaging Spectroradiometer) Aqua and Terra respectively, with a spatial resolution of 4km.These data were used as inputs for individual-based population dynamics model, following a procedure previously 40 applied -see e.g.Filgueira et al. (2013).Due to the lack of long term time series of data, an average POC concentration had to be imposed, 0.1 mg L -1 , on the basis of a time series of monthly data collected at a nearby farm between 2006 and 2007 Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-206Manuscript under review for journal Biogeosciences Discussion started: 7 June 2017 c Author(s) 2017.CC BY 3.0 License.(Brigolin et al., 2009).The POM/TSM ratio was fixed on the basis of the time series collected within the same work (Brigolin et al., 2009;Rampazzo et al., 2013), providing an average AE of 0.6.The median daily value of faeces and pseudofaeces fluxes from the 10 simulations was used as an input for the deposition model.Modelling deposition requires an input time series of water velocity at an hourly time step.These data were provided on the basis of a current meter deployment carried out between March and September 2010 at a station located approximately 500 5 m from the NE edge of the farm (Boldrin A. pers.comm).Current meter data were first processed by means of a classical harmonic analysis (Pawlowicz et al., 2002).On the basis of the procedure proposed by Jusup et al. (2007), the residual currents were therefore edited randomly for short periods of time in order to reproduce the variability recorded from current meter measurements during extreme events (i.e.storms).According to these authors, the above mentioned procedure should be preferred for forecasting simulations (while for hind-casting the current meter time series can be directly loaded into the 10 model).
A steady-state EDM simulation was carried out at station EST1 located out of the mussel farm.The model was applied by using the same parameterization adopted for the North Western Adriatic shelf by Brigolin et al. (2011) -see Tables A3-A4

Sediment traps measurements
As part of this study, sedimentation fluxes were measured by means of sediment traps positioned at the bottom.Two 48 hours sediment trap deployments were carried out.The first experiment was performed between 29/08/2014 and 31/08/2014, at the end of the annual mussel rearing cycle.The second, was carried out between 11/09/2015 and 13/09/2015, at the very 25 beginning of the annual mussel rearing cycle.Three PVC sediment traps were deployed at each station: cylindrical shape; aspect ratio 5:1; collecting area of 0.0095 m 2 each (Cromey et al., 2002;Jusup et al., 2009).In the first experiment (see Figure 1), sediment traps were deployed at four locations; two stations IN1 and IN2 located inside the modeled depositional footprint, and two stations outside EST1 and EST2.For the second experiment, one station was located inside IN1 and one outside EST1.Upon collection traps content was filtered through pre-combusted (450°C, 4h) and pre-weighed Whatmann 30 GF/F filters.Filters were dried at 60 •C for 24 h and re-weighed.Organic carbon was determined using a Thermo Elementar Analyzer (Flash -EA 1112), after acidification with HCl for removing carbonates before analysis.

Sediment coring, microelectode measurements, porewater and solid phase analyses
Sediment were sampled at stations IN1 and EST1 in June 2015 (respectively on 23/06 and 24/06).Cores were collected by means of an Uwitech corer (10 cm diameter; 20 cm avg.penetration depth).Water was sampled 2 m above the bottom by 35 means of a Niskin bottle, for dissolved oxygen, salinity and temperature determinations.Cores were immediately brought back to the lab.and prepared for microprofiling which was conducted a few hours after coring.Cores were bubbled with air during measurements to allow aeration and gentle stirring.Microprofiling was conducted with a Unisense motorized microprofiler.Four oxygen microprofiles were performed using 100 µm tip microsensors which were calibrated by a twopoints method: Winkler titration of the overlying water (with a precision of 2 permil) and zero-oxygen signal in the anoxic 40

POC reaching the sediment-water interface: model simulation results and sediment trap experiments 10
Results obtained by means of the population dynamic model are reported in Figures 3a-d.The time series of satellite data, sea surface temperature and chlorophyll-a concentration, used to force the model are shown in Figure A1.Growth trajectories, here expressed in terms of soft tissues dry weight and shell length are in agreement with previous observations and model results obtained for this area (Brigolin et al., 2009).The minimum commercial size of 5 cm is achieved rapidly (<6 months), and the mussels reach the final length of 7 cm within 10-11 months.At the end of the cycle mussels present 15 weights comprised between 1.5 and 2.5 g.Total faeces and pseudofaeces release rates per line are shown in Figure 3c,d.These fluxes are highly variable, both at the seasonal (average CV within-year for faeces=0.67;pseudofaeces=1.47)and inter-annual time scales (average among years CV faeces=0.30;pseudofaeces=1.46).(Mann-Whitney one-tailed; n1=n2=6; p=0.03), while differences among IN1 and EST1 detected in September 2015 at the beginning of the rearing cycle were on the order of 3.3 mmol C m -2 d -1 , and not significant (Mann-Whitney two-tailed; n1=n2=3; p>0.5).

Early diagenesis processes underneath the farm and at the reference station 30
Bottom water temperature and salinity were respectively 22 °C and 36.3 psu at both stations.Oxygen in bottom waters, measured through Winkler titration, was 223.5 ± 0.3 µM at EST1 and 229.8 ± 0.1 µM at IN1.The grain size analysis classified sediments as medium silt at EST1 and very fine sandy-coarse silt at IN1. Porosity profiles (Figure 5a) show a decreasing trend going down-core, with a steep gradient in the upper 2 cm.A discontinuity is visible in IN1 core, between 4 and 5 cm.A total of 7 oxygen profiles were collected at the two stations (Figure 5b).Oxygen shows a quasi-monotonous 35 decrease in concentration.A slight increase at interface (~ 20 µM), most probably due to microphytobenthic production, is visible at station EST1 (bubbles were visible on the core surface after long-term exposure to sunlight).Indeed, results obtained at the two stations suggested a low variability in the oxygen behavior within the same core -profiles were measured by sampling randomly the available portion of the core in which no shell debris and macrobenthos were visible at the surface.Oxygen is consumed within the first mm, showing a higher penetration depth at EST1 (2.3 mm on average) with   2).
Figure 8 shows the partitioning among mineralization pathways, indicated by the electron acceptors, of the three pools of organic matter.At EST1, appoximately 27% of the total OM (3.0 mmol C m -2 d -1 ) undergoes oxic degradation, fraction similar to that observed at IN1 (30%), where the total OM flux, but also bioirrigation rate, are higher.As reported in Figure 7, 41% of the oxygen is consumed at EST1 by organic matter mineralization, and the remaining 59% is required for re-25 oxidizing reduced compounds (55% by reduced compounds oxidation, and 4% by nitrification).At IN1, 51% of the O2 is consumed by organic matter mineralization, 8% by nitrification and 41% by reduced compounds oxidation.
The extent of the depositional area obtained in this study (on average 50 m from the edge of the farm) is comparable with the 35 results obtained in previous studies.In an exposed site, Weise et al. (2009) (Cascapedia Bay, Canada; 20 m depth; mean current velocity of 10 cm s -1 ), constrained the area of higher organic enrichment within 90 m from the edge of the farm.
Dispersal area reported by Hatstein and Stevens (2005) was smaller, extending with a radius of approximately 30-40 m from the edge of the farm (20-30 m depth; mean current velocity of 3.4-4.0cm s -1 ).These differences in extent of the dispersal areas seem to be primarily associated to the action of currents and wave energy inducing resuspension of biodeposits 40 accumulated on the seabed (see Cromey et al., 2002) (5.82-6.56%).Corg content in deposited organic matter found in the area of study (4.23-6.34%) falls within the range reported by Giani et al. (2001) in different locations of the Northern Adriatic Sea (1.05-21.81%).Also, background total mass fluxes measured in the present work are comparable with data measured with sediment traps by the same authors, who reported fluxes of 5.8 g m -2 d -1 (TSS) at an off-shore station located in the Northern Adriatic Sea, and mean fluxes of 15 approximately 30.0 g m -2 d -1 in prodelta areas of Po and Adige rivers, with high annual variability (range 0.08 -240 g m -2 d -1 ). Relatively low background values obtained in our study can be associated to lower annual discharge (10.9 m 3 s -1 ) of Sile River than of Adige and Brenta ones (respectively, 139.5 and 85.9 m 3 s -1 -1994-2008 averages from Cozzi and Giani (2011)), being characterized by a particularly low value of the discharge rate in August and September (ARPAV, 2010) with respect to other months.Summarizing, the predicted fluxes agree well with sediment traps data, and mass fluxes measured 20 by traps are close to values reported for similar environments, being markedly lower than those recorded in lower energy environments.

Influences of perturbed POC fluxes on organic matter mineralization in sediments
Absolute values of the POC fluxes modelled and measured in this work were corroborated by the ED model calibration.The value of 11.6 mmol C m -2 d -1 estimated through the inverse use of the ED model at EST1 is an independent verification of 25 the results obtained from the sediment trap experiments.This value accounts for approximately 50 % of the flux measured from the traps (22 mmol C m -2 d -1 on average).This difference can be primarily related to the fact that the ED model accounts only for the reactive C, while the sediment trap captures all types of C, including the inert fraction.On the top of this, traps provided only a snapshot of the deposition in August and September -disregarding the seasonal variability in this deposition, and the presence of pulse deposition events.Finally, it is worth remarking that resuspension is a common 30 mechanism in shelf trap measurements: in the northern Adriatic Giani et al. (2001)  records not available from the study).In our case, even if the relative change of POC flux between the mussel farm and the reference is large, it does not affect significantly the porewater composition at depth.We cannot rule out spatial heterogeneity within each station, which would smooth the difference between IN1 and EST1.However, the fact that O 2 fluxes reflect the increased input by the mussels, points towards internal mechanisms that regulate porewater composition.
The shape of the DIC-NH 4 + profiles indicates bio-irrigation (Meile et al, 2001;Canavan et al., 2006), although the deeper 15 increase is not visible in our data profiles due to limited penetration of the cores.Indeed, the very limited increase in concentration profiles in the first centimeters can only be linked to input of bottom water with lower DIC and NH 4 + by irrigation, given the large recycling intensity in surface sediments as exemplified by O 2 profile.Therefore, one way to explain the lack of differences between the IN1 and EST1 porewater composition is to relate it to changes in the bioirrigation profiles at the two stations.We remark here that α 0 and xirr 1 were the only two parameters calibrated at IN1, and they 20 suggest a higher infauna activity, shifted towards the surface at this site.This feature was independently confirmed by a set of macrobenthos samples collected at the two stations as a part of a complementary study (Colla, 2017).Macrobenthos samples showed a higher diversity (48 vs 31 taxa recorded) and abundance (on average 1900 vs 1000 ind.m -2 ) at IN1 with respect to EST1, accompanied by the presence of larger organisms (0.065 g ind. - at IN1 versus 0.034 g ind. - at EST1).This is in agreement with the expected influence of biodeposition from mussel culture (McKindsey et al., 2011).Calibrated values 25 of α 0 (20 y -1 at EST1 and 40 y -1 at IN1) are close to those estimated by Meile et al. (2001) for Buzzards Bay site, 30-60 y -1 (water depth 15m), and slightly higher than those estimated by Canavan et al. (2006), through data fitting -10 y -1 at a coastal freshwater lake (Haringvliet, in the Netherlands).
According to model estimations, total mineralization at EST1 accounts for 96 % of deposited OM (OM1+OM2), with 0.5 mmol C m -2 d -1 escaping mineralization through burial.This fraction decreases to 94 % at IN1 (1.1 mmol C m -2 d -1 of non-30 mineralized POC), where the model is run under transient conditions.The relative contribution of different mineralization pathways to the total OM degradation is comparable to that reported by Pastor et al. (2011) for stations with comparable OM fluxes in the Rhône River prodelta and shelf area.At stations undergoing organic carbon deposition ranging between 7.3 and 16.2 mmol C m -2 d -1 (their stations L, I, C, J, F), these authors found oxic mineralization ranging between 44 % and 67 %, nitrification between 1 % and 6 %, and anoxic mineralization between 27 % and 51 %.The largest fraction of oxygen, 67 %-35 87 % was consumed by OM degradation, nitrification consumed up to the 31 % of the O 2 , while the re-oxidation of anaerobic products accounted for only the 2.1 % of oxygen consumption.This marks a difference with respect to the results obtained in the present work, since in our model a larger fraction of oxygen is consumed for reduced substances oxidation (54 % at EST1 and 41 % at IN1).This is also visible in Figure 6, where relevant drops in SO The higher influx of oxygen and enhanced bioirrigation at the mussel site (IN1) reflected on a substantial change in the pathways of oxygen consumption, with an increase of oxic degradation of OM2, and a relative decrease of oxygen demanded for reduced substances re-oxidation -this is clearly linked to the input of fresh OM from the mussels which is mineralized aerobically.Nonetheless, due to the higher absolute fluxes of OM, the oxygen consumed by reduced substances oxidation is higher at IN1 7.4 mmol O 2 m -2 d -1 , than at EST1, 6.2 mmol O2 m -2 d -1 .Higher NH 4 + concentration predicted at station EST1 5 with respect to field data could be explained with a higher rate of nitrification at this station.However, in the calibration performed within this work, the kinetic constant for nitrification was kept at its original value (see Table A4), due to the lack of data concerning NO 3 -.

Integrated model features
The pelagic deposition model allowed simulating the extent of the deposition area, and its variability with time.Being 10 integrated with a daily time step, the Mediterranean mussel population dynamic model allows to combine instantaneously the non-linear effects of the different environmental variables and physiological processes acting on deposition (i.e.water temperature, chlorophyll-a concentration, allometric dependence of the clearance rate on body size), and integrate these effects along the time of the farming cycle.The combination of the bioenergetics-based population model, which allows estimating organic matter production from the lines, with the deposition model, accounting for particles dispersion, 15 represents a novel aspect of the present work with respect to previous modeling studies on mussel deposition.Hartstein and Stevens (2005) modeling study applied a sensitivity approach to study organic deposition from Perna canaliculus in New Zealand, comparing sites characterized by different hydrodynamic exposure, and assuming an arbitrary particle-release rate.Weise et al. (2009) modeled mussel biodeposition at different sites in the eastern coast of Canada, imposing organic wastage from the farm lines as a model input, on the basis of site-specific field measurements (Callier et al., 2006), and extrapolation 20 from other sites.It is worth remarking here that the integration of growth and deposition models can represent a resource, allowing to apply the model at different sites in which environmental variables are known, without the need of performing in situ estimations of biodeposis production.On the top of this, the model could be used to explore the effect of climatechange-induced long term trends of variation in water temperature and particulate organic matter concentrations, which are expected to have an influence on mussel growth performances (Cochrane et al., 2009;Rosa et al., 2012).25

Conclusions
The combined application of an early diagenesis model and of a model of POC production and deposition from shellfish filter-feeders, allowed to study quantitatively the differences induced on sediment biogeochemistry by a local perturbation of the natural POC downward flux.This is one of the few existing attempts to couple pelagic mussel production models with early diagenetic models in order to investigate the effects of a local gradient of disturbance on coastal sediments 30 biogeochemistry and benthic-pelagic coupling (see the review by Paraska et al., 2014) Reaction-transport equations:

5
Rate laws  DB was supported by an individual mobility grant from the University of Venice Ca' Foscari (IRIDE-DHAMACO project).
Model runs were performed on SCSCF (www.dais.unive.it/scscf),a multiprocessor cluster system owned by Ca' Foscari University of Venice running under GNU/Linux.We gratefully acknowledge G. Monvoisin from GEOPS France for the assistance in the analyses, A. Boldrin from CNR for giving access to the current meter data, and Nautica Dal Vì for providing space for setting up lab facilities.10 . This model 40 requires in input time series of water velocity and fluxes of faeces and pseudofaeces.The latter two time series are provided by the individual-based population dynamic model with a daily resolution.Fluxes associated with the metabolic activities of Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-206Manuscript under review for journal Biogeosciences Discussion started: 7 June 2017 c Author(s) 2017.CC BY 3.0 License. the entire farmed population are computed by integrating individual fluxes over the size structure of the population, and over time, see Eq.( and by calibrating the POC downward flux, Φ OM1+OM2 and the parameters, α 0 and x irr1 , defining respectively intensity and depth of bioirrigation.Initial values of POC for the calibration were set on the basis of the results of the sediment trap 15 experiment (see section 2.4 below).This steady-state model calibration was carried out by fitting the O 2 , DIC, NH 4 + and SO 4 2-profiles observed at station EST1.The transient simulation carried out at station IN1 was performed by imposing as initial conditions model outputs obtained at station EST1.The model was run for 20 years (time of activity of the farm) by taking into account the average OM 3 flux (faeces+pseudofaeces) predicted by the deposition model at station IN1.α 0 and x irr1 parameters were calibrated by fitting the O 2 , DIC, NH 4 + and SO 4 2-profiles observed at station IN1.Diffusive oxygen uptake 20 was calculated from profiles (both model and data, see section 2.5 below) by means of the 1-D Fick's first law of diffusion.
layer below the oxic zone.Porewaters were extracted and sampled under N 2 within 4 hours after coring.Samples were Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-206Manuscript under review for journal Biogeosciences Discussion started: 7 June 2017 c Author(s) 2017.CC BY 3.0 License.preserved using HgCl 2 saturated for DIC and Total Alkalinity (TA) analysis, by freezing for NH 4 + determination and by acidifying for SO 4 2-analysis.Measurements were performed in the laboratory: DIC was analysed on a DIC analyzer (Apollo SciTech®) using 1mL sample volume with 4 to 6 replicates which provided a standard deviation of 0.5%.TA concentrations were measured in a potentiometric opencell titration on 3mL sample volume (Rassmann et al., 2016) with an uncertainty of 0.5%.Ammonium concentration was measured spectrophotometrically following Grasshof et al. (1983) with an uncertainty 5 of about 5%.Sulfate were measured by HPLC on a Dionex Ionic Chromatography System (ICS 1000) using a Dionex IonPac AS9-HC Carbonate Eluent Anion-Exchange Column (4 x 50 mm) and a Dionex IonPac AG9-HC Guard Column (4 x 50 mm) with an uncertainty of 1%.

Figure 4
Figure 4 shows the map of the model-predicted fluxes of organic C induced by the mussels reaching the sediment at days 1, 120, 240 and 360 of simulation (Sep 10; Dec 29; Apr 28; Aug 26).As can be seen, magnitude of the fluxes increases with 20 mussel growth.Maximum organic carbon fluxes predicted at day 10 and 360 are 2.5 and 13.3 mmol C m -2 d -1 , respectively.Footprint induced by the presence of the lines is clearly visible at days 120 and 360.At day 240, a remarkable displacement of deposition towards SW (approximately 200 m) was detected.In the other cases (Figs.4a,b,d), the maximum footprint distance from the edge of the farm is 50m.TSS, OC%, and POC fluxes measured in August 2014 and September 2015 are reported in Table1.In August 2014,POC  25    fluxes at the end of the rearing cycle showed significantly higher values in stationsIN1 and IN2 than in EST1 and EST2 Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-206Manuscript under review for journal Biogeosciences Discussion started: 7 June 2017 c Author(s) 2017.CC BY 3.0 License.respect to IN1 (1.4 mm on average).DIC, NH 4 + and SO 4 2-data are shown in Figures 5c-e.DIC concentration profiles at the two stations are comparable, although at depths >10 cm they stabilize at values around 3.6 mM at EST1 and 2.8 mM at IN1.A similar pattern is visible for ammonium, which average concentration below 10 cm depth reaches values of 50 µM at IN1 and 125 µM at EST1.The effect of bioirrigation is visible within the upper 10 cm, and, is particularly marked at station IN1, where DIC shows a decrease in concentration starting from 4 cm depth.SO 4 2-concentrations are similar between the two 5 stations, and do not present marked variations going down-core.Model-predicted profiles (EDM) calibration are given in Figure 6, and compared with the measured O2, DIC, NH 4 + and SO 4 2-field data.Porosity parameters used in the model, and given in Table 2, were estimated by independently fitting the two porosity profiles shown in Figure 5a.A metabolisable OC flux of 11.6 mmol C m -2 d -1 was estimated by calibrating the model at steady-state at station EST1.At IN1, an additional flux of labile organic matter (8.2 mmol C m -2 d -1 ) was imposed 10 for 20 years, based on the median value predicted for POC flux by the pelagic deposition model.In the EDM, DIC and NH 4 + profiles are both characterized by a concentration enhancement within the first cm, controlled by the degradation of the labile organic matter (OM1 and OM3), and subsequently modulated by the action of bioirrigation, which causes a down-core decrease.In general, simulated profiles are in good agreement with field data, although, at station EST1, predicted NH 4 + exceeds observed concentrations.The depth distributions of bioirrigation coefficients were estimated independently at the 15 two stations, obtaining values reported in Table2, which indicate that infauna activity is higher and more concentrated in the superficial layer at IN1 than at EST1.
Figure 7 compares these fluxes to the ones estimated by means of microelectrode profiles at the same stations.As can be seen, model results are in good agreement with the ones calculated from microelectrode profiles.Oxygen fluxes predicted by the model profiles are, respectively, 11 and 18 mmol O2 m -2 d -1 at EST1 and IN1, while micro-electrode ones range between 10 and 12 mmol O2 m -2 d -1 , at EST1, and 19 and 31 at IN1. 20 Total mineralization calculated by the model is 11.1 mmol C m -2 d -1 at EST1, and 18.7 mmol C m -2 d -1 at IN1 (see Table . Magnitudes of OC fluxes predicted by the model were corroborated by Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-206Manuscript under review for journal Biogeosciences Discussion started: 7 June 2017 c Author(s) 2017.CC BY 3.0 License. the two sediment traps.Simulated biodeposition fluxes of 2.5 mmol C m -2 d -1 on September 10, at the beginning of the cycle, agree well with the very limited and not statistically valid POC flux difference (IN-EST) between average trap measurements at the same time (3.3 mmol C m -2 d -1 ).The 13.3 mmol C m -2 d -1 value predicted on August 26 is very close to the maximum difference between IN1 and EST1 recorded by sediment traps in August 2014, 14.2 mmol C m -2 d -1 .These values are lower than those reported for less exposed sites (lagoons and bays).Hartstein and Stevens (2005) recorded 5 increases in mass fluxes on the order of 60 g m -2 day -1 at a mussel farm in Catherine Cove and Elayne Bay (Malborough Sounds, New Zealand), and Jaramillo et al. (1992) observed differences of 153 g m -2 d -1 in the Queule river estuary (southernChile).These are one to two orders of magnitude higher than the difference of 2 g m -2 d -1 recorded between IN1 and EST1 in August 2014.Values lower and closer to those presented here, on the order of 10 g m -2 d -1 , were reported byWeise et al. (2009) for a longline mussel farm located at 20 m depth in a high energy environment (Cascapedia Bay, Canada).Regarding10Corg contents, values found in the present study are in good agreement with those reported byHartstein and Stevens (2005) estimated the contribution of resuspension to the gross flux sedimented in traps reaching 43.7 %.Dedieu et al. (2007)  compared seasonal cycles of C, N and O inside and outside a shellfish farming area in a lagoon located in Southern France, by combining a steady state diagenetic model(Soetaert et al., 1996) with a comprehensive set of experimental data.Model results showed that organic matter accumulation at the farming area enhanced the anaerobic 35 metabolism.Oxygen microprofiles recorded byDedieu et al. (2007)  inside and outside the mussel farm presented differences which are comparable to those recorded in the present work, with a decrease of 50% in oxygen penetration depth, and an increase up to 3 times of diffusive O 2 fluxes (30 mmol O 2 m -2 d -1 versus 90 mmol O 2 m -2 d -1 ).InDedieu et al. (2007)  this was accompanied by a remarkable enhancement in NH 4 + concentrations, which was not visible in the field data reported in the present study.This can certainly be related to the difference in the rate of biodeposit accumulation, although other factors 40 can also play a role (changes of local macrobenthic activity, and of coupled nitrification-denitrification rates).Mean OC fluxes estimated by Dedieu et al. (2007) by calibrating the diagenesis model at the mussel farm site, were 38.4,108.0, and Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-206Manuscript under review for journal Biogeosciences Discussion started: 7 June 2017 c Author(s) 2017.CC BY 3.0 License.108.0 mmol C m -2 d -1 , respectively in winter, spring and summer, with increases by 26.4,53.4, and 52.4 mmol C m -2 d -1 with respect to the reference station in the respective seasons.In our work, the background flux (OM1+OM2 flux at EST1) was lower, 11.6 mmol C m -2 d -1 with an increase of 8.2 mmol C m -2 underneath the farm, at IN1.The relative increase, with respect to the background flux, was of 71 %, which is comparable to the 69 % increase found by Dedieu et al. (2007) in winter.Nonetheless, in absolute terms, Dedieu et al. (2007) reported a difference of approximately 50 mmol C m -2 d -1 in 5 summer, in a system already characterized by large fluxes of organic matter, and hence dominated by SO 4 2-reduction, whereas in this study, the increase is 5 times less -10 mmol C m -2 d -1 , over a system which has mostly oxic and suboxic diagenesis.Difference in absolute values of OC fluxes can be related to a higher mussel stocking density (production in the Thau lagoon is 2333 tons km -2 versus 1450 tons km -2 at the farm in Jesolo), lower depth (7 m Thau Lagoon, 14 m Jesolo), and hydrodynamic regime of the site (lagoon environment with low energetics in term of hydrodynamics -current meter 10 4 2-concentration are predicted by the model below the bioirrigated layer -approximately 2 mM at EST1 and 5 mM at IN1.The different behavior may be 40 due to the processes controlling H 2 S and Fe in Pastor et al.'s study where most H 2 S is precipitated as FeS 2 thus escaping reoxidation by sulfides.A more precise estimation of the fate of this oxygen could be obtained by introducing in the model FeS precipitation, for which at least Fe 2+ measurements in pore waters would be required.Biogeosciences Discuss., https://doi.org/10.5194/bg-2017-206Manuscript under review for journal Biogeosciences Discussion started: 7 June 2017 c Author(s) 2017.CC BY 3.0 License.

Figure 1 :
Figure 1: Study site and stations sampled in this study.

5 Figure 2 :
Figure 2: Information flow within this study.

Figure 3 : 10 Figure 4 :
Figure 3: a,b) Mussel soft tissues dry weight and shell length; c,d) faeces and pseudofaeces release fluxes per mussel line per day.

Figure 8 : 25 30
Figure 8: Model estimations: relative importance of different pathways on the total organic matter at EST1 (a,b), and IN1 (d,e,f); ratios of different redox pathways on total oxygen consumption at EST1 (c), and IN1 (g).25