The diagenetic Si model was modified from Khalil et al. (2007) and implemented in the R software (http://cran.r-project.org). A short description of equations, processes and parameters is given below.

The general equation of the reactive-transport model affecting solids and solutes (Eq. 1) represent diffusive and advective transport (first and second terms, respectively) and reactive processes (third term):

All fluxes are in mmol m^-3 d^-1 and are detailed below.

The originality of the work of Khalil et al. (2007) is to implement two variables of aSiO₂ with different reactivity. The model has thus three sets of variables: less reactive aSiO₂ concentrations (C_aSiS , µmol l^-1), highly reactive aSiO₂ concentrations (C_aSiF , µmol l^-1) and dissolved Si(OH)₄ concentrations (C_dSi , µmol l^-1) implemented in Eqs. 2, 3, 4 and 5.

The transport of the two solid fractions - less and highly reactive aSiO₂ - is controlled by the biodiffusion rates (D_b , m² d^-1) linked to sediment mixing by benthic organisms, and by the advection rates (w, m d^-1) linked to the accumulation of newly deposited particles and steady-state compaction.

The biodiffusion rate is depth-dependent and calculated at each depth as follows:

where D_b0 is the surface biodiffusion coefficient (m² d^-1), z_b is the biodiffusion depth for solids (m) and coeff_b is the exponential decrease constant for biodiffusion (m).

The aSiO₂ dissolution term (R_dissS or R_dissF ; Eqs. 4 and 5) is parametrized by different dissolution rates (k_aSiS and k_aSiF , d^-1, dependent on in situ temperature and salinity) and equilibrium concentrations (C_dSieqS and C_dSieqF , µmol l^-1) for each aSiO₂ phase.

Note that the porosity ϕ is depth-dependent and defined at each depth as follows:

Where ϕ 0 is the porosity at sediment-water interface, ϕ ∞ the asymptotic porosity and z p o r the exponential decrease constant for porosity.

The reactive-transport of solutes - Si(OH)₄ (noted dSi in equations) - is modified from Khalil et al. (2007) and described in Eq. 6. Advection is neglected because of the dominance of the diffusive process (Peclet number >> 1; McManus et al., 1995). The processes of molecular diffusion, dissolution (R_diss = R_dissF + R_dissS ) and reprecipitation (R_precip = K_p (C_dSi – C_dSieqp)) are kept as described by Khalil et al. (2007). Secondly, a bioirrigation process (R_irr ) is added to take into account the transport of Si(OH)₄ with solutes by bioturbation activities.

Bioirrigation is represented with a non-local term described in Eq. 7 for depth z ≤ z_irr (Emerson et al., 1984; Boudreau, 1994).

The bioirrigation rate parameter α (d^-1) is constant over the bioirrigation depth (m) and exponentially decreases below that depth. C d S i 0 is the concentration of Si(OH)₄ at the sediment-water interface (z=0).

At the sediment-water interface (z=0), diffusive fluxes of Si(OH)₄ are prescribed by Fick's first law of diffusion (Eq. 8).

Where D_dSi0 is the diffusion coefficient for Si(OH)₄ (m² d^-1) and ∂ C d S i ∂ z is the concentration gradient of Si(OH)₄ (mmol m^-3 m^-1) at sediment-water interface.

The upper boundary of the solid fraction is prescribed as a net deposition flux (including resuspension and benthic primary production). At the lower boundary, zero-gradients are imposed for both solid and solute fractions.

The model is solved to steady-state, using methods implemented in R-package rootSolve (Soetaert, 2009). Mass balance is calculated in order to ensure the internal validity of the mathematical steady-state solution. At all times, the difference between the deposition flux of aSiO₂ and the sum of benthic fluxes of Si(OH)₄, burial and reprecipitation fluxes of aSiO₂ was < 10^-12.

Vertical profiles of sediment porosity: Sediment porosity in each sediment slide over depth was obtained after drying wet sediment of known volume for 5 days, after which the loss of weight was determined (Berner, 1980). This data was used to fit vertical porosity decreasing profiles.

Pelagic and benthic aSiO₂ contents: Pelagic aSiO₂ concentrations were determined using the sequential alkaline digestion method of Ragueneau et al. (2005b) and benthic aSiO₂ contents were quantified as in DeMaster (1981). Both methods allow to correct amorphous silica concentrations from lithogenic silica interference which is essential in environments rich in aluminosilicates - e.g. estuaries. Benthic aSiO₂ concentrations are expressed as % to refer to % g gDW^-1.

As aSiO₂ concentrations can vary from values less than 1% to more than 50% depending on study sites, we used aSiO₂ profiles measured in 2008 and surface aSiO₂ concentrations in 2009 to constrain modelled aSiO₂ profiles at each station. These profiles were used to estimate the range of aSiO₂ contents ( ± 1-3%) rather than to represent the fine scale vertical discontinuities (that cannot be captured by steady-state modeling). As no aSiO₂ profile was available for July, an annually averaged aSiO₂ profile of February, April and November was calculated for each station.

Vertical profiles of Si(OH)₄ concentrations: Si(OH)₄ concentrations in porewaters of each sediment slice were measured with an Auto Analyser III (Bran+Luebbe^®) using the method of Tréguer and Le Corre (1975). The triplicate Si(OH)₄ profiles measured at each station and season in 2009 were averaged and used for fitting.

Chl a and Phae concentrations: were performed on surface water samples by using the SOMLIT protocol (http://somlit.epoc.u-bordeaux1.fr) and on surface sediment (0.5 cm) by using a method adapted from Lorenzen (1966) and Riaux-Gobin and Klein (1993). For surface sediment, 10 ml of 90% acetone was added to each sample that was stored in the dark under constant agitation at 4°C for approximately 18 h. Chl a and Phae were respectively measured in the supernatant before and after acidification with a KONTRON fluorimeter (Kontron Instruments).

Proportion of highly reactive aSiO₂ in sediments: Dissolution experiments were carried out on 3 sediment layers (0-1, 2-3 and 6-7 cm), at all stations sampled in February 2008. The temporal increase of Si(OH)₄ concentrations was monitored during 21 days in batches containing the sediments and artificial seawater close to the measured in situ conditions: in situ salinity, pH of 8, constant temperature and in the dark to avoid any production of biogenic silica due to benthic organisms. The concentrations were normalized to the initial introduced Si concentrations in order to determine the dissolution rates of highly and less reactive aSiO₂ phases (the first and second slope of the curve). The proportion of highly reactive aSiO₂ was determined by statistical inverse modeling (Model 2 described in Moriceau et al., 2009). Note that we used the proportion of highly reactive aSiO₂ but not the dissolution rates to constrain the benthic Si model. Indeed, discrepancies between predicted and measured dissolution rates are expected in response to 1) the high detrital-to-amorphous opal ratio in our sediments (> 30%; Ragueneau and Tréguer, 1994), which increases aluminium concentrations and decreases the dissolution rate and solubility (Dixit et al., 2001); 2) the transport of fresh aSiO₂ to deep sediment layers by nonlocal bioturbation (Gallinari et al., 2008); 3) the absence of feedback between reprecipitation and dissolution in the model formulation (Khalil et al., 2007); 4) the cleaning of aSiO₂ particles from their alumino-silicate coating at the beginning of batch experiments (Khalil et al., 2007); and 5) agitation during experiments, which increases dissolution rates compared to stationary sediments (Fabre et al., 2019).

Biodiffusion coefficients: Luminophores were visualized and counted at each sediment layer by image processing (Michaud, 2006). The vertical profiles of luminophores were adjusted to a reaction diffusion type model in order to quantify biological sediment transport (Duport et al., 2006; Oleszczuk et al., 2019) and more specifically the biodiffusion-like transport coefficient (Db; m² d^-1). The best fit between the observed and modelled tracer distribution is estimated by the least-squares method and produces the best Db coefficients. Sediment biodiffusion rates and depths were quantified at each station in February, July and October.

As detailed above, the model contains a large number of parameters that could potentially lead to multiple solutions. In order to reduce uncertainties and ensure the uniqueness of the solution calculated by the model, the strategy was to constrain parameter values with a high number of direct observations and experimental data obtained in this study, as done in (Khalil et al., 2007). All model parameters were determined from measurements and/or inverse modeling and summarized in Tables 1 , 2 .

Direct observations: Temperature T, salinity S and Si(OH)₄ concentrations in overlying water C d S i 0 were determined from direct measurements and summarized in Table 1 . Parameters ϕ 0 , ϕ ∞ and c o e f ϕ were adjusted to fit vertical porosity profile measurements ( Table 2 ).

Assumptions based on experimental data and other measurements ( Table 2 ): Most parameters were obtained from measurements detailed in sections 2.2 and 2.3. The biodiffusion rates D_b0 and depths z_b were determined from ex situ experiments. Equilibrium concentrations for the two phases were assumed to be identical and equal to deep asymptotic Si(OH)₄ concentrations obtained from vertical profiles of Si(OH)₄. Based on the dissolution experiments, the proportion of highly reactive aSiO₂ in the deposition flux was set to 0.5 as a compromise between the low values ~ 0.3-0.5 - resulting from statistical inverse modeling of dissolution experiments performed on sediments collected in February 2008 (with the model 2 of Moriceau et al., 2009) - and the generally higher values ~ 0.5-0.9 from Khalil et al. (2007). The bioirrigation depths z_irr were visually adjusted to the vertically constant portions of Si(OH)₄ profiles where fauna was observed. Accumulation rates w were estimated from radionuclide measurements (~ 0.005 cm d^-1; Khalil et al., 2018). Reprecipitation rates K_p were set in the range 10^-5-10^-4 d^-1 because of the high detrital sediment content known to enhance reprecipitation rates (Khalil et al., 2007).

Inverse modeling ( Table 2 ): was performed to estimate the values of undetermined parameters by using the package FME (Soetaert and Petzoldt, 2010). The non-linear fitting procedure using the Levenberg-Marquardt algorithm aims to minimize the sum of squared residuals of model outputs (aSiO₂ and Si(OH)₄ profiles) versus data and reproduce the curvature of these profiles. The following parameters were estimated: aSiO₂ deposition flux, dissolution rates of highly and less reactive aSiO₂ and non-local bioirrigation rate. When necessary, final adjustments were performed on the proportion of highly reactive aSiO₂ and reprecipitation rates. During the fitting process, upper and lower bounds in the range of expected likely values were imposed for each fitted parameter.

Model outputs (deposition fluxes of less and highly reactive aSiO₂, benthic fluxes of Si(OH)₄ including diffusion and bioirrigation, burial fluxes, reprecipitation fluxes) were used to build benthic Si budgets and to estimate the seasonal and spatial variation of fluxes at each station of the two estuaries (in mmol m^-2 d^-1).

Fluxes of the benthic Si budgets were then extrapolated to the entire surface of each estuary and compared to river Si fluxes and pelagic production. All fluxes are presented in a conceptual scheme ( Figure 3 ).

River Si fluxes were determined as the sum of Si(OH)₄ and aSiO₂ fluxes. Si(OH)₄ fluxes were estimated as the product of weekly Si(OH)₄ concentrations (from the citizen-science network ECOFLUX; Abbott et al., 2018) by weekly river flow (from Banque Hydro) at the outlet of Elorn and Aulne rivers. aSiO₂ fluxes were estimated as the product of seasonal aSiO₂ concentrations (this study) by the same weekly river flow. Weekly fluxes were summed to obtain seasonal and annual fluxes.

Deposition, burial and benthic Si(OH)₄ fluxes were estimated by multiplying modelled fluxes at the 3 stations of the 2 estuaries (section 2.5) by the estuary surface of each estuarine section in order to estimate the retention of Si in the estuaries (in kmol) at the seasonal and annual scales. The surface of spatial integration for the 3 stations in each estuary was determined using a GIS-based approach (Khalil et al., 2018).

Burial fluxes in tidal saltmarshes were estimated for saltmarshes invaded and non-invaded by Spartina alterniflora in the Elorn estuary (Querné, 2011). The ratio of intertidal burial over subtidal burial was calculated for the Elorn estuary (0.69) and used to estimate intertidal burial in the Aulne estuary.

Export was the difference between river fluxes and burial (in subtidal sediments and in intertidal saltmarches).

Pelagic primary production was estimated as the sum of pelagic production rates (measured in this study at salinity 0, 5, 10, 15, 20, 25 and 30) vertically integrated on the seven estuarine boxes in each estuary for each season. Incubations were performed at in situ conditions (temperature, attenuated light) for 24h with ¹⁴C tracer following the method of Le Bouteiller et al. (2003). The range of pelagic BSi production rates was estimated by multiplying carbon primary production rates by the Si:C factor of 0.03 and 0.13. The factor 0.13 is characteristic of 100% diatoms and preservation, while 0.03 is used to consider that only 25% is related to diatoms and/or preserved. The volume of spatial integration was described with a power function (Eq. 9; Soetaert et al., 2006) and validated with morphological observations (Bassoulet, 1979; Google Map^©; this study).

with x is the estuary volume. x_riv and x_sea are the values at the freshwater and marine water end-members of the estuary, respectively. The exponent a regulates the steepness of the relationship.

Temporal changes of Si(OH)₄ loads at the freshwater end-member

The Si(OH)₄ concentrations at the freshwater end-member (station E1 and A1) were similar in Elorn and Aulne from December to April (140 ± 29 and 140 ± 22 µmol l^-1, respectively), but were higher in Elorn from May to September (203 ± 15 and 108 ± 31 µmol l^-1, respectively; Figure 2D ).

The Si(OH)₄ fluxes at the freshwater end-member were higher 1) in winter when the river flow increased, and 2) in the Aulne River from December to April (42 ± 30 t d^-1 for the Aulne River, versus 7 ± 5 t d^-1 for the Elorn River; Figure 2E ). The fluxes were lower and more similar from May to October (5 ± 4 t d^-1 in the Aulne River versus 2 ± 1 t d^-1 in the Elorn River), but they were occasionally higher at the mouth of the Elorn River in July and October ( Figure 2E ). The relative contribution of the Aulne River to the Si(OH)₄ fluxes in the Bay of Brest ranged from a minimum of 40% in summer to a maximum of 80% in winter.

Seasonal distributions of Si(OH)₄, aSiO₂, and Chl a concentrations along the estuaries

The Si(OH)₄ concentrations behaved almost conservatively and continuously decreased along the salinity gradients of both the Elorn and Aulne estuaries at all seasons except July ( Figure 4A ). In the Elorn Estuary, the Si(OH)₄ concentrations were similar at these three seasons. In the Aulne Estuary, the Si(OH)₄ concentrations at the freshwater end-member decreased from fall/winter to spring/summer, leading to seasonal variations of the slope. In summer, the curve was no longer conservative, which suggested Si(OH)₄ consumption in both estuaries. Note that the Si(OH)₄ (and aSiO₂; Figure 4B ) concentrations at the freshwater end-member in the Elorn Estuary were lower than at salinity 5 in February. This non conservative profile resulted probably from the transient increase of river flow ( Figure 2C ) as explained by Regnier et al. (1998).

In both estuaries, the aSiO₂ concentrations decreased with increasing salinity at all seasons, except July ( Figure 4B ). In the Aulne Estuary, a steep decrease of aSiO₂ concentrations occurred between salinity 0 and 5. The highest aSiO₂ concentrations were observed in July and reached a maximum at the salinity range of 5-15 (data not available for the Elorn Estuary).

In spring and summer, maximal concentrations of Chl a were generally observed between salinity 5 and 20, as observed for aSiO₂ concentrations in July. In November, Chl a concentrations were low and constant with salinity, which completely differed from the high and decreasing aSiO₂ concentrations. In both estuaries, the Chl a:(Chl a+Phae) ratios were higher in spring/summer than in winter and reached an annually constant marine end-member ratio in the range 0.6-0.7 ( Figure 4D ). The aSiO₂:Chl a ratios were low in spring/summer and very high in winter, especially in upper estuaries ( Figure 4E ).

The aSiO₂ concentrations were in the range 0.5-3% in the surface sediments ( Figure 5A ). At each season, they decreased from upstream to downstream in the Elorn Estuary, and they were minimal at the intermediate station in the Aulne Estuary ( Figure 5A ). At all stations, surface Chl a concentrations increased from February to July, and then decreased in November ( Figure 5B ). The Chl a:(Chl a+Phae) ratios were always > 0.4 ( Figure 5C ). The aSiO₂:Chl a ratios decreased from February to July regardless of station, and from upstream to downstream in winter ( Figure 5D ).

Observed and simulated aSiO₂ and Si(OH)₄ profiles

The observed aSiO₂ profiles were vertically constant with small discontinuities along the sediment cores. The model represented well the range of aSiO₂ concentrations but failed to reproduce the small vertical heterogeneities. Although variable between stations, Si(OH)₄ concentrations in general increased with depth ( Figure 6 ; black circles). A regular hyperbolic profile, commonly observed in stable or oceanic sediments, was observed at station A1, while stations E1, E2, E3, and A3 showed concentrations that became stable over the first centimeters (of variable thickness) which then increased again at greater depths. In general the model provided a good fit to the observed Si(OH)₄ trends, with the exception of the subsurface maxima of Si(OH)₄ concentrations observed at station A2 in February and May ( Figure 6 ).

Estimated and calibrated parameters

The biodiffusion rates of sediments (obtained from core incubations with luminophores) increased from February to July ( Table 2 ). These rates were higher in the Elorn (1-4 10^-7 m² d^-1) than in the Aulne Estuary (< 1 10^-7 m² d^-1), and they were almost null at stations A1 and A2. The bioirrigation rates varied relatively similarly to biodiffusion rates ( Table 2 ). The calibrated dissolution rates of the less reactive aSiO₂ were in the range 10^-5- 2 10^-4 d^-1 ( Table 2 ). The lowest rates were observed at stations E1 and A2. The calibrated dissolution rates of the highly reactive aSiO₂ were in the range 10^-3- 2 10^-2 d^-1 ( Table 2 ), and generally increased from February to July or October; however, at stations A1 and A2, the highest rates were found in February. Reprecipitation rates generally varied between 10^-5 and 10^-4 d^-1.

Simulated fluxes

The simulated deposition fluxes of aSiO₂ were in the range 2-4.5 mmol m^-2 d^-1, regardless of season and station ( Figure 7A ). In both estuaries, maximal deposition fluxes were observed upstream of the estuaries in winter, and midstream/downstream during other seasons. The burial fluxes of the less and highly reactive aSiO₂ were 25-80% of the deposition flux of aSiO₂ ( Figures 7B, C , 8). In the Elorn Estuary and at station A3, the proportion of burial decreased from winter to summer. This decrease was due to an increase of sediment-water Si(OH)₄ effluxes by diffusion and non-local bioirrigation, with a generally higher contribution of bioirrigation ( Figures 7E, F ). The decrease of burial fluxes in summer is not observed at stations A1 and A2, where bioirrigation is low. Reprecipitation fluxes were generally < 5% regardless of station and season ( Figures 7D, 8 ).

The Aulne river-estuary had a larger contribution than the Elorn river-estuary to all fluxes because of the larger size of its watershed and estuary. The Table 3 summarizes the estuarine Si budget and retention for the Elorn, Aulne and the two estuaries, at seasonal and annual scales.

As previously observed ( Figure 2E ), total Si river fluxes are two to three orders of magnitude higher in winter (134074 kmol per season) than in summer (8417 kmol per season). The Si loads are dominated by Si(OH)₄ in both estuaries (64 to 90% from winter to summer). Compared to the Elorn river, the Aulne river brings 66 to 90% of river Si fluxes to the Bay of Brest from summer to winter.

Deposition fluxes are similar for all seasons (8570-10784 kmol per season) but strongly vary from 6% of river loads in winter to more than 100% during summer.

Recycling leads to benthic fluxes of Si(OH)₄ ranging from 3585 kmol in winter to 6409 kmol in summer. Bioirrigation accounts for 46% of benthic fluxes (15-68%), with stronger bioirrigation in the Elorn estuary (64-68%) than in the Aulne estuary (15-34%). During summer, benthic recycling (6409 kmol) provides a quantity of Si(OH)₄ close to river fluxes (7348 kmol Si(OH)₄; 87% of 8417 kmol total Si).

The reprecipitation flux is low and negligible with 168 kmol per year for the two rivers.

Burial in subtidal sediments accounts for 10-39% of river Si loads in the Elorn estuary and 3-38% in the Aulne estuary, from winter to summer. Burial in tidal marshes is estimated to ~2/3 of subtidal burial (based on data on the Elorn estuary; Querné, 2011). Including burial in tidal marshes leads to an estimated burial of 17-41% of river Si loads in the Elorn estuary and 16-64% in the Aulne estuary.

Finally, export is high during winter and fall (125 649 kmol; 93% of river flux; 51% of annual export flux), but lower during summer (2 779 kmol; 33% of river flux; 1% of annual export flux).

In this study, steady-state diagenetic modeling is used to quantify the averaged benthic processes at seasonal scale along two estuaries. Modeling is particularly efficient for determining the deposition fluxes, which otherwise are very difficult to assess by direct measurements (Ridd et al., 2001). Regardless of the station and season, the deposition fluxes in the Elorn and Aulne estuaries (2-4.5 mmol m^-2 d^-1; Figure 7A ) are in the same order of magnitude than in the seaward Bay of Brest (0.6-3 mmol m^-2 d^-1; Ragueneau et al., 2005a). These high deposition fluxes are consistent with the high accumulation rates measured through radionuclide measurements that constrain the model (Khalil et al., 2018). High deposition and accumulation of aSiO₂ are generally observed seaward in large rivers (e.g., the Amazon, Congo, and Yangtze rivers; DeMaster et al., 1985; Shiller, 1996; Raimonet et al., 2015) that have high river flow rates. In contrast, in our study, the aSiO₂ deposition fluxes are high in the estuary due to lower river flow rate and the semi-enclosed shape of the system. These deposition fluxes calculated in the macrotidal Aulne and Elorn estuaries are in the range of values determined for the macrotidal Scheldt river and estuary: from ~0.028-0.034 mmol m^-2 d^-1 in subtidal sediments in the Scheldt river (Carbonnel et al., 2009), to ~1.6 mmol m^-2 d^-1 in marsh sediments in the Scheldt estuary (Struyf et al., 2006) and a maximum of 8 mmol m^-2 d^-1 in subtidal sediments in the Scheldt river/estuary (Arndt and Regnier, 2007).

In this study, maximal deposition fluxes are observed upstream of the estuaries during winter, which is mostly related to higher aSiO₂ river concentrations and fluxes. High aSiO₂ concentrations and fluxes brought to the estuary are related to enhanced soil weathering and river flows ( Figure 2C ), as previously reported by Smis et al. (2011). These winter deposition fluxes are associated with detrital aSiO₂, as shown by high aSiO₂:Chl a and low Chl a:(Chl a+Phae) benthic ratios ( Figure 5 ).

Deposition is observed more downstream and all along the estuary in spring and summer. The matter deposited then is also more reactive because of a higher contribution of riverine and estuarine, pelagic and benthic, primary production, stimulated by increasing light and temperature, often observed in estuaries (Heip et al., 1995).

The increase in river primary production (either benthic and/or pelagic) is highlighted by a decrease in Si(OH)₄ concentrations (consumption) and an increase in Chl a concentrations at the freshwater end-members from winter to summer ( Figures 2 , 4 ). This riverine material generally settles before salinity 5 (Anderson, 1986) and fuels deposition in the upper estuary.

The contribution of estuarine pelagic primary production to aSiO₂ deposition in the two estuaries is highlighted by an increase in Chl a and aSiO₂ concentrations in the water column at salinity 5-20, and decreasing aSiO₂: Chl a ratios ( Figures 4 , 5 ), as well as production estimates along estuaries (estimated to be 57-225 10³ mol d^-1 from incubations performed in this study). Benthic primary production (subtidal or intertidal by lateral transport) are suggested by the Chl a concentrations in surface sediments at stations E1 and E2 that are 4-15 times higher than in Aulne Estuary in summer, whereas Chl a concentrations in the water column are similar in both estuaries. This pelagic and benthic production must have settled inside the estuary depending on tidal range and river flow.

The relatively low increase of deposition fluxes during summer compared to winter ( Figure 7A ), as well as lower benthic than pelagic Chl a:(Chl a+Phae) ratios in summer ( Figures 4D , 5C ), also suggest the enhancement of pelagic dissolution in summer; dissolution increases the degradation state of deposited aSiO₂ and decreases the quantity of aSiO₂ at the sediment-water interface. Pelagic dissolution has been estimated to account for ~ 50% of pelagic aSiO₂ production in the Bay of Brest (Beucher et al., 2004) and > 80% in estuaries reaching the Chesapeake Bay (Anderson, 1986), and is expected to be enhanced by bacterial degradation in estuaries (Roubeix et al., 2008a). Moreover, the constant benthic Chl a:(Chl a+Phae) ratios observed from winter to spring, which decrease during the summer and co-occur with the highest pelagic ratios, suggests enhanced benthic macrofauna grazing (Cariou-Le Gall and Blanchard, 1995), alterations caused by light (Nelson, 1993) and redox conditions (Sun et al., 1993). All of these results confirm the estimation of 50% highly reactive aSiO₂ calculated through an inverse statistical method (Moriceau et al., 2009) and used in our diagenetic model of benthic Si cycle.

As suggested in the previous section, deposition fluxes are tightly related to seasonal changes of river fluxes, detrital loads, benthic and pelagic primary production, and/or lateral transport. Once aSiO₂ is deposited at the sediment-water interface, benthic processes lead either to the return of Si to the water column through dissolution and/or bioirrigation, or to its sequestration in sediments through burial and/or reprecipitation.

The formation of Si(OH)₄ through benthic recycling is known to potentially suffer reverse weathering in benthic sediments (Michalopoulos and Aller, 2004). Indeed, very high reprecipitation has even been reported in high-detrital sediments of the subtropical Amazon delta, in which ~ 90% of initial benthic aSiO₂ has been converted to clay (Michalopoulos and Aller, 2004), or in the Mississippi River Delta, where it accounts for ~ 40% of Si storage (Presti and Michalopoulos, 2008) or in the Barentz Sea (37% at station B13; Ward et al., 2022). In our study, the reprecipitation processes are limited, and reprecipitation fluxes estimated by the model are < 5% of deposition fluxes regardless of station ( Figure 8 ). This is consistent with the values found in the Scheldt Estuary (Rebreanu, 2009). Laboratory experiments are however needed to confirm our model estimates.

The high burial fluxes determined in this study may be linked to 1) high deposition fluxes related to high aSiO₂ fluxes coming from rivers (high concentrations during high river flow; Figures 2C , 4B ), 2) high detrital and aluminium contents, which increase aSiO₂ preservation (Van Cappellen et al., 2002), 3) the macrotidal regime and associated resuspension events, which increase recycling in the water column (Gehlen and Van Raaphorst, 2002), which decreases aSiO₂ reactivity, and 4) bioturbation and more specifically sediment mixing, which increases the transfer of newly settled aSiO₂ deeper in the sediment and increases its preservation, as observed in organic matter (Aller and Mackin, 1984). Such high burial rates have already been observed in highly accumulating zones, e.g., Antarctic sediments in which one third of aSiO₂ deposited accumulates in sediments (Pondaven et al., 2000; DeMaster, 2002). On the contrary, estuaries of large rivers (e.g., Amazon, Congo) have low retention due to export to the coastal margins where deposition and accumulation take place (Michalopoulos and Aller, 2004; Dürr et al., 2011; Raimonet et al., 2015).

In summer, the estimated fluxes are similar to those observed within the Bay of Brest during the productive period, where burial fluxes were estimated to be ~ 32% (Ragueneau et al., 2005a). However, distinct processes cause these similar proportions. In the Bay of Brest, the absence of total recycling is explained by the presence of an invasive filter feeder species, Crepidula fornicata, that increases aSiO₂ preservation in sediments through 1) its incorporation in feces covered by organic matter, and 2) the presence of aluminium-rich sediments (Ragueneau et al., 2005a). In the Elorn and Aulne estuaries, invasive species are absent at the sediment-water interface, and the high burial fluxes are instead explained by the enhanced preservation due to high aluminium and detrital contents (Van Cappellen et al., 2002), high deposition rates (Aller and Mackin, 1984; Conley and Johnstone, 1995), and burrowing depth reaching more than 20 cm in depth.

During the productive period, particularly during summer, the increased proportion of Si(OH)₄ fluxes associated to dissolution and bioirrigation (> 60%; Figure 8 ) compared to burial fluxes is explained by the enhancement of dissolution rates by temperature (Kamatani, 1982; Ragueneau et al., 2002a) as well as deposition of more reactive (autochthonous) material, and bioirrigation (Green et al., 2004).

First, our study confirms the importance to include aSiO₂ when investigating Si retention in estuaries (Carbonnel et al., 2013) or building global silica budgets (Tréguer et al., 2021).

Regardless of the method used, estuarine budgets do not generally explicitly account for the benthic ecosystem in estuaries, due to the assumed small contribution of benthic fluxes compared to that of river fluxes (Arndt et al., 2009). However, benthic sediments may be subject to non-local transport associated with bioirrigation, which is known to strongly increase benthic Si(OH)₄ fluxes at the sediment-water interface (Marinelli, 1992; Forja and Gómez-Parra, 1998), and this is also observed in our study (~50%). Although Arndt and Regnier (2007) have recently analytically resolved benthic processes in a reactive-transport silica model which couples benthic and pelagic processes, methods used to estimate estuarine Si budgets do not vertically resolve the nonlocal benthic processes, e.g., bioirrigation. Contrary to reactive-transport modeling of water column processes, which provides high spatio-temporal resolution, the goal of the present study is to precisely investigate benthic processes in Elorn and Aulne estuary sediments in order to build Si budgets. Our results show the usefulness to incorporate irrigation effects in reactive-transport modeling along estuaries.

During winter floods, the high river flows increase the aSiO₂ and Si(OH)₄ river fluxes to both estuaries as well as the export to the coastal waters ( Table 3 ). This transient and high export of particulate matter, mostly detritic (see section 4.1), has previously been observed in this ecosystem for organic matter (Savoye, 2001). The increase in river flow also decreases the contribution of benthic Si(OH)₄ fluxes (2.7%), aSiO₂ burial (3.7%) and aSiO₂ deposition fluxes (6.4%) compared to aSiO₂ river fluxes. These results confirm that generally low benthic-pelagic coupling occurs during high river flow conditions (Eyre and Ferguson, 2006).

From winter to summer, the drop in the river flow reduces aSiO₂ and Si(OH)₄ fluxes and increases water residence times throughout the estuary. The aSiO₂ deposition becomes higher than river fluxes ( Table 3 ), which highlights that other sources become significant in fueling the sediment-water interface: not only river loads, but also pelagic primary production in the estuary, lateral transport from intertidal sediments or from mixing of marine waters.

First, the aSiO₂ deposition fluxes are enhanced by the increase in pelagic primary production in summer. This is highlighted by the high Chl a concentrations measured in May and even more in July (from 6 to 35 µg l^-1 along the salinity gradient; Figure 4C ). Note that these concentrations are higher than those in the Bay of Brest (generally < 10 µg l^-1 even during blooms, except a maximum of 21 µg l^-1 observed in 1981; Chauvaud et al., 2000), and that even higher concentrations have already been measured in these estuaries during summer (88 and 46 µg l^-1; Savoye, 2001). In this study, the extrapolation of the pelagic primary production measurements in the Aulne Estuary allows us to roughly estimate a total pelagic primary production of ~ 57-225 10³ mol d^-1. This estimate is slightly lower than the production of 448 10³ mol d^-1 estimated in the Bay of Brest (Ragueneau et al., 2002a). However, this estimation is close to the deposition flux (107 10³ mol d^-1; Table 3 ), which consequently appears to be related to pelagic production. Since dissolution also occurs in estuarine waters (accounting for 20-80% of pelagic aSiO₂ recycling; DeMaster et al., 1983; Ragueneau et al., 2002b; Beucher et al., 2004), it follows that other sources must also contribute to the deposition flux of aSiO₂. This exercise highlights the importance of directly quantifying the pelagic silica dissolution rates in estuaries to better constrain their budgets. In the Bay of Brest, pelagic dissolution rates were shown to recycle ~ 50% of aSiO₂ production annually (Beucher et al., 2004), but these rates could be even higher throughout estuaries because of enhanced aSiO₂ dissolution by bacterial degradation of organic coatings (Roubeix et al., 2008a).

As suggested in section 4.1., other possible sources of deposited aSiO₂ are benthic primary production, lateral resuspension and redistribution of particles from tidal mudflats and small water outlets. The high deposition fluxes of estuarine materials are particularly highlighted by benthic Chl a concentrations (140 µg gDW^-1; Figure 5B ) that are much higher than those already reported in the Bay of Brest (Ragueneau et al., 1994; Sagan and Thouzeau, 1998; Ni Longphuirt et al., 2006). Such high Chl a (or fucoxantin) concentrations have previously been reported to mirror the higher pelagic and benthic production taken place in estuaries, the lateral resuspension and redistribution of particles, with the higher deposition rates related to a shallower water column (Sun et al., 1994; Mangalaa et al., 2017; Wallington et al., 2023). Marine loads might be low due to the absence of a connection between estuarine and coastal waters observed in the Elorn and Aulne estuaries (Savoye, 2001), except in the downstream parts of these estuaries where marine materials are more significant (Khalil et al., 2013). Marine aSiO₂ loads can however be much higher in some estuaries; for instance, it reaches 25% in the Scheldt estuary (Carbonnel et al., 2013).

The relative contribution of benthic Si(OH)₄ fluxes compared to river fluxes is consequently enhanced ( Figure 8 ), which highlights the important role of benthic recycling to the pelagic ecosystem, especially in summer. Even if benthic fluxes are generally lower compared to river fluxes in large estuaries (Arndt et al., 2009), their contribution can become significant in small and shallow estuaries (Anderson, 1986), such as in the Elorn and Aulne, as well as during the low-flow season (this study). The contribution of benthic fluxes to the pelagic ecosystem in summer is higher in the Aulne than in the Elorn Estuary, not just due to the lower river flow, but also due to the higher benthic surface of the Aulne Estuary and its meandering shape (Raimonet et al., 2013a).

Finally, this work bring new data on Si retention in small macrotidal estuaries, at seasonal and annual scales. While retention is ≤ 5% of river Si flux in winter, it increases to 38% in summer, and even 67% when accounting for burial in intertidal mudflats ( Table 3 ). The annual burial of Si in the Elorn and Aulne estuaries is estimated to be 7% in subtidal sediments, and 12% when accounting for intertidal mudflat burial. These estimates are similar for the two estuaries, but lower than a ten times bigger estuary, the Scheldt estuary, where annual retention of Si(OH)₄ and aSiO₂ are estimated to attain 28% and 64%, respectively (Carbonnel et al., 2013). Such difference is expected because of (1) higher water residence time in the Scheldt estuary, which increases retention and limits flushing, and (2) intense dredging activity, which reduce sediment export to the coastal zone (Carbonnel et al., 2013). In this study, we gathered and measured a large number of data to better constrain and limit uncertainties on our estimates. Our values could however be refined, e.g. by coupling our results with a hydro-sedimentary model e.g. Grasso et al. (2018) in order to quantify retention at a fine spatial and temporal scale. This could help in quantifying the spatial heterogeneity and the transient regime of erosion and export of estuarine sediments.