The evolution of early diagenetic processes at the Mozambique margin during the last glacial-interglacial transition

The Mozambique continental margin experienced large variations in sedimentation rates, primarily due to re-routing of sediment deposition from the Zambezi River during the last glacial-Holocene transition. As changes in sediment accumulation and organic matter deposition impose a strong control on the formation of authigenic minerals in the sediment, the distribution of these minerals may reﬂect the regional paleoenvironmental and paleoclimatic evolution. Combining geochemical analyses of porewaters and sediments with a reactive transport modeling approach, we reconstruct the depositional history and its eﬀect on pyrite formation and other biogeochemical transformations at a site on the Mozambique margin over the past 27 kyr. Fitting the model to match the observed geochemical patterns, most importantly authigenic pyrite, allowed for the recon-struction of past sulfate-methane transition zone depth, which migrated in response to changes in the sediment accumulation and organic matter deposition. Changes in sediment deposition quickly aﬀected organoclastic sulfate reduction and associated pyrite formation, but the eﬀect on anaerobic methane oxidation and subsequent pyrite formation occurred with a lag on the order of thousands of years. Model results reveal a transition from high diagenetic reaction rates representative of near-shore depositional environments during the late glacial maximum, to a setting typical of oﬀshore sediments with low reaction rates at the present day. Notably, the remnants of methane and dissolved iron pools produced in the past still shape the diagenetic processes at and below the sulfate-methane transition zone today. Since deglacial shelf-ﬂooding and corresponding changes in sediment deposition occurred along continental margins worldwide, our analysis highlights the important role of non-steady state diagenesis in continental margin sediments and its relevance for paleoceanographic interpretation of sediment cores experiencing strong variations in sediment input. (cid:1) Published by Elsevier Ltd.


INTRODUCTION
The siliciclastic and organic material delivered by rivers to the ocean accumulates at deposition centers near river https://doi.org/10.1016/j.gca.2021.02.024 0016-7037/Ó 2021 Published by Elsevier Ltd. mouths. The large influx of organic matter (OM) leads to fast biogeochemical turnover rates in these environments, marked by high sediment O 2 uptake (Glud, 2008;Pastor et al., 2011;Pastor et al., 2018). High sediment accumulation rates (SAR) can increase the export of relatively reactive organic matter into deeper sediments (Mü ller and Suess, 1979;Betts and Holland, 1991;Canfield, 1994), stimulating anaerobic turnover rates of OM. This includes the dissimilatory reduction of metal oxides (e.g., Aller et al., 1986;Taillefert et al., 2017;Pastor et al., 2018), sulfate reduction (Canfield, 1989), and methanogenesis in deeper sediments (Martens and Klump, 1984;Zhuang et al., 2018). Anaerobic oxidation of upward-diffusing methane (AOM) coupled to the reduction of downwarddiffusing sulfate (SO 4 2À ) creates a sulfate-methane transition zone (SMTZ; Martens and Berner, 1977;Reeburgh, 1980;Niewö hner et al., 1998;Boetius et al., 2000) where a number of biogeochemical transformations take place, including the release of alkalinity and free sulfide (HS À ) into the porewater. The free sulfide released during AOM can react with iron (oxyhydr)oxides or dissolved Fe 2+ , leading to the formation of iron sulfide minerals (Berner, 1970;Berner, 1984). In systems with significant rates of sulfate reduction coupled to AOM and where reactive iron persists at the depth of the SMTZ, the resulting accumulation of such minerals (mainly in the form of pyrite, FeS 2 ) can serve as a proxy for the depth and migration of the SMTZ in sediment archives (März et al., 2008;Borowski et al., 2013;Roberts, 2015). The position of the SMTZ in marine sediments tends to shift in response to changes in OM loading and SAR (e.g., Hensen et al., 2003;Kasten et al., 2003) as a higher methane generation at depth can push the SMTZ upwards, while lower OM loadings and sedimentation can cause a deepening of the SMTZ (Contreras et al., 2013). Through their connection to the depth of the SMTZ, iron sulfide minerals and their isotopic signature observed in sediment cores can thus provide an important window into environmental conditions of the past. A major challenge, however, is to disentangle the impact of early diagenetic alterations from the climatic factors that are imprinted in the observed geochemical signal.
Non-steady state diagenesis in response to changes in paleoclimate and sedimentation patterns has been identified at several continental margins with major rivers like the Amazon or Rio de la Plata Riedinger et al., 2005) and at the Mozambique margin off the Zambezi (März et al., 2008;März et al., 2018). In this study, we address the impact of early diagenetic alteration and quantify the effect of depositional history on the paleoclimatic signatures within a 33 m long sediment core taken off the shelf-break on the Mozambique continental slope. We fit an early diagenetic model that is subject to transient depositional forcing to the measured geochemical profiles. The model reconstructs the early diagenetic history over the last 27 kyr on the continental slope and evaluates the effect of changes in OM input and SAR, evoked by global paleoclimatic transitions, on the depth of the SMTZ and the signature of authigenic mineral accumulations.

Site description and sediment coring
Site MOZ4-CS17 was cored on the Mozambique margin off the shelf-break $85 km off the Zambezi river mouth at 19°07 0 40.8 00 S and 37°01 0 43.6 00 E. The site is situated within the deposition center of the Zambezi during Pleistocene low sea-level stand, but is not reached by the bulk of sediments deposited in modern times of high sea-level stands (van der Lubbe et al., 2014;Jouet and Deville, 2015).
The sediment core MOZ4-CS17 was collected with a Calypso piston corer aboard the R/V ''Pourquoi pas?" at 550 m water depth during the PAMELA-MOZ4 cruise on December 1, 2015 ( Fig. 1; Jouet and Deville, 2015). The 33 m long core was cut into 1 m sections immediately after recovery on deck.

Sampling
Immediately after core cutting, 3 cm 3 of sediment was collected at the bottom of each section using a pre-cut syringe. The samples were directly transferred into 20 mL glass vials with 5 mL 1 M NaOH, sealed and stored head-down at 4°C until measurement of CH 4 concentrations. The presence of H 2 S was measured at the bottom of each section using a Unisense sulfide electrode mounted with a steel needle with a detection limit of 0.3 mM. The electrode was calibrated with standards up to 200 mM. Each core section was then transferred to a cold room (at 8°C). Porewaters ($0.5 to 7 mL) were extracted with Rhizon soil moisture sampler devices (10 cm length) inserted directly into small holes drilled into the core liner (3 per core section) down to the zone where no water could be extracted anymore (around 17.5 m sediment depth). Depths of all samples in the core liner were corrected using the CINEMA2 software (Woerther et al., 2012), which allows estimating sediment disturbances by monitoring the coring kinematics.

Porewater and gas analyses
Methane concentrations were measured on headspace gas by gas chromatography with a PR2100 gas chromatograph equipped with a flame ionization detector (GC/FID Perichrom, France) connected to a headspace injector (dani HSS 86.50; Sarradin and Caprais, 1996). Methane concentrations of the headspace gas (in ppmV) were converted to porewater concentration (in mM) using sample volume and measured porosity. The relative standard deviation for this method is 3% and the limit of quantification is 0.1 mM. The stable carbon isotopic ratio of methane (d 13 C-CH 4 ) was determined with a G2201-i CRDS analyzer from Picarro. Isotopic composition is given in permil (‰) in standard delta notation relative to the Vienna Pee-Dee Belemnite (V-PDB) composition. The global analytical precision was 0.1 permil with 2 replications per sample. To calibrate the d 13 CH 4 measurements, four different isotopic standards were used (Isometric Instruments, Inc., product numbers L-iso1, B-iso1, T-iso1, and H-iso1 at respectively À66.5‰, À54.5‰, À38.3‰ and À23.9‰).
A porewater split was acidified (to a final concentration of $0.4% HNO 3 ) and stored in acid-washed polypropylene vials at 4°C before SO 4 2À and Fe 2+ analysis onshore. Sulfate was measured with an ion-exchange chromatograph Dionex ICS-5000 from Thermo Scientific Ò , equipped with an Ionpac AG 22-SC precolumn, an Ionpac AS22-SC column, and an electrical conductivity detector. The eluent was Na 2 CO 3 (5 mM) / NaHCO 3 (1.75 mM). IAPSO standard seawater was used for calibration and as a certified reference. Analyses were performed in triplicates and precision was better than 2%. Quantification limit was 0.2 mM. Concentrations of dissolved Fe 2+ were analyzed by high resolution inductively coupled plasma mass spectrometry (HR-ICP-MS Element XR, Thermo Scientific). The calibration was based on the standard addition method of known amounts of each element in the IAPSO standard seawater to avoid matrix effects. Indium (2 ppm) was added in the dilution solution of 2% HNO 3 to correct for instrument drift.

Solid phase characterization
Sediment samples for geochemical analyses were freezedried, weighed pre-and post-drying to determine porosity and homogenized with an agate mortar and pestle.
For total organic carbon content (TOC), around 250 mg of samples was first acidified with successive additions of HCl 1 N at 40°C until complete carbonate dissolution and rinsed twice with milli-Q water to remove residual acid. After freeze-drying, TOC was determined using a LECO TruMac Ò CNS auto-analyzer while d 13 C of the sedimentary OM was determined using Combustion Module-Cavity Ring Down Spectroscopy (CM-CRDS -Picarro; Balslev-Clausen et al., 2013). For the latter, calibration was performed using International Atomic Energy Agency (IAEA) reference materials: calcite (NBS-18), sucrose (CH-6) and lithium carbonate (LSVEC). An acetalinide standard (Costech) was inserted every ten samples to correct for the drift, with a measured value of À33.50 ± 0.02‰ (6 replicates). Precision was typically within 0.03‰ for a triplicate.
Pyrite was quantified with a Bragg-Brentano (D8 Advance BRUKER) X-ray Diffractometer equipped with a Cu X-ray tube and a scanning detector with a Cu-kb filter in Ni (Vantec, BRUKER model). Measurements were made in 0.01°steps at an angle of 5-70°. Each step lasted 1 s and was performed at 40 kV and 30 mA. Pyrite was also extracted following the protocol by Canfield et al. (1986) for S isotope measurements. Approximately 0.5 g of dried sediment was boiled with 12 mL 50% HCl under a flow of ; arrows indicate direction of sediment transport during Pleistocene low sea-level (light blue); and during modern high sea-level (red) (Schulz et al., 2011); Thickness (increasing from yellow = 0-400 m to orange = 400-800 m to gray > 800 m shaded areas) of Pleistocene lowstand sediment deposition center centers inferred from the seismic investigations (Walford et al., 2005); modern sediment deposition (green mud-belt interpreted by van der Lubbe et al., 2014); Paleo-Zambezi incised valley (yellow dotted area; Walford et al., 2005). The 100 m isobath approximately indicates dry shelf area during Pleistocene sea-level lowstand $20 kyr BP (brown hues). Bathymetric map created with Ocean Data View (Schlitzer, 2015). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) nitrogen gas to dissolve the operationally defined acid volatile sulfur (AVS) fraction which consists of Fe-monosulfides (predominantly mackinawite) and greigite (Fe 2 S 3 ). After AVS extraction, 6 mL CrCl 2 were added to the boiling mixture to dissolve pyrite (FeS 2 ). In both steps, the released H 2 S was trapped with AgNO 3 and the precipitated Ag 2 S was stoichiometrically converted into moles S and sulfurbound Fe. Standard error of this method was within 10% for three samples analyzed in triplicates and an inhouse standard (WHIT; Alcott et al., 2020). However, the standard error exceeded 30% for a fourth sample with a very low pyrite content. Sulfur isotopes ( 32 S and 34 S) of pyrite were analyzed in the Ag 2 S precipitates. Around 200 lg of Ag 2 S were analyzed using an Elementar Pyrocube coupled to a GC Isoprime stable isotope ratio mass spectrometer. Results are given in the standard delta notation as d 34 S relative to the Vienna Canyon Diablo Troilite (V-CDT) international reference material and each sample was measured in duplicate. Samples were calibrated to the V-CDT scale using a lab seawater sulfate standard (SWS-3), an interlab chalcopyrite standard (CP1), and the international standard IAEA-S-3 with assigned values of +20.3, À4.56 and À32.06‰ respectively. SWS-3 was calibrated against the international standards (assigned values in brackets) NBS-127 (+20.3‰), NBS-123 (+17.01‰), IAEA S-1 (À0.30‰) and IAEA S-3 (À32.06‰). Repeat analyses of a barium sulfate check standard produced a standard deviation of <0.3‰.
Reactive iron, i.e., mainly poorly crystalline and crystalline Fe-(oxyhydr)oxides, was semi-quantitatively determined by dithionite extraction (Kostka and Luther, 1994;Anschutz et al., 2000), followed by ICP-AES analysis of the extraction solution using a Horiba Jobin Yvon Ò Ultima 2 spectrometer. A certified sediment reference material (MESS-4, NRC-CNRC) was used to check the repeatability of the extractions through time. The global variability of MESS-4 was around 5% (standard error).

Chronology and sediment accumulation rates
The chronology of the core MOZ4-CS17 is based on 14 C dating in 16 depth horizons. To this end, sediments were sampled at regular intervals along a continuous and nonreworked typical hemipelagic sequence composed of fine clayey-silt typical for deltaic margins. The disturbancecorrected depth below seafloor was used to estimate SAR.
Sediment was dated using Accelerator mass spectrometer (AMS) standard radiocarbon methods on individual marine mollusk shells and the bulk assemblage of planktonic foraminifera. Analyses were performed at the Beta Analytic Laboratory (Florida, USA) and all 14 C ages were calibrated using the MARINE13.14C calibration curve (Reimer et al., 2013). A local marine reservoir correction of mean DR = 158 ± 42 yrs (Reimer and Reimer, 2001) inferred from radiocarbon measurements in prebomb known-age shells and corals from the tropical SW Indian Ocean (Southon et al., 2002) was applied.
Sediment ages were interpolated in order to build an accurate age model using the age-depth modeling software CLAM 2.2 (Blaauw, 2010) in R 3.4.3. (R Core Team, 2017).

Model structure
The model simulates early diagenetic processes over the past 27.4 kyr. The governing equation accounting for diffusion, burial, and reaction of solutes is and that of solids is where u is the sediment porosity, C i and C j are the concentrations of dissolved and solid species, respectively, z is the depth in the sediment, h is the tortuosity, D b is the biodiffusion coefficient, v and w are the burial velocities of solutes and solids, respectively, and s i k and s j k are the stoichiometric coefficients of the reaction rates R k (Boudreau, 1997). The biodiffusion coefficient is described as an exponentially decaying function where D b0 and q are given in Table 1. The ionic/molecular diffusion coefficients were corrected for in-situ pressure, temperature, and salinity conditions (Boudreau, 1997), and corrected for tortuosity, which was described as a function of porosity (Boudreau, 1996): Assuming steady-state compaction, the porosity is constant over time and described as an exponentially decreasing function with depth where u 0 and u 1 are the porosities at the sediment-water interface and at depth, respectively, and c is a coefficient fitted to measured porosity profiles (Table 1). The effect of steady-state compaction on burial velocities of solids (wðz; tÞ) and solutes (vðz; tÞ) is described by : reflecting that solutes move upward relative to the surrounding sediment matrix due to compaction (Boudreau, 1997). Below the depth of compaction (in the consolidated sediment layer) the burial velocities of solids and solutes equal SAR , which allows the velocities of solutes and solids to be related through v 1 ¼ w 1 . SAR corrected for compaction can be determined from the radiocarbon ages (Table 4) by rearranging equation (6), inserting equation (5) for the porosity, and solving where ðt a ; z a Þ and ðt b ; z b Þ denote the age and depth from two adjacent data points from the dated sediment sequence, and w 1 ðt a < t < t b Þ is SAR for a given time-interval; thus yielding piecewise constant w 1 values between measured ages. The solids in the model include organic matter, Fe-(oxyhydr)oxides and pyrite. OM is available in three fractions, where OM a is the most labile fraction, OM b is less reactive, and OM c is also less reactive and only reacts with O 2 (Tables 1, 2). Fe-(oxyhydr)oxides are represented by Fe 2 O 3 and are divided into two fractions, representing varying degrees of crystallinity. The second fraction Fe 2 O 3 b does not react with OM in the model in order to allow for Fe 2 O 3 to persist through and below the zones of organoclastic sulfate reduction and methanogenesis (Table 2). Authigenic Fe sulfides are represented by pyrite (FeS 2 ). The model also includes dissolved O 2 , SO 4 2À , CH 4 , HS À and Fe 2+ . This set of state variables was selected to represent the main aspects of C, O, Fe and S cycling at the study location.
The reaction network (Table 2) comprises 12 reaction pathways occurring during sedimentary organic matter degradation and authigenic mineral formation. OM miner-alization can be coupled to O 2 (R1), Fe 2 O 3 (R2) and SO 4 2À (R3) reduction, and methanogenesis (R4). Following the observations compiled by Tromp et al. (1995), aerobic respiration (R1) was modeled with a 100-fold faster rate constant than R2-4 (Table 2). Methane can be oxidized aerobically (R5) and anaerobically (coupled to SO 4 2À reduction, R6). The model accounts for the reoxidation of reduced metabolites (HS À , R7; Fe 2+ , R8) and sulfidic reduction of Fe 2 O 3 (R9). During authigenic pyrite formation (R10) H 2 can be produced (Rickard, 2012), which is highly reactive and therefore does not accumulate in porewaters . Similar to the commonly used model formulation of organic matter mineralization, in which H 2 /acetate production and subsequent reoxidation steps are implemented as single reaction pathways, the H 2 produced during pyrite formation is in our model directly coupled to the reduction of Fe 2 O 3 (R10a), SO 4 2À (R10b) or methanogenesis (R10c). Pyrite can react with oxygen (R11). More reactive, less crystalline Fe 2 O 3 a can convert into less reactive, more crystalline Fe 2 O 3 b (R12). As upper boundary conditions, the influxes of OM and Fe 2 O 3 at the sediment-water interface are modeled as linear functions of SAR Parameters A and B are generally kept constant (Table 1), except for OM a and OM b where values are increased by 75% during two periods, i.e., from 21.43 to 20.43 and from 12.54 to 10.14 kyr BP (Table 3; see Section 2.3.2 for an explanation). All solid Fe is initially deposited in the (less crystalline) alpha-phase (Fe 2 O 3 a ). The influx of pyrite is zero. For dissolved species, concentrations are fixed at the sediment-water interface, representative of conditions in the bottom water. The upper boundary concentration of O 2 is 0.225 mM, the concentration of SO 4 2À is 28 mM. All other solutes have an upper boundary concentration of 0. For all solids and solutes, zero-gradient boundary conditions are imposed at the bottom of the model domain. All CH 4 is produced in situ.
Initial concentrations for all dissolved and solid species were set to 0. The model was first run to steady-state with a constant SAR of 1 m/kyr. Since we cannot constrain the conditions during the spin-up (e.g., SAR, OM deposition), the first 5 kyr after the spin-up (27.4-22.4 kyr BP) are treated with caution and we only discuss model results produced for sediments younger than 22.4 kyr BP, which corresponds to the upper 20 m of the sediment core.
The size of the grid cells in the model increases exponentially from 1 mm at the top to 10 cm at 4.7 mbsf and remains constant from there on to the bottom of the domain at 29.8 mbsf. The time steps are variable, but limited to a maximum of 1 year. The model implementation is an improved version of the early diagenetic modeling framework presented in Rooze et al. (2020), which discretizes partial differential equations (1) and (2) with the finite volume method and couples transport to reaction rates with the sequential iterative operator-splitting approach (Steefel and MacQuarrie, 1996). The code is written in R and C, it uses CVODE from the SUNDIALS package (Hindmarsh et al., 2005) to resolve reaction rates, and 1.82 * 10 À3 m 2 yr À1 b q 2 * 10 À2 m b k a 3.3 * 10 À2 yr À1 a k b 2.01 * 10 À4 yr À1 a k c 2.01 * 10 À4 yr À1 a K mO2 8 * 10 À3 mol m À3 c K mSO4 1 * 10 À1 mol m À3 c k 5 1 * 10 7 m 3 mol À1 yr À1 c k 6 1 * 10 4 yr À1 a k 7 1 * 10 7 m 3 mol À1 yr À1 d K mSO4AOM 1 * 10 À1 mol m À3 c K mFeSolid 1.68 * 810 2 mol m À3 c k 8 1.4 * 10 6 m 3 mol À1 yr À1 c k 9 1 * 10 À1 m 3 mol À1 yr À1 a k 10 1 * 10 3 m 3 mol À1 yr À1 e k 11 1 m 3 mol À1 yr À1 a k 12  Table 3).

Transient modeling approach
The initial model parameterization was developed based on literature values (see Table 1), combined with adjustments for local environmental conditions, the measured porosity profile, and a preliminary fit of the total organic carbon (TOC) profile. Subsequently, model parameters were adjusted to reproduce the observed contemporary predominantly ferruginous state of the porewaters, to reasonably fit the measured solid phase profiles, and to match the present depth of the SMTZ.
The influx OM and Fe 2 O 3 depend linearly on SAR. The OM a phase depends both on SAR and a constant input (nonzero value for B, Eq. (9); Table 1), which allows for sufficiently high OM deposition to prevent too deep O 2 penetration during the Holocene, when SAR is very low. The Fe 2 O 3 a phase also has a non-zero constant component (Table 1), causing an increase in Fe 2 O 3 concentrations near the sediment-water interface during the Holocene.
When the OM loading strictly follows the linear relationship (Eq. (9)) only one pyrite peak is formed (see Section 4). To fit two additional peaks observed in the data, the OM a and OM b deposition from 21.43 to 20.43 and from 12.54 to 10.14 kyr BP were increased by 75%.

Geochemistry
Sediments at Site MOZ4-CS17 are mainly homogenous, featureless, hemipelagic mud of dark greenish to brownish gray color. The top 2 mbsf are of coarser grain size than the remainder of the core. From the top of the core, SO 4 2À concentrations drop nearly linearly from seawater values (28 mM) to below detection limit at 12 mbsf ( Fig. 2A). Throughout the core, no free sulfide was detected. While CH 4 has been detected in the topmost sample where it reaches 0.3 mM at 0.9 mbsf, it stays below 1 mM throughout nearly the entire sulfate-containing zone and only starts to rise again below 10 mbsf. It reaches 1.6 mM at 14 mbsf and stays elevated below that at around 1 mM, close to the methane saturation under shipboard conditions. An SMTZ is identified between $10 and $12 mbsf ( Fig. 2A, B). Dissolved Fe concentrations are highly variable between 0 and 120 mM (Fig. 2C). Dissolved Fe is prevalent in two separated zones above and below the SMTZ. Within the SMTZ, dissolved Fe persists at low val- Table 2 Reaction network. No.

Stoichiometry
Kinetic Rate Subscripts 'a' and 'b' denote the C:N ratio in organic matter; a = 106, b = 16 (Redfield, 1958). Superscripts 'a' and 'b' indicate reactive and less reactive fractions of organic matter and iron oxides, whereas the 'c' fraction only reacts with oxygen. detected in the topmost sample exhibits the least negative measured d 13 C of À75‰ (Fig. 2D). Sediment porosity gradually declines from 0.7 at the top of the core to 0.6 at the base (Fig. 2E). Pyrite is more abundant in the upper part of the record (Fig. 2F). The XRDderived pyrite profile exhibits three maxima: around 3.5 mbsf (1.3 wt%); within the current SMTZ around 10 mbsf (2.1 wt%); and around 16 mbsf (1.8 wt%). The chromium reduction derived pyrite profile slightly differs from this pattern. The peak in the current SMTZ is not very pronounced and only reaches 1.2 wt% whereas the peak at 16 mbsf is a bit higher with 2.1 wt%. The differences can partly arise from different analysis times. XRD analysis was undertaken directly after the cruise and chromium reduction has been performed approximately 1.5 years later. Thus, postsampling oxidation could have affected the CRS profile.
The XRD derived pyrite is therefore considered more reliable. Pyrite d 34 S is most negative at the sediment-water interface (À40‰) and gets steadily heavier downcore into the SMTZ where values of À25‰ are reached (Fig. 2G). The pyrite peak at 16 mbsf exhibits the least negative d 34 S of À10‰. Such values are again reached at $25 mbsf, without a corresponding pyrite peak. Reactive Fe ranges between 1.0 and 1.8 wt% (Fig. 2H). The highest values occur near the top and bottom of the core, and the lowest values between 15 and 20 mbsf. The measured TOC content exhibits an increase from 0.6 to 1.1 wt% from 1.2 to 3.5 mbsf, shows a slight decrease between 4 and 6 mbsf, and remains relatively constant at 1.14 ± 0.07% below that (Fig. 2I).
The deepest cored interval at 33 mbsf has an age of 27.33 kyr. SAR over the entire period of deposition varies by an order of magnitude. They reached values between 200 and 450 cm kyr À1 during the Late Glacial Maximum (LGM 26.5-20 kyr BP; Clark et al., 2009), dropped to 100 cm , CH 4 , Fe 2+ , pyrite, dithionite extractable total Fe and total OM, red and blue lines in I) indicate simulated OM fraction OM a (red) and OM b (blue); dashed blue vertical line in B) indicates shipboard CH 4 saturation (Duan and Mao, 2006) (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) kyr À1 between 18 and 16.5 kyr BP, rose to 250 cm kyr À1 again at 16.5 kyr BP, and dropped stepwise until the top of the record. After 7500 kyr BP, SAR remained lower than 20 cm kyr À1 (Table 4; Fig. 3).

Modeling
In a first step, we performed a transient model simulation where OM deposition followed a fixed linear relation with sediment accumulation through time (Fig. 3). This model simulation did not reproduce the three-peak shaped pyrite profile, but showed only one peak (see Section 4 below). We thus tuned the model by varying the OM loading during two specific time periods (Fig. 3), until a good fit was achieved to the observed geochemical profiles. The geochemical depth profiles that are predicted for the present day situation by this ''best fit" transient model simulation are displayed in Fig. 2. The model reproduces the main features as seen in the data, yielding a contemporary SMTZ positioned at 11.2 mbsf ( Fig. 2A, B) and a pyrite profile exhibiting three distinct peaks at 3.3, 8.3 and 17.5 mbsf (Fig. 2F). The simulated Fe 2 O 3 profile captures the main trends in the data (Fig. 2H). Porewaters are ferruginous over most of the sediment column (Fig. 2C) except for near the sediment water interface where oxic conditions prevail (Fig. 2J). Oxygen penetrates $1 cm into the sediment (Fig. 2J). No HS À accumulates. The model provides reaction rate estimates for the present day situation. In the topmost sediment layer, most OM is processed by aerobic degradation, which accounts for 0.57 mol C m À2 yr À1 (92% of the total OM decomposition rate). Organoclastic sulfate reduction accounts for 0.045 mol C m À2 yr À1 (7%), while organoclastic iron reduction and methanogenesis are responsible for 3.4 * 10 À4 mol C m À2 yr À1 (<0.1%) and 0.004 mol C m À2 yr À1 (0.6%), respectively. Methane producing processes (R4 and R10c) are currently generating 0.00315 mol CH 4 m À2 yr À1 while AOM (R6) is consum-ing 0.0126 mol CH 4 m À2 yr À1 . Thus CH 4 consumption (exceeding production in the model domain) occurs at a net rate of 0.0095 mol CH 4 m À2 yr À1 . The maximum rate of organoclastic sulfate reduction (0.5 mol C m À3 yr) is found at 1.8 cmbsf and AOM peaks (0.07 mol C m À3 yr) just above the SMTZ at 11.1 mbsf. Methanogenesis spikes at 1.8 cmbsf (2 * 10 À3 mol C m À3 ) and starts again in the SMTZ (4 * 10 À4 mol C m À3 ) from where it remains elevated over depth. Pyrite is formed in the zones of organoclastic sulfate reduction (7 * 10 À2 mol Fe m À3 yr À1 , 6 cmbsf) and AOM (3 * 10 À2 mol Fe m À3 yr À1 , 11.1 mbsf). Of the depth-integrated pyrite formation rate (1.42 * 10 À2 mol FeS 2 m À2 yr À1 ), 58% is formed in the organoclastic sulfate reduction zone and 42% in the AOM zone.
The ''best fit" transient model simulation suggests that geochemical depth profiles changed appreciably throughout the considered time period. The SMTZ migrated upwards between 21.6 and 17.8 kyr BP from 7.0 to 5.9 mbsf ( Fig. 3; Fig. 4G). From 16.5 to 15.4 kyr BP the SMTZ migrated downwards to a depth of 6.5 mbsf and stayed there till 10.8 kyr BP. Since then the SMTZ has continually moved downwards, and it is currently located at 11.2 mbsf.
The pyrite peak currently located at 17.5 mbsf (Fig. 3E) first emerged 21.4 kyr BP in the zone of organoclastic sulfate reduction, was buried and grew while it moved through the zone of AOM (18.9-16.4 kyr BP, SMTZ located at $6 mbsf), and was further preserved and buried (Fig. 3A-E). On top of pyrite produced in the zone of organoclastic sulfate reduction, a second pyrite peak (currently located 8.3 mbsf) was formed in the zone of AOM from $13 to 8 kyr BP (Fig. 3C, D). During this period the highest concentrations were located at the depth of the SMTZ (where new pyrite is formed). Because the burial velocity exceeded the downward movement of the SMTZ, the formation of new pyrite in this zone resulted in a skewed peak (Fig. 3D). The peak grew rapidly from 7.5 to 4 kyr BP, a period during which the downward movement of the SMTZ exceeded Table 4 Radiocarbon ages performed at the Beta Analytic Laboratory (Florida, USA) and calibrated age ranges at the 95% confidence interval used for the depth-age model. The MARINE13.14C curve (Reimer et al., 2013) was used to calibrate the ages after applying a marine reservoir correction (Reimer and Reimer, 2001;Southon et al., 2002).
burial of the pyrite peak and moved past it. The pyrite concentrations in the lower part of the peak keep increasing to the present day (Fig. 3D, E). The third peak currently located at 3.3 mbsf grew rapidly from 12.5 to 10.2 kyr BP and keeps growing in the zone of organoclastic sulfate reduction ( Fig. 3C-E). Reaction rates integrated over the model domain were higher than today throughout the past 27 kyr (Fig. 4). Organoclastic iron reduction, organoclastic sulfate reduction, methanogenesis methanogenesis and pyrite formation all exhibit maximum rates between 22 and 20 kyr BP, coinciding with the highest SARs and organic matter deposition. At their highest, the reaction rates were were 0.046 mol C m À2 yr À1 (organoclastic iron reduction), 0.75 mol C m À2 yr À1 (organoclastic sulfate reduction), 0.043 mol CH 4 m À2 yr À1 (methanogenesis) and 0.2 mol pyrite m À2 yr À1 (pyrite precipitation). AOM peaks later around 17 kyr BP with values of 0.03 mol CH 4 m À2 yr À1 .

Present day early diagenesis
Sulfate concentrations decrease from 28 mM at the sediment surface to 0 in the SMTZ, with model results suggesting that organoclastic sulfate reduction mostly takes place in the upper $5 m of the sediment column. Deeper, towards the SMTZ, AOM predominates SO 4 2À consumption explaining the nearly linear SO 4 2À profile above the SMTZ (Borowski et al., 1996). The observed SO 4 2À profile is slightly S-shaped which may be caused by changes in lithology and porositya feature that is not captured by the model. Pyrite above the SMTZ (0.9-1.3 wt%) has more negative d 34 S values (À45 to À25‰), indicating that the sulfide was produced during organoclastic sulfate reduction and formed closer to the sediment-water interface (Hartmann and Nielsen, 1968;Jørgensen, 1979). The significant contribution of organoclastic sulfate reduction to the isotopic signal of the pyrite is supported by the model outcomes (see Section 4.3).
Organoclastic iron reduction is a minor contributor to organic matter degradation over the entire simulated period. Yet, the porewaters are mostly in a ferruginous state with abundant dissolved Fe 2+ but no free HS À . In our model, most of the Fe 2+ is produced by sulfidic reduction of Fe 2 O 3 . To titrate out all HS À , iron (oxyhydr)oxides must be able to persist beyond the depth of OM oxidation into the zone of AOM, which is consistent with our measurements (Fig. 2H) and accomplished in the model by the Fe 2 -O 3 b pool that is unreactive towards OM (Reed et al., 2011;Rooze et al., 2016). This emphasizes the importance of iron cycling in the sediment at depths >10 m, and contrasts with findings at the Congo lobes where Taillefert et al. (2017) documented strong organoclastic iron reduction in the top meter of sediment. Fig. 3. Temporal evolution of modeled profiles (in small plots A-F) of sulfate (black), methane (red) and pyrite (blue) in response to sediment accumulation rate (blue), OM deposition (green) and paleo-sea-level (red; Camoin et al., 2004) in the upper plot. Roman numerals indicate phases of OM deposition (see Table 3). Blue shaded background indicates phases of normal OM deposition; white shaded phases indicate phases of increased OM deposition. Gray stippled vertical lines indicate corresponding depth in sediment core MOZ4-CS17. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) The high pyrite and low dissolved Fe 2+ concentrations at the depth of the current SMTZ (11 mbsf, Fig. 2A, B) indicate active pyrite formation, which is supported by the model. The amount of pyrite in this interval relative to deeper in the core suggests that for a longer period of time this SMTZ has been within the same sediment horizon, i.e., that the distance between the SMTZ and the sediment-water interface increases approximately at the speed of SAR. This is possible under the low sedimentation conditions dominating the last 10 kyr. The CH 4 budget shows net consumption, implying that the CH 4 currently diffusing upwards into the SMTZ has been produced in the past and the CH 4 pool is now progressively being depleted. As such, the SMTZ is expected to continue to migrate deeper into the sediment as sulfate diffusion into the sediment outpaces methane production. The d 13 C of CH 4 indicates that it has been produced by biogenic methanogenesis alone with no thermogenic contribution (Schoell, 1980;Whiticar, 1999). This has also been confirmed by Deville et al. (2020) and justifies the model assumption of no methane influx from below.

Variations in sediment accumulation rate and organic matter input
Sea-level has risen by 120 m following the LGM (Fairbanks, 1989;Camoin et al., 2004;Peltier and Fairbanks, 2006), leading to flooding of the previously dry Mozambique continental shelf and a marine transgression shifting the coastline 50-100 km inland (Beiersdorf et al., 1980). Periods of rapid post-glacial and deglacial sea-level rise have been detected as melt water pulses in paleo-sea-level records (Fairbanks, 1989;Bard et al., 1990;Deschamps et al., 2012). At the Mozambique margin sea-level changes have had major impacts on the distance from the Zambezi river mouth to the study site and the near-shore current system (Fig. 1). Changes in climate also affected rainfall, land erosion, vegetation along the river, and connectivity between different parts of the Zambezi catchment (Schefuß et al., 2011;Just et al., 2014). These factors can all have had an impact on the amount and characteristics of deposited OM and the delivery of sediments to our study site. During the last glacial period, sediments from the Zambezi were transported directly to deeper waters of the continental slope and displaced southward by the prevailing ocean currents (van der Lubbe et al., 2014). At Site MOZ4-CS17 highest SARs are determined between 25 and 24 kyr BP, potentially related to the LGM sea-level lowstand prior to the onset of the largescale global sea-level rise (Ramsay and Cooper, 2002).
The meltwater pulse (MWP) 19 kyr ago (Clark et al., 2004) is associated with 15 m sea-level rise from 120 to 105 m below present level, which was not sufficient to raise the sea-level over the Mozambique shelf break which today is between 90 and 100 m water depth (Beiersdorf et al., 1980). Hence, it did not cause a significant shoreline retreat. Thus, it also had a limited effect on SAR at Site MOZ4-CS17 (Fig. 3).
MWP 1A (14.6 kyr BP, during Phase IV, Table 3; Deschamps et al., 2012) is the major deglacial MWP associated with a significant shoreline retreat at the Mozambique margin (Wenau et al., 2020). This event also appears to have had a weak effect on the SAR. The apparent weak imprint of MWP 19 kyr and MWP 1A on SARs determined at Site MOZ4-CS17 (Fig. 3) may in part be due to the resolution of our age model (and potentially also the sea level reconstructions), wherein SARs is averaged  (Table 4). Their effect may also be obscured by high spatiotemporal variability in SARs (Schulz et al., 2011) and deflections of the river plume. Previous studies showed that wetter periods in South-East Africa, i.e., Heinrich Stadial 1 (17-15 kyr BP) and the Younger Dryas (13-11.4 kyr BP), affected the predominant source area and composition of the material transported by the Zambezi (Schefuß et al., 2011;Just et al., 2014). However, these events appear to have neither strongly affected SAR nor the OM loading at MOZ4-CS17.
At the beginning of the Holocene, MWP 1B at 11.3 kyr BP (Bard et al., 1990) led to flooding of a large part of the Mozambique shelf (van der Lubbe et al., 2014). From this time onwards, Zambezi sediments are transported by a strong countercurrent flowing on the shelf along the coast in northeasterly direction (Fig. 1), not reaching the main Pleistocene deposition centers beyond the shelf break anymore (van der Lubbe et al., 2014). This deflection of the main sediment source is reflected in much lower SAR at Site MOZ4-CS17 during the last 10 kyr.

Effect of transient depositional environment on early diagenesis
Simulated pyrite formation in the zone of organoclastic sulfate reduction, near the sediment-water interface, responds almost immediately to increases in OM deposition (Fig. 4C), i.e., FeS 2 peaks are discernible within 100 years of increased OM deposition during Phases II and IV. In contrast, there is a considerable lag between changes in SAR/OM deposition and their effect on AOM-driven pyrite formation in the SMTZ. The SMTZ is located relatively stable, at 5.9-6.5 mbsf depth between 15.4 and 10.8 kyr BP (Fig. 4), despite large changes in SAR during this period. Significant decreases in SAR from 15.3 to 12.7 kyr BP led only considerably later (around 10.8 kyr BP) to the downward movement of the SMTZ and lower rates of AOM (Fig. 4). Also, methanogenesis rates dropped substantially in the past 10.8 kyr and nearly all CH 4 consumed by AOM today is thousands of years old (Fig. 4). The availability of previously accumulated CH 4 makes AOM less sensitive to decreases in OM deposition than organoclastic sulfate reduction and led to an increasing relative contribution of AOM-derived HS À to the total FeS 2 formation rate from $20% to $40% in the last 11 kyr.
To highlight the impact of transient early diagenetic conditions, the outcome of the ''best fit" transient simulation is compared to simulation results arising from simplified scenarios (Fig. 5). When the OM loading is strictly linearly dependent upon the sedimentation rate (Eq. (9)), i.e., without phases with increased OM deposition relative to SAR, the model only reproduces the largest pyrite peak at $8.3 mbsf (compare ''no imposed OM variations" with ''best fit scenario", Fig. 5B and A, respectively). In this scenario (''no imposed OM variations" in Fig. 5B), the SMTZ is located slightly deeper due to lower OM input during the period 12.54 to 11.14 kyr BP (Phase IV, Table 3); and the Fe 2 O 3 profile does not exhibit the troughs at depths corresponding to the pyrite peaks (3.3 and 17.5 mbsf) formed during Phases II and IV (Table 3), which are also observed in the data (Fig. 2H). Increases in OM loading relative to SAR (as in the best fit scenario) give a straightforward explanation for the formation of these pyrite peaks and Fe 2 O 3 troughs. However, as the deeper buried OM at the Zambezi margin only contains the unreactive fraction of the initially deposited OM (Deville et al., 2020) and we lack other direct constraints on OM deposition, we cannot rule out alternative scenarios.
When ignoring the time dependence of sediment input and imposing the average SAR over the past 27.4 kyr, which is on the order of magnitude of glacial SAR, the simulations reach a steady state (''steady state scenario", Fig. 5D). It yields an SMTZ 25.7 mbsf and low (<0.1 mM) methane concentrations below that. Pyrite formation is almost exclusively associated with organoclastic sulfate reduction, which leads to a gradual increase of pyrite and decrease of Fe 2 O 3 with depth. When the organic matter pulses during Phase II and IV (Table 3) are added (''constant SAR with imposed OM variations", Fig. 5C) two pyrite peaks are formed in the zone of organoclastic sulfate reduction and buried to 14.9 and 26.0 mbsf. Phase IV pushes the SMTZ considerably upwards compared to the steady-state simulation, and both pyrite peaks also contain AOM-derived sulfur. The movement of the SMTZ is much greater than in the best fit scenario, since this scenario does not accurately reflect the relatively low SARs (and hence lower OM loading) evident from the age model during this period. In the simulations with constant SAR (Fig. 5C and D), the OM is highest at the sediment-water interface and decreases with depth (not shown). The opposing trend in the data and best fit simulation (Fig. 2I) is a transient feature caused by lower SAR in the last $7.5 kyr allowing for aerobic respiration to consume a larger fraction of the more refractory organic matter.
The contributions of both organoclastic sulfate reduction and AOM to the formation of pyrite are supported by the measured isotopic composition of the pyrite peaks (Fig. 2G). The pyrite peak at 3.3 mbsf formed -according to the model simulations -in the zone of organoclastic sulfate reduction. It exhibits a more negative d 34 S isotope signature of À40‰ (Fig. 2G), which is consistent with open system conditions and a constant supply of fresh SO 4 2À (Hartmann and Nielsen, 1968;Jørgensen, 1979). The simulated pyrite peaks at 8.3 and 17.5 mbsf contain sulfur derived from both organoclastic sulfate reduction and AOM, which is qualitatively in agreement with the intermediate d 34 S values of À25‰ and À10‰, respectively. The decreased apparent isotope effects suggest pyrite formation under conditions with increasingly limited SO 4 2À supply (Borowski et al., 2013), but the values have a considerable offset from seawater SO 4 2À ($21‰), which is at odds with the total SO 4 2À depletion in the SMTZ and indicates that some HS À produced by organoclastic sulfate reduction (without substrate depletion) during an earlier depositional period contributes to this peak as well.

SYNTHESIS
Along ocean margins, SAR, OM and Fe oxide loading have varied substantially since the LGM with overall decreasing trends as sea-level rise increased the distance of previous depositional centers from the river mouth. On the Mozambique continental slope, abrupt decreases in SAR were observed at $12.7 and 8 kyr BP, which were likely to be caused by subsequent changes in the current system on the shelf as it became gradually flooded (See Section 4.2; Beiersdorf et al., 1980). The re-routing of sediments in response to deglacial sea-level rise led to a transition from coastal conditions to a more open marine setting at our study site, characterized by lower depthintegrated rates of organic mineralization including organ-oclastic sulfate reduction, AOM, and pyrite precipitation (Fig. 4). The simulations show that the methane being consumed today was produced in the past, causing the downward migration of the SMTZ starting 6500 years ago and still ongoing today (Fig. 3D-F; Fig. 4G).
Our work shows that the sediment geochemical profiles are the result of repeated changes in depositional forcings driven by a global transition out of a glacial maximum climate state. The long-term evolution of the diagenetic system at Site MOZ4-CS17 can be described as follows (Fig. 6): Fig. 6. Schematic illustration for diagenetic evolution at the Mozambique margin (not to scale). Panels 1-3 correspond to the three-step conceptual model discussed in the text.
(1) Before 20 kyr BP: High SAR are associated with high OM deposition and high rates of methanogenesis. During this glacial period, high SAR dilute pyrite concentrations. In our simulations, this period is less well constrained, as it is influenced by environmental processes that occurred prior to the onset of our record.
(2) 20-10 kyr BP: The post-glacial and deglacial are characterized by varying SAR and variable but overall low OM loading causing (with considerable lag) the depth of the SMTZ to fluctuate and temporarily stabilize in distinct sediment intervals. Here, focused pyrite peaks form. Pyrite peak formation is also affected by fluctuations in OM deposition leading to stronger organoclastic SO 4 2À reduction and HS À formation.
(3) 10 kyr BP -Present: Continued low SAR, low OM input and decreasing CH 4 production at depth during the Holocene lead to gradual deepening of the SMTZ at approximately the speed of the burial velocity. The CH 4 diffusing into the SMTZ was mostly produced during earlier times. Despite simultaneously decreasing pyrite formation rates, pyrite continues to accumulate in the SMTZ as the lower SAR allows for overall higher pyrite accumulation.
The most intense accumulation of diagenetic minerals at Site MOZ4-CS17 took place at the transition between the coastal/shelf-like high SAR and the deep-sea-like low SAR regimes. During periods of high SAR, the SMTZ was driven to shallower depth due to strong upward CH 4 diffusion, preventing the accumulation of pyrite in a single sediment interval. This can still be seen in the observed sedimentary pyrite profile below 20 mbsf where the pyrite concentrations are relatively uniform. During these time periods, total sulfate consumption was $0.4 mol S m À2 yr À1 , typical for inner continental shelf sites (e.g., Bowles et al. 2014). When SAR dropped, the system was initially charged with pre-formed CH 4 and dissolved Fe 2+ , which caused high overall SO 4 2À consumption rates. The decreasing SAR, however, lowered total sulfate consumption to 0.0375 mol m À2 yr À1 under modern conditions, representative for slope or deep-sea environments, eventually causing a downward movement of the SMTZ in relation to the sediment surface, keeping the zone of pyrite formation within the same sediment interval. This led to stronger pyrite accumulations over the last 10 kyr in the respective sediment interval. Provided that SAR and OM deposition remain low, methanogenesis will further decrease and the SMTZ will continue to move deeper into the sediment, AOM rates will decrease and new pyrite formation will stop. Eventually, the pyrite peaks formed in a paleo-SMTZ will be located above the position of the actual SMTZ (Fig. 3F). The stark difference in sulfate consumption rates stated above illustrate how fundamentally the depositional environment and the resulting diagenetic processes at Site MOZ4-CS17 have changed over the last deglacial period.
Our analysis and quantitative interpretation of the sediment record on the Mozambique margin documents that global climate and associated sea level change is recorded in a massive shift in biogeochemical processes, which in turn affects and in part overprints the primary climatic signal recorded in the sedimentary record. As the sea-level related changes in SAR occurred along continental margins worldwide, this mechanistic model is applicable to other settings as well. Therefore, our findings can help quantifying marine sedimentary C turnover along continental margins over the last glacial-interglacial transition. They can further provide a guideline for the interpretation of pyrite accumulations in marine sediments as indicators for changes in the depositional environment.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.