Benthic Dissolved Silicon and Iron Cycling at Glaciated Patagonian Fjord Heads

Abstract Glacier meltwater supplies silicon (Si) and iron (Fe) sourced from weathered bedrock to downstream ecosystems. However, the extent to which these nutrients reach the ocean is regulated by the nature of the benthic cycling of dissolved Si and Fe within fjord systems, given the rapid deposition of reactive particulate fractions at fjord heads. Here, we examine the benthic cycling of the two nutrients at four Patagonian fjord heads through geochemical analyses of sediment pore waters, including Si and Fe isotopes (δ30Si and δ56Fe), and reaction‐transport modeling for Si. A high diffusive flux of dissolved Fe from the fjord sediments (up to 0.02 mmol m−2 day−1) compared to open ocean sediments (typically <0.001 mmol m−2 day−1) is supported by both reductive and non‐reductive dissolution of glacially‐sourced reactive Fe phases, as reflected by the range of pore water δ56Fe (−2.7 to +0.8‰). In contrast, the diffusive flux of dissolved Si from the fjord sediments (0.02–0.05 mmol m−2 day−1) is relatively low (typical ocean values are >0.1 mmol m−2 day−1). High pore water δ30Si (up to +3.3‰) observed near the Fe(II)‐Fe(III) redox boundary is likely associated with the removal of dissolved Si by Fe(III) mineral phases, which, together with high sedimentation rates, contribute to the low diffusive flux of Si at the sampled sites. Our results suggest that early diagenesis promotes the release of dissolved Fe, yet suppresses the release of dissolved Si at glaciated fjord heads, which has significant implications for understanding the downstream transport of these nutrients along fjord systems.

Early diagenetic processes in sediments affect the solubility of reactive particulate silicon (BSi, ASi) and iron (pFe) phases. In general, BSi and ASi are subject to undersaturation and dissolution in marine pore waters (Aller, 2014). The solubility of Si can be regulated by the incorporation of aluminum (Al) and Fe into the amorphous silica structure (Aller, 2014;Van Cappellen et al., 2002), and reverse weathering-the formation of authigenic aluminosilicates or clay minerals (Ehlert et al., 2016;Michalopoulos & Aller, 2004). In addition, sediment cycling of Fe can influence pore water DSi through adsorption onto/desorption from Fe (oxy) hydroxides and co-precipitation of amorphous Fe-Si minerals (Delstanche et al., 2009;Geilert, Grasse, Doering, et al., 2020;Schulz et al., 2022;Zheng et al., 2016). Redox state is the main control of Fe availability in pore waters (Aller, 2014). Dissolved Fe 2+ is released into pore waters through microbially mediated and abiotic reduction of Fe(III)-rich minerals (Herbert et al., 2021;Laufer-Meiser et al., 2021). Precipitation of the dissolved Fe 2+ can also occur through oxidative formation of Fe (oxy)hydroxides in the sediments (Herbert et al., 2021;Laufer-Meiser et al., 2021), and formation of pyrite in anoxic, sulfidic sediments (Aller, 2014).
The balance between competing diagenetic processes in fjord sediments governs not only elemental solubility but also pore water gradients, which drive the flux of these nutrients back to the overlying water column (benthic flux) via diffusion, pore water advection, bioturbation, and other sediment mixing processes (Schulz & Zabel, 2006). DSi and dFe released from the sediments may then be transported to sustain fjord, coastal and marine primary production and ecosystems via a combination of upwelling, diffusion, bottom-water flushing, and shelf-fjord exchange (Bianchi et al., 2020;Hopwood et al., 2020). Determining the relative magnitudes of these benthic nutrient fluxes is also important to evaluate the modification of downstream nutrient ratios and phytoplankton composition (Hopwood et al., 2020).
Here, we present a range of geochemical analyses, including dissolved δ 30 Si and δ 56 Fe of pore water samples from four sediment cores that were acquired from the head of Baker-Martínez Fjord Complex (BMFC) in Chilean

Study Sites
The BMFC is located between the Northern and Southern Patagonian Icefields (NPI and SPI) and is connected to the Gulf of Penas and the Pacific Ocean to the west (Figure 1). The BMFC generally has a two-layered water mass structure: the upper layer (top 20-30 m depth) consists of estuarine water that is formed from a mix of continental freshwater and seawater, while the lower layer is composed of Subantarctic water that flows in from the Pacific Ocean through the Gulf of Penas (Quiroga et al., 2016). The major freshwater and terrestrial inputs into the fjord complex are the Baker River, the Huemules River, the Pascua River, and Jorge Montt Glacier (Troch et al., 2021). Glacier meltwater is a significant contributor to the runoff of these rivers (Amann et al., 2022;. Colored symbols on the main map show the locations of the sediment core sites at the fjord heads, which are marked by the colored polygons. The bold red box on the inset map shows the extent of the main map. The figure was made using QGIS software (QGIS, 2020) and ESRI satellite imagery (Esri, 2020 The NPI sits atop the Patagonian batholith-a series of intrusive igneous rocks with a range of magmatic compositions from gabbro to granite. These igneous rocks formed from the subduction of the Antarctic plate beneath the South American plate between the Late Jurassic and the Neogene (Hervé et al., 2007). The SPI is predominantly underlain by the Eastern Andes Metamorphic Complex, which is mainly composed of Paleozoic greenschist-facies metasediments (Hervé et al., 2008). Subglacial weathering of these igneous and metamorphic bedrocks produces high concentrations of colloidal and reactive particulate Si and Fe phases , which are delivered into the BMFC.

Field Sampling and Sediment Cores
Bottom water temperature, salinity, and dissolved oxygen concentration measurements at or near the core sites were acquired using an RBR maestro conductivity-temperature-depth (CTD) sensor. A gravity corer fitted with 1.5 m long transparent core liners was used to acquire sediment cores from four sites in the BMFC (Figure 1) on the research vessel Sur-Austral in February 2017. The equipment available on the research vessel was suited for relatively shallow coring operations (<300 m depth), and hence sediment cores were only collected from the fjord heads. The four core sites (Table S1 in Supporting Information S1): SJ (106 m depth), SH (203 m depth), SP (248 m depth) and SB (151 m depth) were situated close to the outflows of the Jorge Montt Glacier, Huemules River, Pascua River, and Baker River, respectively. One set of sediment cores collected from these four sites was used for sampling filtered core-top water, followed by sampling of filtered pore water, using syringes attached with Rhizon pore water samplers (0.15 μm; Rhizosphere) that were inserted through holes pre-drilled on the core liners (Ng et al., 2020). The collected water samples were acidified by adding 1% v/v 10 M lab distilled HCl prior to storage.
Another set of cores was obtained from the four fjord sites for sediment analysis (Table S1 in Supporting Information S1). Sediment physical properties such as porosity and bulk density were measured on this set of cores using Gamma ray attenuation on a Geotek Multi-Sensor Core Logger, prior to subsampling of solid sediments for subsequent analyses. Sedimentation rate(s) of core SJ was based on annual variations in ice-rafted debris deposition (De Wilde, 2016), cores SH and SP were based on 137 Cs (Piret et al., 2021(Piret et al., , 2022, and core SB was based on 210 Pb (Vandekerkhove et al., 2021).

Dissolved Nutrient and Element Concentration Measurements
Pore water and core-top water dissolved Si and nitrate concentrations, [DSi] and [NO 3 − ] respectively, were analyzed in the LOWTEX facility at the University of Bristol. Sample [DSi] was measured with the molybdate-blue method (DeMaster, 1981) using the LaChat 8500 Series-2 flow injection analyzer. [NO 3 − ] was measured using a Thermo Gallery discrete analyzer after samples were reduced with vanadium chloride and reacted with the Griess reagent (Wang et al., 2016). The analytical errors of the methods were ∼2.5% for [DSi] and ∼0.5% for [NO 3 − ], derived from 2 SD (standard deviation) of replicate measurements of samples and TraceCERT ® standards.
Concentrations of dissolved (<0.15 μm) iron [dFe] and manganese [dMn] were measured using a Thermo Scientific X-Series 2 quadrupole inductively coupled plasma-mass spectrometer (ICP-MS) with a collision/reaction cell at the National Oceanography Centre, Southampton, after a ×30 gravimetric dilution in 2% v/v HNO 3 . Beryllium, indium and rhenium were used as internal standards to correct for drift and matrix effects. Multi-element external calibration solutions were gravimetrically made to match the concentration range observed in samples. All measurements were above the instrument's limit of detection (∼3 nM for Fe, ∼0.1 nM for Mn). Measurement accuracy was determined by analyzing the gravimetric multi-element external calibration solutions, and precision was determined by repeat measurements on the gravimetric calibration solutions of intermediate concentrations.
Analytical precision for the elements analyzed ranged from ±0.7 to ±8.9%, and accuracy was better than ±10%.

Fe Isotope Analysis
Aliquots of the pore water and core-top water samples were used for Fe isotope analysis at the University of South Florida (USF). All work was carried out under ISO-5 clean conditions, all water was 18.2 MΩ. cm from a Genpure Barnstead System and all acids were single-distilled from trace metal grade acids via Savillex™ DST-1000 stills. A 57 Fe-58 Fe double spike was added to samples in a sample-to-spike ratio of 1:2. Samples were Global Biogeochemical Cycles NG ET AL.
10.1029/2022GB007493 5 of 22 dried down, redissolved and refluxed in concentrated HNO 3 with 10% H 2 O 2 (Fisher Ultrapure grade) to digest any organic material. Samples were then processed using an AGMP-1 resin microcolumn purification scheme following Conway et al. (2013), as modified by Sieber et al. (2021). Following purification, samples were dried down and then reconstituted in 2% v/v HNO 3 for instrumental analysis. The method has a total procedural blank of <0.4 ng (Sieber et al., 2021).
Iron isotope ratios were measured using a Thermo Neptune MC-ICP-MS in "high" resolution mode at the USF Tampa Bay Plasma Facility, using a ∼100 μL min −1 PFA nebulizer and an ESI Apex Ω introduction system. Analysis proceeded as described in Sieber et al. (2021); beams were corrected for background using an on-peak measurement of 2% v/v HNO 3 immediately prior to each sample, the abundances of 53 Cr and 60 Ni were measured to correct 54 Cr and 58 Ni isobaric interferences on 54 Fe and 58 Fe, and correction for instrumental mass bias was carried out using the double spike method (Siebert et al., 2001). Sample δ 56 Fe are expressed relative to the mean of two IRMM014-double spike bracketing standards (per block of 6 samples). At USF the long-term 2 SD of analysis (±0.04‰) of NIST-3126a is used as a conservative estimate of analytical precision (Sieber et al., 2021). We therefore express uncertainty on δ 56 Fe measurements as 0.04‰, or the 2 SE (standard error) of an individual measurement when >0.04‰.

Analysis of Sediment Reactive Si Phases
Operationally-defined reactive Si phases in the sediments analyzed here are (a) mild acid-leachable silica (Si-HCl) associated with Fe, Mn and Al oxides (Michalopoulos & Aller, 2004;Pickering et al., 2020), and (b) mild alkaline-leachable amorphous silica (Si-Alk) that includes BSi and glacier-derived ASi (Ng et al., 2020;. Analyses of these phases were carried out using sequential extraction, with procedures modified from Michalopoulos and Aller (2004), Pickering et al. (2020), and Ng et al. (2020).
Si-HCl was extracted from ∼50 to 100 mg of air-dried sediments by treating the samples with 36 mL of 0.1 M HCl (lab distilled) and shaking for 18 hr. Samples were centrifuged, and the supernatant was filtered through 0.2 μm Sartorius™ polyethersulfone (PES) syringe filters. The pH of the Si-HCl leachates was adjusted to ∼5 using 4 M NaOH (Honeywell Fluka™; brown Fe precipitates would form at pH 7) before the leachates were analyzed for [DSi] and the Si isotopic composition. The remaining sediments were rinsed with 18.2 MΩ.cm Milli-Q water, treated with 5 mL of 10% H 2 O 2 (Romil SpA™) for 0.5 hr to remove organic matter, and rinsed again with Milli-Q water.
Si-Alk was extracted from the HCl-and H 2 O 2 -treated sediments by digesting the sediment in 40 mL of 0.1 M Na 2 CO 3 (ACROS Organics™) in a water bath at 85°C (DeMaster, 1981). A portion of the leachate was taken 20 min into the digestion for Si isotope analysis before a significant amount of lithogenic material started to dissolve. Aliquots were also taken at 2, 3, and 5 hr to quantitatively distinguish the rapidly dissolving BSi and ASi (Si-Alk) from the slowly released lithogenic silica. Leachates were filtered through 0.2 μm syringe filters (as above) and neutralized with HCl prior to instrumental analysis. A linear regression was fitted through the measured leachate Si content at 2, 3, 5 hr of digestion, and the intercept at 0 hr was used to quantify the Si-Alk concentration.
Si content in the acid and alkaline leachates were measured with the molybdate-blue method (DeMaster, 1981) using a V-1200 Vis spectrophotometer (Ng et al., 2020). Procedural blanks were below the detection limit (<0.01 μM). Uncertainties in the sequential extraction procedures were calculated using average standard deviations of repeat sediment leaches and measurements of Si content in the Si-HCl and Si-Alk phases, as ±15% and ±25%, respectively.

Si Isotope Analysis
For the pore water and core-top water samples, pre-concentration of Si was carried out using a Mg-induced co-precipitation method (de Souza et al., 2012). The method involved precipitation of Mg(OH) 2 with 1 M NaOH (Titripur ® Reag. Ph Eur grade), followed by rinses of 0.001 M NaOH to quantitatively remove excess major ions: Na + , Cl − , SO 4 2− , Ca 2+ , and K + , after which the Si-containing precipitates were dissolved with HCl (lab distilled) to a smaller volume. Using the spectrophotometer (Section 2.5), aliquots of the concentrated samples were measured for [DSi] to monitor the yields of the co-precipitation method, which were >95%.

Global Biogeochemical Cycles
NG ET AL.
10.1029/2022GB007493 6 of 22 The pre-concentrated pore water, core-top water samples and the sediment leachate samples were further purified with cation exchange column chemistry following the procedures of Georg et al. (2006). The column chemistry consisted of polypropylene columns containing 1.8 mL Bio-Rad AG50W-X12, 200-400 mesh cation exchange resin in H + form, which was pre-rinsed with HCl and Milli-Q before the samples were loaded onto the columns and eluted with ultra-pure water.

Reaction-Transport Modeling for Si
The concentration and isotope profiles of dissolved and solid Si phases in the fjord sediments were simulated using the Biogeochemical Reaction Network Simulator (BRNS), an adaptive simulation environment for mixed kinetic-equilibrium reaction networks (Aguilera et al., 2005;Regnier et al., 2002). Our BRNS design followed the model setup detailed in Cassarino et al. (2020) and Ward et al. (2022a). We assumed that Si diagenetic processes were at steady state. Although fjord systems are subject to transient episodes of glacial runoff, phytoplankton blooms and sediment deposition (Bianchi et al., 2020), the assumption of steady state allows us to define the long-term dynamics within fjord sediments and offers a first-order estimate of the diagenetic processes governing the benthic Si cycle at fjord heads.
The model is based on a vertically-resolved reaction-transport mass conservation equation that computes concentration changes of dissolved and solid phases in porous media (Berner, 1980;Boudreau, 1997): t, z, C Si are time, depth, and concentration of Si respectively. ω denotes the measured sedimentation rate (Table S2 in Supporting Information S1). For the dissolved phase: σ = porosity or ϕ, while for the solid phase: σ = (1ϕ).
The model assumed an exponential decrease in porosity due to sediment compaction (e.g., Freitas et al., 2020), which was calculated based on the measured core-top porosity, ϕ 0 , measured bottom porosity, ϕ x , and a fitted porosity attenuation coefficient, β ϕ (Table S2 in Supporting Information S1). Sediment bioturbation was modeled as a diffusive-like process, constrained with the bioturbation coefficient, D bio and bioturbation depth, z bio (Table  S2 in Supporting Information S1). Values of D bio and z bio were inferred from a benthic macrofauna study in fjords (Oleszczuk et al., 2019). The coefficient of molecular diffusion of DSi, D m_Si (D m_Si = 0 for solids) was derived based on an empirical relationship (Rebreanu et al., 2008) with measured temperature, T and measured salinity, S (Table S2 in Supporting Information S1). Bioirrigation rate, α irr (α irr = 0 for solids) was computed as a function of the bioirrigation coefficient, α 0 and bioirrigation attenuation depth, x irr (Table S2 in Supporting Information S1). α 0 was inferred from the benthic macrofauna study in fjords (Oleszczuk et al., 2019), while a global mean value Ward et al., 2022a) was assumed for x irr (Table S2 in Supporting Information S1).
The final term of Equation 1 represents the sum of reactions j affecting the Si phases, with reaction stoichiometric coefficient λ Si,j and rate R j . We implemented a reaction network that resolved the dissolution of an undersaturated Si-Alk phase (BSi, glacial ASi) and the precipitation of authigenic aluminosilicate clay minerals (AuSi)

Global Biogeochemical Cycles
NG ET AL.
10.1029/2022GB007493 7 of 22 associated with reverse weathering. The kinetic rate laws associated with the dissolution and precipitation reactions, R d and R p respectively, are expressed as: k diss and k precip are fitted rate constants. ad, bd, ap, and bp are fitted coefficients defining the exponential decreases of the dissolution and precipitation reactions Ehlert et al., 2016;Ward et al., 2022a).
[ASi] sat and [AuSi] sat are the saturation concentrations of the amorphous silica and the authigenic aluminosilicates, respectively, based on values derived from previous studies ( Table S2 in Supporting Information S1). We assume an exponential decrease of R d (Equation 2) to account for the rapid decline of amorphous silica solubility due to the incorporation of Al and Fe in the sediments (Aller, 2014;Van Cappellen et al., 2002). Similarly, we assume that R p also decreases exponentially (Equation 3), following observations that authigenic clay precipitation occurs predominantly in the upper sediments (Michalopoulos & Aller, 2004), potentially due to the greater availability of reactive Al and Fe phases.
The core-top concentrations and isotope compositions of DSi, Si-Alk, and Si-HCl measured in Section 2.5 and 2.6 were utilized to define the upper boundary conditions of the BRNS-Si model. The δ 30 Si measurements from Section 2.6 were converted to explicit 28 Si and 30 Si isotope concentrations for the model simulation through a series of isotope mass balance equations (Equation S1-S4 in Supporting Information S1; Cassarino et al., 2020;Ward et al., 2022a). The model assumes that the total abundance of Si is represented by the sum of 28 Si and 30 Si ( 29 Si not included), and so there is an error of under 5% associated with the simulated concentration of a given Si phase. Potential isotopic fractionation associated with the dissolution and precipitation reactions are accounted for by the fractionation factors α d and α p respectively (Equation S5-S8 in Supporting Information S1), based on values derived from previous studies (Table S2 in Supporting Information S1).
Our understanding of coupled Si-Fe cycling in marine sediments is relatively new. Detailed mathematical formulation of the associated reactions requires future experimental studies and more observational data, especially on the fraction of Si adsorbed onto or co-precipitated with Fe(III) phases in seawater and pore water composition.
Here, simplified constant reaction rates (Ward et al., 2022a) were used to simulate Si-Fe reactions. R up28 and R up30 are model-optimized constant rates of uptake of pore water 28 Si and 30 Si by reactive Fe phases through adsorption and co-precipitation in the relatively oxic upper sediment layer. In contrast, R re28 and R re30 are model-optimized constant rates of release of 28 Si and 30 Si into pore water from reductive dissolution of Fe in the ferruginous zone (Table S2 in Supporting Information S1). The upper and lower depth limits of the oxidizing layer, z ox-up and z ox-low , and those of the ferruginous zone, z fe-up and z fe-low (Table S2 in Supporting Information S1), were inferred from measured pore water dFe concentration profiles. The sediment Si-HCl phase is thought to be partly composed of Si associated with metal oxides including Fe oxides . Therefore, in the model we allowed sediment Si-HCl content and its isotope composition to vary following the modeled R up28 , R up30 , R re28 , and R re30 (Equation S9-S12 in Supporting Information S1).
The length of the model domain, L M , was set based on the maximum depth where pore water DSi data were available (Table S2 in Supporting Information S1). The model run time was set to allow for two full sediment deposition cycles (runtime = 2 × L M ÷ ω), thus ensuring that steady state was achieved. The model was run on an uneven depth resolution grid, which allowed for a fine resolution of the upper, most active sediment layers (e.g., Freitas et al., 2020). The depth intervals were resolved at 0.1 cm interval from 0 to 3 cm depth, at 0.2 cm interval from 3 to 5 cm depth, at 0.4 cm interval from 5 to 10 cm depth, at 0.6 cm interval from 10 to 20 cm depth, and at 1 cm interval for sediment depths greater than 20 cm.
We ran model ensembles until a best-fit solution was found. The best-fit solution was determined by minimizing the misfit (specifically, the root-mean-square error or RMSE: Equation S13 and S14 in Supporting Information S1) between the measured and modeled depth profiles of Si concentration and isotope composition (Table  S3 in Supporting Information S1). Our model's outputs include rates of reaction R d and R p , as well as the concentrations and isotope compositions of pore water DSi, sediment Si-Alk and Si-HCl. As such, the model allows us 10.1029/2022GB007493 8 of 22 to quantify parameters defining the rates of Si-Alk dissolution (k diss , ad and bd), AuSi precipitation (k precip , ap and bp), and coupled Si-Fe cycling (R up28 , R up30 , R re28 , and R re30 ).

Redox Conditions and Pore Water dFe in BMFC
Bottom water dissolved oxygen concentrations at the fjord heads ( (Figures 2 and 3a). In particular, the ferruginous zones at all the study sites ( Figure 2) are deeper, broader, and with higher maxima in pore water [dFe] (550-770 mmol m −3 ) than other coastal marine and open ocean sediments observed, but similar to selected sites from glaciated fjords in Svalbard (Herbert et al., 2020(Herbert et al., , 2021(Herbert et al., , 2022Laufer-Meiser et al., 2021;Wehrmann et al., 2014Wehrmann et al., , 2017. The broader ferruginous zones and [dFe] maxima at our BMFC sites are most likely due to a high supply of reactive Fe phases from glacier-fed rivers , land-terminating and marine-terminating glaciers (Henkel et al., 2018;Wehrmann et al., 2014).
In all BMFC cores, pore water δ 56 Fe at depths below 0.5 m cluster around −1 to −0.2‰ (Figure 3b), which are lower than the average crustal δ 56 Fe value (∼+0.1‰; Beard et al., 2003). These lower δ 56 Fe values most  (Beard et al., 2003). Lines in the [DSi] plots are least squares exponential fitting of the data carried out using SigmaPlot (Systat Software Inc.). Error bars represent 95% confidence intervals of δ 56 Fe and δ 30 Si measurements, and some are smaller than the symbols. NG ET AL.
10.1029/2022GB007493 10 of 22 likely reflect an imprint of microbially mediated dissimilatory Fe reduction within the ferruginous zone (Crosby et al., 2007;Homoky et al., 2009;Severmann et al., 2006). In contrast, there is a notable difference in the shallower pore water δ 56 Fe profiles between the BMFC sites (Figure 3b), which is likely the result of multiple processes that involve lithogenic colloids/nanoparticles, organic ligands, and redox-driven dissolution and re-precipitation of Fe. Fe isotope signatures close to or heavier than crustal compositions (+0.3 to +0.8‰) are observed in both the shallowest, oxic pore waters (<0.1 m depth) and the overlying core-top waters of SH and SP sites. The heavier Fe isotope signatures potentially reflect isotopic fractionation by "non-reductive" formation of lithogenic Fe(III) colloids/nanoparticles (<0.15 μm; Radic et al., 2011;Homoky et al., 2021), or by the binding of Fe(III) to organic ligands (Dideriksen et al., 2008;Morgan et al., 2010), that mainly occur in oxygenated waters.

Spatial Differences in Fe Redox Cycling
Sites SJ and SH show pore water δ 56 Fe minima (−2-−3‰) at the Fe(II)-Fe(III) redox boundary (just above the ferruginous zone, Figure 2) around 0.1-0.3 m depth (Figure 3b). These isotopically-light Fe compositions can be explained by intense Fe redox cycling that involves: [i] reductive dissolution of Fe yielding isotopically light Fe 2+ (Crosby et al., 2007) in the ferruginous zone, [ii] upward advection and diffusion of Fe 2+ to the relatively oxic sediment layer, [iii] oxidation of Fe 2+ to amorphous Fe(III) (oxy)hydroxides (Homoky et al., 2009) driving the remaining Fe 2+ even lighter (Severmann et al., 2006). In contrast, the lack of an evident pore water δ 56 Fe minimum around the Fe(II)-Fe(III) redox boundary (∼0.1-0.3 m depth) in the other two cores (Figure 3b) suggests that the Fe redox cycling described above occurs to a far lesser extent at the SP and SB sites.
Here, we examine the potential factors that could be responsible for the observed spatial differences in the intensity of Fe redox cycling. The intensity of Fe(III) reduction (step [i]) is partly linked to the supply of easily reducible Fe (oxy)hydroxides, such as ferrihydrites, which differ between the fjord sites. The Huemules river upstream of the SH site ( Figure 1) has been found to yield 10 times higher amounts of highly reactive Fe than the Baker and Pascua rivers upstream of the SB and SP sites respectively (Table S4 in Supporting Information S1; . The higher Fe input to SH has been linked to the greater extent of glacial cover in the Huemules catchment ). An even greater yield of highly reactive Fe is expected at the SJ site given the supply of glacial sediments through ice rafting and subglacial discharge without proglacial lake processing . The greater rate of supply of highly reactive Fe available at SJ and SH is thought to sustain increased microbial Fe reduction, which could contribute to the lighter pore water δ 56 Fe signals observed at these two sites (Figure 3b).
A dFe uptake mechanism in the relatively oxic upper sediment layer involves an upward transport of Fe 2+ from the underlying ferruginous zone and subsequent Fe 2+ oxidation (steps [ii] and [iii] of Fe redox cycling). The intensity of the upward transport (step [ii]) is a function of pore water advection, which can differ substantially between the study sites depending on sediment physical properties and depositional environment. For example, given its proximity to the marine-terminating Jorge Montt Glacier (Figure 1), the SJ site has the highest sedimentation rate and the highest sediment porosity (Table S1 in Supporting Information S1), which would maintain a significantly larger compaction-driven, upward pore water advection (Santos et al., 2012;Schulz & Zabel, 2006) compared to the other study sites. Bottom currents flowing over the uneven surface sediments is another common mechanism maintaining pore water advection through the top several tens of centimeters of permeable sediments (Santos et al., 2012;Schulz & Zabel, 2006). The bottom current-induced pore water advection most likely occurs at the SJ and SH sites where the upper sediments have permeable sandy layers (Table S1 in Supporting Information S1), likely sourced from direct glacial sediment inputs from Jorge Montt Glacier, and from periodic glacial lake outburst floods (GLOF) from Steffen Glacier (Piret et al., 2021), respectively. Meanwhile, bottom current-induced pore water advection is most likely impeded at the SP and SB sites, where the upper sediments (<0.5 m, Table S1 in Supporting Information S1) are less permeable muds and silts (Santos et al., 2012;Schulz & Zabel, 2006) mainly deposited from riverine outflows (Piret et al., 2022;Vandekerkhove et al., 2021).
The spatial difference in sediment physical properties between the fjord sites likely indicates decreasing intensity of pore water advection from SJ, SH, SP to SB. Therefore, upward flow of Fe 2+ and subsequent Fe 2+ oxidation (step [ii] and [iii] of Fe redox cycling) are likely most intense at SJ, followed by SH, SP, and SB. The inferred trend is consistent with the pore water δ 56 Fe observation at around the Fe(II)-Fe(III) redox boundary, where the lowest δ 56 Fe minimum is observed at SJ, followed by SH, and the lack of such δ 56 Fe minima at SP and SB (Figure 3b). In short, the intensity of Fe redox cycling at the fjord sites appears to be related to sediment physical properties, depositional environment, and the rate of supply of highly reactive Fe such as amorphous Fe (oxy) hydroxides.

Processes Influencing Sediment Pore Water Si in BMFC
The core-top water [DSi] ranges between 21 and 26 mmol m −3 . Pore water [DSi] at all sites exhibit asymptotic exponential profiles (Figure 3c), where [DSi] increases with increasing core depths at ever-decreasing rates toward values ranging between ∼200 and 260 mmol m −3 . The increase in pore water [DSi] with depth indicates a net supply of Si from the dissolution of undersaturated sediment silica phases, which are dominantly represented by Si-Alk (BSi, glacial ASi) in the sediments (Ng et al., 2020;. For sites SJ, SH and SP, the δ 30 Si of the Si-Alk phases fall between 0 and −0.5 ‰ (Figure 4, Figure S1 in Supporting Information S1), which is within the range of previously published glacially-derived ASi isotope endmembers Hawkings et al., 2018). Lower Si-Alk δ 30 Si is observed at the SB site (Figure 4, Figure S1 in Supporting Information S1), which suggests a greater relative abundance of other Si-Alk phases that could have a lighter isotope composition than glacial ASi.
If Si-Alk dissolution is the sole process governing benthic Si cycling, pore water [DSi] and δ 30 Si should lie on the mixing line between the core-top water endmember: [DSi] ctw and δ 30 Si ctw , and the Si-Alk endmember. Theoretical mixing of DSi concentration, [DSi] mix , and its isotope composition, δ 30 Si mix , between the two endmembers can be calculated with the following model (Geilert, Grasse, Doering, et al., 2020): [ASi] sat denotes the saturation concentration of amorphous silica (Table S2 in Supporting Information S1). δ 30 Si-Alk represents the average δ 30 Si of the Si-Alk phase in a given core. f is mixing fraction, with f = 0 signifying 100% fluids derived from Si-Alk dissolution, and f = 1 indicating 100% core-top water.
All pore water data lie well above the linear mixing models computed for each of the sediment cores (Figure 4), providing clear evidence that benthic Si cycling at the fjord sites is regulated by other processes, in addition to Si-Alk dissolution. Two major diagenetic processes that can account for the observation above are: (a) reverse weathering-precipitation of AuSi (Ehlert et al., 2016;Ng et al., 2020) including Fe-rich aluminosilicates (Geilert, Grasse, Doering, et al., 2020;Isson & Planavsky, 2018); (b) interaction with Fe, including adsorption onto Fe (oxy)hydroxides and co-precipitation of amorphous Fe-Si minerals (Delstanche et al., 2009;Geilert, Grasse, Doering, et al., 2020;Zheng et al., 2016). These diagenetic processes preferentially remove isotopically light Si isotopes from solution and lead to high pore water δ 30 Si values, which are far above the mixing line between the core-top water and Si-Alk endmembers (Figure 4). AuSi precipitation and coupled Si-Fe reactions are further examined in Section 3.4 and 3.5.
A significant drop of pore water δ 30 Si ( Figure 3d) and an evident transition of sediment Si-Alk content ( Figure S1 in Supporting Information S1) are observed at 0.55 m in core SP. These changes coincide with a sharp decrease in sediment grain size above 0.55 m, which is attributed to a hydrological reorganization of the Pascua River following a GLOF event that occurred in the 1940s (Piret et al., 2022). Our data suggest that the GLOF and resultant turbidity current may have transported a significant amount of isotopically light Si phases, which explains the anomalous pore water δ 30 Si (Figure 3d) at 0.55 m in core SP. The particularly low Si-Alk values ( Figure S1 in Supporting Information S1) in the coarser sediments below 0.55 m (Piret et al., 2022) may suggest a diminished presence of glacial-derived ASi and BSi, supporting the hypothesized change in depositional environment post-GLOF (Piret et al., 2022).

Redox-Driven Coupling of Si and Fe in Pore Water
Pore water δ 30 Si of the fjord sites ( Figure 3d) display strikingly opposing trends to δ 56 Fe (Figure 3b) (Figure 3).
The opposing trends of δ 30 Si and δ 56 Fe between the fjord sites suggests a tight coupling of benthic Si and Fe cycling. This benthic Si-Fe coupling likely involves DSi adsorption onto freshly precipitated amorphous Fe(III) (oxy) hydroxides (Delstanche et al., 2009) and co-precipitation of an amorphous Fe(III)-Si phases (Zheng et al., 2016) that preferentially enriches pore waters with heavier Si isotopes (Geilert, Grasse, Doering, et al., 2020).
Intense pore water advection sustains higher rates of upward Fe 2+ transport (Section 3.2) and the formation of amorphous Fe(III) minerals in the relatively oxic upper sediment layer. This, in turn, drives additional DSi adsorption and/or co-precipitation that progressively increases the δ 30 Si of pore waters. The trend of decreasing pore water δ 56 Fe ( Figure 3b) and increase of pore water δ 30 Si ( Figure 3d) from SB, SP, SH to SJ therefore likely reflects increasing rates of amorphous Fe(III) precipitation driven by pore water advection, and [DSi] adsorption and Si-Fe co-precipitation across the fjord sites.

Uncovering the Benthic Coupling of Si and Fe Through Early Diagenetic Modeling
The BRNS-Si model is used to deconvolve the competing processes in the early diagenetic cycling of Si and to examine the potential significance of benthic Si-Fe coupling in regulating pore water [DSi] at the fjord sites. Model simulations are performed for the two sites with the highest contrasts in depositional setting and pore water

Global Biogeochemical Cycles
NG ET AL.
The model's best-fit of the pore water [DSi] and δ 30 Si profiles (lowest RMSE, Table S3 in Supporting Information S1) can only be achieved with a combination of Si-Alk dissolution, AuSi precipitation, and redox-sensitive coupled Si-Fe reactions (Figures 5 and 6), highlighting the importance of these three diagenetic processes in regulating benthic Si cycling at fjord heads. The modeled dissolution rates of Si-Alk at SJ are more than twice as high as SB (Figures 5g and 6g), implying that ASi delivered from Jorge Montt (the marine-terminating glacier) to SJ may be more soluble compared to ASi that have traveled a long distance along the Baker River to reach SB. Model simulations show relatively invariant downcore profiles of sediment Si-Alk content (Figures 5c and 6c), indicating that only a small fraction of the Si-Alk pool has dissolved. This finding is consistent with the rapid decline of amorphous silica solubility caused by the incorporation of Al and Fe in the sediments (Aller, 2014;Van Cappellen et al., 2002). Observational data show some downcore variation of sediment Si-Alk content at SB (Figure 6c) and δ 30 Si-Alk at SJ (Figure 5d) that are not accounted by the model, which might be due to transient changes in the depositional flux of Si-Alk and its composition, or other poorly understood diagenetic activities.
AuSi precipitation mainly occurs after pore water [DSi] reaches the modeled AuSi saturations, which is below 0.25 m sediment depth at the SJ and SB sites (Figures 5 and 6). The increase of pore water δ 30 Si above 0.25 m depth may therefore be attributed mainly to coupled Si-Fe reactions under oxidizing conditions. Following the argument above, isotopic fractionations (Δ 30 Si solid-aqueous ) of up to −3.8‰ and −2.2‰ are potentially required for the coupled Si-Fe reactions at SJ and SB respectively to balance the isotopically light DSi fluid sourced from Si-Alk dissolution. Such large Si isotopic fractionation may be caused by a combination of DSi adsorption onto Fe(III) (oxy)hydroxides (up to −0.8‰; Delstanche et al., 2009), co-precipitation of amorphous Fe(III)-Si phases (up to −3.2‰; Zheng et al., 2016), and successive cycles of DSi-Fe(III) precipitation and dissolution near a seasonally migrating Fe(II)-Fe(III) redox boundary (Bianchi et al., 2020;Herbert et al., 2022).
The modeled R up28 and R up30 at the SJ site are higher than those at the SB site (Table S2 in Supporting Information S1), consistent with greater DSi adsorption and co-precipitation with Fe at SJ, as inferred from pore water δ 30 Si observations (Figure 3d; Section 3.4). The modeled data also indicates relatively invariant downcore profiles of sediment Si-HCl content and δ 30 Si-HCl, in agreement with observations (Figures 5e, 5f, 6e, and 6f).
These modeling results indicate that the addition of Fe-bound Si to the Si-HCl pool in the relatively oxic sediment layer and corresponding solubilization in the ferruginous zone result in only subtle variations of the Si-HCl pool, which could easily be masked by other less reactive Si species that make up the bulk of Si-HCl phase. These other Si-HCl species may be Si associated with other less soluble metal oxides  or other unidentified secondary weathering products. In addition, the δ 30 Si-HCl values observed in BMFC sediments show relatively large variation between the study sites (from −0.8‰ to −2.1‰; Figure S1 in Supporting Information S1), and are markedly different from δ 30 Si-HCl observed in the Gulf of Mexico river plume sediments and 10.1029/2022GB007493 14 of 22 the Barents Sea shelf sediments (around −2.9‰; Pickering et al., 2020;Ward et al., 2022b). These observations further highlight the heterogeneous nature of the Si-HCl phase in the BMFC sediments, which is likely not dominantly composed of Fe (oxy)hydroxide-bound Si.  sediment layer, which is incorporated in the equations below (Raiswell & Anderson, 2005;Wehrmann et al., 2014).
Given the asymptotic exponential pore water [DSi] profiles observed at the fjord sites (Figure 3c), diffusive benthic fluxes of DSi, J_Si diffusion (mmol m −2 day −1 ), are calculated by Fick's First Law of Diffusion using the following equations (McManus et al., 1995;Ng et al., 2020): D m_Si (m 2 day −1 ) is DSi molecular diffusion coefficient in free solution. β DSi (m −1 ) is the parameter that describes the curvature of the exponential profiles (Table S6 in Supporting Information S1).
[DSi] 0 , [DSi] d , and [DSi] z (mmol m −3 ) represent the DSi concentration in core-top water, asymptotic concentration in pore water, and DSi concentration in pore water at a given sediment depth, z (m), respectively.
J_Fe diffusion and J_Si diffusion derived for the BMFC sites range from 0.0002 to 0.0243 mmol m −2 day −1 and 0.015-0.047 mmol m −2 day −1 , respectively (Table 1). Our four Patagonian sites provide the first estimation of the diffusive benthic DSi flux in a fjord setting, as well as the first diffusive benthic dFe flux outside Arctic (Svalbard) fjords. The Patagonian fjords' diffusive benthic dFe fluxes are lower than those from Svalbard fjords (0.01-0.13 mmol m −2 day −1 ; Herbert et al., 2021Herbert et al., , 2022Wehrmann et al., 2014). However, diffusive benthic fluxes of dFe in these fjords are larger than the open ocean average and of comparable magnitude to the coastal margin average (Figure 7).
In contrast, diffusive benthic DSi fluxes at the fjord sites are small relative to other marine settings (Figure 7). Given the knowledge of sediment Si-HCl and Si-Alk content (mmol kg −1 ), sediment bulk density, ρ b (kg m −3 ), and sedimentation rates, S (m day −1 ) ( Table S7 in Supporting Information S1), the sediment burial rates of the reactive silica phases, J_RSi burial at the fjord sites are estimated to range from 14 to 264 mmol m −2 day −1 (Equation 10).
These values suggest that diffusive benthic DSi fluxes are only returning 0.006%-0.25% of the sediment's reactive Si back to the water column, which is very small compared to the global marine average J_Si diffusion : J_RSi burial (>500%) (Tréguer et al., 2021). The small J_Si diffusion : J_RSi burial (0.006%-0.25%) of the study area is consistent with the high sedimentation rates at the fjord sites (>0.7 cm yr −1 , Table S1 in Supporting Information S1) compared to the ocean (Bianchi et al., 2020), which supports reactive Si burial in these fjord sediments.

The Role of the Benthic Flux in Fjord Nutrient Cycling
Our findings support active benthic cycling that sustains a significant flux of dFe, but a smaller than expected flux of DSi from the sediments at the fjord heads to the overlying fjord waters. These results have important implications for the ongoing debate on the glacier-to-ocean transport of the two nutrients, and for the modulation of nutrient ratios and phytoplankton composition in fjords and the adjacent coastal ocean (Hopwood et al., 2020). The elevated sediment dFe efflux relative to the open ocean average (Figure 7), despite high sedimentation rates, highlights the importance of benthic cycling in processing glacially-sourced solid Fe phases to dissolved forms Figure 7. Global data compilation of diffusive benthic Fe and Si fluxes derived from pore water profiles. The data set is presented using logarithmic vertical scales in a combination of box, whisker and violin plots, showing respectively the median and quartiles (Q 1 or 25th percentile, Q 3 or 75th percentile), the data spread (1.5 × [Q 3 −Q 1 ]), and the data density smoothed by a kernel density estimator. Colored annotations indicate the number of data points available for the three marine settings: fjord, coastal margin, and open ocean. Data are plotted on log scales to handle skewness and to visualize the small values. The plots are created with the R program (R, 2022) using the ggplot2 data visualization package (Wickham, 2016). Thirteen J_Fe diffusion data points (eight from the coastal margin and five from the open ocean) and three J_Si diffusion data points (from the open ocean) have zero values, which are not included in the log scale plots. Our four Patagonian sites represent the first estimations of J_Si diffusion in a fjord setting, while fjordic J_Fe diffusion compiled here also includes estimations from Svalbard fjord sites (Herbert et al., 2021(Herbert et al., , 2022Wehrmann et al., 2014). No recalculation has been performed on the published benthic fluxes. Previous pore water Fe and Si studies that have not evaluated benthic fluxes are not included in the compilation. Full references for the compiled data are listed in Table S8 and S9 in Supporting Information S1.

10.1029/2022GB007493
18 of 22 through intense redox cycling. These findings support the notion that fjord sediments recycle significant fluxes of potentially-bioavailable Fe to the overlying water column (Herbert et al., 2021;Laufer-Meiser et al., 2021;Wehrmann et al., 2014), which could be supplied to downstream ecosystems via a combination of upwelling, estuarine circulation and fjord-shelf exchange (Bianchi et al., 2020;Hopwood et al., 2020).
A low sediment DSi efflux at the fjord sites relative to other marine settings (Figure 7) is likely the combined result of high sedimentation rates at the fjord heads, low pore water [DSi] asymptotes, and the coupling of benthic Si-Fe cycling. Indeed, when Si adsorption and co-precipitation with Fe (and co-dissolution) are "turned off" in the BRNS-Si model (i.e., R up28 , R up30 , R re28 , R re30 are zero), the diffusive benthic DSi fluxes calculated from the new pore water [DSi] profiles ( Figure S2 in Supporting Information S1) are 0.043 mmol m −2 day −1 for SJ and 0.080 mmol m −2 day −1 for SB, which are 2-3 times larger than the original observations ( Table 1). The experiment above supports the role of coupled Si-Fe reactions in suppressing the diffusive benthic DSi flux. This coupling is thought to be particularly pronounced at fjord heads, where there are heightened supplies of reactive Fe from the melting glaciers (Henkel et al., 2018;Wehrmann et al., 2014) and oxidative precipitation of amorphous Fe(III) minerals in the upper sediment layer. Whilst there is a ready supply of dissolvable ASi from the melting glaciers , efficient burial of this phase deposited beneath the fjord heads and the low DSi return flux out of these sediments may limit the amount of bioavailable Si that gets transported to downstream ecosystems. While this study focuses on sites at the Patagonian fjord heads, benthic Si and Fe cycling at sites further away from glacier calving front and riverine outflows may exhibit different trends, given the lower sedimentation rates, changes in primary productivity, organic matter supply and sedimentary Fe(III) reduction rates down the fjords (Herbert et al., 2021;Laufer-Meiser et al., 2021). Further constraining the evolution of benthic Si-Fe coupling along the fjords requires future dedicated sampling campaigns that will obtain sediment cores and pore water samples from middle and outer fjords, in addition to fjord heads.
The sediment nutrient effluxes derived here are likely conservative estimates, as the benthic flux is not solely driven by diffusion. Other processes that contribute to benthic nutrient fluxes at fjords include benthic macrofaunal activities such as bioturbation and bio-irrigation, pore water release driven by bottom currents and physical disturbance such as iceberg scouring (Schulz & Zabel, 2006;Wehrmann et al., 2014). These processes will likely result in a net sediment efflux that is higher than the flux maintained by diffusion alone (Ng et al., 2020;Wehrmann et al., 2014). Our study provides important insights into the primary geochemical processes that govern the return fluxes of nutrient Si and Fe deposited in the sediments back to fjord waters. These benthic geochemical processes should be considered when evaluating nutrient cycling in fjords and the future evolution of nutrient delivery to downstream coastal areas amid rapidly retreating glacier cover.

Conclusions
Analyses of pore water and sequential sediment chemical extractions were performed on four sediment cores acquired from the head of the BMFC in Chilean Patagonia. The fjord sediments exhibit some of the highest published concentrations of dFe in pore waters (up to 800 mmol m −3 ), which are associated with an elevated supply of highly reactive Fe (oxy)hydroxides and other Fe-rich reactive minerals from upstream glaciers. Delivery of Fe as colloids/nanoparticulate aggregates may also play a role in influencing the isotopic composition of Fe in pore waters and its transport away from sediments, as has recently been observed in open ocean sediments.
Our stable isotope analyses and model simulation suggest that pore water DSi at the fjord sites are governed by the dissolution of glacially-derived ASi and BSi, authigenic precipitation of clays, and a close coupling between the benthic cycling of Si and Fe. Fjord sites that show lower pore water δ 56 Fe, indicative of more intense Fe redox cycling and in-situ production of Fe (oxy)hydroxides, is concurrent with elevated pore water δ 30 Si. This coupling of Si and Fe isotopes can be explained by Si adsorption onto the newly produced Fe (oxy)hydroxides from upwardly diffusing Fe 2+ rich pore water and co-precipitation of amorphous Fe-Si minerals.
The tight coupling of benthic Si and Fe cycling has an impact on the diffusive fluxes of the two nutrients across the sediment-water interface. Sediment effluxes of dFe at the Patagonian fjord sites, maintained by the high pore water concentrations, are relatively large (up to 0.02 mmol m −2 day −1 ) compared to open ocean sites (typically <0.001 mmol m −2 day −1 ). In contrast, sediment effluxes of DSi in fjord sediments (0.02-0.05 mmol m −2 day −1 ) are relatively low (typical ocean values are >0.1 mmol m −2 day −1 ). We hypothesize that upwardly diffusing Si