# Two dimensional mapping of iron release in marine sediments at submillimetre scale

Thibault De Chanvalon A<fnref n="1" />,<fnref n="2" />,*, Metzger E<fnref n="1" />, Mouret A<fnref n="1" />, Knoery Joel<fnref n="2" />, Geslin E<fnref n="1" />, Meysman Fjr<fnref n="3" />

## Abstract :

Coastal and shelf sediments are considered as an important source of dissolved iron to the ocean. Here, we present a new numerical approach to estimate geochemical fluxes and production rates in an estuarine sediment at sub-millimetre resolution. This approach is based on application of Savitsky-Golay filter (SGF) procedure to two-dimensional concentration distributions of dissolved iron. We verified the procedure by applying it to artificial data of known production rates, and analysed the resulting uncertainty on production rates and fluxes across the water-sediment interface. This SGF procedure was applied to data from an intertidal mudflat that is densely inhabited by macrofauna (e.g. 630 ind m⁻<fnref n="2" /> of Hediste diversicolor, I. Métais, pers.com.). Our analysis reveals an apparent recycling rate of 3780 ± 1399 µmol m⁻<fnref n="2" /> d⁻<fnref n="1" /> and a mean residence time of iron in the dissolved phase of 2.3 days. Visual identification of burrows permitted to calculate separately the diffusive flux across the sediment-water interface (104 ± 20 µmol m⁻<fnref n="2" /> d⁻<fnref n="1" />) and the bio-irrigational flux (410 ± 213 µmol m⁻<fnref n="2" /> d⁻<fnref n="1" />). Reactive iron particles will undergo on average 7.4 cycles of dissolution/precipitation before being released to the water column. These results show that estuarine sediments support intensive iron recycling that has probably a large impact on terrigeneous particles before being released into the ocean.

## Highlights

* New method for benthic flux and production rate calculation in 2 dimensions
* Bio-irrigation accounts for more than 80% of iron benthic flux in estuarine mudflat.
* Dissolved iron time residence within estuarine sediment of 2.3 days
* Reactive iron particles undergo 7.4 cycles of dissolution/precipitation before being released to the water column.

**Keywords** : Coastal sediment, Benthic flux, Bioturbation, 2D-DET

<footnote n="1">¹ Université d'Angers, LPG-BIAF, UMR CNRS 6112, 49045 Angers Cedex, France</footnote>

<footnote n="2">² Ifremer, LBCM, Rue de l'Ile d'Yeu, 44300 Nantes, France</footnote>

<footnote n="3">³ Department of Ecosystem Studies, The Netherlands Institute of Sea Research (NIOZ), Korringaweg 7, 4401 NT Yerseke, The Netherlands</footnote>

<footnote>* Corresponding author : A. Thibault De Chanvelon, email address : athibaultdc@gmail.com</footnote>
# 1. Introduction

Dissolved iron acts as a limiting micronutrient controlling oceanic primary production, and in this way, is considered to play an important role in the carbon cycle in the ocean (Martin and Fitzwater, 1988). However, the sources of dissolved iron to the ocean and the mechanisms controlling iron speciation and bioavailability are still weakly constrained, due to the important temporal variability in the iron inputs (Johnson et al., 1999), the spatial heterogeneity of the different sources (Boyd and Ellwood, 2010; Conway and John, 2014; Raiswell, 2006) and the complex speciation and biogeochemical interactions that take place once iron reaches the ocean (Chever et al., 2015; Homoky et al., 2012; Rose and David Waite, 2007).

Coastal and shelf sediments are considered as an important source of dissolved iron to the ocean, with a source strength comparable to riverine inputs and hydrothermal vent emissions (Chever et al., 2015; Raiswell and Canfield, 2012). However the uncertainty on the efflux from near-shore sediments is large, because experimental quantification of the dissolved iron released from sediments is challenging. One important factor that complicates iron flux measurements in coastal and shelf sediments is the presence of benthic fauna. In these sediments, the efflux of a solute out of the sediment is no longer solely driven by diffusion across the sediment-water interface, but results from a large part from bio-irrigation, i.e., the enhanced sediment-water exchange due to burrow pumping activities by infauna (Aller, 2001; Meysman et al., 2006). As burrows have complex, three dimensional geometries, this bio-irrigational flux component cannot be quantified with classical methods such as one-dimensional pore water profiling. Flux measurements based on benthic chambers (e. g. Sundby et al., 1986; McManus et al., 1997; Noffke et al., 2012) can be biased by iron reactivity after release in the water column and therefore can underestimate the efflux reduced iron due to oxidative precipitation in the chamber (Pakhomova et al., 2007; Severmann et al., 2010). Furthermore, the enclosure of sediments can disturb the natural bio-irrigation activity
of fauna, thus influencing the bio-irrigational fluxes of dissolved iron. Finally, sedimentary iron release can also be constrained by the analysis of vertical distributions of iron concentrations and isotopes in the water column. However, interpretation of such isotopic data is often challenging, as it is dependent on numerous and still poorly understood transformations in the water column, which can have a large impact on isotopic signature (Chever et al., 2015; Homoky et al., 2013).

Coastal and shelf sediments are also characterized by a large internal cycling of iron compounds between reduced and oxidized forms (Raiswell and Canfield, 2012). These redox oscillations are likely to control the speciation of iron (i.e. the relative abundance of Fe(II) and Fe(III) in both pore water and solid phase) and therefore, they can impact the bioavailability of the iron that is released to the coastal ocean. Traditionally the behaviour of freshly precipitated iron oxide is described by Ostwald's rule (Steefel and Van Cappellen, 1990), which implies an increase in cristallinity and an associated decrease in bioavailability as the iron minerals age (Bligh and Waite, 2011; Kuma et al., 1992; Yoshida et al., 2006). However, this unidirectional “aging” may not be representative for coastal sediments, where intensive redox cycling takes place. As a result of the reductive dissolution and subsequent oxidative precipitation, smaller and more reactive (nano)particles are formed (Magnuson et al., 2001; Taillefert et al., 2000; Von Gunten and Schneider, 1991). Redox cycling hence can have an important impact on the composition and reactivity of iron minerals in the sediment. Raiswell (2011) assumes that each redox cycle simply rejuvenates iron oxides and thus increases the bioavailability, whereas Thompson et al. (2006) reports an increase of iron oxide crystallinity after redox cycling due to the preservation of the most crystalline phase during each dissolution cycle. Moreover, new authigenic minerals, such as Fe-P precipitates (Egger et al., 2015; Hyacinthe and Van Cappellen, 2004; Senn et al., 2015), can be formed after reductive dissolution, and so the redox cycling of iron may also regulate the fluxes, transformation and burial of other elements (e.g. phosphorus, trace metals) in the coastal marine environment.
Yet, at present, there has been no direct estimation of the recycling intensity (i.e. the number of reduction/oxidation cycles and Fe atom experiences before being buried) in coastal sediments. This is because these redox transformations occur at millimetre scale, and are often localized near ventilated burrows of fauna.

In the present study, we propose a new technique to estimate the flux of iron across the sediment-water interface as well as its internal recycling intensity. The technique is based on the numerical analysis of 2D distributions of dissolved iron concentrations in coastal sediments, which are obtained at high resolution (200µm) by the technique of Dissolved Equilibrium in Thin film (DET; Davison et al., 1991), combined with planar colorimetry (2D-DET; Jézéquel et al., 2007; Robertson et al., 2008; Shuttleworth et al., 1999; Zhu and Aller, 2012).
## 2. Material and methods

### 2.1. Two-dimensional sampling of dissolved Fe

Pore water data were obtained from the intertidal mudflat Les Brillantes in the Loire estuary (France; 47°16′56.00″N 2° 3′47.00″W). Data collection (02/05/2013) occurred during a period of high river discharge, when fresh water penetrates far downstream into the estuary, and hence, water with a low salinity (range 2-6) inundates the mudflat at high tide. Two-dimensional (2D) distributions of total dissolved iron ($Fe_d$) were obtained by insertion of a 2D-DET probe into the sediment. The 2D-DET technique quantifies the distribution of a solute in a vertical plane by diffusive equilibration of the pore water with the gel. These $Fe_d$ distributions were complemented by photographs of the sediment section that faced the 2D-DET probes, which then allowed to correlate the pore water data with the presence of burrows and sediment structures. These sediment photographs were obtained by introduction of a so-called “jaw device” into the sediment, which recovers the 2D-DET probe from the sediment, together with the sediment slice that faces the gel probe (Fig. 1a); a detailed description of this instrument is given in Thibault de Chanvalon et al., (2015). Overall, the sampling procedure provides two separate images: (1) a normal bright field image of the face plane of the sampled sediment, which allows the visual identification of burrows and other sediment structures (Fig. 1b), and (2) a false-coloured image generated from the 2D-DET gel probe, which reveals the distribution of dissolved iron in the pore water in two dimensions and at high resolution (Cesbron et al., 2014; Jézéquel et al., 2007).

The 2D-DET probe is composed out of a polycarbonate plate with a central depression of 1 mm filled with a polyacrylamide gel (pore size ~2 nm; Zhang and Davison, 1999) that was prepared according to Jézéquel et al., (2007). The plate is subsequently covered by a protective porous PVDF membrane (pore size 0.2 µm) that protects the gel and prevents it from falling out the central depression. The 2D-DET probe was prepared with one week of sampling and was conserved in the dark at 4°C in a moist plastic bag (humidity prevents the
gel from drying out). Before deployment, the plastic container bag was deoxygenated by N₂ bubbling for 6 hours. After insertion into the sediment at the beginning of low tide, the 2D-DET probe was deployed for 5 hours during emersion, and recovered before the flood water came in. All deployments occurred during day time. The polycarbonate frame of the 2D-DET probe supports two vertical rails that can guide a stainless steel plate to slide the jaw device into the sediment. This way, at the end of the deployment, the 2D DET probe could be retrieved along with the 2 cm thick sediment layer that was in contact with the gel probe (Fig. 1a)

Upon retrieval, the 2D-DET probe was carefully separated from the sediment, covered with a plastic-coated aluminium plate (to ensure a homogenous freezing process), and stored in a cooling box with dry ice pellets for transport. A picture of the face plane of the sampled sediment was taken with a compact digital camera (Olympus Stylus 7000). Upon arrival in the laboratory, the 2D-DET probe was stored in a freezer (-18°C) until further analysis. This preservation procedure allows storage for up to three months without concentration gradient smoothing due to diffusion within the gel (Cesbron et al., 2014). During analysis, the 2D gel probe was thawed at ambient temperature, and covered by a reactive gel impregnated with specific colorimetric reagents for total dissolved iron (Fe<sub>d</sub>). The preparation of the reactive gel was modified from Cesbron et al. (2014) and Jézéquel et al. (2007), and consists of a 0.46 mm thick polyacrylamide gel incubated during 1 hour in a reactive solution (ascorbic acid 30 mM, sulfuric acid 0.558 M, potassium antimony(III) tartrate hydrate 0.40 mM, ammonium molybdate tetrahydrate 18.5 mM and ferrozine 12.2 mM). This reactive solution is sensitive to both dissolved reactive phosphorus (DRP; results not used in this study) and total dissolved iron (Fe<sub>d</sub>; Fig. 1c). Iron complexation by the colorimetric reagent is achieved in less than 5 minutes, while DRP coloration takes about 20 minutes (Cesbron et al., 2014). After 25 minutes of contact time, an image of the superposed gels (Fe/P reactive gel + sample gel) was taken with a hyperspectral camera (HySpex VNIR 1600) thus producing a two-dimensional
reflectance map. A subsequent spectral analysis of the data allows to discriminate the DRP
signal from the Fe<sub>d</sub> signal (Cesbron et al., 2014). Eight standard solutions of iron (II) were
prepared (from 0 to 200 µM) and exposed to an identical sampling protocol (diffusion in a
polyacrylamide gel followed by colorimetric reaction, hyperspectral imaging and spectral
analysis). This calibration dataset was used to convert the acquired two-dimensional
reflectance map into a Fe<sub>d</sub> concentration map (Fig 1d). The resulting uncertainty on the Fe<sub>d</sub>
concentration has been estimated as < 10% (Thibault de Chanvalon et al., 2015). The spatial
resolution of the Fe<sub>d</sub> analysis is determined by the pixel size of the hyperspectral camera (200
µm), the ability of the sediment to buffer the pore water concentration (Harper et al., 1997)
and the relaxation effect induced by diffusion, once the DET probe boundary is no longer
constrained by the sediment (~3 millimetres in 5 minutes; Davison et al., 1994).

Temperature (T) and salinity (S) were measured in the first centimetre of the sediment by
insertion of a WTW Series 3110 conductivity meter. For the purpose of this study, depth
variations in S and T were not considered. A porosity depth profile was determined on a
sediment core (plastic core liner; 3 cm inner diameter) that was retrieved 1 m from the
location of the 2D-DET probe deployment. This sediment core was sliced at 2 mm intervals
over the first 2 cm, and using 5 mm intervals up to 5 cm depth. For each depth layer, the
porosity was calculated from the mass difference before and after freeze drying, accounting
for the salt content of the pore water.
## 2.2. Production rate estimates

Both fluxes and production rates of Fe<sub>d</sub> were determined from the recorded 2D distributions based on a two-dimensional reactive transport analysis. From a theoretical point of view, the 2D distribution of a solute in an aquatic sediment can be described by the mass balance equation (Berner, 1980; Boudreau, 1997; Meysman et al., 2005b):

$$ \Phi \frac{\partial C}{\partial t} = -\operatorname{div} \mathbf{J} + \mathbf{R}(y, z) \quad (1) $$

In this expression, C denotes the concentration, Φ is the porosity, and R represents the production rate, i.e., the net balance of all production and consumption due to biogeochemical reactions. The term div J = ∂J<sub>y</sub>/∂y + ∂J<sub>z</sub>/∂z represents the divergence of the flux vector J = (J<sub>y</sub> / J<sub>z</sub>).

Concentrations, fluxes and rates are dependent on the z-coordinate representing depth into the sediment, and the y-coordinate representing the lateral distance. This flux vector can be specified as (Meysman et al., 2005a):

$$ J = -\Phi D_s \operatorname{grad} C + \Phi C v \quad (2) $$

$$ \Leftrightarrow \left\{ \begin{array}{l} J_y(y, z) = \Phi C v_y - \Phi D_s \frac{\partial C}{\partial y} \\ J_z(y, z) = \Phi C v_z - \Phi D_s \frac{\partial C}{\partial z} \end{array} \right. $$

where D<sub>s</sub> is the effective diffusion coefficient of Fe<sub>d</sub> in the pore water, and v = (v<sub>y</sub> / v<sub>z</sub>) is the advective velocity vector of the pore water. If we assume that (1) advective transport of pore water is negligible (see discussion below) and (2) that the solute distribution is at steady state, we obtain the following expression for the production rate:

$$ R(y, z) = R_y + R_z $$

$$ R(y, z) = \left[ -\Phi D_s \frac{\partial^2 C}{\partial y^2} - \left( \Phi \frac{\partial D_s}{\partial y} + D_s \frac{\partial \Phi}{\partial y} \right) \cdot \frac{\partial C}{\partial y} \right] + \left[ -\Phi D_s \frac{\partial^2 C}{\partial z^2} - \left( \Phi \frac{\partial D_s}{\partial z} + D_s \frac{\partial \Phi}{\partial z} \right) \cdot \frac{\partial C}{\partial z} \right] \quad (3) $$

Equation (3) reveals that it is possible to estimate the total production rate R(y, z) at each point as the sum of two separate production terms: the first term R<sub>z</sub> is based on the vertical
variation of the concentration along the z-axis, and the second term $R_y$ is solely based on the
horizontal variation of the concentration along the y-axis.

The neglect of the advective terms in Eq. (3) requires some consideration. The
sediment at the sampling location is a bioturbated, cohesive sediment with a high porosity ($> 0.8$). Accordingly, the permeability is low, and advective flows in the pore water due to tidal
pumping and wave action will not be important (Berner, 1980; Glud et al., 1996; Huettel et
al., 2014). Yet, another important type of advection is pore flow inside the burrows. Due to
burrow flushing by macrofauna, an advective flow will be created inside the burrow. These
burrow flows can be highly variable in time, as the burrow irrigation intensity changes on the
order of minutes (Kristensen, 1989; Riisgård and Larsen, 2005; Volkenborn et al., 2012,
2010). The dominant species on Les Brillantes mudflat is the polychaet *Nereis diversicolor*,
which shows highly frequent (5-10 minutes) intermittent ventilation (Kristensen, 2001, 1981;
Pischedda et al., 2012). As our DET deployments lasted for 5 hours, they will provide a time-
averaged picture of this burrow ventilation activity (Harper et al. 1997). Over this time-scale,
the mean transport of $Fe_d$ typically implies that dissolved iron diffuses into the burrow across
the burrow wall, and is subsequently flushed out by burrow ventilation. By ignoring the
advective terms in Eq. (3), the effect of $Fe_d$ removal by burrow flushing will show up as a
reactive term R. Therefore, inside actively flushed burrows, the net production rate R will
likely be negative, and should be not be interpreted in terms of geochemical reactions, but in
terms of bio-irrigation. So overall, the neglect of advection in our reactive-transport analysis
can be justified, provided that the net production term R is correctly interpreted. Outside the
burrows, R represents the effect of biogeochemical reactions. Inside the burrows, R
predominantly represents the removal of solutes due to burrow ventilation.

The assumption that the solute distribution is at steady state is more difficult to justify
*a priori*. There are various processes acting in intertidal mudflats that can drive the pore water
into a transient state. These processes also work on different time scales, e.g.,
microphytobenthos activity (diurnal scale), the deposition of a new sediment layer during floods (month scale) or variations in mineralization intensity (seasonal). Processes that act on longer time scale than the 2D-DET measurement (5 hours) will not provide large contribution; the pore water iron distribution will attain a quasi-steady state. Still, important transient concentration changes can occur on shorter time scales in bioturbated sediments, due to variations in the intensity of burrow flushing. Ferrous iron concentrations (and other reduced species) may increase when burrow ventilation is stopped, and rapidly decrease, when burrow ventilation restarts. So especially near burrows, we cannot rule out concentration change terms $dC/dt$, and in our analysis, these terms will be lumped into the production rate $R$. For this reason, the quantity $R$ should be regarded as an “apparent production rate”, which hence may also include transient effects. The assessment of the importance of such concentration change terms $dC/dt$ is beyond the present study, and should be the focus of future work.

Equation (3) reveals what information is needed to calculate the local production rate of dissolved iron: the spatial distributions profiles of the porosity $\phi$, the effective diffusion coefficient $D_s$, and the concentration $C$, as well as the spatial derivatives of these quantities. Typically, the effective diffusion coefficient $D_s$ is calculated from the relation $D_s=D_{mol}/\theta^2$, where the molecular diffusion coefficient $D_{mol}$ is calculated as an empirical function of temperature and salinity (Boudreau, 1997). Here, we assume that salinity and temperature do not vary in the sediment, and hence, the molecular diffusion coefficient $D_{mol}$ will not be spatially dependent. The molecular diffusion coefficient of $Fe_d$ was calculated from the R package "marelac" ($D_{mol} = 5.16 \ 10^{-10} \ m^2 \ s^{-1}$ at T= 12.7 °C and salinity 5.6). The tortuosity $\theta^2$ is typically taken as a function of porosity, as in the modified Weissberg relation, $\theta^2=1-2\ln(\phi)$ (Boudreau, 1996), and so, once the spatial distribution of the porosity is known, the associated tortuosity can be calculated analytically. The porosity in the vertical is usually assumed to be mainly controlled by compaction, while the porosity can vary in the lateral around burrows linings. So both vertical and lateral gradients in porosity and tortuosity can
exist. Presently however, such variation of the porosity at sub-millimetre scale cannot be quantified experimentally. Accordingly and to a first approximation, the lateral variation in porosity is neglected here, and only the vertical gradients are accounted for. The porosity is described by an exponential decreasing function of depth $\phi = (\phi_0 - \phi_{inf}) e^{-\alpha z} + \phi_{inf}$, which was fitted to the measured porosity depth profile and the derivatives $\frac{\partial\phi}{\partial y}$ and $\frac{\partial D_s}{\partial z}$ were subsequently calculated analytically. As a result, the only unknowns that remain in equation (3) are the spatial derivatives of the concentration (i.e. $\frac{\partial c}{\partial y}$, $\frac{\partial c}{\partial z}$, $\frac{\partial^2 c}{\partial y^2}$ and $\frac{\partial^2 c}{\partial z^2}$). We will now investigate how these derivatives can be derived from the available dataset.
2.3. The Savitzky Golay filter procedure

The Savitzky Golay filter (SGF) procedure (Savitzky and Golay, 1964) is a widely
used method for smoothing noisy experimental data and calculating associated derivatives.
Briefly, for each point within a data series, the SGF procedure generates a polynomial
function that is fitted by the method of least-squares to a central point and its n closest
neighbours on either side (i.e., on 2n+1 points, the window size). Then the SGF procedure
estimates the value (and its derivatives) at the central point by means of evaluating and
differentiating the fitted polynomial function. The standard implementation of the SGF filter
is described in Gorry, (1990) and the algorithm is given in Press et al., (1992). Here, we have
implemented the routines in the open-source programming framework R, which use the
function "sgolayfilt" in the extension package "signal" (see Meysman and Burdorf (in prep)
for details).

Two parameters have to be fixed in the SGF procedure: the order (p) of the fitted
polynomial and the filter window size (2n+1). The polynomial order is typically set at p = 3,
as is done here. The larger the filter window size, the more noise is suppressed, but also the
more the data signal is smoothed. Accordingly, an optimal window size suppresses most of
the (undesired) experimental noise, while it removes as little signal as possible. The
automated (numerical) determination of this optimal window size still forms an active field of
research. Meysman and Burdorf (in prep) have proposed an automated procedure based on
fractal analysis which calculates the fractal dimension D as a function of the window size
using the R function "fd.estimate" from the extension package "fractaldim", selecting the
"madogram" method to calculate the fractal dimension (Gneiting et al., 2012). As the window
size increases, the fractal dimension D decreases, and there is a clear break in the slope of D
when it reaches the value of 1. This break point is selected as the optimal window size. More
details on this automated procedure to determine the optimal window size are given in
Meysman and Burdorf (in prep).
Here we implemented the one-dimensional SGF procedure as proposed by Meysman and Burdorf (in prep) (i.e. with polynomial order p = 3 and optimal window size selection) to both the rows and columns of the two-dimensional concentration distribution of Fe<sub>d</sub>. Application to the columns (z-coordinates) estimates the concentration derivatives in the vertical direction, thus yielding the vertical flux J<sub>z</sub>(y, z) and the production rate R<sub>z</sub>(y,z) in each point of the concentration matrix. Similarly, application to the rows (y-direction or horizontal) yields the horizontal flux J<sub>y</sub>(y, z) and the associated production rate R<sub>y</sub>(y,z). As derived above, the total production rate in any given location is then simply obtained as the summation of both production rates, i.e., R(y,z) = R<sub>y</sub>(y,z) + R<sub>z</sub>(y,z). Overall, the end product of our reactive transport analysis is a two-dimensional matrix of production rates, which is then visualized as a false-colour image at the same high resolution as the original concentration dataset.
## 2.4. Model validation procedure

To verify the performance of our SGF procedure, we applied it to a synthetic dataset, which was first created by a forward numerical simulation of the reactive transport model. In a first step, various types of “idealized” or “true” two-dimensional patterns of production rate were defined, i.e., $R^{ideal}(y,z)$. The most basic production rate pattern consisted of a localized production spot combined with a uniform first order consumption (R= - k C) throughout the whole model domain. The localized production was modelled by a symmetrical two-dimensional Gaussian function

$$f(y, z) = R_{max} \exp \left( - \left( \frac{(y - y_0)^2}{2 \sigma^2} + \frac{(z - z_0)^2}{2 \sigma^2} \right) \right) \quad (4)$$

where the coefficient $R_{max}$ is the amplitude, $y_0$ and $z_0$ are the coordinates of the centre, and $\sigma$ is the spreading parameter which is alternatively quantified by the full width at half maximum, i.e. FWHM = 2.355 $\sigma$. More complex production patterns were constructed by combining multiple of these localized Gaussian production sites (with different amplitudes and spreads) with a uniform first order consumption. The latter generated a distribution that was roughly similar to the natural production rate pattern observed in the field data.

In a second step, the “idealized” 2D production patterns $R^{ideal}(y, z)$ were used as input into a reactive transport model, and the resulting simulations thus generated the “idealized” 2D concentration distribution $C^{ideal}(y,z)$ and the fluxes $J_x^{ideal}(y, z)$ and $J_y^{ideal}(y, z)$. To this aim, the governing 2D mass balance equation (Eq. 1) was numerically solved for steady state using reactive transport model procedure as fully detailed in Soetaert et al. (2010) and Soetaert and Meysman (2012). This numerical solution procedure was implemented in the open-source programming language R.

In a third step, random white noise (i.e. with a Gaussian distribution) was added to these idealized concentration distributions to mimic the random noise generated during experimental data collection. This then provided to obtain suitable “test” concentration
distributions. The applicable noise level was determined by analysis of the 2D-DET calibration standards. Around 5000 concentration data points from each of the eight 2D-DET calibration standards were analyzed and plotted. The resulting distribution was fitted by a normal distribution N (0, var) (Gaussian distribution) with a zero mean. The variance of this normal distribution then determined the applicable noise level. This noise level was quantified by two parameters: the relative standard deviation of the highest concentration standard (rsd<sub>high</sub> = 1%) and the absolute standard deviation of the lowest concentration standard (sd<sub>low</sub> = 0.4 µmol L<sup>-1</sup>). The "test" dataset was subsequently calculated according to:

$$C_{test} = C_{ideal} (1 + C_{high}^{err}) + C_{low}^{err} \quad (5)$$

where $C_{low}^{err} \sim N(0, sd_{low})$ and $C_{high}^{err} \sim N(0, rsd_{high})$

In a fourth and final step, the SGF filter procedure was applied to $C_{test}$. The synthetic datasets were treated in the exactly same manner as the real datasets. The result is an estimated 2D concentration distribution $C^{estimated}(y, z)$, a vertical and a horizontal flux called $J_z^{estimated}(y, z)$ and $J_y^{estimated}(y, z)$, and an estimated 2D production pattern called $R^{estimated}(y, z)$.

The *reconstruction bias* was calculated as the difference between idealized and reconstructed quantities, i.e.:

$$\Delta C(y, z) = C^{ideal}(y, z) - C^{estimated}(y, z) \quad (6)$$

$$\Delta J_x(y, z) = J_x^{ideal}(y, z) - J_x^{estimated}(y, z) \quad (7)$$

$$\Delta R(y, z) = R^{ideal}(y, z) - R^{estimated}(y, z) \quad (8)$$

These reconstruction bias values are calculated at every individual pixel. To summarize the overall accuracy of the model reconstruction, a single Mean Reconstruction Bias (MRB) is introduced
$$ \text{MRB (\%)} = \frac{\text{sd}(\Delta C(y, z))}{\text{sd}(C^{\text{ideal}}(y, z))} * 100 \quad (9) $$

This indicator requires that we calculate the standard deviation $sd(C^{\text{ideal}}(y, z))$ of the sampling distribution of the concentrations over the domain, as well as the standard deviation $sd(\Delta C(y, z))$ of the sampling distribution of the concentration reconstruction bias. As a result, the Mean Reconstruction Bias provides a single indicator for a whole 2D domain, summarizing the misfit between the “true” concentration values and the “reconstructed” concentration values. A similar MRB indicator is defined in an analogous way for the 2D reconstruction of fluxes and production rates.
# 3. Results

## 3.1. Model analysis of synthetic datasets

To verify the performance of the reactive transport analysis based on the SGF procedure, it was applied to three synthetic datasets. Each dataset has the same model domain (size = 8.7 cm x 17 cm), but a different production pattern of Fe<sub>d</sub>: (1) a single production zone with a broad spreading (FWHM = 2.36 cm, pixel size = 200 µm, R<sub>max</sub> = 33 µmol L<sup>-1</sup> d<sup>-1</sup>, y<sub>0</sub>= 4.35 cm z<sub>0</sub>= 8.5 cm), (2) a single and highly localized production zone (FWHM = 0.59 cm, pixel size = 200 µm, R<sub>max</sub> = 286 µmol L<sup>-1</sup> d<sup>-1</sup>, y<sub>0</sub>= 4.35 cm z<sub>0</sub>= 8.5 cm), and (3) a superposition of 25 localized production zones (FWHM range = 0.12 to 5.89 cm, pixel size = 200 µm, R<sub>max</sub> range = 2.6 to 13 µmol L<sup>-1</sup> d<sup>-1</sup>) thus providing a patchy consumption distribution. To enable a suitable comparison, the idealized production rate R<sup>ideal</sup>(y,z) was scaled in such a way in all three cases, that the resulting 2D concentration distribution showed a maximum concentration of 100 µmol L<sup>-1</sup>.

Figure 2 shows the results for the concentration, i.e., (a) the idealized 2D concentration distribution C<sup>ideal</sup>(y, z) as generated by the reactive transport model (Fig. 2a) and (b) the 2D distribution of the “estimated” concentration C<sup>estimated</sup>(x, y) as recovered by the SGF procedure from the purposely noise-compromised data (Fig. 2b). The difference between “idealized” and “estimated” concentration (ΔC(y, z); Fig. 2c) is vanishingly small for all three synthetic datasets, thus showing that the SGF procedure can accurately reproduce the idealized concentration patterns. The maximum reconstruction bias of 2.5 µmol L<sup>-1</sup> occurs at the maximum concentration, i.e., at the peak of the localized production zone. When comparing the concentration profiles extracted from the middle section of the model domain, there are no visible differences between the two vertical profiles of C<sup>ideal</sup> and C<sup>estimated</sup> (Fig. 2d). The Mean Reconstruction Bias (MRB) ranges from 0.6 to 2.2 % over the three datasets (Table 1, first column).
Figure 3 shows the results of the vertical flux estimation (horizontal fluxes lead to same conclusions and are not shown). Figure 3a shows the original fluxes $J_z^{\text{ideal}}(y, z)$, where the blue and red colours represent upward and downward fluxes respectively. As expected for datasets 1 and 2, there are two localized zones of high fluxes on either side of the production maximum, confirming that dissolved iron is diffusing out the localized production zone in opposite directions. Because fluxes are based on first order derivatives of the "noisy" concentration, these flux distributions are *a priori* more difficult to reconstruct than concentration distributions (the higher the order of the derivatives, the more sensitive to the noise). Figure 3b depicts $J_z^{\text{estimated}}(y, z)$ and the reconstruction bias $\Delta J_z(y, z)$ is shown in the Fig. 3c. Vertical profiles (Fig. 3d) show small differences between $J_z^{\text{ideal}}(y, z)$ and $J_z^{\text{estimated}}(y, z)$. The Mean Reconstruction Bias Indicator ranges from 3.8% to 7.2% (Table 1, second column).

The 2D distribution of the “idealized” production rate $R^{\text{ideal}}(y, z)$ is shown in Fig. 4a for the three synthetic datasets. An accurate reconstruction of $R^{\text{estimated}}(y, z)$ (Fig. 4b) is more challenging than for fluxes, as it is based on the second-order spatial derivatives of a "noisy" concentration profile. The higher the order of the derivatives, the more susceptible their numerical reconstruction is to noise. The MRB (Table 1, third column) is 12.2% for the broad production patch (dataset n°1) and 31.9% for the narrow production patch (dataset n°2; table 1). These results illustrate the difficulty to simulate narrow localized production zones that generate steep gradients in concentrations. The third dataset presents a higher MRB of 34.8%. Yet overall, the reconstructed production rate pattern (Fig. 4b) closely resembles the original pattern (Fig. 4a). The SGF procedure is able to correctly localize all production patches, even those production zones that are not immediately evident from the concentration dataset. The narrow production zone in synthetic dataset n°3 showed the highest reconstruction error ($\Delta R(y, z)$; Fig. 4c). The vertical transect in Fig. 4d showed a good
agreement between R<sup>ideal</sup>(y, z) and R<sup>estimated</sup>(y, z), but also revealed noise that is not
completely smoothed out by the SGF procedure.

<figure><img src="image_1.png" /></figure>
## 3.2. Model analysis of the field dataset

The measured two-dimensional distribution of dissolved iron in the pore water of the Loire estuary mudflat is shown in Fig. 5, alongside the one-dimensional depth profile in porosity. Overall, the sediment is highly porous (mean porosity of 0.86 in the top 5 cm), indicating that fine-grained cohesive mud is accumulating at the field site. However, this porosity depth profile shows a clear subsurface minimum at 0.9 cm depth (Fig. 5a), which is most likely the result of a recent sedimentation event, where a layer of more coarse grained material has been deposited in between two layers of fine-grained sediment. Below this subsurface negative peak, the porosity depth profile gradually declines with depth, as expected from a compacting cohesive sediment. To keep the analysis tractable, we did not attempt to reconstruct the subsurface maximum of porosity, but instead, we captured the overall porosity decrease by fitting an exponential function to the measured porosity depth profile (Fig. 5a; $\alpha=0.39$, $\phi_0=0.93$ and $\phi_{inf}=0.81$).

Figure 5b shows the 2D $Fe_d$ distribution that was experimentally retrieved at the field site. This 2D profile was cropped at the left, right borders (by 1 cm) and at the bottom (by 0.8 cm) to remove boundary effects, such as lateral diffusion from the pore water (Harper et al., 1997) or partial desiccation of the gel during experimental handling. The sediment-water interface was estimated from graduations drawn on the DET-2D probe frame during deployment and the overlying water section of the 2D profile was also removed (the final model domain was 7.7 cm x 16.2 cm). Overall, the 2D distribution of $Fe_d$ concentration was highly heterogeneous, with clear zones of iron accumulation interspersed with zones of iron depletion. Apart from the 3 first millimetres near the sediment-water interface, which formed a layer uniformly depleted in iron, no vertical stratification was observable. The dominant features were two vertically elongated "ridges" (width ~2 cm; length ~10 cm) showing high $Fe_d$ concentrations (maximum 150 -200 µmol L⁻¹) that extended between 5 and 15 cm depth. Alongside, there were narrow, vertically elongated zones of low $Fe_d$ concentrations (0-20
µmol L⁻¹), which penetrated as iron-depleted “fingers” into deeper sediment (top left and bottom right of Fig. 5b).

Fig 5c shows the reconstructed 2D flux pattern in the vertical direction (z-axis) as obtained by the SGF procedure. In an idealized 1D-stratified marine sediment, dissolved iron produced by dissimilatory iron reduction typically precipitates near the sediment-water interface upon contact with O₂ and also in the deeper layers upon contact with H₂S. Accordingly, one would expect a bimodal flux distribution, with upward fluxes in the upper layer of the sediment and downward fluxes in deeper sediment horizons. Such a pattern was however not observed. Instead, a very heterogeneous “crest and trough” pattern was recorded, where localized zones with a downward flux (yellow areas) were interspersed with localized zones having an upward flux (blue areas - Fig. 5c). Only near the sediment-water interface, there was a horizontal layer between 0.2 and 0.7 cm depth, where the flux was uniformly directed upwards (a mean upward flux of 193 µmol m⁻² d⁻¹ and maximum flux of 736 µmol m⁻² d⁻¹). The fluxes in the upper 0.2 cm of the sediment are however lower (mean upward flux of 101 µmol m⁻² d⁻¹). This transition matches the oxygen penetration depth of 0.2 ± 0.02 cm (n=18; Thibault de Chanvalon et al., 2015), thus suggesting an oxidative removal of dissolved iron near the oxygen penetration depth, which hence reduces the efflux of dissolved iron across the sediment-water interface.

Figure 5d shows the reconstructed production rates, which were also characterized by a heterogeneous spatial distribution. The reconstructed rates span a broad range between consumption and production (99% percentile range = -223 to 274 µmol L⁻¹ d⁻¹). Below the top 1 cm, we observed a fragmented pattern, where localized “crests” (small production zones elongated in the vertical direction) were running in parallel to localized “troughs” (narrow consumption zones stretched in the vertical). The production crests matched the Fe₉ maxima in the concentration plots, while the production troughs corresponded to the edges of zones with high Fe₉ concentrations. Conversely, along the zones of low Fe₉ concentrations, the
position of the consumption troughs matched the Fe<sub>d</sub> minima, while the production crests corresponded to the edges of the zones with low Fe<sub>d</sub> concentrations. Between 0.2 and 0.4 cm depth, Fig. 5d indicates a horizontal layer of iron consumption (minimum -380 µmol L<sup>-1</sup> d<sup>-1</sup>; mean -30 µmol L<sup>-1</sup> d<sup>-1</sup>), which is followed by a deeper layer of production between 0.6 and 0.9 cm (maximum 586 µmol L<sup>-1</sup> d<sup>-1</sup> ; mean 59 µmol L<sup>-1</sup> d<sup>-1</sup>). This alternation matches the flux pattern discussed above.
### 3.3. Diffusive flux across the sediment-water interface

As coastal and shelf sediments are important sources of Fe<sub>d</sub> to overlying waters, it is crucial to be able to reliably estimate the diffusive flux across the sediment water interface. Typically, this is done by applying Fick's first law to a one-dimensional pore water depth profile, obtained by pore water extraction from vertical slices of sediment. The concentration gradient required in Fick's first law is then determined by linear regression on the data point immediately below the sediment water interface. Here, this pore water extraction procedure was simulated by suitably averaging the 2D concentration distribution (equivalent to 0.5 cm depth slices; Fig. 6a).

The diffusive flux of Fe<sub>d</sub> obtained by applying Fick's first law to the resulting 1D profile is 395 µmol m<sup>-2</sup> d<sup>-1</sup> (Fig. 6a). This flux value can be compared to the average of the J<sub>z</sub> values obtained by the SGF procedure at each pixel point along the sediment-water interface (Figs. 6b and 6c), which provides a flux value of 97 µmol m<sup>-2</sup> d<sup>-1</sup> (sd= 99 µmol m<sup>-2</sup> d<sup>-1</sup>, n=369). The 2D method calculates the mean flux across a straight horizontal line, and hence, it neglects any surface topography, which may impact the diffusive flux. However, as shown by Røy et al. (2005) for the diffusive flux of oxygen, this topography effect leads to relatively small deviations on the order of ~10%. Instead, the correct vertical positioning of the SWI has a more important impact. The topography of the sediment, and possibly also the sediment sampling and the handling of the 2D-DET gel, generates an uncertainty of at least 1 mm on the vertical positioning of the sediment-water interface (SWI). To assess this uncertainty, we calculated the flux at 1 mm above and below the original interface, which provided a flux of 110 µmol m<sup>-2</sup> d<sup>-1</sup> and 105 µmol m<sup>-2</sup> d<sup>-1</sup> respectively (Fig. 6b and 6c). These 2D fluxes are hence highly comparable (sd=7 µmol m<sup>-2</sup> d<sup>-1</sup>; n=3) despite the large lateral variability in flux values generated by the small scale of our sampling. Overall, the 2D flux obtained (mean of three sections 104 µmol m<sup>-2</sup> d<sup>-1</sup>) is 4 times lower than the value obtained by the corresponding 1D analysis (395 µmol m<sup>-2</sup> d<sup>-1</sup>). These differences in fluxes are likely due to difference in
resolution between the conventional 1D approach (resolution of 5 mm) and the high-resolution 2D approach implemented here (resolution of 0.2 mm). Because the oxygen penetration is shallow (~ 2 mm), intense iron oxidation and precipitation occurs within the upper millimetres, which is beyond the resolution of the conventional 1D approach. Oxidative consumption of Fe<sub>d</sub> in the 0.2-0.5cm layer (Fig 5d) likely decreases the iron release to the overlying water column. This highly localized zone of aerobic iron oxidation is not suitably accounted for by the 1D procedure, thus overestimating the iron efflux from the sediment.
### 3.4. Uncertainty assessment

The SGF procedure estimates the fluxes and rates within the pore water, and to this end, it requires a set of input parameters (salinity, temperature, porosity). These parameters are always determined with some experimental uncertainty, and hence, one could ask how sensitive the reconstructed fluxes and rates are to this experimental uncertainty. At the same time, the SGF procedure is also crucially dependent on the choice of the size of the filter window (the smoothing parameter *n*), which is in the approach here, automatically determined by fractal analysis (see methods). Accordingly, one can also ask how sensitive the reconstructed fluxes and rates are to the numerical selection of the filter window.

To analyse the sensitivity of our reconstruction approach to these input parameters (salinity, temperature, porosity, and filter window), we used the field dataset and varied a single input parameter over a suitable range (each time repeating the SGF procedure and keeping the other input parameters at their baseline value). Salinity variation did not have an important effect (results not shown). This is because the effect of salinity variations on the diffusion coefficient is one hundred times less important than temperature. Figure 7 shows the importance of the variation of surface porosity, temperature and filter window size, on a characteristic production rate (the standard deviation of the estimated production rate distribution) and on a characteristic flux (the mean estimated flux across the sediment-water interface). For each parameter, the SGF procedure is applied for a range of values and the corresponding values of the characteristic production rate and the characteristic flux are calculated.

The horizontal red lines in Fig. 7 indicate the experimental uncertainty on the values of surface porosity and temperature, as estimated for the sediment at the field site at the time of sampling (surface porosity range from 0.88 to 0.96; temperature range from 11°C to 15°C). The uncertainty on the filter window size was determined by verifying the uncertainty on the automated detection of the breakpoint in the slope of the fractal dimension *D* (see methods).
The associated uncertainty in the characteristic production rate and characteristic flux are indicated by the vertical red lines in Fig. 7, and are also summarized in Table 2. Overall, the reconstructed production rate is particularly sensitive to the filter window size (Fig. 7c) which induces an uncertainty of ± 9% while the reconstructed flux is particularly sensitive to surface porosity variations which induces an uncertainty of ± 11% (Table 2).

Figure 7c shows that an increase of the filter window size (n), decreases the characteristic production rate. This is due to a smoothing of reconstructed concentration data that averages localized production and consumption zones (Fig. 7c). In contrast, Fig. 7f shows a non-linear relation between the characteristic flux and the window size. For larger window sizes (more important smoothing), the characteristic flux first decreases and then increases again. This specific response arises because of the consumption of Fe<sub>d</sub> near the oxygen penetration depth, so the Fe<sub>d</sub> concentration quickly increases below the oxygen penetration depth (thus increasing the flux for large filter windows). The minimum flux estimated for n=24 is induced by a short plateau on Fe<sub>d</sub> at 4mm depth that generates a decrease of the slope (thus of the estimated flux) for intermediate filter windows.

Finally, the quadratic sum of each input parameter uncertainty (temperature, porosity and windows size) generates a "parameters sensitivity uncertainty" of about 11% for the reaction rate and 13% for the flux across the SWI (Table 2). The uncertainty of the whole SGF procedure is calculated by the quadratic sum of (1) the maximal error of numerical reconstruction estimated from the synthetic dataset (§3.1; Table 1, bolt values), (2) the "parameters sensitivity uncertainty" (see above; Table 2), (3) the impact of non-linear topography and (4) the error for surface positioning (§3.3). Combination by quadratic sum is summarized in the Table 3 and leads to an error of 37% on the characteristic production rate and of 20% on the characteristic flux. These values are comparable with benthic chamber measurement (e.g. Severmann et al., 2010).
# 4. Discussion

## 4.1. Impact of burrow networks on iron cycling

Figure 8 illustrates the close relation between sedimentary features and Fe<sub>d</sub> production. Macrofauna renew the water of their burrows with oxygen-rich and Fe<sub>d</sub>-poor water obtained from the overlying water column (e.g. Meysman et al., 2006). Accordingly, one expects low Fe<sub>d</sub> inside actively flushed burrows, as dissolved iron is either transported out of the burrow or aerobically oxidized. This bio-irrigation effect is clearly observable in the top 4 cm of the sediment, where the lumens of macrofauna burrows show a strong depletion in dissolved iron (from 150 µmol L<sup>-1</sup> in the surrounding area to detection limit in the burrow lumen; Fig 8b; burrows in upper left part). These burrow linings are generally also zones of dissolved iron consumption, surrounded by “linings” of iron production (Fig. 8c). Below 4 cm depth, the relation between burrow systems and the associated Fe<sub>d</sub> concentration distribution becomes more complex to interpret. For example, the right hand side of the investigated sediment domain shows a large and well-defined burrow system, of which the deeper parts are associated with elevated iron concentrations (Fig. 8b; red zone at z=5-14 cm, y=4-6.5 cm). From the concentration image alone, one would associate the deeper burrow part with iron production. Similar burrow enrichment is reported by Cesbron et al. (2014) and Zhu and Aller (2012) and is interpreted as abandoned burrows favourable for iron transport. However, when closely scrutinizing the image of the reconstructed production rate (Fig. 8c), one clearly observes a long stretch of consumption penetrating downwards, right where the burrow is, which is flanked by zones of intense production. Overall, the “vertical fingering pattern” with elongated zones of alternating production and consumption zones (Fig 8c) is striking, and corresponds to a large extent to the presence of vertically oriented burrows. Accordingly, it is clear that iron mobilization takes place in the immediate vicinity of the burrow walls, confirming that burrow construction and ventilation have a considerable impact on sedimentary iron cycling (Aller, 1982; Cesbron et al., 2014; Zhu and Aller, 2012).
To better document the iron cycling near the burrows, Figure 9 displays a detailed
zoomed-in image of three separate burrows (pink rectangles in Fig. 8a). The Fe<sub>d</sub> concentration
profiles and corresponding production rates are shown along transects across the burrow
(white line segments in Fig. 8a). Dissolved iron concentrations typically decrease towards the
burrow centre, designating iron removal inside the burrow. The production rate transects
reveal more detail and allow four separate zones to be distinguished: (1) the burrow lumen
displaying strong consumption (2) an inner burrow wall zone with moderate consumption
rate, less than inside the burrow, (3) an outer burrow wall zone showing intense production,
and (4) the far field away from the burrow that shows either low iron production or
consumption.

The removal of Fe<sub>d</sub> within the burrow lumen (the grey shaded area in Fig. 9b) is most
likely the result of burrow flushing, and the oxidation state of the iron that is exported into the
overlying water will largely depend on the oxygen level in the ventilated burrow water.
Recently, planar optode studies have demonstrated that burrow flushing can result in plumes
of anoxic water that are expelled from the sediment (Volkenborn et al., 2012). If this is the
case, ferrous iron will be transported out of the burrow lumen. However, if the burrow water
is sufficiently oxygenated, ferrous iron will be oxidized as nanoparticles in suspension
(Raiswell et al., 2008; Taillefert et al., 2000), colloïds (Homoky et al., 2011; Huerta-Diaz et
al., 2007; Wells and Goldberg, 1994) or form stable dissolved organic-Fe(III) complexes
(Gledhill and Buck, 2012; Liu and Millero, 2002). The removal of Fe<sub>d</sub> within the first 1-4 mm
of the burrow wall (zone 2) is likely the result of oxygen diffusing from the burrow into
burrow wall, and reacting with reduced iron (Aller, 1982). This process will produce the
accumulation of iron (hydr)oxides within the burrow wall, which are sometimes visible as
such orange mineral enrichments that line the burrow. In a transient situation (e.g. when the
O<sub>2</sub> level in the burrow decreases), these iron (hydr)oxides may again give rise to ferrous iron
mobilization, when used for organic matter degradation through dissimilatory iron reduction
(Lovley, 1991) or abiotic reduction with hydrogen sulphide (Canfield, 1989; Kostka and Luther, 1995; Theberge and Luther, 1997). This may be one explanation for the observed production in the subsequent outer burrow wall zone. Such reduction of freshly precipitated of Fe(III) may be particularly enhanced once the burrows are abandoned. However, alternative processes for Fe<sub>d</sub> formation could be FeS dissolution due to low pH around burrows (Hulth et al., 2002; Rickard, 2006; Zhu et al., 2006). Finally, the iron production/consumption in the far field away from the burrow will be determined by the local availability of geochemical reagents (e.g. sulphide, labile organic matter).

## 4.2. Bio-irrigational flux of dissolved iron.

In section 3.3, we estimated the diffusive flux of Fe<sub>d</sub> across the sediment-water interface ($J_{\square_{dif}} = 104 \pm 21 \mu\text{mol Fe m}^{-2} \text{ d}^{-1}$; uncertainty of 20%, see Table 3). However, in coastal sediments subject to bio-irrigation, this only forms one component of the sediment-water solute exchange. One has to additionally account for the bio-irrigational flux across the sediment-water interface $J_{\square_{irr}}$, which is generated by the flushing of ventilated burrow networks (Aller, 2001; Meile et al., 2001; Wenzhofer and Glud, 2004). To this end, we can calculate the amount of Fe<sub>d</sub> removed from each individual burrow ($F_{burrow}$), expressed as mass per unit of time (e.g. mmol Fe d<sup>-1</sup>). The bio-irrigational flux is then given by the summation over all burrows and normalized per unit of sediment surface area:

$$ \bar{J}_{irr} = (\sum F_{burrow})/(d_{sed}L_{sed}) \quad (10) $$

The sediment area normalization is done by the width $L_{sed}$ of the studied sediment domain (width of the 2D-DET gel analysed: 7.7 cm) and the thickness $d_{sed}$ of the analysed sediment domain (note that $d_{sed}$ drops out further, and that there is no need for knowing the sediment thickness sampled by 2D-DET in front of the gel).
The burrow removal rate ($F_{burrow}$) can be calculated in two separate ways. In a first procedure, one assumes that $F_{burrow}$ must be equal to the mass of Fe$_d$ that diffuses across the burrow lining into the burrow. The SGF procedure is able to provide the horizontal and vertical fluxes in each point (i.e. the two component of the flux vector – see Fig. 5c). So, in theory, $F_{burrow}$ can be calculated as the line integral of the diffusive flux ($J_{dif}$) over the burrow lining (described by the infinitesimal distance $dl$ oriented by its normal vector $n$):

$$F_{burrow} = d_{sed} \int_{burrow lining} J_{dif} \cdot n \, dl \quad (11)$$

Note that the “burrow lining” is not a closed surface area, but the line contour of the burrow in the 2D image of the sediment (see Fig 8c). Because of our planar 2D approach, the line integral is multiplied by the thickness $d_{sed}$ of the analysed sediment domain (note that $d_{sed}$ drops out of the equation when applying Eq. 11 to calculate the bio-irrigation flux).

This procedure hence requires (1) the visual identification of all burrow linings from the sediment photograph. Only burrows that showed clear iron depletion were considered as actively flushed burrows and were taken into account, (2) the superimposition of these burrow linings onto the 2D false coloured image of the reconstructed flux, and (3) the numerical calculation of the line integral. According to the lining shown in black in the Fig. 8c, this procedure generates a bio-irrigation flux of $J_{□}_{irr} = 250 \pm 50 \mu mol m^{-2} d^{-1}$. The uncertainty corresponds to the general uncertainty (20%) derived above for flux reconstruction by the SGF procedure (Table 3). However, we also found that this procedure was very sensitive to the correct delineation of the burrow linings, and the uncertainty is likely much larger. For example, small spatial offsets can occur between the sediment image and the 2D-DET gel, due to a slight modification of the sediment position during retrieval and transport, or the transformation of the gel geometry during freezing and thawing. For this reason, this first approach for calculating the bio-irrigation flux was not pursued any further.
In a second procedure, we still assume that $F_{burrow}$ represents the mass of Fe$_d$ that diffuses across the burrow lining into the burrow. By means of mass balance, and assuming steady state, the burrow removal rate must be equal to the integral of the production rate over the burrow lumen (Eq. 13), with $dA$ the infinitesimal surface of the burrow lumen in contact with the 2D-DET gel.

$$F_{burrow} = - \iint_{burrow \ lumen} R \ dA \qquad (12)$$

Again we made a distinction between actively flushed and non-flushed burrows, but this time, this was done based on the 2D production image (only burrows that showed a clear consumption of dissolved iron were considered as actives). So, based on a combination of the sediment image and consumption areas, we identified actively flushed burrows. This approach avoids any spatial offsets between the sediment image and the 2D-DET gel. Subsequently we delineated the areas of dissolved iron consumption around burrow lumens on the 2D production image. These areas were constrained by the “zero production contour”, i.e., the line where dissolved iron consumption (negative R) changed into iron production (positive R; boundaries indicated by white lines in Fig. 8c). Application of the above surface integral over these zones provided a value for the bio-irrigational flux of 745 ± 276 µmol m⁻² d⁻¹. The uncertainty corresponds to the general uncertainty (37%) derived above for production rate reconstruction by the SGF procedure (Table 3). However, this is an overestimate of the true bio-irrigational flux, as the delineated areas of dissolved iron consumption around burrow lumen likely include parts of the burrow walls. As already discussed, there is also consumption within the sediment of the inner burrow wall, which occurs through re-oxidation, and hence, this does not contribute to the bio-irrigational flux. Future studies hence should concentrate on ways to better delineate burrow linings and burrow lumens in 2D-DET images. For now, we compensated for the overestimation of the bio-irrigational flux in a simplified manner. For three burrow cross-section profiles (Fig. 9), the consumption that occurs inside the burrow lumen represents 55 ± 20% of the total consumption occurring in the
burrow surrounding. When rescaling the value obtained above (745 x 0.55), we arrive at a bio-irrigational flux of 410 ± 213 µmol m⁻²d⁻¹. The uncertainty of 52% is obtained by standard error propagation (quadratic sum) of parameters sensitivity (11%), numerical reconstruction (35%) and burrow lumen positioning (37%).

The total flux across the interface (sum of bio-irrigational flux and diffusive flux) is hence 514 ± 233 µmol m⁻² d⁻¹. Since the 2D-DET procedure provides a measurement of the average pore water concentration on the time scale of a half to a few hours, this value hence represents the average flux over the 5 hour period of measurement. The bio-irrigational flux would account for 80% of the total efflux of dissolved iron, which is in agreement with the 68-99% as estimated by diagenetic model simulation for a “standard” shelf sediment (Raiswell and Canfield, 2012). These values illustrate the importance of bio-irrigation for dissolved iron fluxes in coastal environments. These total iron fluxes are also higher than those based on benthic chamber measurements (Elrod, 2004; Severmann et al., 2010), which report total Fe₉ fluxes less than 30 µmol d⁻¹ m⁻² for environments where oxygen levels in the overlying water exceed 100 µmol L⁻¹, while higher values up to 568 µmol d⁻¹ m⁻² are reported for hypoxic environments. Although intertidal, the overlying water during high tide in the estuarine setting studied here is fully oxygenated, and hence, the total Fe₉ flux reported here is substantially larger than previously reported values for well oxygenated water columns. First, this can be related to our model approach which assumes a complete transfer to the water column of all removed iron, while in reality, some of this iron may freshly precipitate inside the burrows and stick to the burrow lining (Homoky et al., 2012; Raiswell, 2011). Secondly, an important part of dissolved iron release in benthic chamber is likely to precipitate before being measured (Pakhomova et al., 2007; Severmann et al., 2010). Thirdly, our study site is intertidal and has a high burrow density and active bio-irrigating fauna (e.g. Kristensen et al. (2014) for a low salinity environment compared to Wheatcroft (2006) for continental shelf). This likely induces higher bio-irrigational fluxes than the environments where dissolved iron
fluxes have been measured by benthic chambers. A fourth important factor governing benthic iron recycling is the relative weakness of sulphur cycling ($SO_4^{2-}/Na^+$ constant within the top 12 cm; unpublished data) compared to the coastal environments where previous flux measurements have been performed (Elrod, 2004; Severmann et al., 2010). Consequently, sulphate reduction is likely substantially lower, and hence less sulphide is produced that would react with iron and immobilize it. Important iron fluxes under oxygenated water column conditions were recently reported by Pakhomova et al. (2007) for a coastal, low salinity environment (Vistula Lagoon; Baltic Sea), reporting $Fe_d$ fluxes between 40 and 270 µmol $d^{-1} m^{-2}$. Similarly, an Fe(III) flux between 8 and 227 µmol $d^{-1} m^{-2}$ was estimated by Jones et al. (2011) for estuarine sediment at low salinity (Satilla river estuary, Georgia, USA). These studies, in combination with the results obtained here, indicate that coastal sediments can release substantial amounts of dissolved iron under an oxygenated water column.

### 4.3. An iron budget of a coastal sediment

Figure 10 presents the tentative iron budget calculated for one square meter of sediment at the Loire field site. The "apparent recycling rate" (ARR) is defined as the transfer of iron from the dissolved phase back to the solid phase of the sediment, expressed for 1 $m^2$ of mudflat. The ARR is calculated by (1) integrating the production rate for all pixels that show a negative production rate (i.e. iron consumption) over the whole sediment domain (7.7 cm x 16.2 cm), (2) dividing this quantity by the domain width ($L_{sed} = 7.7$ cm), and (3) subtracting the bio-irrigational flux (this part of the "consumption" is not truly recycled, but actually leaves the sediment). The word "apparent" is chosen to highlight that the calculation assumes (1) steady state and (2) that within a given pixel only consumption occurs (and no production). This may not always be the case, as within a given localized spot, crypting Fe cycling could occur within the pore water (i.e. production and consumption occur at the same time, so the actual consumption is larger). Moreover, it should be noted that if the steady state
assumption is not satisfied, the ARR also includes temporal variations (Eq. 1). Applied to the 2D-DET data, we arrive at an ARR of 3780 ± 1510 µmol m⁻² d⁻¹. An equivalent estimation can be done based on the production rate (i.e. only accounting for areas where production is positive) and this provides a value for the “apparent production rate” (APR) of 4240 ± 1700 µmol m⁻² d⁻¹. Note that in a steady state the APR should match the sum of the ARR (3780 µmol m⁻² d⁻¹), the diffusive flux across the sediment-water interface (104 µmol m⁻² d⁻¹) and the bio-irrigational flux (410 µmol m⁻² d⁻¹). Accordingly, the iron budget estimated here shows only a small deficit of 1.4% (54 µmol m⁻² d⁻¹), which remains within the uncertainty of the estimated fluxes (Fig. 10). These high ARR values imply that iron recycling is likely more important than previously thought in some coastal sediments. In comparison, the model used by Raiswell and Canfield (2012) for a “typical” shelf sediment calculated an ARR of 574 µmol m⁻² d⁻¹, which is almost one order of magnitude below our estimate obtained here.

As the "apparent recycling rate" corresponds to the fluxes between the interstitial water and the mineral particles, it characterizes the rejuvenation rate of iron oxide in the sediment and therefore controls the age of iron particles exported to the ocean through resuspension events. The dissolved iron inventory in the pore water of the mudflat (9680 µmol m⁻²) can be divided by the total production rate of dissolved Fe (4290 µmol m⁻² d⁻¹; Fig. 10) to arrive at the mean residence time of iron in the dissolved phase (here ~2.3 days). This residence time characterizes the rate of turnover of the dissolved iron pool. Furthermore, in a steady-state situation, the total efflux of dissolved iron across the sediment-water interface (514 µmol m⁻² d⁻¹) should be matched by an equal input of iron (hydr)oxides into the sediment (Fig. 10). Each iron atom that enters the sediment goes to a number of cycles in which it is first reductively dissolved and subsequently oxidatively precipitated again. The ratio of the apparent recycling rate over the input of iron (hydr)oxides (or equally the total Feₐ efflux) reveals that each iron atom goes through n = 7.4 redox cycles.
As already noted above, this number does not account for crypting Fe cycling within the pore water. For example, in the water column of seasonal hypoxic systems, it has been shown that suspended Fe(III) (oxy)hydroxides can be reduced by H₂S to Fe(II) at the oxic–anoxic interface, and subsequently, the Fe(II) is re-oxidized to Fe(III) by O₂ above the interface (Ma et al., 2006; Findlay et al., 2014.). These transformations form a catalytic redox cycle for iron, which fully occurs in suspension. One could imagine that a similar catalytic redox cycle could occur within the lumen of burrows or within the pore water of burrow linings. The 2D-DET probes employed here do not discriminate between dissolved Fe(II) and Fe(III), and hence, such internal cycling of Fe within the dissolved phase would remain undetected. According the calculated number of n = 7.4 redox cycles should be considered as a conservative estimate.

**5. Conclusion**

The SGF procedure was successfully applied on a two dimensional distribution of total dissolved iron concentrations in an estuarine sediment. The dissolved iron that is produced inside the sediment can follow three ways: (1) diffusion across SWI (2.3%), (2) removal by bio-irrigational flushing in macrofauna burrows (9.6%) and (3) reoxidation and precipitation in the sediment (88.1%). These results hence lead to two main conclusions. Firstly, macrofaunal bio-irrigation substantially increases the efflux, and therefore controls export of dissolved iron from sediments to the coastal ocean. Secondly, our results imply that iron has a high chance of being recycled, and will undergo, on average, to at least 7.4 cycles of dissolution/precipitation before being released in the water column. This intense recycling is likely to decrease the age of crystalline iron in the sediment and therefore controls the age of particles exported to the ocean through resuspension events. Overall, the 2D SGF procedure developed here has enabled an improved quantitative insight into the geochemical cycling of iron in dynamic benthic environments such as estuarine or coastal mudflats. A
better understanding of the spatial dynamics, local geochemical cycling and complex interactions could be potentially achieved by applying the 2D SGF procedure to other chemical compounds.

<figure><img src="image_2.png" /></figure>
Acknowledgements

The research leading to these results has received funding from the Région Pays de la Loire
(France) through the project RS2E-OSUNA and from the European Research Council under
the European Union's Seventh Framework Programme (FP/2007-2013) through ERC Grant
306933 (FJRM). The authors want to thank the associate editor and reviewers for their
valuable suggestions to improve this manuscript.
References

Aller, R.C., 2001. Transport and reactions in the bioirrigated zone. Benthic Bound. Layer 5, 269-301.

Aller, R.C., 1982. The Effects of Macrobenthos on Chemical Properties of Marine Sediment and Overlying Water, in: McCall, P.L., Tevesz, M.J.S. (Eds.), Animal-Sediment Relations, Topics in Geobiology. Springer US, pp. 53-102.

Berner, R.A., 1980. Early diagenesis : a theoretical approach. Princeton : Princeton University Press.

Bligh, M.W., Waite, T.D., 2011. Formation, reactivity, and aging of ferric oxide particles formed from Fe(II) and Fe(III) sources: Implications for iron bioavailability in the marine environment. Geochim. Cosmochim. Acta 75, 7741-7758.
doi:10.1016/j.gca.2011.10.013

Boudreau, B.P., 1997. Diagenetic Models and Their Implementation Modelling Transport and Reactions in Aquatic Sediments. Springer Berlin Heidelberg, Berlin, Heidelberg.

Boyd, P.W., Ellwood, M.J., 2010. The biogeochemical cycle of iron in the ocean. Nat. Geosci. 3, 675-682. doi:10.1038/ngeo964

Canfield, D.E., 1989. Reactive iron in marine sediments. Geochim. Cosmochim. Acta 53, 619-632.

Cesbron, F., Metzger, E., Launeau, P., Deflandre, B., Delgard, M.-L., Thibault de Chanvalon, A., Geslin, E., Anschutz, P., Jézéquel, D., 2014. Simultaneous 2D Imaging of Dissolved Iron and Reactive Phosphorus in Sediment Porewaters by Thin-Film and Hyperspectral Methods. Environ. Sci. Technol. 48, 2816-2826.
doi:10.1021/es404724r

Chever, F., Rouxel, O.J., Croot, P.L., Ponzevera, E., Wuttig, K., Auro, M., 2015. Total dissolvable and dissolved iron isotopes in the water column of the Peru upwelling regime. Geochim. Cosmochim. Acta 162, 66-82. doi:10.1016/j.gca.2015.04.031

Conway, T.M., John, S.G., 2014. Quantification of dissolved iron sources to the North Atlantic Ocean. Nature 511, 212-215.

Davison, W., Grime, G.W., Morgan, J. a. W., Clarke, K., 1991. Distribution of dissolved iron in sediment pore waters at submillimetre resolution. Nature 352, 323-325.
doi:10.1038/352323a0

Davison, W., Zhang, H., Grime, G.W., 1994. Performance Characteristics of Gel Probes Used For Measuring the Chemistry of Pore Waters. Environ. Sci. Technol. 28, 1623-1632.
doi:10.1021/es00058a015

Egger, M., Jilbert, T., Behrends, T., Rivard, C., Slomp, C.P., 2015. Vivianite is a major sink for phosphorus in methanogenic coastal surface sediments. Geochim. Cosmochim. Acta 169, 217-235. doi:10.1016/j.gca.2015.09.012

Elrod, V.A., 2004. The flux of iron from continental shelf sediments: A missing source for global budgets. Geophys. Res. Lett. 31. doi:10.1029/2004GL020216

Gledhill, M., Buck, K.N., 2012. The Organic Complexation of Iron in the Marine Environment: A Review. Front. Microbiol. 3. doi:10.3389/fmicb.2012.00069
Glud, R.N., Forster, S., Huettel, M., 1996. Influence of radial pressure gradients on solute exchange in stirred benthic chambers. Mar. Ecol. Prog. Ser. 141, 303–311.

Gorry, P.A., 1990. General least-squares smoothing and differentiation by the convolution (Savitzky-Golay) method. Anal. Chem. 62, 570–573. doi:10.1021/ac00205a007

Harper, M.P., Davison, W., Tych, W., 1997. Temporal, spatial, and resolution constraints for in situ sampling devices using diffusional equilibration: Dialysis and DET. Environ. Sci. Technol. 31, 3110–3119.

Homoky, W.B., Hembury, D.J., Hepburn, L.E., Mills, R.A., Statham, P.J., Fones, G.R., Palmer, M.R., 2011. Iron and manganese diagenesis in deep sea volcanogenic sediments and the origins of pore water colloids. Geochim. Cosmochim. Acta 75, 5032–5048. doi:10.1016/j.gca.2011.06.019

Homoky, W.B., John, S.G., Conway, T.M., Mills, R.A., 2013. Distinct iron isotopic signatures and supply from marine sediment dissolution. Nat. Commun. 4. doi:10.1038/ncomms3143

Homoky, W.B., Severmann, S., McManus, J., Berelson, W.M., Riedel, T.E., Statham, P.J., Mills, R.A., 2012. Dissolved oxygen and suspended particles regulate the benthic flux of iron from continental margins. Mar. Chem. 134–135, 59–70. doi:10.1016/j.marchem.2012.03.003

Huerta-Diaz, M.A., Rivera-Duarte, I., Sañudo-Wilhelmy, S.A., Flegal, A.R., 2007. Comparative distributions of size fractionated metals in pore waters sampled by in situ dialysis and whole-core sediment squeezing: Implications for diffusive flux calculations. Appl. Geochem. 22, 2509–2525. doi:10.1016/j.apgeochem.2007.07.001

Huettel, M., Berg, P., Kostka, J.E., 2014. Benthic Exchange and Biogeochemical Cycling in Permeable Sediments. Annu. Rev. Mar. Sci. 6, 23–51. doi:10.1146/annurev-marine-051413-012706

Hulth, S., Aller, R.C., Engstrom, P., Selander, E., 2002. A pH plate fluorosensor (optode) for early diagenetic studies of marine sediments. Limnol. Oceanogr. 47, 212–220.

Hyacinthe, C., Van Cappellen, P., 2004. An authigenic iron phosphate phase in estuarine sediments: composition, formation and chemical reactivity. Mar. Chem. 91, 227–251. doi:10.1016/j.marchem.2004.04.006

Jézéquel, D., Brayner, R., Metzger, E., Viollier, E., Prévot, F., Fiévet, F., 2007. Two-dimensional determination of dissolved iron and sulfur species in marine sediment pore-waters by thin-film based imaging. Thau lagoon (France). Estuar. Coast. Shelf Sci. 72, 420–431. doi:10.1016/j.ecss.2006.11.031

Johnson, K.S., Chavez, F.P., Friederich, G.E., 1999. Continental-shelf sediment as a primary source of iron for coastal phytoplankton. Nature 398, 697–700.

Jones, M.E., Beckler, J.S., Taillefert, M., 2011. The flux of soluble organic-iron(III) complexes from sediments represents a source of stable iron(III) to estuarine waters and to the continental shelf. Limnol. Oceanogr. 56, 1811–1823. doi:10.4319/lo.2011.56.5.1811

Kostka, J.E., Luther, G.W., 1995. Seasonal cycling of Fe in saltmarsh sediments. Biogeochemistry 29, 159–181.

Kristensen, E., 2001. Impact of polychaetes (Nereis spp. and Arenicola marina) on carbon biogeochemistry in coastal marine sedimentsPresented during the ACS Division of Geochemistry symposium ?Biogeochemical Consequences of Dynamic Interactions
Between Benthic Fauna, Microbes and Aquatic Sediments?, San Diego, April 2001.
Geochem. Trans. 2, 92. doi:10.1039/b108114d

Kristensen, E., 1989. Oxygen and carbon dioxide exchange in the polychaete Nereis virens: influence of ventilation activity and starvation. Mar. Biol. 101, 381-388.

Kristensen, E., 1981. Direct measurement of ventilation and oxygen uptake in three species of tubicolous polychaetes (Nereis spp.). J. Comp. Physiol. 145, 45-50.

Kristensen, E., Delefosse, M., Quintana, C.O., Flindt, M.R., Valdemarsen, T., 2014. Influence of benthic macrofauna community shifts on ecosystem functioning in shallow estuaries. Mar. Ecosyst. Ecol. 1, 41. doi:10.3389/fmars.2014.00041

Kuma, K., Nakabayashi, S., Suzuki, Y., Matsunaga, K., 1992. Dissolution rate and solubility of colloidal hydrous ferric oxide in seawater. Mar. Chem. 38, 133-143. doi:10.1016/0304-4203(92)90072-I

Liu, X., Millero, F.J., 2002. The solubility of iron in seawater. Mar. Chem. 77, 43-54.

Magnuson, M.L., Lytle, D.A., Frietch, C.M., Kelty, C.A., 2001. Characterization of Submicrometer Aqueous Iron(III) Colloids Formed in the Presence of Phosphate by Sedimentation Field Flow Fractionation with Multiangle Laser Light Scattering Detection. Anal. Chem. 73, 4815-4820. doi:10.1021/ac010702m

Martin, J.H., Fitzwater, S.E., 1988. Iron deficiency limits phytoplankton growth in the north-east Pacific subarctic. Nature 331, 341-343. doi:10.1038/331341a0

McManus, J., Berelson, W.M., Coale, K.H., Johnson, K.S., Kilgore, T.E., 1997. Phosphorus regeneration in continental margin sediments. Geochim. Cosmochim. Acta 61, 2891-2907. doi:10.1016/S0016-7037(97)00138-5

Meile, C., Koretsky, C.M., Cappellen, P.V., 2001. Quantifying bioirrigation in aquatic sediments: an inverse modeling approach. Limnol. Oceanogr. 46, 164-177.

Meysman, F.J., Burdorf, L., in prep. The estimation of geochemical rates and fluxes in aquatic sediments.

Meysman, F.J.R., Boudreau, B.P., Middelburg, J.J., 2005a. Modeling reactive transport in sediments subject to bioturbation and compaction. Geochim. Cosmochim. Acta 69, 3601-3617. doi:10.1016/j.gca.2005.01.004

Meysman, F.J.R., Galaktionov, O.S., Madani, S., Middelburg, J.J., 2005b. Modelling Biological Interactions in Aquatic Sediments as Coupled Reactive Transport, in: Kristensen, E., Haese, R.R., Kostka, J.E. (Eds.), Interactions Between Macro- and Microorganisms in Marine Sediments. American Geophysical Union, pp. 359-388.

Noffke, A., Hensen, C., Sommer, S., Scholz, F., Bohlen, L., Mosch, T., Graco, M., Wallmann, K., 2012. Benthic iron and phosphorus fluxes across the Peruvian oxygen minimum zone. Limnol. Oceanogr. 57, 851-867. doi:10.4319/lo.2012.57.3.0851

Pakhomova, S.V., Hall, P.O.J., Kononets, M.Y., Rozanov, A.G., Tengberg, A., Vershinin, A.V., 2007. Fluxes of iron and manganese across the sediment-water interface under various redox conditions. Mar. Chem. 107, 319-331. doi:10.1016/j.marchem.2007.06.001

Pischedda, L., Cuny, P., Esteves, J.L., Poggiale, J.-C., Gilbert, F., 2012. Spatial oxygen heterogeneity in a Hediste diversicolor irrigated burrow. Hydrobiologia 680, 109-124. doi:10.1007/s10750-011-0907-x

Press, W.H., Flannery, B.P., Teukolsky, S.A., Vetterling, W.T., 1992. Numeric recipes in C: the art of scientific computing. Camb. Univ. Press Camb.
Raiswell, R., 2011. Iron Transport from the Continents to the Open Ocean: The Aging-Rejuvenation Cycle. Elements 7, 101–106. doi:10.2113/gselements.7.2.101

Raiswell, R., 2006. Towards a global highly reactive iron cycle. J. Geochem. Explor., Extended Abstracts presented at the 7th Symp. on the Geochemistry of the Earth's Surface (GES-7) 7th Symp. on the Geochemistry of the Earth's Surface 88, 436–439. doi:10.1016/j.gexplo.2005.08.098

Raiswell, R., Benning, L.G., Tranter, M., Tulaczyk, S., 2008. Bioavailable iron in the Southern Ocean: the significance of the iceberg conveyor belt. Geochem. Trans. 9, 7. doi:10.1186/1467-4866-9-7

Raiswell, R., Canfield, D.E., 2012. The Iron Biogeochemical Cycle Past and Present. Geochem. Perspect. 1, 1–222.

Rickard, D., 2006. The solubility of FeS. Geochim. Cosmochim. Acta, A Special Issue Dedicated to Robert A. Berner 70, 5779–5789. doi:10.1016/j.gca.2006.02.029

Riisgård, H.U., Larsen, P.S., 2005. Water pumping and analysis of flow in burrowing zoobenthos: an overview. Aquat. Ecol. 39, 237–258. doi:10.1007/s10452-004-1916-x

Robertson, D., Teasdale, P.R., Welsh, D.T., 2008. A novel gel-based technique for the high resolution, two-dimensional determination of iron (II) and sulfide in sediment. Limnol Ocean. Methods 6, 502–512.

Rose, A.L., David Waite, T., 2007. Reconciling kinetic and equilibrium observations of iron(III) solubility in aqueous solutions with a polymer-based model. Geochim. Cosmochim. Acta, Physical Chemistry of Soils and Aquifers: A Special Issue in Honor of Garrison Sposito 71, 5605–5619. doi:10.1016/j.gca.2007.02.024

Røy, H., Huettel, M., Jørgensen, B.B., 2005. The influence of topography on the functional exchange surface of marine soft sediments, assessed from sediment topography measured in situ. Limnol. Oceanogr. 50, 106–112.

Savitzky, A., Golay, M.J.E., 1964. Smoothing and Differentiation of Data by Simplified Least Squares Procedures. Anal. Chem. 36, 1627–1639. doi:10.1021/ac60214a047

Senn, A.-C., Kaegi, R., Hug, S.J., Hering, J.G., Mangold, S., Voegelin, A., 2015. Composition and structure of Fe(III)-precipitates formed by Fe(II) oxidation in water at near-neutral pH: Interdependent effects of phosphate, silicate and Ca. Geochim. Cosmochim. Acta 162, 220–246. doi:10.1016/j.gca.2015.04.032

Severmann, S., McManus, J., Berelson, W.M., Hammond, D.E., 2010. The continental shelf benthic iron flux and its isotope composition. Geochim. Cosmochim. Acta 74, 3984–4004. doi:10.1016/j.gca.2010.04.022

Shuttleworth, S.M., Davison, W., Hamilton-Taylor, J., 1999. Two-Dimensional and Fine Structure in the Concentrations of Iron and Manganese in Sediment Pore-Waters. Environ. Sci. Technol. 33, 4169–4175. doi:10.1021/es990184l

Soetaert, K.E.R., Petzoldt, T., Setzer, R.W., 2010. Solving differential equations in R: package deSolve. J. Stat. Softw. 33.

Soetaert, K., Meysman, F., 2012. Reactive transport in aquatic ecosystems: Rapid model prototyping in the open source software R. Environ. Model. Softw. 32, 49–60. doi:10.1016/j.envsoft.2011.08.011

Steefel, C.I., Van Cappellen, P., 1990. A new kinetic approach to modeling water-rock interaction: The role of nucleation, precursors, and Ostwald ripening. Geochim. Cosmochim. Acta 54, 2657–2677. doi:10.1016/0016-7037(90)90003-4
Sundby, B., Anderson, L.G., Hall, P.O., Iverfeldt, V.Ake, Westerlund, S.F., 1986. The effect of oxygen on release and uptake of cobalt, manganese, iron and phosphate at the sediment-water interface. Geochim. Cosmochim. Acta 50, 1281–1288.

Taillefert, M., Bono, A.B., Luther, G.W., 2000. Reactivity of Freshly Formed Fe(III) in Synthetic Solutions and (Pore)Waters: Voltammetric Evidence of an Aging Process. Environ. Sci. Technol. 34, 2169–2177. doi:10.1021/es990120a

Theberge, S.M., Luther, G.W., 1997. Determination of the electrochemical properties of a soluble aqueous FeS species present in sulfidic solutions. Aquat. Geochem. 3, 191–211.

Thibault de Chanvalon, A., Metzger, E., Mouret, A., Cesbron, F., Knoery, J., Rozuel, E., Launeau, P., Nardelli, M.P., Jorissen, F.J., Geslin, E., 2015. Two-dimensional distribution of living benthic foraminifera in anoxic sediment layers of an estuarine mudflat (Loire estuary, France). Biogeosciences 12, 6219–6234. doi:10.5194/bg-12-6219-2015

Thompson, A., Chadwick, O.A., Rancourt, D.G., Chorover, J., 2006. Iron-oxide crystallinity increases during soil redox oscillations. Geochim. Cosmochim. Acta 70, 1710–1727. doi:10.1016/j.gca.2005.12.005

Volkenborn, N., Meile, C., Polerecky, L., Pilditch, C.A., Norkko, A., Norkko, J., Hewitt, J.E., Thrush, S.F., Wethey, D.S., Woodin, S.A., 2012. Intermittent bioirrigation and oxygen dynamics in permeable sediments: An experimental and modeling study of three tellinid bivalves. J. Mar. Res. 70, 794–823.

Volkenborn, N., Polerecky, L., Wethey, D.S., Woodin, S.A., 2010. Oscillatory porewater bioadvection in marine sediments induced by hydraulic activities of *Arenicola marina*. Limnol. Oceanogr. 55, 1231–1247. doi:10.4319/lo.2010.55.3.1231

Von Gunten, U., Schneider, W., 1991. Primary products of the oxygenation of iron (II) at an oxic–anoxic boundary: Nucleation, aggregation, and aging. J. Colloid Interface Sci. 145, 127–139.

Wells, M.L., Goldberg, E.D., 1994. The distribution of colloids in the North Atlantic and Southern Oceans. Limnol. Oceanogr. 39, 286–302.

Wenzhofer, F., Glud, R.N., 2004. Small-scale spatial and temporal variability in coastal benthic O2 dynamics: Effects of fauna activity. Limnol. Oceanogr. 49, 1471–1481.

Wheatcroft, R.A., 2006. Time-series measurements of macrobenthos abundance and sediment bioturbation intensity on a flood-dominated shelf. Prog. Oceanogr. 71, 88–122. doi:10.1016/j.pocean.2006.06.002

Yoshida, M., Kuma, K., Iwade, S., Isoda, Y., Takata, H., Yamada, M., 2006. Effect of aging time on the availability of freshly precipitated ferric hydroxide to coastal marine diatoms. Mar. Biol. 149, 379–392. doi:10.1007/s00227-005-0187-y

Zhang, H., Davison, W., 1999. Diffusional characteristics of hydrogels used in DGT and DET techniques. Anal. Chim. Acta 398, 329–340.

Zhu, Q., Aller, R.C., 2012. Two-dimensional dissolved ferrous iron distributions in marine sediments as revealed by a novel planar optical sensor. Mar. Chem. 136-137, 14–23. doi:10.1016/j.marchem.2012.04.002

Zhu, Q., Aller, R.C., Fan, Y., 2006. Two-dimensional pH distributions and dynamics in bioturbated marine sediments. Geochim. Cosmochim. Acta 70, 4933–4949. doi:10.1016/j.gca.2006.07.033
<figure><img src="image_3.png" /><figcaption>Figure 1: (a) schematic view of the "jaw device" for simultaneous sampling sediment and pore water, (b) photograph of the sediment after retrieval, (c) photograph of the DET-2D gel after colorimetric reactions: dissolved iron appears in dark pink and dissolved phosphorus in turquoise, (d) concentration of dissolved iron in false colour. The analysis of the correspondences between iron and sediment is performed in Fig. 8.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 2: Savitsky-Golay filter (SGF) procedure applied on 2D synthetic datasets (a) idealized concentration distribution $C_{x,y}^{\text{ideal}}$ (white crosses on dataset 3 indicate the position of most important production patches), (b) “estimated” concentration $C_{x,y}^{\text{estimated}}$, (c) error on the concentration estimation $\Delta C_{x,y}$ and (d) vertical profile extracted from the middle section of $C_{x,y}^{\text{ideal}}$ (red) and $C_{x,y}^{\text{estimated}}$ (black).</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 3: SGF procedure applied on 2D synthetic datasets: (a) idealized flux distribution $J_x^{\text{ideal}}$, (b) “estimated” flux $J_x^{\text{estimated}}$, (c) error on the downward flux estimation $\Delta J_{xx}$, and (d) vertical profile extracted from the middle section of $J_x^{\text{ideal}}$, (red) and $J_x^{\text{estimated}}$, (black).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 4: SGF procedure applied on 2D synthetic datasets: (a) idealized reaction rate distribution $R^{\text{ideal}}_{x,y}$, (b) “estimated” reaction rate $R^{\text{estimated}}_{x,y}$, (c) error on the reaction rate estimation $\Delta R_{x,y}$ and (d) vertical profile extracted from the middle section of $R^{\text{ideal}}_{x,y}$ (red) and $R^{\text{estimated}}_{x,y}$ (black).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 5: Two dimensional Savitzky-Golay procedure applied on experimental dataset: (a) porosity profile (triangle) and modelled exponential functions of porosity (plain line); the uncertainty on the porosity used to estimate parameter sensitivity is represented by the dotted lines, (b) experimental dataset of dissolved iron concentration in false colour, (c) estimated downward flux from SGF (downward flux in red and upward flux in blue) and (d) estimated reaction rate from SGF.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 6: Sediment-water interface estimated flux: (a) 1D iron profile averaged over 0.5 cm thick layers extracted from the 2D iron distribution, grey area represent ±1 standard deviation; the arrow indicates the chemical gradient extracted and the corresponding flux, (b) position of the sediment water interface (SWI) in red, and the uncertainty of ±1 mm on its position (in black and blue); the corresponding averaged dissolved iron upward fluxes estimated from the 2D SGF procedure are indicated in µmol m⁻² d⁻¹ and (c) lateral variability of benthic fluxes according to the SWI position.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 7: Sensitivity of the model to parameters uncertainty: (a-b-c) standard deviation of the production rate distribution, (d-e-f) Flux across water-sediment interface (black line) with uncertainty due to interface positioning (± 1 mm; grey area); red crosses indicate experimental parameters uncertainties (horizontal branch) and corresponding impact on modelling uncertainty (vertical branch).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 8: (a) Photograph of the sediment plate with a centimetre grid (white dashed lines); continuous lines within pink rectangles correspond to the 3 examples of cross-sections studied (Fig. 9), (b) dissolved iron concentration in false colour with schematic representation of burrows and (c) Fe<sub>d</sub> production rate estimated by SGF, with delimitations of active burrows' surfaces; in black, direct identification by superimposition on the sediment photograph and in white, delimitation according to the interpretation of consumption.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 9: (a) Zoom of the photograph on the three red rectangles from the Fig. 8a: in white, the selected transect across burrow with vertical white bars identifying burrow linings, (b) Fe<sub>d</sub> concentration (dotted line) and production rate (black line) along the selected transects with grey areas identifying burrow lumen; percentages of iron consumption inside the burrow bore relatively to the total consumption estimated on the transect are indicated.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 10: Small scale iron recycling in a suboxic iron-rich sediment in µmol m⁻² d⁻¹; diffusive flux (100±26) and bio-irrigational flux (410±185) across the water sediment interface are differentiated; values underlined and in italic are calculated assuming steady state.</figcaption></figure>
<table><caption>Table 1: Estimation of uncertainties from numerical reconstruction by the Mean Reconstruction Bias (MRB; Eq. 9) for the misfit between the “true” concentration, diffusive flux and production rate values and the “reconstructed” concentration, flux and rate values. Maximum values for diffusive flux and production rate (in bold) defined the uncertainty due to numerical reconstruction (Table 3).</caption><thead><tr><th></th><th>Concentration</th><th>Downward flux</th><th>Production rate</th></tr></thead><tbody><tr><td>dataset n°1</td><td>0.58</td><td>3.75</td><td>12.2</td></tr><tr><td>dataset n°2</td><td>2.17</td><td>7.2</td><td>31.9</td></tr><tr><td>dataset n°3</td><td>1.6</td><td>5.52</td><td>34.8</td></tr><tr><td>Environmental dataset</td><td>1.9</td><td>-</td><td>-</td></tr></tbody></table>
<table><caption>Table 2: Estimation of uncertainties from parameters sensitivity. Parameter uncertainties (1st column) are represented by the vertical branch of red crosses from Fig.7. The uncertainties they induce on the characteristic flux and on the characteristic production rate correspond to the vertical branch of the Fig.7 and are indicated in the 2nd and 3rd column. The overall uncertainties from parameters sensitivity are calculated by quadratic sum and are reported in the Table 3.</caption><thead><tr><th rowspan="2"></th><th rowspan="2">Range of parameter uncertainty</th><th colspan="2">Model uncertainty</th></tr><tr><th>Characteristic flux (µmol m<sup>-2</sup> d<sup>-1</sup>)</th><th>Characteristic production rate (µmol L<sup>-1</sup> d<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Temperature</td><td>12-16 °C</td><td>100,5-112,7 (± 6%)</td><td>72,7-81,5 (± 6%)</td></tr><tr><td>Salinity</td><td>2-4</td><td>-</td><td>-</td></tr><tr><td>Surface porosity</td><td>0,88-0,96</td><td>92,5-116,1 (± 11%)</td><td>73,8-79,3 (± 3%)</td></tr><tr><td>Window size</td><td>10-14 n</td><td>101,5-106,6 (± 2%)</td><td>84,3-70,1 (± 19%)</td></tr><tr><td>Parameters sensitivity</td><td></td><td>± 13%</td><td>± 11%</td></tr></tbody></table>
Table 3: Overall estimation of the SGF procedure uncertainties. Uncertainty from numerical reconstruction and parameter sensitivity are from Table 1 and Table 2 respectively. Uncertainty for surface topography is from Røy et al. (2005) and uncertainty for surface positioning corresponds to the relative standard deviation of the 3 transects used to calculate the diffusive flux across the SWI. The overall uncertainties are calculated by quadratic sum.

<table><thead><tr><th></th><th>Diffusive flux (µmol m<sup>-2</sup>d<sup>-1</sup>)</th><th>Production rate (µmol L<sup>-1</sup>d<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Parameters sensitivity</td><td>± 13%</td><td>± 11%</td></tr><tr><td>Numerical reconstruction</td><td>± 7%</td><td>± 35%</td></tr><tr><td>Surface topography</td><td>± 10%</td><td>-</td></tr><tr><td>Surface positioning</td><td>± 7%</td><td>-</td></tr><tr><td>(consumption inside the burrow lumen)</td><td>-</td><td>(± 37%)</td></tr><tr><td>SGF uncertainty</td><td>± 19%</td><td>± 37% (± 52%)</td></tr></tbody></table>