The rain rate of particulate organic carbon (POC) to the sediment surface and its degradation in the sediment drives all the biogeochemical processes taking place in the model. The depth and rate where POC mineralization occurs will affect the distribution of solutes in sediment porewaters, their turnover rates and the flux of solutes between the sediment and the ocean. It is therefore critical to describe POC degradation as accurately as possible.

Bulk POC can be envisaged as a continuum of reactive particulate fractions, each having its own discrete first-order rate constant (Middelburg, 1989; Boudreau and Ruddick, 1991). In the steady state version of our model, the reactivity of bulk POC pool is constrained directly from the ammonium (NH4+) profile and decreases smoothly with depth in the sediment, that is, as a reactive continuum (Bohlen et al., 2011; Dale et al., 2016). However, this approach is not suitable for non-steady state applications because of complications in defining the boundary conditions for a single POC pool whose reactivity decreases with sediment depth (i.e., over time). It is thus necessary to define several POC fractions, or “G,” that capture the spectrum of our known bulk reactivity on the Peruvian margin; otherwise known as the multi-G model (Jørgensen, 1978).

Our aim is to parameterize a multi-G model consisting of several POC fractions, i, each being degraded with a rate constant k_Gi. The fluxes of each fraction to the seafloor and the rate constants were constrained using the known distribution of POC degradation rate with depth in sediments on the Peruvian margin (RPOC). In the first step, the RPOC profile was obtained from the average model-derived RPOC profile from 3 stations (128, 142, and 195 m) on the outer shelf at 12°S which neighbor the time-series station (Dale et al., 2016). POC degradation at each of these stations was tightly constrained in that study using porewater data, POC content and benthic fluxes of dissolved inorganic carbon and nitrogen. The average RPOC profile for these stations is shown in Figure 2A (dotted red curve). RPOC decreases by three orders of magnitude over the upper 20 cm, reflecting rapid respiration of the labile organic components.

The second step is to use this rate profile to predict the corresponding POC content. The POC content and the rate profile are both needed in the subsequent step to fully constrain the parameters of the multi-G model. To achieve this, the RPOC profile was used as a forcing function in a steady state 1-D reaction-transport model that considers POC degradation as well as burial with compaction and bioturbation in the upper 20 cm. The model simply calculates the POC content for the imposed RPOC profile, POC rain rate and transport processes. The rain rate was taken as 17 mmol m⁻² d⁻¹ (204 mg m⁻² d⁻¹), which is the average rain rate determined for the outer shelf stations (Dale et al., 2015). From the same study, the sediment accumulation rate and bioturbation coefficient were determined by radiotracer measurements to be 0.15 cm year⁻¹ and 2 cm² year⁻¹, respectively. The resulting simulated surface POC content of 5.7 wt.% (Figure 2B, dotted red curve) is close to the measured value on the outer shelf of 7.5 ± 3.2 wt.% (Dale et al., 2015). The carbon burial efficiency is 50%, which is well within the range of 24–74% determined for these stations based on measured POC accumulation rates and in situ flux measurements of dissolved inorganic carbon.

Finally, in step 3, the multi-G model was fit to the RPOC and POC profiles using a reaction-transport model with several discrete POC fractions, employing identical sediment accumulation and bioturbation rates as in step 2. We found that a model with 4 POC fractions was adequate to give optimal fits to both profiles over the sediment depth range (Figure 2, dashed lines and Table 1). Fewer fractions were insufficient, and more fractions introduce additional and unwarranted parameter uncertainty. As an initial model constraint, we chose a value of k_G1 = 8 year⁻¹ and k_G2 = 0.8 year⁻¹ for the two most reactive fractions (POC1 and POC2) based on degradation experiments of fresh phytodetritus under anoxic conditions (Westrich and Berner, 1984). These appear to be satisfyingly appropriate for this setting (Figure 2). The two remaining rate constants were constrained solely from the RPOC and POC profiles. The results show that around 35% of the POC pool is highly reactive (POC1) with a further 10% attributed to POC2. The corresponding lifetimes (1/k_Gi) of these fractions are 0.13 and 1.3 year. Around half of the organic material is barely degradable on the time-scales of interest to this study (POC4), with a lifetime of 1700 years, with a smaller fraction (6%) attributed to POC3 with a lifetime of 20 years. Only the least reactive fraction survives below ca. 10 cm, and is thus mainly responsible for the 50% carbon burial efficiency in the OMZ (Dale et al., 2015).

It is important to note that the rate constants and POC fractions derived from this procedure bear additional uncertainties since they have been derived using RPOC profiles that are assumed to be in steady state but which clearly are not (see Figure 4 in Results). The reader should therefore be aware that our POC kinetics, especially of the more reactive fractions, probably better reflect the summer situation from which the RPOC profiles were derived, rather than winter. We would require data on the temporal variability in benthic mineralization to minimize this uncertainty.

The current model contains a more complete N cycle than the non-steady state version developed previously (Dale et al., 2016). Mass balances, model calculations and biogeochemical observations suggest that sulfide (H₂S) oxidation to sulfate (SO42−) coupled to nitrate (NO3−) reduction by vacuolated Thioploca spp. is a major process in Peruvian shelf sediments (Henrichs and Farrington, 1984; Fossing et al., 1995; Huettel et al., 1996; Otte et al., 1999; Bohlen et al., 2011). This process, DNRA, can be described as (Zopfi et al., 2001):

where “BNO3−” denotes NO3− stored within vacuolated bacteria. Being motile, Thioploca can directly access the bottom water NO3− reservoir by extending their filaments through the benthic boundary layer (Huettel et al., 1996). Enhanced transport of NO3− into the sediment by Thioploca is simulated in the model (Dale et al., 2016), and this leads to an increase in the intracellular biological NO3− pool (see Supplementary Material). Reduction of BNO3− during DNRA is simulated as a first-order loss term, where the rate constant was set to 30 year⁻¹, equal to a mean lifetime of around 2 weeks. This value was determined from the observed decrease in the intracellular NO3− pool during a period when Thioploca were surviving on their internal reserves (Dale et al., 2016). Simulated BNO3− concentrations are reported on a porewater basis, that is, the equivalent concentration of the internal BNO3− pool if it were distributed evenly in the porewater. The surface BNO3− content is used as a boundary condition in the model and was scaled to bottom water NO3− concentrations as described in the Supplementary Material.

On the Peruvian shelf, DNRA may produce more NH4+ in the sediment than ammonification:

where CH₂O represents POC of zero oxidation state and r_NC is the molar N-to-C ratio that links organic N, defined above as a NH₃ moiety, to POC degradation. Organic phosphorus and phosphorus dynamics in general are not considered. The role of Thioploca on P dynamics in these settings is not well understood (Lomnitz et al., 2016), and will be addressed in a forthcoming manuscript. Mineralization of POC by an oxidant thus produces NH4+ in addition to reducing equivalents and dissolved inorganic carbon (DIC) and a net change in protons (ΣH⁺). Acid-base speciation is not dealt with explicitly in this model. POC mineralization pathways include aerobic respiration, NO3− reduction to nitrite (NO2−), denitrification of NO2− to N₂, and sulfate reduction. Methanogenesis and metal oxide reduction are unimportant in surface sediments on the Peruvian shelf (Dale et al., 2016).

Initial simulations at the time-series station using the above-derived rate constant for DNRA showed that this process alone cannot completely diminish the benthic H₂S flux because the rates of sulfate reduction are so high. Given that no benthic H₂S flux was detected on the shelf at 150 m (Sommer et al., 2016), an additional H₂S removal process was required in the model to prevent its escape to the water column. An immediate increase in H₂S concentration inside benthic chambers when NO3− and NO2− fell below 0.5 μM (Dale et al., 2016) indicates that H₂S removal may be coupled to the reduction of NO3− and/or NO2−, that is, chemoautotrophic denitrification, CD (Lam and Kuypers, 2011):

We ruled out NH4+ as the end product of this reaction since the NH4+ budget is already well constrained by benthic chamber fluxes and porewater concentrations (Dale et al., 2016). The microorganisms that mediate CD are assumed to be non-vacuolated and non-motile and so rely on NO3− and NO2− that diffuses into the sediment or produced in situ by nitrification. These aspects are discussed further later on.

Heterotrophic denitrification (HD) produces N₂ during the oxidation of POC using NO2−:

The rate of HD is determined directly from the rate of POC degradation (see Supplementary Material).

Finally, anammox is the other reaction in the model that leads to fixed N loss to N₂:

Other processes included are aerobic nitrification of NH4+ and NO2− and aerobic oxidation of H₂S. We do not consider precipitation of particulate S in the current model. For typical particulate S content on the shelf of 1.5 wt.% and sediment mass accumulation rate of 1 g m⁻² d⁻¹ (Suits and Arthur, 2000; Scholz et al., 2011; Dale et al., 2016), S burial amounts to ca. 0.4 mmol m⁻² d⁻¹, equivalent to around 15% of sulfate reduction (see Results).

Model results are shown for the 12 year period between 1 January 1998 and 31 December 2009 when bottom water chemistry and rain rate data are simultaneously available. Prior to this, the model was spun-up to generate a realistic initial condition. The mean RRPOC data for the first 3 years (1998–2000, Figure 3C) was used to derive seasonally-cycling profiles of the POC fractions. At the same time, constant bottom water concentrations equivalent to 3 μM of O₂, 14 μM of NO3−, and 2 μM of NO2− were imposed. It should be noted that porewater concentrations nonetheless vary over time in this initialization step due to the seasonal variability in rain rate.

The initial condition is only an approximation of reality since time-series rain rate estimates prior to 1998 are unavailable. However, the uncertainty created at the start of the simulation (1 January 1998) will be limited by the fact that the lifetime of almost half of the POC is around 1 year or less, and the other half is mostly refractive on the time scales relevant to the simulated sediment column (Table 1). The reactive fractions will mainly determine the reactivity of the surface sediments and the solute flux across the sediment-water interface. Specifying the initial condition for O₂, NO3− and NO2− is less critical since their distribution is limited to the surface millimeters or centimeters and their residence time in the porewater is on the order of hours. All initial concentration profiles for 1 January 1998 were taken from the final day (31 December) of the spin-up simulation.

Selected geochemical profiles from the final year of the simulation are shown to illustrate the intra-annual dynamicity of the sediments (Figure 4). The model results agree well with observations made close to the time-series station in summer 2013 (Dale et al., 2016). The reactive POC fractions (1 and 2) show clear seasonal variability due to their high reactivity. With increasingly refractivity, less variability is observed. The variability in POC content is linked to the rate of POC degradation which varies by an order of magnitude in the surface layers (Figure 4B). Dissolved O₂, NO3− and NO2− show pronounced variability in the upper millimeters caused both by the changing bottom water concentrations and the input of labile organic matter (Figures 4C–E). NO3− and NO2− concentrations also occasionally show subsurface peaks. Further model simulations with constant O₂ boundary conditions (not shown) indicate that these peaks are driven by nitrification when O₂ is available in the bottom water, rather than changes in bottom water NO3− and NO2−. The results demonstrate that down-core measurements of NO3− using microsensors (e.g., Glud et al., 2009) from a single time point could not be used to accurately extrapolate NO3− fluxes over an annual cycle at this site.

Pronounced temporal variability is even seen in the solutes that have constant boundary conditions (Figures 4F–H). Yet, whereas the variability in DIC concentrations is only due to the seasonality in rain rate and POC degradation, the profiles of NH4+ show more variability since they are modified by biological NO3− transport and reduction by Thioploca (Figure 4I). The effect of anammox and nitrification on NH4+ concentrations is limited to the surface 0.5 cm where NO3− and NO2− are available.

Finally, the model predicts a lack of H₂S in the upper 3 cm except for a small subsurface peak in H₂S of ca. 20 μM that is also apparent in porewater measurements made at 101 and 128 m water depth (Dale et al., 2016). This reflects extremely high rates of POC respiration by sulfate reduction and, at the same time, very efficient microbial H₂S oxidation, but not quite enough to completely eliminate H₂S in the top centimeter.

Solute fluxes across the sediment surface demonstrate high temporal variability (Figures 5A–F and Table 3). O₂ and NO3− fluxes show no obvious trends because their rate of sediment uptake is largely determined by their bottom water concentrations which are highly variable due to displacement and exchange of shelf waters (Graco et al., 2016). O₂ flux was highest in 1998 (−4.0 mmol m⁻² d⁻¹) when bottom waters became ventilated during a strong El Niño event. It should be noted that during our fieldwork in 2008 and 2013 the O₂ flux was essentially zero on the shelf because bottom water O₂ levels were below detection limit (Table 3, Dale et al., 2016). NO2− uptake by the sediment tends to be higher in summer due to elevated POC flux (−0.3 mmol m⁻² d⁻¹) in combination with abundant NO2− in the bottom water. During winter, the rain rate is diminished and NO2− fluxes are more directed toward the water column (0.2 mmol m⁻² d⁻¹). Fluxes of NH4+ and DIC show the same seasonal trends as POC rain rate with highest values in summer and lowest in winter.

The model predicts negligible efflux of H₂S (Figure 5F) which is in line with observations on the shelf in summer (Bohlen et al., 2011; Sommer et al., 2016). Notable H₂S release events up to 1.0 mmol m⁻² d⁻¹ were predicted during early summer in other years (e.g., 2006) and at the end of summer in 1999. These emissions coincide with periods of elevated carbon rain rate, demonstrating a direct link between carbon input and H₂S release.

Organic matter mineralization in the sediment shows regular seasonal patterns driven by organic matter input (Figures 5H–N). Mean POC degradation is 5.7 mmol m⁻² d⁻¹, and ranges between a minimum of 2.6 mmol m⁻² d⁻¹ in August 1998 during the strong El Niño event to 14.0 mmol m⁻² d⁻¹ in April 2004. Summer and winter averages for the whole period are 7.4 and 3.9 mmol m⁻² d⁻¹, respectively (Table 3). On average, 83% of the POC is respired by sulfate reduction, with a mean rate of 2.4 mmol m⁻² d⁻¹ of S (Figure 5I, Table 3). Furthermore, 81% of the H₂S is consumed by DNRA, with the remainder either escaping the sediment by diffusion or oxidized by chemoautotrophic denitrification (CD). As for DNRA, CD acts as an important barrier for H₂S emissions and a sink for fixed N (NO3− and NO2−). The maximum CD rates occurred in 2004 to 2006 when rain rate and sulfate reduction were highest. Heterotrophic denitrification is equivalent to around 15% of CD. Both these processes contribute to N loss, which increases from 1.1 mmol m⁻² d⁻¹ in winter to 1.8 mmol m⁻² d⁻¹ in summer, with a mean of 1.5 mmol m⁻² d⁻¹ (Figure 5N; Table 3). These are similar to rates in sediments in the OMZ offshore Chile measured by Farías et al. (2004) and to those determined by mass balance at Peru (Bohlen et al., 2011; Dale et al., 2016; Table 3). In contrast, DNRA conserves fixed N in the system and produces twice as much NH4+ as ammonification (Figure 5M). Given that anammox consumes around one-third of NH4+ produced by ammonification and DNRA (Table 3), most of the NH4+ diffuses out of the sediment and enters the water column. This is important for the system-wide N balance, because NH4+ production by DNRA can be an important driver for N loss in the water column by anammox (Kalvelage et al., 2013).

The intra-annual biogeochemistry of the Peruvian upwelling system is well-recognized, resulting from strong seasonality in the intensity of upwelling, primary production and hydrodynamic factors (Pennington et al., 2006; Echevin et al., 2008, 2014; Gutiérrez et al., 2008; Graco et al., 2016). Yet, the corresponding variability in sediment biogeochemistry on the Peruvian shelf is only just beginning to be appreciated. This has been demonstrated here by coupling water column time-series data to an empirical diagenetic model. Time-series data provide invaluable boundary conditions for sediment models that allow the strength and direction of benthic-pelagic feedbacks to be quantified (e.g., Fossing et al., 2004; Dale et al., 2013). We have shown that the sediments on the Peruvian shelf are not static repositories that are independent of changes taking place in the water column. Rather, benthic reaction rates and fluxes vary several-fold over the course of the year in concert with upwelling and organic carbon rain rate. Considerable temporal variability in benthic feedbacks has been observed in the adjacent upwelling zone further south in Conception Bay (Chile); a site also characterized by microbial mat communities (Graco et al., 2001; Farías et al., 2004). A well-defined coupling between the benthos and water column is probably true for continental shelves in general and deeper waters too (Martin and Bender, 1988; Soetaert et al., 1996).

It is interesting to note that the seasonality in benthic biogeochemical reactions is greatly attenuated with regard to the input of organic matter to the sediment. This is illustrated by the rate of POC degradation, which is around two times higher in summer than winter whereas the rain rate is five to six times higher (Table 3). Benthic DIC fluxes are similarly attenuated. The maximum DIC flux also slightly lags behind the maximum rain rate by around 3–4 weeks (vertical lines, Figure 5E). This lag is also apparent in the NH4+ and H₂S flux as well as in the POC degradation rate (Figure 5H). Clearly, then, the POC rain rate to the seafloor does not exclusively determine the benthic solute fluxes at any given time.

These temporal lags have implications for the benthic N balance too. Under steady state conditions, the net sum of the dissolved N fluxes and ammonification across the sediment surface should equal zero. This mass balance, plotted as ΣN in Figure 5G, demonstrates that an imbalance of a few mmol m⁻² d⁻¹ on a month-by-month basis is the norm rather than the exception. There is no obvious seasonal trend, with positive and negative ΣN occurring in both summer and winter. Maximum deviations from steady state occurred in summer 1999 when ΣN increased to ca. 4 mmol m⁻² d⁻¹. Similar N imbalances based on single time point measurements have been observed on the outer Peruvian shelf (Sommer et al., 2016) and in the California Borderland basins (Chong et al., 2012).

The chief causes of the N imbalance are ammonification and DNRA. The former process directly depends on the rate of organic matter mineralization whereas the latter indirectly depends on organic matter mineralization via the rate of H₂S production by sulfate reduction. Organic matter that falls to the seafloor is not instantaneously mineralized, but decays over time according to the derived rate constants in Table 1. The seasonality in rain rate thus leads to temporal storage of solids and solutes in the sediment and, consequently, biogeochemical reactions (Martin and Bender, 1988). Similar deviations from steady state have been reported for sites with transient development of bacterial mats of filamentous Beggiatoa spp. and subsequent release of NH4+ associated with DNRA (Graco et al., 2001).

More generally, the model provides a rough idea of the uncertainty in benthic turnover rates that are derived from single time point measurements. For instance, if benthic DIC fluxes collected in summer are used to determine total benthic respiration rates using mass balances or modeling (e.g., Bohlen et al., 2011; Dale et al., 2016), Figure 5E and Table 3 inform us that the derived rates will have a 30% margin of error (deviation from annual mean). The error using data from winter rises to around 40%. This is greater than the typical uncertainty in measured fluxes due to sediment heterogeneity and experimental artifacts (Hammond et al., 2004; Sommer et al., 2016).

Perhaps one of the most important results from the entire modeling exercise is that the sediments must exert a strong feedback on the nutrient cycles in the water column. For example, benthic N loss at the time-series site (1.5 mmol m⁻² d⁻¹; Table 3) is up to 20% of water column N loss by anammox reported for 140 m water depth at a site close to the time-series station (6.2 mmol m⁻² d⁻¹; Hamersley et al., 2007). Furthermore, NH4+ produced by DNRA can easily cover the NH4+ required by anammox in the water column (Dale et al., 2016). Sediments thus likely play an important direct and indirect role in the pelagic fixed N budget. These basic feedbacks need to be considered in regional models of the Peruvian system. Later on, we suggest ways in which this can be achieved.

Reports of benthic emissions of H₂S on the Peruvian margin are rare (Schunck et al., 2013). Benthic chambers deployed so far on the shelf have only measured H₂S release where the sediments receive massive amounts of phytodetritus (>60 mmol m⁻² d⁻¹ of C) and where bottom waters are depleted in electron acceptors (Sommer et al., 2016). At a station close to the time-series station managed by IMARPE and at a similar water depth (150 m), no H₂S emissions were detected (Sommer et al., 2016). To simulate this in our model, an efficient sink of H₂S in addition to DNRA was required because the rate constant of DNRA was pre-determined based on observed temporal decreases in the BNO3− pool (see section Model). An immediate increase in H₂S concentrations inside benthic chambers following depletion of NO3− and NO2− to below 0.5 μM further alludes to second pathway of H₂S removal coupled to the reduction of NO3− and/or NO2− (Dale et al., 2016). Such an immediate increase may not be expected of NO3−-storing microbial communities that can continue to respire H₂S in the short-term absence of electron acceptors.

In the present model, this sink was assumed to be chemotrophic denitrification (CD) coupled to H₂S oxidation by non-vacuolated bacteria [Equation (3) and (4)]. It would seem reasonable that chemoautotrophic denitrifying bacteria are present and able to compete with Thioploca for H₂S in the organic-rich anoxic Peruvian muds displaying high rates of sulfate reduction (Fossing et al., 2004). These bacteria have been isolated in OMZs and the reaction is highly favorable energetically (Lam and Kuypers, 2011). Other H₂S removal processes do not greatly help to close the mass balance. Aerobic H₂S oxidation is negligible (Table 3), and burial of particulate sulfur (e.g., pyrite) was not considered since, as mentioned (Section Benthic N Cycle), it can only sequester around 15% of total H₂S production. Although, this is roughly the same as the fraction of H₂S consumed by CD, we infer from the immediate increase in H₂S concentrations inside benthic chambers when NO3−/NO2− concentrations fall below 0.5 μM that the sink must be occurring directly at the sediment-water interface (Dale et al., 2016). It is questionable whether the precipitation of iron sulfides will be thermodynamically favorable at the sediment-water interface due to low reactant concentrations and the presence of potential oxidants for reduced sulfur such as O₂ and NO3− (Rickard, 1997). Although, the role of sedimentary CD has not yet been quantified in the field, the model has identified an additional process in the benthic S cycle that may turn out to be just as important as DNRA in maintaining a H₂S-free water column.

Nonetheless, as mentioned, benthic microorganisms are unable to consume H₂S completely when rain rates are exceptional, as seen by the peaks in H₂S emissions (Figure 5F). The peak in summer 2009 coincides with the largest sulfidic event ever measured in the water column on the Peruvian margin (Schunck et al., 2013). These authors further attributed the H₂S plume to benthic H₂S release, and estimated a benthic H₂S source of around 2 mmol m⁻² d⁻¹ on the middle shelf at 100 m water depth. This is in agreement with the present model prediction (0.4 mmol m⁻² d⁻¹) when one considers that organic matter mineralization at 100 m is a factor of 3–4 higher than at 150 m (Dale et al., 2015). The model results suggest that large-scale H₂S emission events may occur more regularly than currently presumed and may have gone undetected in 2004, 2005, and 2006 (Figure 5F). Furthermore, given that H₂S release appears at times when primary production exceeds ca. 20 mmol m⁻² d⁻¹ (c.f. Figure 3A), satellite-derived data could be used as an early warning system to predict sulfidic events on the Peruvian margin. The feasibility of this idea is certainly worth further investigation, especially given the importance of the fishery sector to the Peruvian economy.

Denitrification in continental margin sediments constitutes the largest sink of fixed N in the ocean today (Gruber, 2008). Although, sediments underlying oxygen deficient waters account for less than 1% of the total seafloor, globally they account for 10% of benthic N loss (Bohlen et al., 2012). On average, 1 mole of fixed N is transformed to N₂ for every 4 moles of POC oxidized at the time-series station (Table 3), compared to a ratio of 1:10 in oxic shelf settings (Middelburg et al., 1996; Bohlen et al., 2012). As explained above, benthic denitrification must modulate the biogeochemistry of the water column on the Peruvian shelf. Given that N loss limits productivity, regional and global biogeochemical models that simulate oxygen deficient waters need to be able to predict benthic N loss accurately.

Several approaches of varying complexity exist to couple benthic and pelagic models (reviewed by Soetaert et al., 2000). At the one end, time-dependent, vertically-resolved diagenetic models, such as the one presented here, can be coupled to a water column biogeochemistry model (e.g., Archer et al., 2002). In practice, this approach is seldom used due to the excessive computational demands and parameterization required. At the other end, over-simplified approaches that include the specification of fixed or reflective sea floor boundary conditions in ocean models inadequately capture the diagenetic transformations of organic matter and benthic sources and sinks. Intermediate complexity approaches, such as vertically-integrated sediment models, are a pragmatic alternative. These allow sediment fluxes to be parameterized using zero-dimensional time-dependent approaches (i.e., a single box) at the lower ocean boundary, or by the implementation of transfer functions to predict the immediate response of sediments to changes occurring at the sediment-seawater boundary. These models can be derived from statistical analysis of vertically-resolved benthic models or be purely empirical in nature (Middelburg et al., 1996; Soetaert et al., 2000; Capet et al., 2016).

An example of an empirical transfer function for predicting fixed N loss in global models has been developed by Bohlen et al. (2012). Here, denitrification (all N₂ production pathways in mmol m⁻² d⁻¹ of N) can be predicted from the POC rain rate (mmol m⁻² d⁻¹ of C) and the bottom water concentrations of O₂ and NO3− in μM:

where a (0.060), b (0.19) and c (0.99) are constants determined from a global database of measured benthic O₂ and NO3− fluxes and rain rates. The function has been coupled to Earth System Models to predict benthic denitrification in the modern ocean (e.g., Somes et al., 2013; Yang and Gruber, 2016). The function is also dynamic in the sense that N loss responds instantaneously to changes in the boundary conditions, but does not consider the storage of solutes and particulates in sediments.

For the application of Equation (7) to the Peruvian margin, RRPOC has to be adjusted to account for the elevated carbon burial efficiency (CBE) in sediments under anoxic waters. The more carbon buried, the less degradation and thus less N₂ loss. Previously, the CBE was shown to be 50% for sediments within the OMZ (section POC Degradation and Dale et al., 2015). RRPOC was thus corrected by multiplying by the factor 1− CBE/100:

Fixed N loss predicted using Equation (8) is shown alongside the N loss from the numerical model in Figure 6 jointly with the resultant residuals. The two approaches show a very good agreement; only during periods when the POC rain rate rises above ca. 25 mmol m⁻² d⁻¹ does Equation (8) over-predict the N loss by more than two standard deviations (2σ) of the residuals, such as in 2004 and 2005. This is probably because extreme settings such as the Peruvian shelf that display high rates of DNRA at the expense of denitrification are not well represented in the database used to derive the transfer function. Nonetheless, over the whole time-series, the models agree to within 30%, which is within the empirical uncertainty of the transfer function (Bohlen et al., 2012). The model could be further generalized by making CBE a function of, for example, rain rate (Dunne et al., 2007).

It can nonetheless be seen that the yearly maxima in N loss predicted by Equation (8) precede the model result by 1 to 2 months (Figure 6). This is due to the time lag between POC deposition and organic matter mineralization. An improved transfer function can be formulated that considers these lags by filtering RRPOC using a low pass moving average:

where n denotes the number of months, i, over which the moving average was applied. Best results were found for n = 4 (Figure 6, blue curve), for which the mean and standard deviation of the N loss (1.3 ± 0.8 mmol m⁻² d⁻¹) are within excellent agreement of the model result (1.5 ± 0.5 mmol m⁻² d⁻¹). The 4 month filter captures much of the variability induced by the reactivity of the POC1 and POC2 fractions (Table 1) as well as the time required for molecules to diffuse over the upper 10 cm in response to the changes in rain rate.

On the one hand, these are important results since they show that these functions can be coupled to dynamic regional biogeochemical models of the Peru system such as the Regional Ocean Modeling System (e.g., Vergara et al., 2016). On the other hand, it is perhaps surprising that the N loss predicted by the vertically–resolved and—integrated models agree so well. Nitrogen loss in the numerical model depends on the temporally-variable vertical distribution of organic matter and dissolved species and hence their rates of advective and diffusive transport in the sediment column. The history of the previous boundary conditions is thus entrained in the numerical model. The transfer functions, in contrast, retain no memory of past conditions; only the instantaneous rain rate and bottom water concentrations are important. This is partly compensated for in Equation (9) by the low pass filter on the rain rate. Storage of POC seems not to greatly affect the predictive capability of the functions because the most labile carbon pool, with a lifetime of a few weeks, constitutes a large fraction of bulk POC (Table 1). Furthermore, almost all denitrification takes place within the upper 5 mm of the sediment (Figure 4C). The time scale for a NO3− molecule to diffuse across this distance is only a few hours, such that the relatively gradual changes of solute concentrations in the bottom water allow benthic denitrification to proceed in quasi-steady state. The accuracy of the transfer functions would begin to break down in situations where bottom water chemistry fluctuates on hourly time scales due to e.g., tides.

Almost all of the N loss from the water column in the Eastern Tropical South Pacific occurs by anammox rather than heterotrophic denitrification (Hamersley et al., 2007; Kalvelage et al., 2013). Mass balance calculations suggest that around 40% of the N loss takes place within the coastal OMZ where water depths are <600 m (Kalvelage et al., 2013). Moreover, it has been suggested that NH4+ production by Thioploca could contribute as much as 50% of the total NH4+ oxidized by anammox (Dale et al., 2016). This potentially important feedback linking the benthic and pelagic N cycle is thus focused on shallower waters. Its impact on the wider biogeochemistry of the Peruvian OMZ requires investigation using a regional model that accounts for N transformations in both the water column and sediments.

Analogous to Equation (9), a simple vertically-integrated model for benthic DNRA (mmol m⁻² d⁻¹ of N) can be proposed based on the model findings:

where r_SC is the sulfate-to-carbon ratio during POC respiration by sulfate reduction (0.5) and f is the mean fraction of H₂S produced by sulfate reduction that is oxidized by DNRA (81%, Table 3). The coefficient 0.83 accounts for the fraction of POC respired by sulfate reduction at the time-series station (Figure 5I).

Results with this function, again with n = 4, show a very good agreement with the modeled DNRA rates (Figure 7). Bottom water NO3− is not a variable in this function. DNRA is more sensitive to RRPOC than changes in bottom water NO3− because Thioploca, along with some species of Beggiatoaceae and Thiomargarita, are able to store NO3− within their vacuoles, and thus bridge periods where NO3− availability is low or absent (Schulz et al., 1996; Preisler et al., 2007). The equation would be straightforward to implement as a lower water column boundary condition in a regional model. Applying the function from the shelf down to 400 m, where Thioploca are distributed (Mosch et al., 2012; Sommer et al., 2016), would allow a first-order estimate of benthic NO3− conversion to NH4+ by Thioploca for the entire Peruvian margin to be determined, as well as the impact of benthic DNRA on water column N loss. Ideally, such an approach would also have to consider that the sediments north of 10°S and south of 15°S are sandy due to sediment winnowing by strong bottom currents (Suess et al., 1987), and thus unlikely to support Thioploca communities.

One drawback of this approach is that the ecology of Thioploca is treated in a very simplified manner. Time-series analysis of Thioploca biomass at a nearby station (94 m water depth) between 1993 and 2006 showed that the bacteria are mostly present when O₂ levels are <20 μM but not anoxic (Gutiérrez et al., 2008). These conditions characterize the period of the current model application, with the exception of the strong El Niño in 1998 when bottom waters became well ventilated for several months due to onshore intrusion of Subtropical Surface Water (Figure 3E; Graco et al., 2016). Following this event, Thioploca biomass at the other site (Gutiérrez et al., 2008) was diminished for 2–3 years during which time the benthic biota became dominated by amphipods and polychaetes. Although, Thioploca biomass does not correlate statistically with O₂ concentrations (Gutiérrez et al., 2008), the presence or absence of bacterial mats very likely depends on the occurrence of coastal trapped waves that depress the OMZ downwards and ventilate shelf sediments. The frequency of these waves increases during El Niño, and this sustained ventilation is believed to allow burrowing macrofaunal communities to be established at the expense of Thioploca (Gutiérrez et al., 2008). Furthermore, Thioploca populations can be diminished during La Niña events when bottom waters stagnate and NO3− is exhausted, or following the sinking of massive plankton blooms that smother the seafloor biota (Graco et al., 2001; Gutiérrez et al., 2008). These extreme low frequency (e.g., ENSO) and high frequency forcings (e.g., bottom water stagnation) lead to complex dynamics of Thioploca biomass and DNRA that are not dealt with in the current model.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.