FESDIA (v1.0): exploring temporal variations of sediment biogeochemistry under the influence of flood events using numerical modelling

Episodic events of flood deposit in coastal environments are characterized by deposition of large quantities of sediment containing reactive organic matter within short periods of time. While steady-state modelling is common in sediment biogeochemical modelling, the inclusion of these events in current early diagenesis models has yet to be demonstrated. We adapted an existing model of early diagenetic processes to include the ability to mimic an immediate organic carbon deposition. The new model version (FESDIA) written in Fortran and R programming language was able to reproduce the basic trends from field sediment porewater data affected by the November 2008 flood event in the Rhône River prodelta. Simulation experiments on two end-member scenarios of sediment characteristics dictated by field observation (1–high thickness deposit, with low TOC (total organic carbon) and 2–low thickness, with high TOC), reveal contrasting evolutions of post-depositional profiles. A first-order approximation of the differences between subsequent profiles was used to characterize the timing of recovery (i.e. relaxation time) from this alteration. Our results indicate a longer relaxation time of approximately 4 months for ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ and 5 months for DIC (dissolved inorganic carbon) in the first scenario, and less than 3 months for the second scenario which agreed with timescale observed in the field. A sensitivity analysis across a spectrum of these end-member cases for the organic carbon content (described as the enrichment factor α) and for sediment thickness indicates that the relaxation time for oxygen, sulfate, and DIC decreases with increasing organic enrichment for a sediment deposition that is less than 5 cm. However, for larger deposits (>14 cm), the relaxation time for oxygen, sulfate, and DIC increases with α. This can be related to the depth-dependent availability of oxidant and the diffusion of species. This study emphasizes the significance of these sediment characteristics in determining the sediment's short-term response in the presence of an episodic event. Furthermore, the model described here provides a useful tool to better understand the magnitude and dynamics of flooding event on biogeochemical reactions on the seafloor.

1 Introduction

Coastal margins play a crucial role in the global marine systems in terms of carbon and nutrient cycling (Wollast, 1993; Rabouille et al., 2001b; Cai, 2011; Regnier et al., 2013; Bauer et al., 2013; Gruber, 2015). Due to their relatively shallow depth, sedimentary early diagenetic processes are critical for the recycling of a variety of biogeochemical elements, which are influenced by organic matter (OM) inputs, particularly carbon (Middelburg et al., 1993; Arndt et al., 2013). Furthermore, these processes have the potential to contribute to the nutrient source that fuels the primary productivity of the marine system. In river-dominated ocean margins (RiOmar; McKee et al., 2004), organic matter input can also be enhanced by flood events which provide a significant fraction of the particulate carbon (POC) delivered to depocenters (Antonelli et al., 2008). Organic matter derived from riverine input to sediment has biogeochemical significance in coastal marine systems (Cai, 2011). As a result, the coastal environment serves as both a sink for particulate organic carbon and nutrients, and an active site of carbon and nutrient remineralization (Burdige, 2005; McKee et al., 2004; Sundby, 2006).

In the context of early diagenetic modelling, numerical models with time-dependent capability are well-established (Lasaga and Holland, 1976; Rabouille and Gaillard, 1991; Boudreau, 1996; Soetaert et al., 1996a, b; Rabouille et al., 2001a; Archer et al., 2002; Couture et al., 2010; Yakushev et al., 2017; Munhoven et al., 2021; Sulpis et al., 2022), and they are used in many coastal and deep-sea studies. However, because of the scarcity of observations and their unpredictability, the role of massive deposit of sediment in these early diagenesis models has frequently been overlooked (Tesi et al., 2012). As these rare extreme events are being currently documented in various locations, there is a growing appreciation for their impact on the coastal margin (Deflandre et al., 2002; Cathalot et al., 2010; Tesi et al., 2012).

Attempts to use mathematical models to understand perturbation-induced events such as sudden erosion/resuspension event, bottom trawling, and turbidity-driven sediment deposition on early diagenetic processes have resulted in a variety of approaches that incorporate this type of phenomenon. As an example, Katsev et al. (2006) demonstrated that the position of the redox boundary (depth zone beneath the sediment–water interface that separates the stability fields of the oxidized and reduced species of a given redox couple) in organic-poor marine sediment can undergo massive shifts due to the flux of new organic matter on a seasonal basis, whereas on a longer timescale (e.g. decadal), redox fluctuation linked to organic matter deposition can induce the redistribution of solid-phase manganese with multiple peaks (due to depth-wise oxidation reduction of Mn). Another study in a coastal system revealed that coastal sediments change as a result of an anthropogenic perturbation in the context of bottom dredging and trawling (van de Velde et al., 2018). More recently, using a similar model, De Borger et al. (2021) highlighted that perturbation events such as trawling can possibly decrease total OM mineralization.

In river-dominated ocean margins, episodic flood events can deliver sediment with varying characteristics depending on its source origin, frequency, and intensity (Cathalot et al., 2013). Therefore, the flood characteristics have direct impact on the deposited sediment's characteristics such as scale/thickness of the deposited layer, composition (mineralogy and grain size), OM content, and so on. For example, In the Rhône prodelta, a single flash flood can deliver up to 30 cm of new sediment material in a matter of days (Cathalot et al., 2010; Pastor et al., 2018). Despite the large amount of sediment introduced by this episodic loading, vertical distribution of porewater species like oxygen (O₂) can be restored after a few days (Cathalot et al., 2010). It has also been noticed (Rassmann et al., 2020) that spring and summer porewater compositions measured for several years following fall and winter floods show quasi-steady-state profiles for sulfate and DIC (dissolved inorganic carbon). Similar massive deposition was also reported in the Saguenay fjord (Quebec, Canada) (Deflandre et al., 2002; Mucci and Edenborn, 1992). The recovery timescale from this perturbation has only been roughly estimated for species with short relaxation time such as oxygen, but this is not always the case for sulfate (${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$), dissolved inorganic carbon (DIC), or other redox species. Furthermore, due to the limitation in temporal resolution of the observations, the short-term post-depositional dynamics in the aftermath of this flood deposition event are scarcely described, making it difficult to discern how the system responds after the event. While experimental approaches (Chaillou et al., 2007) can provide useful insight into how they work, they lack the ability to provide continuous system dynamics and are often difficult to set up. A modelling approach can assist in addressing these issues, providing useful feedback in terms of the scale and response of the sediment to this type of event.

The goal of this study is to better understand episodic events in the context of flood-driven sediment deposition and their impact on benthic biogeochemistry, post-flood evolution dynamics, and relaxation timescale. As the relaxation dynamics represent a gap in our understanding of how coastal systems respond to external drivers, we characterize the timescale of the recovery of sediment porewater profiles using a first-order approximation. To accomplish this, we developed an early diagenetic model called FESDIA. The ability to explicitly simulate non-steady early diagenesis processes in systems subject to perturbation events such as massive flood or storm deposition is a novel contribution of FESDIA to early diagenetic models. In the following ways, FESDIA differs therefore from the OMEXDIA model (Soetaert et al., 1996a) by implementing:

• an explicit description of the anoxic diagenesis including (i) iron and sulfur dynamics, and (ii) methane production and consumption. In comparison, OMEXDIA has a single state variable (ODU: oxygen demand unit) to describe reduced species;

• the possibility to include sediment perturbation events such abrupt deposition of sediment.

In this paper, we only discuss part of the FESDIA model concerned with implementation of a perturbation event as it relates to some biogeochemical indicators. The model is implemented in Fortran (for speed) and linked to R (for flexibility). We demonstrate the model's utility in describing data collected from a flood event in November/December 2008 (Pastor et al., 2018) and numerically investigate the impact of varying degrees of flood type characteristics on the system's relaxation dynamics. This work is a foundation for a more in-depth investigation of the model–data biogeochemistry of the porewater and solid phase components of core samples from Pastor et al. (2018), and it provides a useful baseline for understanding the spatiotemporal dynamics of coastal marine systems subject to event-driven organic matter pulses.

2 Materials and methods

2.1 Site and events description

The Rhône prodelta serves as a case study for the development of the model used to evaluate sediment perturbation dynamics. This particular coastal area acts as the transitory zone between the inland river channel and the continental shelf (Gulf of Lion) of the Mediterranean Sea. The Rhône River with a drainage basin of 97 800 km² and mean water discharge of 1700 m³ s⁻¹ delivers up to 1.6×10¹⁰ moles C of particulate carbon (POC) annually (Sempéré et al., 2000) to the pro-deltaic part (i.e. where the river meets the sea). The Rhône prodelta covers an area of approximately 65 km² with depth ranging from 2 to 60 m (Lansard et al., 2009) and is characterized by high sedimentation rates reaching up to 41 cm yr⁻¹ in the proximal zone (Rassmann et al., 2016; lat: 43^∘18.680^′ N, long: 4^∘51.038^′ E and average depth of 21 m) (Radakovitch et al., 1999; Miralles et al., 2005). The organic matter delivered to the depocenter typically reflects the different compositional materials derived from the terrestrial domain (Pastor et al., 2018), whereas the magnitude of material transported and the quantity of organic carbon transferred laterally varies according to seasons and the period of massive instantaneous deposition (Lansard et al., 2008; Cathalot et al., 2013).

Figure 1 Map showing the locations of sampling sites off the Rhône River mouth.

Relating to the episodic pulse of organic matter, numerous studies have documented instances of flood-driven deposition from the Rhône River from a hydrographic perspective (Boudet et al., 2017; Hensel et al., 1998; Pont et al., 2017). Pastor et al. (2018) go beyond sedimentology and hydrographic characteristics to provide a concise description of the various flood types, their diagenetic signatures, and biogeochemical implications. Furthermore, published porewater chemistry and solid-phase data have highlighted sediment characteristics following such an event (Cathalot et al., 2010; Toussaint et al., 2013; Cathalot et al., 2013; Pastor et al., 2018).

2.2 Model development and implementation

Following the description of the Rhône River flood types and the composition of the flood deposit (mainly in terms of organic carbon) at the proximal station A (Pastor et al., 2018), we proceed to describe the model developed to explore the observed data and their diagenetic implications in terms of relaxation times and their evolution following this transient perturbation.

Table 1 Description of notations, phrases, acronyms, and abbreviations, as used in this paper.

Our model combines the development in the OMEXDIA model (Soetaert et al., 1996a), applied in the Rhône prodelta area (Ait Ballagh et al., 2021; Pastor et al., 2011) and which has recently been equipped with event-driven processes (De Borger et al., 2021). In De Borger et al. (2021), the authors specifically addressed the issue of bottom trawling as a mixing and an erosional process that removes an upper layer of sediment and mixes a certain layer below. In addition, the model considers a bulk categorization of reduced substance in a single state variable, ODU (oxidative oxygen unit). For our approach, the event is defined by an addition of a new layer on top of the former sediment–water interface (Table 1). Furthermore, we explicitly modelled pathways involving sulfur and iron. Following this preamble, the following sections go over aspects of the model description and parameterization. Table 1 provide some key glossary of mathematical notations used in the model.

2.2.1 Model state variables

The complete model describes the concentration of labile ($C_{\mathrm{org}}^{\mathrm{fast}}$) and semi-labile ($C_{\mathrm{org}}^{\mathrm{slow}}$) decaying organic matter, oxygen (O₂), nitrate (${\mathrm{NO}}_{\mathrm{3}}^{-}$), ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$), and dissolved inorganic carbon (DIC), following the classic early diagenetic equation of Berner (1980) and Boudreau (1997). In addition to the model from De Borger et al. (2021), our model includes sulfate (${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$), hydrogen sulfide (H₂S), methane (CH₄), and iron species (Fe²⁺ and Fe(OH)₃) (Table 2).

Table 2 State variables described in the model.

In some coastal settings, oxidation via sulfate reduction has been highlighted as the primary pathway for organic carbon (OC) mineralization, with minor contributions from manganese and iron oxidation (Burdige and Komada, 2011). In addition, the flux of integrated remineralization products such as DIC has previously been estimated to contribute up to 8 times that of diffusive oxygen uptake (Rassmann et al., 2020) – thus highlighting its importance in describing the amplitude of benthic recycling in coastal water. As such in this paper, we focus our analysis on these proxy variables (O₂, ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, DIC) because they serve as indicators of the integrated effect of the main diagenetic processes.

2.2.2 Biogeochemical reaction

Early diagenesis processes on the seafloor are driven by organic matter deposition. For areas such as the Rhône prodelta, continental organic carbon input is dominant, and it is difficult to identify the fraction of labile fraction responsible for fast OM pool consumption (Pastor et al., 2011). Moreover, observations show that some organic compounds are preferentially degraded and become selectively oxidized (Middelburg et al., 1997; Pozzato et al., 2018). As a result, the model assumed solid phase organic carbon with two reactive modelled fractions with different reactivities and $\mathrm{C}/\mathrm{N}$ ratios (Westrich and Berner, 1984; Soetaert et al., 1996a). The mineralization of OM occurs sequentially, with the labile fraction mineralizing faster than the slow decaying carbon. During the timescales considered here, the refractory organic matter class is not reactive. To compare with the observation, we consider an asymptotic OC constant (TOC_ref) for the inert fraction that scales the model-calculated TOC output to the observation (Pastor et al., 2011) (see Sect. 2.2.8). This organic carbon degradation requires oxidants, and the depth dependency in sequential utilization of terminal electron acceptors assumption first proposed by Froelich (1979) is used here. Oxygen is consumed first, followed by nitrate, iron oxides, sulfate, and finally methanogenesis occurs (Eq. 3). Because the quantity of organic matter and the relative proportions of fast and slow degrading materials decrease with depth, the overall organic matter degradation rate decreases accordingly. In the formulation of the individual biogeochemical processes, we use a similar paradigm as Soetaert et al. (1996a) (Eq. 2).

This rate of carbon mineralization of organic matter (mmol m⁻³ d⁻¹) can be expressed as

$$\begin{array}{}\text{(1)}& \text{Cprod}=\left(\text{rFast}\times C_{\mathrm{org}}^{\mathrm{fast}}+\text{rSlow}\times C_{\mathrm{org}}^{\mathrm{slow}}\right)\times {\displaystyle \frac{\left(\mathrm{1}-\mathit{\varphi }\right)}{\mathit{\varphi }}},\end{array}$$

where the rFast and rSlow are the decay rate constant (d⁻¹) for fast and slow detritus component. ϕ and (1−ϕ) are the volume fraction for both solutes and solid respectively. This process is mediated by microorganisms and oxidant availability. The primary redox reaction includes (1) oxic respiration, (2) denitrification, (3) Fe (III) reduction, (4) sulfate reduction, and (5) methane production:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}\text{OM}& +{\mathrm{O}}_{\mathrm{2}}\to {\mathrm{CO}}_{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\left(\mathrm{C}:\mathrm{N}\right)}}{\mathrm{NH}}_{\mathrm{3}}+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\\ \text{OM}& +\mathrm{0.8}{\mathrm{NO}}_{\mathrm{3}}^{-}+\mathrm{0.8}{\mathrm{H}}^{+}\to {\mathrm{CO}}_{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\left(\mathrm{C}:\mathrm{N}\right)}}{\mathrm{NH}}_{\mathrm{3}}\\ & +\mathrm{0.4}{\mathrm{N}}_{\mathrm{2}}+\mathrm{1.4}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\\ \text{OM}& +\mathrm{4}\mathrm{FeOOH}+\mathrm{8}{\mathrm{H}}^{+}\to {\mathrm{CO}}_{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\left(\mathrm{C}:\mathrm{N}\right)}}{\mathrm{NH}}_{\mathrm{3}}\\ & +\mathrm{4}{\mathrm{Fe}}^{\mathrm{2}+}+\mathrm{7}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\\ \text{OM}& +\mathrm{0.5}{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}+{\mathrm{H}}^{+}\to {\mathrm{CO}}_{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\left(\mathrm{C}:\mathrm{N}\right)}}{\mathrm{NH}}_{\mathrm{3}}\\ & +\mathrm{0.5}{\mathrm{H}}_{\mathrm{2}}\mathrm{S}+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\\ \text{OM}& \to \mathrm{0.5}{\mathrm{CO}}_{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\left(\mathrm{C}:\mathrm{N}\right)}}{\mathrm{NH}}_{\mathrm{3}}+\mathrm{0.5}{\mathrm{CH}}_{\mathrm{4}}.\end{array}\end{array}$$

This reaction can be modelled using a Monod type relationship with each oxidant having a half-saturation constant (k_s[C]) represented as ks* in the model code. The inhibition of mineralization by the presence of other oxidants is also modelled with a hyperbolic term (subtracted from 1), where k_in[C] is concentration at which the rate drops to half of its maximal value. Using these limitation and inhibition functions, a single equation for each component across the model–depth domain can be realized (Rabouille and Gaillard, 1991; Soetaert et al., 1996a; Wang and Van Cappellen, 1996), together with some possible overlap (Froelich et al., 1979; Soetaert et al., 1996a). For a generic species, this can be described mathematically as

$$\begin{array}{}\text{(3)}& \text{lim}={\displaystyle \frac{\left[C\right]}{\left(k_{\mathrm{s}\left[C\right]}+\left[C\right]\right)}}\prod \left(\mathrm{1}-{\displaystyle \frac{\left[C\right]}{\left(k_{\mathrm{in}\left[C\right]}+\left[C\right]\right)}}\right),\end{array}$$

where C is one oxidant. Formulation for individual pathways as well as values of half-saturation and inhibition constants for each oxidant can be found in Appendix A1. With this limitation term, mineralization rate per solute can be estimated using potential carbon produced via OM degradation in (Eq. 1),

$$\begin{array}{}\text{(4)}& {\mathrm{rate}}_{\min }=\text{Cprod}\times \text{lim}\times {\displaystyle \frac{\mathrm{1}}{\sum \text{lim}}},\end{array}$$

with the ∑lim the sum of all limitation terms which normalizes the term in order to always achieve the maximum degradation rate. See Soetaert et al. (1996a) for more details on the derivative of this equation.

Secondary redox reaction includes reoxidation of reduced substances (nitrification, Fe oxidation, H₂S oxidation, methane oxidation) (Eq. 5) and the precipitation of FeS. Anaerobic oxidation of methane occurs in the absence of O₂ following upward diffusion of methane to the sulfate–methane transition zone (SMTZ) (Jørgensen et al., 2019):

$$\begin{array}{}\text{(5)}& \begin{array}{l}{\mathrm{NH}}_{\mathrm{4}}^{+}+\mathrm{2}{\mathrm{O}}_{\mathrm{2}}\to {\mathrm{NO}}_{\mathrm{3}}^{-}+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}+\mathrm{2}{\mathrm{H}}^{+}\\ \mathrm{4}{\mathrm{Fe}}^{\mathrm{2}+}+{\mathrm{O}}_{\mathrm{2}}+\mathrm{6}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\to \mathrm{4}\mathrm{FeOOH}+\mathrm{8}{\mathrm{H}}^{+}\\ {\mathrm{CH}}_{\mathrm{4}}+{\mathrm{O}}_{\mathrm{2}}\to {\mathrm{CO}}_{\mathrm{2}}+\mathrm{2}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\\ {\mathrm{H}}_{\mathrm{2}}\mathrm{S}+{\mathrm{O}}_{\mathrm{2}}\to {\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}+\mathrm{2}{\mathrm{H}}^{+}\\ {\mathrm{CH}}_{\mathrm{4}}+{\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}\to \mathrm{HCO}{\mathrm{3}}^{-}+{\mathrm{HS}}^{-}+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\\ {\mathrm{Fe}}^{\mathrm{2}+}+{\mathrm{H}}_{\mathrm{2}}\mathrm{S}\to \mathrm{FeS}+\mathrm{2}{\mathrm{H}}^{+}.\end{array}\end{array}$$

These reactions are mathematically described using a coupled reaction formulation. Nitrification is limited by the availability of oxygen and the other reactions are described with a first-order term.

$$\begin{array}{}\text{(6)}& \begin{array}{rl}\text{Nitri}& =R_{\mathrm{nit}}\times {\mathrm{NH}}_{\mathrm{4}}\times {\displaystyle \frac{{\mathrm{O}}_{\mathrm{2}}}{\left({\mathrm{O}}_{\mathrm{2}}+k{s}_{\mathrm{nitri}}\right)}}\\ & \left(\mathrm{Nitrification}\right)\\ \text{Feoxid}& =R_{{\mathrm{FeOH}}_{\mathrm{3}}}\times \mathrm{Fe}\times {\mathrm{O}}_{\mathrm{2}}\\ & \text{(Iron oxidation)}\\ \text{H2Soxid}& =R_{{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}\times {\mathrm{H}}_{\mathrm{2}}\mathrm{S}\times {\mathrm{O}}_{\mathrm{2}}\\ & \text{(sulfide oxidation)}\\ \text{CH4oxid}& =R_{{\mathrm{CH}}_{\mathrm{4}}}\times {\mathrm{CH}}_{\mathrm{4}}\times {\mathrm{O}}_{\mathrm{2}}\\ & \text{(Methane oxidation)}\\ \text{AOM}& =R_{{\mathrm{CH}}_{\mathrm{4}}}\times {\mathrm{CH}}_{\mathrm{4}}\times {\mathrm{SO}}_{\mathrm{4}}\\ & \text{(Anaerobic oxidation of methane)},\end{array}\end{array}$$

where R_nit is the rate of nitrification (d⁻¹), $R_{{\mathrm{FeOH}}_{\mathrm{3}}}$, $R_{{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}$, and $R_{{\mathrm{CH}}_{\mathrm{4}}}$ are the maximum rate of oxidation of iron, sulfide, and methane via oxygen, respectively (mmol⁻¹ m³ d⁻¹). Because sulfide precipitation can occur in some coastal sediments, we accounted for this sink process by removing produced sulfide from sulfate reduction as a first-order FeS formation.

$$\begin{array}{}\text{(7)}& \text{FeSprod}=R_{\mathrm{FeSprod}}\times \mathrm{Fe}\times {\mathrm{H}}_{\mathrm{2}}\mathrm{S}\text{(FeS production)},\end{array}$$

with R_FeSprod the rate of production of FeS (mmol⁻¹ m³ d⁻¹).

2.2.3 Transport processes

Transport processes in the model are described by molecular diffusion and bio-irrigation for dissolved species whereas bioturbation is the main process for mixing the solid phase. In addition, advection occurs in both the solid and dissolved species. The model dynamics described as a partial differential equation (PDE) is the general reaction–transport equation (Berner, 1980). We use a similar paradigm and formulations to that of Soetaert et al. (1996a). For substances that are dissolved:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}{\displaystyle \frac{\partial \mathit{\varphi }C}{\partial t}}& =-{\displaystyle \frac{\partial }{\partial z}}\left[-\mathit{\varphi }\times D_{\mathrm{sed}}\times {\displaystyle \frac{\partial C}{\partial z}}+w_{\mathrm{\infty }}\times {\mathit{\varphi }}_{\mathrm{\infty }}\times C\right]\\ & +\sum \mathit{\varphi }\times \text{REAC}.\end{array}\end{array}$$

With special consideration of ammonium adsorption to sediment particles, the governing equation is given by

$$\begin{array}{}\text{(9)}& \begin{array}{rl}{\displaystyle \frac{\partial \mathit{\varphi }C}{\partial t}}& =-{\displaystyle \frac{\partial }{\partial z}}\left[-{\displaystyle \frac{\mathit{\varphi }\times D_{\mathrm{sed}}}{\left(\mathrm{1}+k_{\mathrm{ads}}\right)}}\times {\displaystyle \frac{\partial C}{\partial z}}+w_{\mathrm{\infty }}\times {\mathit{\varphi }}_{\mathrm{\infty }}\times C\right]\\ & +\sum {\displaystyle \frac{\mathit{\varphi }\times \text{REAC}}{\left(\mathrm{1}+k_{\mathrm{ads}}\right)}},\end{array}\end{array}$$

where we assumed that the immobilization of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ is in instantaneous, local equilibrium (i.e. any changes caused by the slow ${\mathrm{NH}}_{\mathrm{4}}^{+}$ removal process results in an immediate adjustment of the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ equilibrium; so, can be modelled with a simple chemical species) and k_ads is the adsorption coefficient. The inclusion of this formulation for the diffusion and reaction term has the effect of slowing down ammonium migration in sediment. Derivation of this formulation is given in Berner (1980) and Soetaert and Herman (2009).

For the solid phase:

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{\displaystyle \frac{\partial \left(\mathrm{1}-\mathit{\varphi }\right)S}{\partial t}}& =-{\displaystyle \frac{\partial }{\partial z}}\left[-\left(\mathrm{1}-\mathit{\varphi }\right)\times D_{\mathrm{b}}\times {\displaystyle \frac{\partial S}{\partial z}}\right.\\ & \left.+w_{\mathrm{\infty }}\times {\left(\mathrm{1}-\mathit{\varphi }\right)}_{\mathrm{\infty }}\times S\right]\\ & +\sum \left(\mathrm{1}-\mathit{\varphi }\right)\times \text{REAC},\end{array}\end{array}$$

where C is the concentration of porewater (unit of mmol m⁻³ liquid) for Eq. (8) and S for solid (unit of mmol m⁻³ solid) Eq. (10). w (cm d⁻¹) and D_sed (cm² d⁻¹) represent the burial/advection and molecular diffusion coefficient in the sediment respectively, and REAC is the source/sink processes linked to biogeochemical reactions in the sediment. This term includes both biological and chemical reaction within the sediment column and non-local bio-irrigation transport term (see next section). D_b is the bioturbation term for solid driven by the activities of benthic organisms. For dynamic simulation, w can change as a function of time but in most cases we assumed a constant value.

Diffusive fluxes of solutes across the sediment–water interface are driven by the concentration gradients between the overlying seawater and the sediment column. Fick's first law is used to describe the solute flux due to molecular diffusion,

$$\begin{array}{}\text{(11)}& J_{\mathrm{d}}=-\mathit{\varphi }{D}_{\mathrm{sed}}{\displaystyle \frac{\partial C}{\partial z}},\end{array}$$

where the D_sed (cm d⁻¹) is the effective diffusion coefficient corrected for tortuosity and given as $D_{\mathrm{sed}}=\frac{D^{\mathrm{sw}}}{{\mathit{\theta }}^{\mathrm{2}}}$, with D^sw the molecular diffusion coefficient of the solute in free solution of sea water, and θ is the tortuosity derived from the formation factor (F) and porosity (ϕ) of a sediment matrix (Berner, 1980; Boudreau, 1997). This molecular diffusion coefficient is calculated as function of temperature and salinity using compiled relation of Boudreau (1997), implemented in the R package Marelac (Soetaert and Petzoldt, 2020).

As a simplifying assumption, material accumulation has no effect on porosity. We further assumed the porosity profile decreased with depth but invariant with time. Although, this assumption is a restrictive as the site of flood deposition can undergo variation in grain size which might affect their porosity (Cathalot et al., 2010), we proceed noting that the fixed parameters which define the porosity curve can be changed when necessary. Thus, using optimized parameters fitted with data in the proximal sites of Rhône prodelta (Ait Ballagh et al., 2021), porosity (ϕ(z)) in Eqs. (10)–(8) is prescribed as an exponential decay,

$$\begin{array}{}\text{(12)}& \mathit{\varphi }\left(z\right)={\mathit{\varphi }}_{\mathrm{\infty }}+\left({\mathit{\varphi }}_{\mathrm{0}}-{\mathit{\varphi }}_{\mathrm{\infty }}\right)e^{\frac{-\left(z-z_{\mathrm{swi}}\right)}{\mathit{\delta }}},\end{array}$$

where ϕ₀ and ϕ_∞ is the porosity surface and at deeper layer, respectively, while Z_swi is the depth of the SWI (sediment–water interface), and δ (cm) is the porosity exponential decay coefficient with depth.

2.2.4 Bioturbation and bio-irrigation

Bioturbation in the model is characterized by the movement and mixing of particles by benthic organisms. This is parameterized as a diffusivity function in space (D(z)) and acts on the concentrations of the different solid species in the sediment. In our model, this bioturbation flux is assumed to be interphase, with porosity ϕ(z) remaining constant over time. Thus, this process is prescribed as

$$\begin{array}{}\text{(13)}& D_{\mathrm{b}}\left(z\right)=\left\{\begin{array}{ll}{D}_{\mathrm{b}}^{\mathrm{0}}& \text{if }Z\le Z_{\mathrm{L}}\\ D_{\mathrm{\infty }}+\left(D_{\mathrm{b}}^{\mathrm{0}}-D_{\mathrm{\infty }}\right)e^{-\frac{\left(Z-Z_{\mathrm{L}}\right)}{{\mathrm{biot}}_{\mathrm{att}}}}& \text{if }Z>Z_{\mathrm{L}}\end{array}\right.,\end{array}$$

where $D_{\mathrm{b}}^{\mathrm{0}}$ is the bio-diffusivity coefficient (cm² d⁻¹) at the SWI and in the mixed layer, Z_L is the depth of the mixed layer (cm), and biot_att is the attenuation coefficient (cm⁻¹) of bioturbation below the mixed layer. D_∞ is the diffusivity at the deeper layer as usually specified as zero. In the model, we did not account for mortality of benthic fauna following the deposition as in De Borger et al. (2021) where they focus on habitat recolonization after trawling.

Bio-irrigation is modelled in an identical manner to that of biodiffusion, and acts as a non-local exchange process between the porewater parcels and the overlying bottom water.

$$\begin{array}{}\text{(14)}& \mathrm{\text{Irr}}\left(z\right)=\left\{\begin{array}{ll}{\text{Irr}}_{\mathrm{0}}& \text{if }Z\le Z_{\mathrm{L}}\\ {\text{Irr}}_{\mathrm{\infty }}+\left({\text{Irr}}_{\mathrm{0}}+{\text{Irr}}_{\mathrm{\infty }}\right)e^{-\frac{\left(Z-Z_{\mathrm{L}}\right)}{{\mathrm{Irr}}_{\mathrm{att}}}}& \text{if }Z>Z_{\mathrm{L}}\end{array}\right.,\end{array}$$

for which Irr₀ is the bio-irrigation rate (d⁻¹) and Irr_att is the attenuation of irrigation (cm) below the depth of the irrigated layer Z_irr (cm). At depth, the bio-irrigation rate (Irr_∞) is generally set to zero.

2.2.5 Model vertical grid

The model is vertically resolved with grid divided into 100 layers (N), of thickness (Δz) increasing geometrically from 0.01 cm at the sediment–water interface to 6 cm at the lower boundary. The result is a 100 cm model domain comprising of a full grid with non-uniform spacing and maximum resolution near the SWI. Depth units are in centimetres. This choice of modelled depth allows for complete carbon degradation. This modelled grid is generated by the grid generation routine of the ReacTran R package (Soetaert and Meysman, 2012), which implements many grid types used in early diagenesis modelling. During the time instance of the event specification, the added grid of new layers (N_pert≈Z_pert) and the current grid (N_cur≈N) is rescaled to the model's common grid of N layer by linear interpolation (see Sect. 2.2.6 and Fig. S1). The concentration of state variables is defined at the layer midpoints, whereas diffusivities, advection (sinking/burial velocities), and resulting transport fluxes are defined at the layer interfaces.

Figure 2 Schematic of model implementation for the deposition event scenario. Profile from previous time step (left) and after addition of new layer over a predefined depth layer (right). For the solid (a), the new layer can be enriched (blue) or depleted (red) relative to the old (average) (black). The dissolved substances (b) are set equal to the bottom water concentration during the deposition. Thereafter, the profile is integrated forward with time. The whole sequence of step occurs dynamically with time capitalizing on the integrator ability to simulate dynamic event process. α is the carbon enrichment factor applied over depth Z_pert (see text for detail).

2.2.6 Deposition event

The inclusion of the deposition event as a separate external routine to modify the sediment properties (i.e. porewater species, C_org) is a fundamental difference between our approach and the other previous early diagenesis model applied in the Rhone Delta, but it bears similarity with De Borger et al. (2021). We assume the event occurred as an instantaneous deposition of organic carbon ($C_{\mathrm{org}}^{\mathrm{fast}}$ and $C_{\mathrm{org}}^{\mathrm{slow}}$) over a depositional layer, Z_pert (Fig. 2).

The event calculation was carried out dynamically within the same simulation time. For the solid species, following the flux of organic carbon via the boundary condition (see Sect. 2.2.7), the portion of organic carbon is split between the fast and slow decaying component using a proportionality constant (pfast) as in Ait Ballagh et al. (2021). pfast varies from 0 to 1 and it is expressed in percentage of carbon flux deposited associated to either fraction (fast and slow). However, at the time when the event is prescribed, the integrated profile of the solid species $C_{\mathrm{org}}^{\mathrm{fast}}$ and $C_{\mathrm{org}}^{\mathrm{slow}}$ from previous time step, defined as (t⁻), was used to create a virtual composite of the deposited layer. This integral calculation was performed over a specified sediment thickness (Z_pert), which corresponded to the vertical extent of the depositional event. This average concentration for the solid, which we define exclusively for the time of deposition as ($C_{\mathrm{org}}^{\mathrm{flood}}$), is scaled with an enrichment factor (α) (see below) and then nudged on top of the old layer which is supposed to be buried beneath after the event. To avoid numerical issues caused by the discontinuity of both layers with different properties, an interpolation of the composite profile was performed over the modelling domain. This smoothes the interface between the deposited layer's base and the current model grid's upper layer. This algorithmic procedure is schematically shown in Fig. 2 and we summarized this process mathematically as

$$\begin{array}{}\text{(15)}& \begin{array}{rl}{C}_{\mathrm{org}}^{\mathrm{flood}}& \approx \mathit{\alpha }{C}_{\mathrm{org}}^i\left(t\right){\mathrm{|}}_{\mathrm{0}}^{z_{\mathrm{pert}}}\\ & =\mathit{\alpha }{\displaystyle \frac{{\int }_{\mathrm{0}}^{z_{\mathrm{pert}}}{C}_{\mathrm{org}}^i\left(t^{-}\right)\mathrm{d}z}{Z_{\mathrm{pert}}}}\left\{\begin{array}{l}\mathit{\alpha }<\mathrm{1},\text{ if }{\text{TOC}}_{z_{\mathrm{pert}}}<{\mathrm{TOC}}_{\mathrm{old}}\\ \mathit{\alpha }>\mathrm{1},\text{ if }{\text{TOC}}_{z_{\mathrm{pert}}}>{\mathrm{TOC}}_{\mathrm{old}}\end{array}\right.;\\ & i=\mathrm{1}\left(\text{fast}\right),\mathrm{2}\left(\text{slow}\right)\end{array},\end{array}$$

where ${\text{TOC}}_{z_{\mathrm{pert}}}$ corresponds to the TOC content introduced by the flood layer and TOC_old represents the TOC in the previous layer prior to the flood deposition. The carbon enrichment factor (α) in the model (confac in the model code) is introduced here in order to scale the deposited OC with those observed from field data. This helps in calibrating the deposited organic matter concentration ($C_{\mathrm{org}}^{\mathrm{fast}}$ and $C_{\mathrm{org}}^{\mathrm{slow}}$) in the new layer relative to the previous sediment fraction, simulating the wide range of TOC content observed in the field. For instance, when the newly deposited organic matter is similar to the former sediment topmost layer (average pre-flood layer concentration over an equivalent Z_pert depth), an (α) value of 1:1 is used. If the new material is lower in organic carbon content compared to what is near the sediment–water interface, then (α)<1, while if the newly deposited material is higher in carbon content than the sediment surface, (α)>1. This flexibility can be used to constrain the simulation to match the corresponding TOC profile from field observation. In the modelling application, this parameter is generally specified by using different value for the magnitude of OC in each fraction depending on the empirical observation of the TOC data. This quantity is therefore tunable and the upper bound of this parameter is dictated by the maximum TOC in the sediment sample.

It is important to note that this parameter differs from pfast. This OC flux partitioning by pfast occurs regardless of the event and it is related to the carbon flux received at the boundary, but the carbon enrichment factor occurs only during the event. The carbon enrichment factor (α) can be viewed as a method of imposing a new initial condition only at the time of the event by using the integral concentration from the previous time. However, using the approach described here, all calculations can be done dynamically without stopping the model.

For the solutes (O₂, ${\mathrm{NO}}_{\mathrm{3}}^{-}$, ${\mathrm{NH}}_{\mathrm{4}}^{+}$, DIC, ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, H₂S, and CH₄), the bottom water concentration is imposed through the perturbed layer at the time of event by assuming this new layer is homogenously mixed.

2.2.7 Boundary Conditions

The boundary conditions for the model are of three types:

• At the sediment–water interface, a Dirichlet concentration condition for most solutes equalling the bottom water concentration was used,

  $$\begin{array}{}\text{(16)}& C{\mathrm{|}}_{z=\mathrm{0}}=C_{\mathrm{bw}}.\end{array}$$

  Both pore water and solid have a zero-flux boundary condition at the bottom of the model,

  $$\begin{array}{}\text{(17)}& {\displaystyle \frac{\mathrm{d}C}{\mathrm{d}z}}{\mathrm{|}}_{z=z_{\mathrm{\infty }}}=\mathrm{0}.\end{array}$$

  For solid, an imposed flux at the upper boundary for most of the year is used,

  $$\begin{array}{}\text{(18)}& {\mathrm{flux}}_{\mathrm{org}}{\mathrm{|}}_{z=\mathrm{0}}=-\left(\mathrm{1}-{\mathit{\varphi }}_{\mathrm{0}}\right)D{b}_{\mathrm{0}}{\displaystyle \frac{\mathrm{d}C}{\mathrm{d}z}}{\mathrm{|}}_{z=\mathrm{0}}+\left(\mathrm{1}-{\mathit{\varphi }}_{\mathrm{\infty }}\right)w_{\mathrm{\infty }}C{\mathrm{|}}_{z=\mathrm{0}}.\end{array}$$

The model also includes the ability to include time-varying organic carbon flux with user-specific time series or a functional representation such as sinusoidal pattern. In the case study presented here, this carbon flux varies over the annual carbon flux (${\stackrel{\mathrm{‾}}{\text{flux}}}_{\mathrm{org}}$) in the region in question. This was expressed mathematically as

$$\begin{array}{}\text{(19)}& {\mathrm{flux}}_{\mathrm{org}}\left(t\right){\mathrm{|}}_{z=\mathrm{0}}=\left[{\stackrel{\mathrm{‾}}{\text{flux}}}_{\mathrm{org}}\times \left(\mathrm{1}+\sin \left({\displaystyle \frac{\mathrm{2}\cdot \mathit{\pi }\cdot t}{\mathrm{365}}}\right)\right)\right].\end{array}$$

At the time of the instantaneous deposition, this deposited carbon is treated as described in Sect. 2.2.6.

2.2.8 Model parameterization and verification

The model parameters in Table 2 (for full model parameter, see Table S1 in Supplement) were derived from previously published model in the Rhône Delta (Pastor et al., 2011; Ait Ballagh et al., 2021). The organic matter stoichiometry for both fractions is represented here by the NC ratio (NCrFdet and NCrSdet) with values of 0.14 and 0.1, respectively. The flux of carbon in the upper boundary of the model was defined using a yearly mean flux (${\stackrel{\mathrm{‾}}{\text{flux}}}_{\mathrm{org}}$) of 150 mmol m⁻² d⁻¹ in Rhône prodelta (Pastor et al., 2011). TOC (in % dw) is estimated from both carbon fractions ($C_{\mathrm{org}}^{\mathrm{fast}}$ and $C_{\mathrm{org}}^{\mathrm{slow}}$) assuming a sediment density (ρ) of 2.5 g cm⁻³ and conversion from the model unit for detrital carbon fraction of mmol m⁻³ to a unit in percent mass. The model diagnostics TOC value is then computed as follows:

$$\begin{array}{}\text{(20)}& \text{TOC}=\left(C_{\mathrm{org}}^{\mathrm{fast}}+C_{\mathrm{org}}^{\mathrm{slow}}\right)\times \mathrm{1200}\times {\displaystyle \frac{{\mathrm{10}}^{-\mathrm{9}}}{\mathrm{2.5}}}+{\mathrm{TOC}}_{\mathrm{ref}},\end{array}$$

with TOC_ref the asymptotic TOC value at deeper layer of the sediment, thus representing concentration of refractory carbon not explicitly modelled. The sedimentation rate used in this modelling study was kept constant at 0.027 cm d⁻¹ (Pastor et al., 2011). The decay rate constant for the labile and semi-labile detritus matter is set as 0.1 and 0.0031 d⁻¹, respectively, with both fractions split equally with a proportionality constant (pfast) of 0.5. Using parameters fitted by the model of Ait Ballagh et al. (2021) to data observed in the Rhône prodelta area, the rate of bioturbation and bio-irrigation is fixed as 0.01 cm² d⁻¹ and 0.23 d⁻¹ with these fauna-induced activities occurring down to a depth of 5 and 7 cm, respectively.

The bottom water temperature was fixed at 20 ^∘C. The bottom water salinity is nearly constant below the Rhône River plume, ranging from 37.8 to 38.2. In the model, the average temperature and salinity is used to calculate the diffusion coefficient for the solute chemical species (Sect. 2.2.3). Bottom water solute concentrations were constrained using previously reported values in previous modelling efforts (Ait Ballagh et al., 2021) and adapted with new data for the time corresponding to the flood deposit event (see Table 3). Porosity decreases exponentially with depth from 0.9 at the sediment water interface to 0.5 at deeper layer with a decay coefficient of 0.3 cm (Lansard et al., 2009).

Table 3 Core parameters used in the model.

For the verification of the model output, data from Pastor et al. (2018) corresponding to the diagenetic situation 26 d after an organic-rich flood were used. We restricted our benchmark to data from the proximal station (station A) near the river mouth, where the impact of this flood discharge is more visible (Fig. 1).

2.2.9 Numerical integration, application, and implementation

Because the procedure is based on OMEXDIA, complete details of the derivation can be found in that paper and is referenced therein (Soetaert et al., 1996a). Here we recap the mathematical formulation of the method of lines (MOLs) algorithm used by FESDIA. Direct differencing of Eqs. (8)–(10) results to

$$\begin{array}{}\text{(21)}& \begin{array}{rl}{\displaystyle \frac{\partial C_i}{\partial t}}& ={\displaystyle \frac{{\mathrm{\Phi }}_{i,i+\mathrm{1}}{D}_{{\mathrm{\Phi }}_{i,i+\mathrm{1}}}\left(C_{i+\mathrm{1}}-C_i\right)}{{\mathrm{\Phi }}_i\mathrm{\Delta }{z}_{i,i+\mathrm{1}}\mathrm{\Delta }{z}_i}}\\ & -w_{\mathrm{\infty }}{\mathrm{\Phi }}_{\mathrm{\infty }}{\displaystyle \frac{{\mathit{\lambda }}_{i,i+\mathrm{1}}{C}_i+\left(\mathrm{1}-{\mathit{\lambda }}_{i,i+\mathrm{1}}\right)C_{i+\mathrm{1}}}{{\mathrm{\Phi }}_i\mathrm{\Delta }{z}_i}}\\ & -{\displaystyle \frac{{\mathrm{\Phi }}_{i-\mathrm{1},i}{D}_{{\mathrm{\Phi }}_{i-\mathrm{1},i}}\left(C_i-C_{i-\mathrm{1}}\right)}{{\mathrm{\Phi }}_i\mathrm{\Delta }{z}_{i-\mathrm{1},i}\mathrm{\Delta }{z}_i}}\\ & +w_{\mathrm{\infty }}{\mathrm{\Phi }}_{\mathrm{\infty }}{\displaystyle \frac{{\mathit{\lambda }}_{i-\mathrm{1},i}{C}_{i-\mathrm{1}}+\left(\mathrm{1}-{\mathit{\lambda }}_{i-\mathrm{1},i}\right)C_i}{{\mathrm{\Phi }}_i\mathrm{\Delta }{z}_i}},\end{array}\end{array}$$

for a generic tracer C with a phase properties index Φ and D_Φ denoting porosity and dispersive mixing term, respectively, for solid or liquid. This equation is calculated such that the variables and parameters are defined both at the centre of each layer z_i and at the interface between layers $z_{i-\mathrm{1},i}{z}_{i,i+\mathrm{1}}$. The position at the centre of the grid is then given as $z_i=\frac{z_{i-\mathrm{1},i}+z_{i+\mathrm{1},i}}{\mathrm{2}}$. Δz_i represents the thickness of the i layer and $\mathrm{\Delta }{z}_{i,i+\mathrm{1}}$ is the distance between two consecutive grid layers. A Fiadeiro scheme (Fiadeiro and Veronis, 1977) based on the model's Peclet number (a dimensionless ratio expressing the relative importance of advective over dispersive processes) is used to set ${\mathit{\lambda }}_{i,i+\mathrm{1}}$, thus providing a weighted difference of the transport terms which reduces numerical dispersion.

Equations (8)–(10) implemented as Eq. (22) are integrated in time using an implicit solver, called lsodes, that is part of the ODEPACK solvers (Hindmarsh, 1983). This solver uses a backward differentiation method (BDF); it has an adaptive time step, and is designed for solving systems of ordinary differential equations where the Jacobian matrix has an arbitrary sparse structure. The model output time and its time step (dt) is set by the user and is generally problem-specific. Because of the aforementioned challenge in observability of the massive flood event deposition, daily resolution is most often used for user dt. However, there is the possibility of obtaining higher resolution by decreasing dt.

The model application starts by estimating the steady-state condition of the model using the high level command FESDIAperturb(). This steady-state condition is calculated using iterative Newton–Raphson method (Press et al., 1992) and is then used as a starting point for a dynamic simulation, with perturbation times as in perturbTimes and depth of perturbation given as perturbDepth in the model function call. As the event can be given as a deposit and mixing process, further specification of the perturbation type (deposit or mix) is provided as an argument to the simulation routine. In our case, we used only the deposit mode. The event algorithm is used at the stated time point to estimate the model porewater and solid properties driven by the instantaneous change in the boundary condition. The concentrations are successively updated by their diagenetic contributions during this time step. Afterward, this modified profile is integrated forward in time. The model is written in Fortran for speed and integrated using the R programming language (R Core Team, 2021) via the “method of lines” approach (Boudreau, 1996). In addition, the model made use of the event-handling capabilities' specific numerical solvers written in the R deSolve package (Soetaert et al., 2010b). The R programming language is used in the preprocessing routine for model grid generation (Soetaert and Meysman, 2012), porewater chemistry parameter (Soetaert et al., 2010a), steady-state calculation (Soetaert, 2014), and time integration (Soetaert et al., 2010b). Further information about the model usage can be found in the model-user vignette found in R-forge page (https://r-forge.r-project.org/R/?group_id=2422, last access: 2 August 2022).

2.2.10 Quantification of sediment diagenetic relaxation timescale

Quasi-relaxation timescale. Given the strong non-linearity and coupled nature of the biogeochemical system in question, we used an approximate approach to define the timescale of relaxation. Recognizing that in a nonlinear system, a perturbed trajectory is frequently arbitrarily divided into a fast, transient phase and a slow, asymptotic stage that closes in on the attractor (i.e. steady-state concentration; Kittel et al., 2017), we proceeded to estimate the relaxation timescale by using the time for which the memory of the perturbed signature disappears. We estimate the relaxation timescale by first calculating the absolute difference (φ(t)) between successive model output after the event, assuming that a slowly evolving state will eventually converge to the pre-perturbed state as time after the disturbance approaches infinity. This point-by-point concentration difference between two successive discretized profiles is then terminated at the point where the sum of absolute differences at each time point is less than the threshold (i.e. given by the median over the entire time duration). The relaxation time, τ, for each porewater profile species is then defined as the first time this threshold is crossed. A similar technique was employed by Rabouille and Gaillard (1990).

$$\begin{array}{}\text{(22)}& \begin{array}{rl}\mathit{\phi }\left(t\right)& ={\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N\mathrm{|}\mathrm{|}{X}_{t+\mathrm{1}}^i-X_t^i\mathrm{|}\mathrm{|}\\ & \text{In the limit of time }\left(t\right):\\ & \mathit{\tau }\left(t\right)\\ & \Rightarrow \mathit{\phi }\left(t\right)\le \mathit{\phi }{\left(t\right)}_{\mathrm{threshold}}\\ & \text{where }\mathit{\phi }{\left(t\right)}_{\mathrm{threshold}}=\stackrel{\mathrm{‾}}{\mathit{\phi }\left(t\right)}\approx \text{seasonal background},\end{array}\end{array}$$

where N is the total number of grid point (i) used to discretize the depth profile (X_t), and X_t+1 is the depth profile at t+1 after the event.

This relaxation timescale calculation based on the disappearance of the perturbed signal (via successive profile similarity) may differ from an approach in which the profile returns to a pre-defined “old profile”. Because the exactness of pre-flood and post-flood profiles is difficult to quantify numerically (Wheatcroft, 1990), and the return to the former is frequently driven by slow dynamics, the approach used here can provide a window of estimate for which a particular signal fades toward the background of a theoretically pre-perturbed signal.

Uncertainty in relaxation timescale estimate. The uncertainty introduced by this technique is quantified using a non-parametric bootstrap of the τ statistics. The objective of bootstrapping is to estimate a parameter based on the data, such as a mean, median, or any scalar or vector statistics but with less restrictive assumptions about the form of the distribution that the observed data came from (Efron, 1992).

Figure 3 The bootstrapping technique used to calculate the uncertainty in the relaxation timescale. The resampled median about a reference provides a replicate over which the standard error estimate is defined. The solid red represents the expected value of the quantity estimated while the vertical red line is the deviation from this expected value.

In this case, we employ a modified bootstrapping technique to estimate the uncertainty in the relaxation timescale by resampling on the cutoff point introduced in Eq. (23) (i.e. median, $\stackrel{\mathrm{‾}}{\mathit{\phi }\left(t\right)}$ of a given reference simulation). This calculation takes advantage of the fact that the time series will be dominated by the slowly varying seasonal cycle over a long time period away from the point of perturbation, with the influence of the perturbation fading to the background. The variation of this reference time series over time reflects the uncertainty in this median threshold point. This variance, along with the reference cutoff value, can be used to generate n random perturbations varying about the normative threshold value. We can proceed to create a histogram of the replicate threshold(s) distribution. The histogram of this distribution is depicted schematically in the left margin of Fig. 3. The relaxation time in each realization of the threshold is calculated ($\widehat{{\mathit{\tau }}_i}$). The median absolute deviation from this ensemble of relaxation times is then used to calculate the level of uncertainty in the statistics of interest (timescale of relaxation – ($\widehat{\mathit{\tau }}$)). Figure 3 depicts this concept schematically. It should be noted that this method eliminates the need to rerun the deterministic model for each iteration, reducing the computational burden of this technique.

The 95 % confidence intervals (${\widehat{\mathit{\sigma }}}_{\mathit{\tau }\left(\text{level}/\mathrm{2}\right)}$ and ${\widehat{\mathit{\sigma }}}_{\mathit{\tau }\left(\mathrm{1}-\text{level}/\mathrm{2}\right)}$) are reported in this paper by calculating the quantiles of this empirical distribution of ${\widehat{\mathit{\sigma }}}_{\mathit{\tau }}$.

2.2.11 Model simulation

The model is initialized as explained in Sect. 2.2.9. Thereafter, for the dynamic simulation, the model is spin-up for a sufficiently long time to attain dynamic equilibrium (≥5 years). A 2 yr run is carried out for the respective model application. The time step (dt) for dynamic simulation is daily in order to match the frequency for which observation of field data is possible. For specific numerical experiment, model configuration required for the simulation will be detailed in Sect. 2.2.11, “end-member type numerical experiment”.

End-member type numerical experiment

For the numerical model experiment, we investigate the sediment's response to two end-member types of deposition that can represent actual field observations in the Rhône prodelta (Pastor et al., 2018).

• Low OC content with high sediment thickness scenario (EM1). In this scenario, we assume that a 30 cm new layer of sediment of degraded sediment was deposited. This scenario can describe old terrestrial material and is similar to the extreme case of flood event of May/June 2008 in the proximal outlet of Rhône River where lateral transfer of low TOC sediment (around 1 %) was deposited on top of the previously deposited sediment (OC around 1.5 %–3 %) (Cathalot et al., 2010). Using the partitioning of the carbon as explained in Sect. 2.2.7, an α value of 0.5 and 0.7 for $C_{\mathrm{org}}^{\mathrm{fast}}$ and $C_{\mathrm{org}}^{\mathrm{slow}}$, respectively, was used to scale the TOC profile in order to mimic this type of trend.

• High OC with low sediment thickness scenario (EM2). For this, we assumed a moderate 10 cm deposition of a new layer enriched in carbon during a flood discharge event. This scenario can correspond to the end-member case of November 2008 flood type with high TOC around 2.5 %, reaching more than 6 % in some sediment cores from the prodelta (Pastor et al., 2018), (most likely composed of freshwater phytoplankton detritus, debris, and freshly dead organisms) overlain on a less labile layer. In order to simulate this type of pattern, an α value of 20 and 10 for $C_{\mathrm{org}}^{\mathrm{fast}}$ and $C_{\mathrm{org}}^{\mathrm{slow}}$, respectively, was used to adjust the TOC profile to such high-deposit OC scenario.

Except for the α and the thickness of the flood deposit, all other parameters were held constant in all numerical experiments. The time of the event occurrence in both scenarios were initialized at a period corresponding to published dates for May and November 2008 flood deposition as reported in Pastor et al. (2018). This helps to provide some realism to this hypothetical case study and appropriate context to the environmental regime when these events occur.

Sensitivity analysis

Lastly, we conducted a sensitivity analysis of the relaxation timescale for oxygen, sulfate, and DIC concentrations in terms of their variation to the thickness of the new sediment layer and the quantity of organic carbon introduced by the deposition.

We assumed a 15 cm average deposit thickness and conducted simulations with a thickness variation ranging from 1 to 30 cm. A 5 cm thickness increment was used for the sensitivity analysis. The α value is calculated in the same way: assuming a 1:1 ratio in the fast and slow OC fractions, and because deposited sediment can be highly refractory in nature, we geometrically conducted simulations with values ranging from 0.3 to 35. This was done only by changing the α corresponding to $C_{\mathrm{org}}^{\mathrm{fast}}$ with the slow fraction fixed as 1. We also made sure that both series are equilateral in length, and that the values were chosen to span the range of values in EM1 and EM2, thus bracketing the normative value for the end-member case. This range encompasses the large spectrum of flood deposits such as those experienced in the Saguenay fjord, Canada (Deflandre et al., 2002; Mucci and Edenborn, 1992), the Rhône prodelta, France (Pastor et al., 2018), and in the Po River, Italy (Tesi et al., 2012).

3 Results

3.1 Qualitative model performance: Cevenol flood in the Rhône prodelta

In order to compare the model evolution to field data, we made a comparison between the simulated profiles 26 d after a flood layer deposition and data collected in the Rhône prodelta in December 2008 (observed data collected 26 d after a Cevenol flood). During this flooding period, riverine discharge delivered 0.4×10⁶t of sediment which amounted to approximately 10 cm of sediment deposited in the site A of depocenter (Fig. 4).

The general pattern of the simulated profile agrees well with the observed data (Fig. 4). The newly introduced organic carbon-rich sediment resulted in rapid oxygen consumption. The data for total organic carbon (TOC) shown in Fig. 4 suggest a good agreement with the model, with high TOC (2.5 wt %–2.0 wt %) deposited at the upper 10 cm. Twenty-six days after the flood, the oxygen concentration dropped from 250 µM at the new sediment interface to nearly zero at 0.2 cm depth, and oxygen may have already returned to pre-flood levels; the simulated porewater profile was within the data's range (Fig. 4). The model diffusive flux of oxygen at this period was 18 mmol m⁻² d⁻¹ while the measured DOU (diffusive oxygen uptake) flux was 16.6±2.9 mmol m⁻² d⁻¹.

Overall, the model–data trend was satisfactory with observed depth distribution of sulfate (${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$) 26 d after the flood event fitted well, without much parameter fine-tuning. Only the sedimentation rate of the sediment was changed from 0.027 to 0.06 cm d⁻¹ in order to match the observed distribution at depth. Sulfate reduction was high in the new layer. However, below the flood layer, the ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentration in the data seem to asymptote to a value of 10 mM at 25 cm, while the model simulates complete sulfate depletion below 20 cm (Fig. 4).

The DIC profile shows a similar trend to the data collected after the flash flood. Within the depth interval of data, the model tends to follow the data. It drifts at lower depths, on the other hand, by overestimating the concentration of DIC observed at deeper layers. Similarly, the modelled ${\mathrm{NH}}_{\mathrm{4}}^{+}$ shows a gradual increase with depth, and the model overestimates the production of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ below 15 cm (Fig. 4).

Figure 4 Model and observation depth profile of TOC (%), SO${}_{\mathrm{4}}^{\mathrm{2}-}$, DIC,NH${}_{\mathrm{4}}^{+}$, and O₂ for November/December event in station A (Rhône prodelta). The green and red dashed lines depict the vertical depth profile of the model before (3 d) and after (1 d) the flood deposit. The blue solid line represents the model result on the day the observations were collected (26 d after the flood, as indicated by the black circle).

3.2 Numerical experiment on end-member scenarios

3.2.1 Low carbon, high thickness scenario (EM1)

With a test case of 30 cm of new material deposited during the event (EM1) in the spring, the sediment changes as thus: prior to the event, the oxygen penetration depth (OPD) was about 0.17 cm. The OPD increases to 1.17 cm after the deposition of these low OC materials. The model showed a gradual return to its previous profile within days, with the OPD shoaling linearly with time (Table 4). By day 5, oxygen has returned to the pre-flood profile with similar gradient to the pre-flood state.

Against a background OM flux following the introduction of the flood layer, the sediment responded quickly. As a result, the perturbation has a significant effect on sulfate penetration depth, with concentration remaining nearly constant within the perturbed depth (≈20 cm). This corresponds to the bottom water concentration (30 mM) trapped within the flood deposit. Within that layer, sulfate reduction rate was low with an estimated integrated rate of 2.14 mmol C m⁻² d⁻¹ from the surface to 30 cm.

Table 4 Model vs. data comparison for oxygen penetration depth (OPD), flux of oxygen, sulfate, and DIC (26 d after deposition).

Figure 5 Scenario 1 (EM1): model evolution for sulfate following deposition. Relative deviation of successive profile with time shown below. Dashed vertical line signifies cutoff point by the median (dashed horizontal line). Inset: histogram of bootstrap estimate of sulfate relaxation timescale for EM1 with 95 % confidence interval.

Figure 6 Scenario 1 (EM1): model evolution for DIC following deposition. Relative deviation of successive profile with time shown below. Dashed vertical line signifies cutoff point by the median (dashed horizontal line). Inset: histogram of bootstrap estimate of DIC relaxation timescale for EM1 with 95 % confidence interval.

Below the interface with the newly deposited layer, the sediment is enriched in OM whose mineralization results in a higher sulfate reduction rate (SRR) at the boundary that delineates the newly deposited layer and the former sediment–water interface. The simulated SRR falls from 437 mmol C m⁻³ d⁻¹ at the former sediment–water interface (now re-located at 26 cm) to 24 mmol C m⁻³ solid d⁻¹. This high interior sulfate consumption at the boundary correlates well with the higher proportion of reactive organic material buried by the new layer containing less reactive material. From day 10, the consumption of this OM stock by sulfate controls the shape of the profile (Fig. 5). This anoxic mineralization via sulfate reduction will continue until the entire stock of carbon is depleted 50 d after deposition. Following that, OM mineralization via sulfate reduction shift becomes more intense at the top layer by day 60 (2 months after the event), when it begins to gradually evolve to the typical depth-decreasing sulfate profile. By day 115 (∼4 months), the profile had almost completely returned to its pre-flood state. We estimate that it took approximately 4 months for sulfate to relax back to within the range of background variability (with lower and upper bootstrap estimate between 92–139 d).

Correspondingly, OC mineralization products (such as DIC) were significantly lower in the upper newly deposited layer, as a consequence of the reduced quantity of OC brought by the flood. This concentration increased with depth to about 80 mM. Starting from the deposition, higher production of DIC below the former SWI led to a distinct boundary in the sediment: a DIC-depleted layer above an increasing DIC with concentrations up to 75 mM trapped in the region below the new–old sediment horizon 20 d after deposition (Fig. 6). This increased DIC production continued despite complete exhaustion of buried labile fraction with mineralization driven by the slow decaying component. The depth gradient caused by the increased DIC production enhances diffusive DIC flux. Following that initial period, DIC began to revert to its previous state. This slow re-organization, mostly driven by diffusion continues, with an estimated recovery time of 5 months (with a 95 % bootstrap confidence interval of 137–147 d respectively), as it temporarily lags behind ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ in its return to the previous pre-flood state.

3.2.2 High carbon, low thickness scenario (EM2)

A flood deposition scenario of 10 cm thick material with enhanced OC content was used for the other end-member case experiment (EM2) in autumn. In this scenario, the modelled sediment exhibits a variety of response characteristics. The newly introduced sediment resulted in rapid oxygen consumption. The OPD decreased to 0.74 cm shortly after the event, according to the model, and stabilized there for days. There was no visible deformation in the shape of oxygen during its recovery trajectory, and total oxygen consumption for organic matter mineralization decreased by 8 % during the first 2 d after the event, from 12 to 11 mmol O₂ m⁻² d⁻¹.

The ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentration that developed as a result of the deposition showed two gradients: a concentration gradient from 30 mM at the “new” sediment water interface to 26 mM in the newly deposited layer (Fig. 7). Accordingly, the DIC in the corresponding depth layer gradually increased up to 20 mM (Fig. 8). An intermittent increase in ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ was simulated below the new interface, at the boundary with the “old” sediment–water interface (SWI), reaching up to 29 mM from 9 to 12 cm (Fig. 7). This layer, which corresponded to the depth horizon where the new layer gradually mixed with the old layer, resulted in less sulfate reduction and DIC production in comparison to the new layer. Porewater ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations decreased monotonically with depth from this interface, with a corresponding increase in DIC. Within 26 d of the event, the sulfate profile appears to be returning to its original shape. By then, 75 % of the newly introduced fraction of OM had been depleted, with OM remineralization in the upper layer fuelled by the small amount of remaining detrital materials. As the temporal memory of the deposition fades, the profile continues to gradually evolve towards the background, fed by the slow decaying OM, up to day 90, when the sulfate profile appears to have reached a similar pre-flood state. In this scenario, the estimated ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ and DIC relaxation timescales were around 3 months (91 d for ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ and 102 d for DIC) (Fig. 7).

Figure 7 Scenario 2 (EM2): model evolution for sulfate following deposition. Relative deviation of successive profile with time shown below. Dashed vertical line signifies cutoff point by the median (dashed horizontal line). Inset: histogram of bootstrap estimate of sulfate relaxation timescale for EM2 with 95 % confidence interval.

Figure 8 Scenario 2 (EM2): model evolution for DIC following deposition. Relative deviation of successive profile with time shown below. Dashed vertical line signify cutoff point by the median (dashed horizontal line). Inset: histogram of bootstrap estimate of DIC relaxation timescale for EM2 with 95 % confidence interval.

Figure 9 Relaxation timescale (τ) in days as function of deposited sediment thickness and enrichment factor (α) for degradable OM.

3.3 Sensitivity of relaxation time to variation in enrichment factor (α) and sediment thickness (z_pert)

We then examine the sensitivity analysis of the relaxation timescale (τ) for oxygen, sulfate, and DIC for variation in sediment deposit thickness (Z_pert) and the concentration factor for $C_{\mathrm{org}}^{\mathrm{fast}}$ enrichment (α) covering values ranging between the two EM scenarios.

Over all runs varying the enrichment factor (α) and the thickness of the flood input layer, relaxation time for oxygen varied from 2 d for a flood-deposited layer consisting of a thin layer of high concentration of labile OC to 9 d for a thicker deposited layer with low concentration of labile OC. In contrast, the relaxation timescales for ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ and DIC were significantly longer than those for O₂ (3 to 4 months). In addition, the relaxation timescale surface structure for ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ appears complex with divergence gradient at mid-depth of 15 cm. For deposited depth layers above 5 cm and at low α value, the relaxation time for ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ varied between 75–100 d (2–3 months). Below 5 cm (bioturbated depth imposed in the model), relaxation time was constant across all α variations (100 d). As organic enrichment (α) and thickness increase, the model estimates a longer relaxation time with a maximum time span of 105 d.

Similar variation of relaxation time for DIC was simulated for different α and sediment deposit thickness. However, unlike ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, relaxation time for DIC varies smoothly across the range of α and thickness combinations with relatively constant relaxation time (100 d) at low thickness and α combinations. The relaxation time increased exponentially as sediment deposit thickness and labile OC concentration increased (α), with maximum recovery time (171 d/6 months) simulated at the extremes of both combinations (Fig. 9).

4 Discussion

In highly dynamic coastal ecosystems, such as RiOmar (river-dominated ocean margins) systems, driven by seasonal variability and meteorologically extreme events, the response of early diagenetic processes to time varying deposition of organic matter is generally non-stationary (Tesi et al., 2012). While dynamic equilibrium as a steady-state condition may be reasonable in the case of seasonal variability, such an assumption may fail in cases of instantaneously event-driven deposition. An intermittent supply of sediment and OC, like those presented here, can cause a change in the system's properties on a short- or long-term basis. Previous works have highlighted excursions in sediment redox boundary (Katsev et al., 2006), flux of solutes at the sediment–water interface (Rabouille and Gaillard, 1990), and modification of other system properties due to depositional flux of organic matter. Thus, the premise of steady-state conditions in early diagenetic processes which often depends on the temporal resolution of the observation might need revisiting especially in areas of episodic sedimentation (Wheatcroft, 1990; Tesi et al., 2012). Here, we discuss the evolution and dynamics of a non-stationary sedimentary system following a singular perturbation.

4.1 Model representation and utility

Non-steady-state models are increasingly being applied in dynamic coastal environments, but they are still primarily based on forcing from smooth varying boundaries that mimic seasonal forcing or long-term variability (Soetaert et al., 1996b; Rabouille et al., 2001a; Zindorf et al., 2021). Explicit consideration of abrupt changes in the upper boundary of the model caused by events such as landslides, flash flooding, turbiditic transfer of materials on a continental slope, and trawling is still relatively uncaptured by these models (but see De Borger et al., 2021, for inclusion of erosion events). In this paper, we adapt OMEXDIA (Soetaert et al., 1996a), a well-known reaction transport model, to investigate the changes in the solid and liquid phases during massive deposition event. Our efforts highlight the algorithm's utility in incorporating this process with minimal numerical issues. The model represented the basic characteristics of the data derived from the November/December 2008 flood event at station A in the Rhône Delta's depocenter (Fig. 4). The simulated flux was also in agreement with the estimate from field data, as diffusive oxygen uptake (DOU) rate sampled 26 d after the event (8 December 2008) was 16.6±2.9 mmol m⁻² d⁻¹ (Cathalot et al., 2010) while the estimate from the model was 18 mmol m⁻² d⁻¹. As the inclusion of such discontinuity in PDE(s) presents numerical challenges in classic solvers, the implementation utilized by our model ensures such difficulties are overcome. This is the result of improved development of solvers adapted to such a problem (Soetaert et al., 2010b). This difference in the approach employed here distinguishes ours from other published models (e.g. Berg et al., 2003; van de Velde et al., 2018) with similar scientific motivation for time-dependent simulation. Overall, the validation of the model output with field observations lends some confidence in using the model in scenarios involving abrupt changes in boundary conditions and investigating biogeochemical changes in the sediment as a result of such an event. This is despite the model underestimation of the amplitude of sulfate and DIC at depth which can be improved with better optimization of some parameters, especially those derived from previous studies that might not be suited for such a flooding regime or with better process resolution relating to these pathways. Nonetheless, there are advantages to this model especially in the case of episodic flood deposit events, where only a snapshot of data is available at any given time. Modelling tools capable of simulating this event with high fidelity can provide continuous information of the system state and help fill in data gaps needed to understand the sediment's response on different timescales.

4.2 Role of end-member flood input OM in the diagenetic relaxation dynamics

Flooding events can transport large amounts of material through the river to transitional coastal environments such as deltas and estuaries. River floods can account for up to 80 % of terrigenous particle inputs (Antonelli et al., 2008; Zebracki et al., 2015), and they can have a significant impact on geomorphology (Meybeck et al., 2007), ecosystem response, and biogeochemical cycles (Mermex Group et al., 2011). If the source materials have a different organic matter composition (Dezzeo et al., 2000; Cathalot et al., 2013), the rapid deposition of these flood materials can alter diagenetic reactions and resulting fluxes.

Furthermore, the relaxation timescale associated with the sediment recovery following this external perturbation can be important in terms of the process affecting the biogeochemistry of solid and solutes species. With a series of numerical experiments ranging in between two end-members of the input spectrum for flood events such as those in the Rhône prodelta (Pastor et al., 2018), our study revealed contrasting sedimentary responses and associated typical timescales at which porewater profiles relax back to undisturbed state. Using a simple metric for estimating relaxation timescale of the perturbation, our calculations for the first end-member scenario (EM1) show that the upper bound of the timescale of relaxation for oxygen is 5±3 d, whereas it was approximately 2±2 d for the second end-member scenario (EM2). This reflects the property of oxygen, which quickly approaches a steady-state situation after an event (Aller, 1998). This viewpoint is supported by an ex situ controlled laboratory setup. In their studies, Chaillou et al. (2007) demonstrated that after gravity-levelled sediment was introduced, oxygen consumption quickly recovered to its first-day level, with a sharp response time of 50 min and gradual shoaling of OPD within 5 d. We conclude that the tiny difference in oxygen relaxation and diagenetic response between the two scenarios can be attributed to the slow kinetic degradability of the refractory carbon deposited in the first scenario versus the labile nature of the deposit in the second scenario. This kinetically driven OM degradation has been extensively studied and provides the basis for the reactive continuum in early diagenesis models (Middelburg, 1989; Jørgensen and Revsbech, 1985; Burdige, 1991).

Other terminal electron acceptors (TEAs), such as ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, relax toward natural variation over a longer timescale than oxygen. For EM1, our simulation predicts a sulfate relaxation time of 117 d with a 95 % confidence interval (CI) estimate between 92 d (lower CI) and 139 d (upper CI) while in the case of EM2, we estimate a sulfate relaxation time of 91 d with comparatively low temporal variability (lower CI–80 and upper CI–103 d). This difference in relaxation time is caused by the differences in sediment characteristics and how their mineralization occurs over the sediment layer. In the first scenario, organic-rich sediment is buried by less reactive new material. The buried sulfate fraction is reduced faster than in the new layer above and controls the short-term recovery. As the buried carbon stock depletes and the physical imprint of the flood deposition fades, the profile begins to revert to its pre-flood shape. The post-flood evolution for the second scenario (EM2), on the other hand, differs in that the OM is consumed in the classical manner, with decreasing sulfate consumption with depth, caused by top-down control of the OM flux that adds OM to the sediment surface.

Such a long time lapse for the recovery of an element with a complex pathway, such as ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, has been reported in the literature (Anschutz et al., 2002; Stumm and Morgan, 2012; Chaillou et al., 2007). Similarly, estimates from our simulation for each end-member scenario indicate that mineralization products such as DIC have a longer relaxation time. This is especially true for the first scenario as opposed to the second, with evidence of slow convergence at depth within the simulation timescale for the first scenario. We estimate that DIC will recover to its pre-deposition state in 5 months for EM1 and in a comparatively shorter time for EM2 (3 months). This lag in DIC recovery could be attributed to the fact that its post-flood dynamics is governed by the slow decaying detrital material that contribute to the already buried refractory carbon. This long-term quasi-static behaviour of the porewater concentration despite such dynamic introduction of flood input can be understood by introducing the concept of a biogeochemical attractor effect – a similar analogy to the Lorenz attractor (Lorenz, 1963). This idea is derived from the mathematical theory that describes chaos in the real world (Strogatz, 2018; Ghil, 2019). The existence of a “biogeochemical attractor” may explain why multiple temporal data sets in the Rhône River prodelta show a similar diagenetic signature from spring to summer (Rassmann et al., 2016; Dumoulin et al., 2018). Our timescale analysis estimates that such rapid system restoration is indeed plausible and of the correct order of magnitude, based on the range of uncertainty reported here.

In addition, our calculations show that the timescale of return to the previous “pre-flood” profile is bracketed by the range of recovery due to purely molecular diffusion, putting an upper bound on our estimate. For example, using the Einstein's approximation, a species such as oxygen with a sediment diffusion coefficient (D_s) of 1.52 cm² d⁻¹ takes approximately 300 d to be transported solely by diffusion through a 30 cm sediment column and approximately 30 d for a 10 cm sediment column. Similar scaling argument could be made for species such as ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ (D_s=0.86 cm² d⁻¹) with >500 d to be transported through 30 cm and ≈60 d for 10 cm. Because our estimates are less than these values, it suggests that processes other than diffusion (thickness effect) may contribute to relaxation control. It emphasizes the importance of biogeochemistry (OM kinetic) in modulating the response after the event. Besides that, any long-term recovery timescale is governed by the solid deposited. In comparison to the timescale of relaxation roughly estimated from field data (Cathalot et al., 2010), our estimate shows the right order of magnitude.

The relaxation time may also vary depending on the diagenetic interaction, and the characteristics of the organic matter available for degradation. This difference in characteristics was partially imposed in our study by assuming variations of α in the new deposit. The empirical observation of sediment characteristics associated with flood input dictates this parametric turning to match the TOC characteristics (Pastor et al., 2018; Deflandre et al., 2002; Mucci and Edenborn, 1992; Tesi et al., 2012; Bourgeois et al., 2011). However, more data from the field and laboratory experiments that resolve the OM composition of flood deposits are required to constrain the choice of this numerical parameter.

4.3 Control of relaxation time by sediment deposit properties

With the sensitivity analysis, we further explore the variation of relaxation timescale under variation of the thickness of layer and enrichment factor of input material given by α in our model. The model's sensitivity analysis reveals that the thickness and concentration of the reactive fraction of TOC control the relaxation time across a wide range of deposited sediment perturbation characteristics (Fig. 9).

In terms of the recovery time as a function of the availability of labile OC, our results revealed a contrasting pattern for oxygen and sulfate. Several factors related to how different oxidants react with sediment matrix disturbances can explain these differences:

• With oxygen that has a high molecular diffusion coefficient, variations in relaxation time depend on the levels of labile OC, with thin sediments containing a high level of labile OC showing a shorter recovery time than thicker sediments with a low OC content. This pattern can be attributed to the higher relative importance of oxygen consumption in OM-poor sediment relative to the OM-rich sediment.

• For low thickness deposits, sulfate and DIC relaxation times were more or less constant. However, a longer relaxation time was simulated for larger deposits and higher labile OC. This can be attributed to the increased distance required for solutes to migrate back after the event. This is clearly the case for sediment thicknesses greater than 14 cm. Such two-way dynamics could be explained by the fact that biological reworking and physical mixing within the surface mixing layer (SML) can improve OC degradation by promoting the replenishment of electron acceptors (i.e. oxygen, sulfate, nitrate, and metal oxyhydroxides) (Aller and Aller, 1992; Aller, 2004), resulting in a shorter recovery time for the porewater profile to reorganize within the SML.

  This critical depth could also be the distance horizon at which the slow diffusion of the profile when retracting back to its pre-flood profile becomes an important factor in controlling the relaxation timescale. This is especially true for DIC, where the connection is more obvious. It has been proposed that when flood deposits extend beyond the sediment bio-mixing depth, the relaxation time for the constituent species is determined by the concentration gradient between the historical and newly deposited layers (Wheatcroft, 1990). In our sensitivity analysis, higher α corresponds to higher $C_{\mathrm{org}}^{\mathrm{fast}}$ concentrations at depth, resulting in a case of enhanced OC degradation (both at the surface and within the sediment matrix). This depletes electron acceptors such as sulfate, which are required for OM mineralization at this depth. The slow diffusion across the displaced distance, on the other hand, cannot quickly compensate for its demands, which may explain the longer relaxation time. In other words, a higher concentration of OC in a region where all oxidants are nearly consumed results in a profile that takes a relatively longer time to recover to its previous state due to the constraints imposed by oxidant availability. This viewpoint is consistent with previous research from the Rhône prodelta area, where a minimum transport distance of 20 cm is suggested for efficient connection with the SWI; above which several processes are decoupled (Rassmann et al., 2020) and other eutrophic systems, where evidence of large accumulation of organic matter in subsurface sediments serves as a constraint on system restoration (Mayer, 1994; Pusceddu et al., 2009). Indeed, more observational and experimental studies are needed to better understand these processes.

4.4 Model limitations and future development

Because it is based on the well-established OMEXDIA model, FESDIA has several capabilities that make it suitable for a wide range of application domains for non-steady-state early diagenetic simulation. However, due to assumptions made during model development, some limitations in model usage must be considered:

• First, we assumed that porosity is time independent. This may not be the case in some coastal systems that receive sediment materials from regions with distributary channels, which contribute particles of varying origin and grain size (Grenz et al., 2003; Cathalot et al., 2010). The composite sediment that is eventually transported to the depocenter by a flood event may differ in porosity, and thus vary temporarily depending on when and where the source materials are derived during the flood event. In this case, model estimates of fluxes in dissolved species may be over/underestimated. The resulting porosity in the new layer is barely predictable and could range between 0.65 and 0.85 in the proximal zone of the prodelta (Grenz et al., 2003; Cathalot et al., 2010), allowing us to justify our assumption.

• Second, in our examples, we assumed that the burial rate and bioturbation were constant. With the introduction of these flood events, such assumptions may be called into question (Tesi et al., 2012). In addition, benthic animals respond to other perturbation event such as trawling in ways that may warrant explicit description of their recovery, which is linked to bioturbation (De Borger et al., 2021; Sciberras et al., 2018). While some coastal sediment burial rates have been shown to vary seasonally (Soetaert et al., 1996b; Boudreau, 1994), in the proximal zone of the Rhône prodelta, approximately 75 % of sediment deposition can occur during the flood (e.g. 30 cm d⁻¹), with the remaining 25 % distributed throughout the year at a low range daily constant rate (0.03 cm d⁻¹). The dominance of flood deposition over non-flood sedimentation, and the low bioturbation rate observed in the Rhône prodelta (Pastor et al., 2011), prompted the use of constant rate in the application shown here. Moreover, we designed the FESDIA model to allow for the use of a temporarily varying rate constant and coefficient for these processes, and the possibility of imposing an observational time series in cases where such data exist.

  The current FESDIA version does not include a diffusive boundary layer, which can be important for material exchange between the overlying bottom water and the sediment. This is critical for calculating fluxes of species such as O₂, where the depth extent of the DBL (diffusive boundary layer) zone is comparable to the depth at which oxygen consumption occurs (Boudreau and Jorgensen, 2001). As a result, the current version of FESDIA may overestimate the flux of O₂. However, because the primary focus of this paper is on the relaxation dynamics of species (SO${}_{\mathrm{4}}^{\mathrm{2}-}$ and DIC), where the DBL has negligible impact on the relaxation time and overall diagenetic processes (Boudreau and Jorgensen, 2001), the simplification presented here is justified. Even for oxygen, the inclusion of DBL which might result to corresponding change in the concentration at the SWI only have a marginal effect on its relaxation time (<2 d – within the range of uncertainty reported here), so the conclusions drawn in the case studies discussed here are still valid.

In terms of future development, we hope to improve the model's diagenetic pathways, particularly for the iron and sulfur cycles. Furthermore, processes such as calcite formation have been shown to affect DIC profile by 10 %–15 % in the proximal sites of Rhône prodelta (Rassmann et al., 2020), thus might necessitate inclusion in future versions of the model. This will enable FESDIA to account for carbonate system dynamics in marine sediment which can play an important role in the coast carbon cycle (Krumins et al., 2013).

4.5 Relaxation time metric: limitation and perspective

While one main focus of this study is on providing a quantitative estimate of relaxation time, the difficulty of objectively defining what relaxation means necessitates some commentary. This difficulty is not unique to marine biogeochemistry, as accurate quantification of recovery time is an open research question in other fields. In the context of a sedimentary system, Wheatcroft (1990) proposed that determining “dissipation time” (analogous to our “relaxation time”) can be subjective when it comes to signal preservation after sediment event layer deposition. The difficulties are exacerbated by previous work on episodic pulse on sediment biogeochemistry (Rabouille and Gaillard, 1990), in which two metrics for estimating relaxation timescale for silica were proposed. Outside of benthic early diagenesis, Kittel et al. (2017) proposed two generic metrics for systems with well-defined asymptotic properties that can be applied to a distance function from a given target (subject to certain mathematical assumptions). Because porewater profiles are inherently nonlinear, and biogeochemical pathways in sediment are tightly coupled, the mathematical suggestion of asymptoticity using such a distance metric for an evolving profile converging toward the target proposed in that paper is frequently not met. This is the case for our investigation. Overall, while we provide a first-order approximation of relaxation time following perturbation for some model-state variables, these studies highlight also some of the challenges associated with defining the timescale at which a signal can be validly assumed to have returned to its prior state. However, our method allows a full discussion of relaxation times for the main biogeochemical pathways.

5 Conclusion

The need to comprehend extreme events and their relationship to marine biogeochemistry prompted the development of novel methods for diagnosing flood-driven organic matter pulses in coastal environments. In this paper, we propose a new model (FESDIA) for characterizing flood deposition events and the biogeochemical changes that result from them. This type of event can have an impact on the benthic communities and the response of the whole ecosystem (Smith et al., 2018; Bissett et al., 2007; Gooday, 2002). Our modelling study shows that the post-depositional sediment response varies depending on the input characteristics of the layer deposit. For instance, we tested the combined effect of enrichment of labile organic carbon and deposit thickness on space–time distribution and relaxation time of key dissolved species (oxygen, sulfate, DIC). This integral timescale of relaxation is constrained by the intrinsic properties of the solutes (diffusion) and the characteristics of the flood input (thickness and concentration of labile organic carbon). In essence, the findings from this study highlight the importance of the quantity and quality of organic carbon in modulating the sediment response following such a singular perturbation, and the role of flood events with heterogeneous quantitative contributions in the coastal ocean.

Appendix A

A1 Biogeochemical reaction

The full model equation explained in Sect. 2.3.2 is described fully below. Organic matter is composed of three fractions: fast degradable organic matter, slow degradable organic matter, and refractory organic matter. Given the long timescale for the degradability of the refractory OM, it is parameterized using TOC_ref as the asymptotic value. For the two other fractions, five mineralization pathways are included: aerobic respiration (AP), denitrification (DE), dissimilatory iron reduction (DIR), sulfate reduction, and methanogenesis (MG).

Degradation of organic matter:

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}\text{Cprod}& =\left(\text{rFast}\times C_{\mathrm{org}}^{\mathrm{fast}}+\text{rSlow}\times C_{\mathrm{org}}^{\mathrm{fast}}\right)\times {\displaystyle \frac{\left(\mathrm{1}-\mathit{\varphi }\right)}{\mathit{\varphi }}}\\ \text{Nprod}& =\left(\text{rFast}\times C_{\mathrm{org}}^{\mathrm{fast}}\times {\text{NCratio}}_{C_{\mathrm{fast}}}+\text{rSlow}\right.\\ & \left.\times C_{\mathrm{org}}^{\mathrm{fast}}\times {\text{NCratio}}_{C_{\mathrm{slow}}}\right)\times {\displaystyle \frac{\left(\mathrm{1}-\mathit{\varphi }\right)}{\mathit{\varphi }}}.\end{array}\end{array}$$

Limitation terms: the limitation of a mineralization pathway by the availability of the oxidant is modelled by a MOND-type hyperbolic limitation function with inhibition of a pathway represented by a reciprocal hyperbolic function.

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}\text{Oxicminlim}& ={\displaystyle \frac{{\mathrm{O}}_{\mathrm{2}}}{{\mathrm{O}}_{\mathrm{2}}+k_{{\mathrm{O}}_{\mathrm{2}}}}}\times {\displaystyle \frac{\mathrm{1}}{\lim }}\\ \text{Denitrificlim}& ={\displaystyle \frac{{\mathrm{NO}}_{\mathrm{3}}}{{\mathrm{NO}}_{\mathrm{3}}+k_{{\mathrm{NO}}_{\mathrm{3}}}}}\times \left(\mathrm{1}-{\displaystyle \frac{{\mathrm{O}}_{\mathrm{2}}}{{\mathrm{O}}_{\mathrm{2}}+k_{\mathrm{in}{\mathrm{O}}_{\mathrm{2}}\mathrm{den}}}}\right)\times {\displaystyle \frac{\mathrm{1}}{\lim }}\\ \text{Feredminlim}& ={\displaystyle \frac{{\mathrm{FeOH}}_{\mathrm{3}}}{{\mathrm{FeOH}}_{\mathrm{3}}+k_{{\mathrm{FeOH}}_{\mathrm{3}}}}}\times \left(\mathrm{1}-{\displaystyle \frac{{\mathrm{NO}}_{\mathrm{3}}}{{\mathrm{NO}}_{\mathrm{3}}+k_{\mathrm{in}{\mathrm{NO}}_{\mathrm{3}}\mathrm{ano}}}}\right)\\ & \times \left(\mathrm{1}-{\displaystyle \frac{{\mathrm{O}}_{\mathrm{2}}}{{\mathrm{O}}_{\mathrm{2}}+k_{\mathrm{in}{\mathrm{O}}_{\mathrm{2}}\mathrm{ano}}}}\right)\times {\displaystyle \frac{\mathrm{1}}{\lim }}\\ \text{BSRminlim}& ={\displaystyle \frac{{\mathrm{SO}}_{\mathrm{4}}}{{\mathrm{SO}}_{\mathrm{4}}+k_{{\mathrm{SO}}_{\mathrm{4}}}}}\\ & \times \left(\mathrm{1}-{\displaystyle \frac{{\mathrm{FeOH}}_{\mathrm{3}}}{{\mathrm{FeOH}}_{\mathrm{3}}+k_{\mathrm{in}{\mathrm{FeOH}}_{\mathrm{3}}\mathrm{ano}}}}\right)\\ & \times \left(\mathrm{1}-{\displaystyle \frac{{\mathrm{NO}}_{\mathrm{3}}}{{\mathrm{NO}}_{\mathrm{3}}+k_{\mathrm{in}{\mathrm{NO}}_{\mathrm{3}}}}}\right)\\ & \times \left(\mathrm{1}-{\displaystyle \frac{{\mathrm{O}}_{\mathrm{2}}}{{\mathrm{O}}_{\mathrm{2}}+k_{\mathrm{in}{\mathrm{O}}_{\mathrm{2}}\mathrm{ano}}}}\right)\times {\displaystyle \frac{\mathrm{1}}{\lim }}\\ \text{Methminlim}& =\left(\mathrm{1}-{\displaystyle \frac{{\mathrm{SO}}_{\mathrm{4}}}{{\mathrm{SO}}_{\mathrm{4}}+k_{\mathrm{in}{\mathrm{SO}}_{\mathrm{4}}\mathrm{ano}}}}\right)\\ & \times \left(\mathrm{1}-{\displaystyle \frac{{\mathrm{FeOH}}_{\mathrm{3}}}{{\mathrm{FeOH}}_{\mathrm{3}}+k_{\mathrm{in}{\mathrm{FeOH}}_{\mathrm{3}}\mathrm{ano}}}}\right)\\ & \times \left(\mathrm{1}-{\displaystyle \frac{{\mathrm{NO}}_{\mathrm{3}}}{{\mathrm{NO}}_{\mathrm{3}}+k_{\mathrm{in}{\mathrm{NO}}_{\mathrm{3}}\mathrm{ano}}}}\right)\\ & \times \left(\mathrm{1}-{\displaystyle \frac{{\mathrm{O}}_{\mathrm{2}}}{{\mathrm{O}}_{\mathrm{2}}+k_{\mathrm{in}{\mathrm{O}}_{\mathrm{2}}\mathrm{ano}}}}\right)\times {\displaystyle \frac{\mathrm{1}}{\lim }}\\ lim& ={\displaystyle \frac{\mathrm{1}}{\begin{array}{c}\mathrm{\text{Oxicminlim}}+\text{Denitrificlim}\\ +\text{Feredminlim}+\text{BSRminlim}+\text{Methminlim}\end{array}}}.\end{array}\end{array}$$

Depth-dependent kinetic reaction: this limitation is used to reconstruct the vertical distribution of the successive mineralization pathways with a rescaling term “lim” to ensure that the sum of the individual pathway equal the total degradation rate.

$$\begin{array}{}\text{(A3)}& \begin{array}{rl}\text{Oxicmin}& =\text{Cprod}\times \text{Oxicminlim}\times \text{lim}\\ \text{Denitrific}& =\text{Cprod}\times \text{Denitrificlim}\times \text{lim}\\ \text{Feredmin}& =\text{Cprod}\times \text{Feredlim}\times \text{lim}\\ \text{BSRmin}& =\text{Cprod}\times \text{BSRlim}\times \text{lim}\\ \text{Methmin}& =\text{Cprod}\times \text{Methminlim}\times \text{lim}.\end{array}\end{array}$$

Secondary reaction: the re-oxidation of reduced substance and other secondary reactions are modelled with a first-order reaction term.

$$\begin{array}{}\text{(A4)}& \begin{array}{rl}\text{Nitri}& =R_{\mathrm{nit}}\times {\mathrm{NH}}_{\mathrm{4}}\times {\displaystyle \frac{{\mathrm{O}}_{\mathrm{2}}}{\left({\mathrm{O}}_{\mathrm{2}}+k{s}_{\mathrm{nitri}}\right)}}\\ & \left(\mathrm{Nitrification}\right)\\ \text{Feoxid}& =R_{{\mathrm{FeOH}}_{\mathrm{3}}}\times \mathrm{Fe}\times {\mathrm{O}}_{\mathrm{2}}\\ & \left(\mathrm{Iron}\mathrm{oxidation}\right)\\ \text{H2Soxid}& =R_{{\mathrm{H}}_{\mathrm{2}}\mathrm{S}}\times {\mathrm{H}}_{\mathrm{2}}\mathrm{S}\times {\mathrm{O}}_{\mathrm{2}}\\ & \left(\mathrm{sulfide}\mathrm{oxidation}\right)\\ \text{CH4oxid}& =R_{{\mathrm{CH}}_{\mathrm{4}}}\times {\mathrm{CH}}_{\mathrm{4}}\times {\mathrm{O}}_{\mathrm{2}}\\ & \left(\mathrm{Methane}\mathrm{oxidation}\right)\\ \text{AOM}& =R_{\mathrm{AOM}}\times {\mathrm{CH}}_{\mathrm{4}}\times {\mathrm{SO}}_{\mathrm{4}}\\ & (\mathrm{Anaerobic}\mathrm{oxidation}\mathrm{of}\mathrm{methane}\\ \text{FeSprod}& =R_{\mathrm{FeSprod}}\times \mathrm{Fe}\times {\mathrm{H}}_{\mathrm{2}}\mathrm{S}\\ & \left(\mathrm{FeS}\mathrm{production}\right)\\ \text{DICprodCH4}& =-\mathrm{0.5}\times \text{Methmin}+\text{CH4oxid}+\text{AOM}\\ & \left(\mathrm{DIC}\mathrm{production}\mathrm{from}\mathrm{methane}\right).\end{array}\end{array}$$

Removal of sulfide via FeS production and oxidation with oxygen:

$$\begin{array}{}\text{(A5)}& \begin{array}{rl}\text{FeSprod}& =R_{\mathrm{FeSprod}}\times \mathrm{Fe}\times {\mathrm{H}}_{\mathrm{2}}\mathrm{S}\\ & \left(\mathrm{FeS}\mathrm{production}\right)\\ {\mathrm{H}}_{\mathrm{2}}{\mathrm{S}}_{\mathrm{oxid}}& ={\mathrm{rH}}_{\mathrm{2}}{\mathrm{S}}_{\mathrm{oxid}}\times {\mathrm{O}}_{\mathrm{2}}\times {\mathrm{H}}_{\mathrm{2}}\mathrm{S}\\ & \left(\mathrm{Sulfide}\mathrm{oxidation}\right).\end{array}\end{array}$$

Rate of change in state variable:

$$\begin{array}{}\text{(A6)}& \begin{array}{rl}{\displaystyle \frac{\partial C_{\mathrm{org}}^{\mathrm{fast}}}{\partial t}}& =\text{transport}+\text{rFast}\times C_{\mathrm{org}}^{\mathrm{fast}}\\ {\displaystyle \frac{\partial C_{\mathrm{org}}^{\mathrm{slow}}}{\partial t}}& =\text{transport}+\text{sFast}\times C_{\mathrm{org}}^{\mathrm{slow}}\\ {\displaystyle \frac{\partial {\mathrm{O}}_{\mathrm{2}}}{\partial t}}& =\text{transport}-\text{Oxicmin}-\mathrm{1.5}\text{Nitri}\\ & -\mathrm{0.25}\text{FeOxid}-\text{2H2Soxid}-\text{2CH4oxid}\\ {\displaystyle \frac{\partial {\mathrm{NH}}_{\mathrm{3}}}{\partial t}}& =\text{transport}+{\displaystyle \frac{\left(\text{Nprod}-\text{Nitri}\right)}{\left(\mathrm{1}+\text{NH3ads}\right)}}\\ {\displaystyle \frac{\partial {\mathrm{NO}}_{\mathrm{3}}}{\partial t}}& =\text{transport}-\mathrm{0.8}\text{Denitrific}+\text{Nitri}\\ {\displaystyle \frac{\partial {\mathrm{CH}}_{\mathrm{4}}}{\partial t}}& =\text{transport}-\text{DICprodCH4}\\ {\displaystyle \frac{\partial \text{DIC}}{\partial t}}& =\text{transport}+\text{Cprod}+{\mathrm{DICprodCH}}_{\mathrm{4}}\\ {\displaystyle \frac{\partial \mathrm{Fe}}{\partial t}}& =\text{transport}+\mathrm{4}\times \text{Feredmin}\\ & -\text{Feoxid}-\text{FeSprod}\\ {\displaystyle \frac{\partial {\mathrm{FeOH}}_{\mathrm{3}}}{\partial t}}& =\text{transport}+\left(\text{Feoxid}-\mathrm{4}\text{Feredmin}\right)\\ & \times {\displaystyle \frac{\mathit{\varphi }}{\left(\mathrm{1}-\mathit{\varphi }\right)}}\\ {\displaystyle \frac{\partial {\mathrm{H}}_{\mathrm{2}}\mathrm{S}}{\partial t}}& =\text{transport}+\mathrm{0.5}\text{BSRmin}-\text{H2Soxid}\\ & -\text{FeSprod}+\text{AOM}\\ {\displaystyle \frac{\partial {\mathrm{SO}}_{\mathrm{4}}}{\partial t}}& =\text{transport}-\mathrm{0.5}\text{BSRmin}\\ & +\text{H2Soxid}-\text{AOM}.\end{array}\end{array}$$