Parameterization of biogeochemical sediment – water fluxes using in situ measurements and a diagenetic model

Diagenetic processes are important drivers of water column biogeochemistry in coastal areas. For example, sediment oxygen consumption can be a significant contributor to oxygen depletion in hypoxic systems, and sediment– water nutrient fluxes support primary productivity in the overlying water column. Moreover, nonlinearities develop between bottom water conditions and sediment–water fluxes due to loss of oxygen-dependent processes in the sediment as oxygen becomes depleted in bottom waters. Yet, sediment– water fluxes of chemical species are often parameterized crudely in coupled physical–biogeochemical models, using simple linear parameterizations that are only poorly constrained by observations. Diagenetic models that represent sediment biogeochemistry are available, but rarely are coupled to water column biogeochemical models because they are computationally expensive. Here, we apply a method that efficiently parameterizes sediment–water fluxes of oxygen, nitrate and ammonium by combining in situ measurements, a diagenetic model and a parameter optimization method. As a proof of concept, we apply this method to the Louisiana Shelf where high primary production, stimulated by excessive nutrient loads from the Mississippi–Atchafalaya River system, promotes the development of hypoxic bottom waters in summer. The parameterized sediment–water fluxes represent nonlinear feedbacks between water column and sediment processes at low bottom water oxygen concentrations, which may persist for long periods (weeks to months) in hypoxic systems such as the Louisiana Shelf. This method can be applied to other systems and is particularly relevant for shallow coastal and estuarine waters where the interaction between sediment and water column is strong and hypoxia is prone to occur due to land-based nutrient loads.


Introduction
Sediment biogeochemistry represents a major component of elemental cycling on continental margins (Middelburg and Soetaert, 2005;Liu et al., 2010).In these shallow, productive areas on average 30 % of photosynthetically produced organic matter is deposited and recycled in the sediment (Wollast, 1998).The recycling of this organic material consumes oxygen (O 2 ) and can result in either a source or a sink of nutrients to the water column (Cowan and Boynton, 1996).For instance, a proportion of the deposited organic matter is remineralized via denitrification, which produces biologically unavailable N 2 gas.Denitrification represents a major removal pathway for nitrogen (N) in coastal areas (Fennel et al., 2009;Bohlen et al., 2012) and buffers the effects of excessive N loads in eutrophic systems (Seitzinger and Nixon, 1985).In this type of environment, high respiration rates in the water column and in the sediment may lead to bottom O 2 depletion under stratified conditions, resulting in bottom water hypoxia (O 2 < 62.5 mmol O 2 m −3 ) or anoxia (absence of O 2 ).Under low O 2 conditions, coupled nitrification-denitrification in the sediment is inhibited, and remineralized N may return entirely to the water column as ammonium (NH + 4 ), readily available to primary producers, which constitutes a positive feedback on eutrophication (Kemp et al., 1990).Conversely, N removal into N 2 may increase due to direct denitrification or due to anammox if a source of nitrate/nitrite is available (Neubacher et al., 2012).O 2 -dependent sediment-water interactions are therefore particularly important in low O 2 environments.
Clearly, the strong benthic-pelagic interaction is a key aspect of coastal biogeochemistry that needs to be represented Published by Copernicus Publications on behalf of the European Geosciences Union.
A. Laurent et al.: Parameterization of sediment-water fluxes on the Louisiana Shelf accurately in biogeochemical models.However, sedimentwater fluxes in models are often difficult to parameterize, being poorly constrained by observations.One of the simplest approaches to parameterizing sediment-water fluxes is using a reflective boundary where fluxes are proportional to particulate organic matter (POM) deposition (e.g., Fennel et al., 2006).Empirical relationships can be used to represent sediment biogeochemical processes, such as denitrification (Fennel et al., 2009) or sediment O 2 consumption (SOC) (Hetland and DiMarco, 2008).An advantage of these firstorder sediment-water flux parameterizations is that they are computationally extremely efficient and can be sufficient depending on the type of environment and the focus of the study (Wilson et al., 2013).However, sediment-water flux parameterizations are a coarse representation of sedimentwater interaction, and typically do not capture nonlinearities in nutrient fluxes that occur under hypoxic/anoxic conditions.Moreover, the choice of parameterization can have a significant effect on model results as shown in Fennel et al. (2013) where different parameterizations of SOC led to dramatically different regions of hypoxia.
Mechanistic models of diagenesis are more realistic representations of sediment biogeochemistry (Rabouille and Gaillard, 1991;Soetaert and Herman, 1995;Soetaert et al., 1996a;DiToro, 2001;Meysman et al., 2003a, b).They are forced by POM deposition and bottom water conditions, and simulate aerobic and anaerobic remineralization pathways including processes such as nitrification, denitrification, the anaerobic production of reduced substances -represented either explicitly or lumped together in O 2 demand units (ODU) -and the resulting flux of O 2 and nutrients across the sediment-water interface.While these models have been useful for studies of sediment biogeochemistry (Middelburg et al., 1996;Soetaert et al., 1996b;Boudreau et al, 1998;Meysman et al., 2003b) and for improving our understanding of sediment-water interactions (Katsev et al., 2007;Reed et al., 2011), their coupling to water column processes in biogeochemical circulation models is often limited or done at the expense of spatial resolution (Eldridge and Roelke, 2010) because of the increased computational cost.Furthermore, the diagenetic model parameter sets are often poorly constrained by observations, and therefore these models do not necessarily perform better than the simple parameterizations (Wilson et al., 2013).
An alternative, computationally more efficient approach is to parameterize sediment-water fluxes from a diagenetic model using a meta-model of diagenetic processes, as recommended by Soetaert et al. (2000).Their mass conservative method is more realistic than the simple reflective boundary and computationally more efficient than a mechanistic model of diagenesis.The method requires addition of a vertically integrated pool of sedimentary particulate organic matter for each horizontal grid cell thus enabling a mass balanced approach, but adding a layer of complexity to the water column model.Here we further simplify the meta-modeling method of Soetaert et al. (2000) by direct metamodeling of sediment-water fluxes.Our method parameterizes sediment-water fluxes of O 2 , NO − 3 and NH + 4 in a coupled biogeochemical-circulation model using in situ measurements, a mechanistic model of early diagenesis and a parameter optimization technique.The method is universal but its application is region-specific due to the local characteristics of the sediment, e.g., sediment quality (POM concentration and lability), type (porosity) and species composition (bioturbation) that influence local sediment biogeochemistry and sediment-water fluxes and are reflected in the choice of diagenetic model parameters.We apply this method to the Louisiana Shelf in the northern Gulf of Mexico, where hypoxia develops annually due to eutrophication (Rabalais et al., 2002).
First, we calibrate the diagenetic model with the help of a genetic optimization algorithm using a set of observations collected on the Louisiana Shelf.We then implement the calibrated model to simulate time-resolved sediment biogeochemistry in the region and use the model results to compute a meta-model parameterization of sediment-water fluxes for O 2 , NH + 4 and NO − 3 .Finally, we compare the fluxes parameterized with the meta-model with previous relationships used for the Louisiana Shelf.

Observations
The data used for optimization of the diagenetic model parameters were collected at two locations along the 20 m isobath on the Louisiana Shelf (Fig. 1) during three cruises in April, June and September 2006 (Murrell et al., 2013).The two locations experience hypoxia in summer but have distinct hydrographic and biological regimes.Station Z02 (see Murrell et al., 2013, for details on sampling design) is located off Terrebonne Bay on the eastern Louisiana Shelf and is influenced by river discharges from the Mississippi Delta with high primary productivity and high POM depositional flux.Station Z03 is located southwest of Atchafalaya Bay on the western Louisiana Shelf with somewhat higher salinity and lower chlorophyll concentrations than station Z02 (Lehrter et al., 2009;2012).The data set includes bottom water properties (temperature, salinity, O 2 and nutrients, Table 1), sediment-water fluxes (O 2 , nutrients) and NH + 4 sediment profiles (Fig. 2).On each date, eight sediment cores were collected at each station (three for O 2 flux, three for nutrient fluxes and two for sediment profiles).O 2 and nutrient fluxes were measured on site from triplicate individual incubations in sediment chambers.Sediment NH + 4 concentration was measured for each 2 cm bin in the duplicate sediment cores.Bottom water temperature and salinity were measured with a CTD, whereas O 2 and nutrient concentrations were measured in the water overlying the sediment cores.Details on the data set are available in Lehrter et al. (2012), Murrell et al. (2013) and Devereux et al. (2015).

Sediment-flux parameterization
The parameterization of sediment-water fluxes was derived using output from a diagenetic model.The diagenetic model was first optimized using the observational data set described in the previous section.The optimized diagenetic model was then run multiple times to derive meta-model parameterizations.

Diagenetic model
The diagenetic model represents the dynamics of the key constituents of the sediment (solids and pore water) involved in early diagenesis, as formulated by Soetaert et al. (1996a, b).The model is vertically resolved, and represents the upper 10 cm of the sediment using 10 layers with increasing resolution toward the surface.The diagenetic model has six state variables: the solid volume of organic carbon (OC), which is split into a labile class (which remineralizes rapidly) and a refractory class (which remineralizes slowly), NH + 4 , NO − 3 , O 2 and ODU.Reduced substances produced by anoxic remineralization are added to the ODU pool rather than being explicitly modeled.Model processes include aerobic remineralization, nitrification, denitrification, anaerobic remineralization and ODU oxidation.Dissimilatory nitrate reduction to ammonium (DNRA) and anaerobic ammonium oxidation (anammox) are not explicitly represented in the model.Vertical transport of solid and pore water constituents depend on sedimentation of POM to the sediment, and on diffusion, bioturbation and permanent burial.The burial of ODU refers to the deposition of ODUs as solids (e.g., pyrite, manganese carbonate) below the bioturbated zone (Soetaert et al., 1996a).The model simulates sediment-water fluxes of pore water constituents, namely, NH + 4 , NO − 3 , O 2 and ODU.We assume that ODUs are oxidized instantaneously in the water column when O 2 is available.Therefore, the net O 2 flux into the sediment is the addition of the direct O 2 flux necessary for nitrification, oxidation of ODUs and of POM in the sediment, termed SOC, plus the O 2 sink in bottom waters necessary to oxidize any ODU efflux from the sediment.
The original model of Soetaert et al. (1996a, b) was modified as follows.A temperature dependency was introduced for the remineralization of the two organic matter pools and the bioturbation of solids following a Q 10 relationship such that where R i (T ) and R T b i (yr −1 ) are the remineralization or bioturbation at ambient temperature (T ; • C) and at the base temperature (T b ; • C; i.e., R T b 1 and R T b 2 for remineralization and Dbio 0 for bioturbation, Table 2) and θ is the Q 10 factor.In the updated model, temperature thus influences the solute diffusivity, the degradability of the two OM pools and bioturbation.This modification allows for the representation of temperature dependence of microbial processes in the sediment (aerobic respiration, denitrification and anaerobic metabolism), which is known to be important in coastal systems (see, e.g., Fig. 5 in Wilson et al., 2013).Nitrification is not temperature dependent in the diagenetic model.It is assumed that O 2 concentration is the main factor limiting nitrification in the Louisiana Shelf sediments.Non-local mixing of pore water constituents due to bioturbation (irrigation) was also introduced and formulated following Boudreau (1997) such that where I (z) (µmol L −1 yr −1 ) is the irrigation at depth z, and C ow and C(z) (µmol L −1 ) are the solute concentrations at the sediment-water interface and at depth z in the sediment, respectively.α(z) is the rate of non-local exchanges at depth z such that α(z) = α 0 × f (z), where α 0 (yr −1 ) is the rate at z = 0 and f (z) is a function representing the decay of α with depth.Here, f (z) is the same function as for the bioturbation of solids (Soetaert et al., 1996a).Bioturbation and non-local mixing of solutes are not dependent on O 2 in the model.Such a dependence could be introduced to account for repeated cycles of eradication/re-establishment of macrofauna due to anoxia.However, given the limited information on the relationship between porewater O 2 , infauna biomass and irriga-tion in this region (Eldridge and Morse, 2008), we assumed that macrobiota does not re-establish itself in the regions affected by recurring severe seasonal hypoxia or anoxia on the Louisiana Shelf, and thus do not expect a strong dependence of bioturbation and bioirrigation on O 2 .The model has a total of 36 parameters (Table 2).Sediment porosity parameters were chosen to obtain a porosity profile that is within the range observed on the Louisiana Shelf.Given a lack of observations, the nitrogen to carbon ratio (N : C; mol N (mol C) −1 ) of the labile and refractory fraction of OC were fixed to constant values following Wilson et al. (2013).The assumption is that N : C follows Redfield (Redfield et al., 1963) in the labile fraction (N : C = 0.15), whereas the proportion of carbon increases in the refractory fraction (N : C = 0.1).Since deposited OC mainly originates from local primary production on the shallow Louisiana Shelf (Redalje et al., 1994;Justić et al., 1996;Rowe and Chapman, 2002), labile OC is assumed to repre-   sent 74 % of total OC in deposited material.This value was used by Soetaert et al. (1996a) to represent the fraction of labile organic matter in surface waters and is in line with previous modeling investigations of the Louisiana Shelf (Justić et al., 1996;Eldridge and Morse, 2008).However, inshore areas adjacent to river discharge may have higher fraction of terrestrial organic matter.The exponential decay coefficient for bioturbation was set as in the original model (Soetaert et al., 1996a).Solute-specific diffusion coefficients (D T i ; cm 2 d −1 ) at ambient temperature T were calculated following Soetaert et al. (1996a) and Li and Gregory (1974) ) is the solutespecific temperature dependency coefficient (Table 2).The 20 remaining parameters of the diagenetic model (Table 2) were optimized to obtain the best match between the observed and simulated sediment profiles and sediment-water fluxes.

Parameter optimization
The diagenetic model parameters were first optimized to match the sediment-water fluxes and sediment NH + 4 concentrations observed in April, June and September 2006 at stations Z02 and Z03.The sampling frequency at these stations did not allow construction of a reasonable time-dependent forcing data set for the diagenetic model (i.e., solute concentrations in overlying water, POM deposition).Thus, we did not run the optimization in a time-dependent mode; instead the model was run for 300 days with constant forcing for  1).Since no observations of POM depositional flux were available, POM depositional fluxes were prescribed using monthly means calculated for stations Z02 and Z03 from a multiyear biogeochemical model simulation (see Sect. 2.2.3).The mean depositional fluxes do not represent short-lived deposition events, which is appropriate for a model with constant forcing.
Optimization of the parameter set was carried out with the help of an evolutionary algorithm.This stochastic technique mimics natural selection by iteratively selecting the "fittest" set of parameters to reproduce the observations.The evolutionary algorithm is a well accepted method for optimization problems (Hibbert, 1993;Fogel, 1994;Chatterjee et al., 1996;Kolda et al., 2003) and has been increasingly used to optimize parameters in biogeochemical models (Kuhn et al., 2015;Robson et al., 2008;Schartau and Oschlies, 2003;Ward et al., 2010).The technique was successfully used for the optimization of parameters of Soetaert et al.'s (1996a) diagenetic model in two independent studies (Wilson et al., 2013;Wood et al., 2013).The advantage of the evolutionary algorithm over traditionally used gradient-descent algorithms is that it explores the parameter space with an element of randomness and therefore is less prone to converging on a local minimum.Each parameter is given a range of variation within which the algorithm will search for the best value to match the observations.Regardless of which minimization technique is used, gradient-descent or an evolutionary algorithm, some parameters may not be identifiable because they are unconstrained by the available observations (Soetaert et al., 1998;Fennel et al., 2001).
The evolutionary algorithm works as follows.Each set of parameters is considered to be a single individual.An initial set of n individuals includes the initial parameter set and n − 1 individuals generated randomly from this initial set of parameters through the addition of log-normally distributed random noise.The diagenetic model is run with the n parameter sets, and the difference between the results and observations is quantified using a cost function, which measures the misfit between the observations and their model counterparts.The fittest n/2 individuals, i.e., those with the lowest cost, become the parent population, and a next generation of n/2 individuals (child population) is created by recombination of the parameters from the fitter half of the population and by mutation, which occurs through the addition of random noise.The model is run again for all the parameter sets of the child population, and the above procedure is repeated for k generations.The fittest individual after k generations is the optimized parameter set.Here, we used n = 30 population members and k = 200 generations.The chosen value of k is large enough to allow the results to converge.
Ideally a single parameter set should capture the temporal and spatial variability of sediment processes throughout the Louisiana Shelf.For this reason, the diagenetic model was run with identical parameters in all six model configurations (three dates, two locations), each corresponding to a set of observed bottom water conditions plus estimated F POM (Table 1).Model results were compared with their corresponding set of sediment observations (NH + 4 porewater concentrations and sediment-water fluxes) using a cost function that includes all model variables at the six locations/times.The smaller the cost, the fitter an individual (i.e., parameter set) during the evolutionary optimization process.The cost function F for the parameter set p was calculated as follows: where s refers to locations Z02 and Z03, t is the sampling date (three in 2006) and i is the observation type: three sediment-water fluxes (SOC, NH + 4 and NO − 3 ) and one sediment profile (NH + 4 ).X obs and X mod represent the observed and simulated variable, respectively; σ 2 s,t,i is the observation standard deviation; and 1/w i represents the weight of each variable in the cost function.The values of w i were calculated for each variable i as the cost of a diagenetic model run using the initial parameter set p 0 such that w i = F i (p 0 ).The weight gives the variables approximately equal influence on the overall cost, at least initially.The weighting approach is common in parameter optimization studies (see, e.g., Friedrichs, 2001;Schartau and Oschlies, 2003;Friedrichs et al., 2007;Kane et al., 2011).To avoid biasing the cost calculation toward the NH + 4 profiles, we computed an average cost per profile.
The sensitivity of the optimized model to parameter changes was assessed by successively varying each parameter by ±50 % and calculating the change in the total cost.Then the influence of observations and forcing data sets on the optimization results was assessed as follows.First, the optimization was carried out for each station individually (to obtain site-specific parameters); then sediment profiles were excluded from the optimization (to obtain site-specific parameters optimized for flux data only); and, finally, POM depositional fluxes were included as additional parameters in the optimization rather than prescribed (to obtain sitespecific parameters and F POM optimized for flux data only).

Meta-modeling procedure
Our meta-modeling procedure parameterizes sedimentwater fluxes by means of a multivariate regression model that relates bottom water conditions and depositional flux to sediment-water fluxes, and was used here to parameterize Louisiana Shelf fluxes at the sediment-water interface.Using a meta-model of sediment-water fluxes is a simplification of the method proposed by Soetaert et al. (2000), who used a meta-model of diagenetic processes (rates) instead.The aim of our technique is to combine the simplicity and efficiency of a sediment-water flux parameterization with the realism of a diagenetic model.It is important to note that our simplified meta-model is not mass conservative; however, as long as the method is used for the system for which it was developed and within the range of conditions that were used for the parameterization, violation of mass conservation should be minor.An advantage of our simplification is that it does not require knowledge of integrated POM concentration in the sediment.
In order to obtain the meta-model parameterization the diagenetic model was run many times in time-varying mode using the single parameter set optimized for the Louisiana Shelf.The diagenetic model was forced with multi-year time series of bottom water conditions obtained from a biogeochemical circulation model of the Louisiana Shelf based on the Regional Ocean Modeling System (ROMS; Fig. 3).The simulation is described in Fennel et al. (2013;case B20clim) and covers the period from 2004 to 2009.The same simulation was used to prescribe POM depositional fluxes during the parameter optimization.For details on the model set up and validation we refer the reader to Fennel et al. (2013).We included only those grid cells on the Louisiana Shelf (z < 50 m) and west of the Mississippi River delta.Each grid cell (3791 in total) provides a time series of bottom water temperature, salinity, NO − 3 , NH + 4 , O 2 and POM depositional flux conditions that was used to run the optimized diagenetic model.We consider 2004 as a spin-up year for the diagenetic model and selected the period 2005-2009 for analysis.Half of the data from each simulation were randomly chosen to derive the meta-model.The multivariate meta-model regressions were then calculated to relate bottom water conditions and depositional flux (model inputs) to the corresponding sediment-water fluxes (model output) using the 3.45 × 10 6 data vectors.To validate the meta-model, we calculated correlation coefficients between the remaining data of each diagenetic model simulation (i.e., at each model grid location) and the corresponding meta-model results.
Each regression model is expressed as follows: where each x i corresponds to an explanatory variable i, and a, b i , c i and d i are the coefficients for the zero-order term, the regular term (x i ), the squared term (x 2 i ) and the cubic term (x 3 i ), respectively.

Other flux parameterizations
The while a fraction of N is lost through denitrification.IR is formulated as follows (Fennel et al., 2006(Fennel et al., , 2009)): The other two parameterizations assume that SOC depends on bottom water O 2 and temperature (T ) only and ignore POM deposition.One, referred to as H&D (Eq.7), is from Hetland and DiMarco (2008) as M&L (Eq.8), is from Murrell and Lehrter (2011), with a temperature dependence added by Fennel et al. (2013).Sediment-water O 2 fluxes are formulated as follows: For each parameterization x, the sediment-water NH + 4 flux is a function of SOC such that with r NH + 4 : SOC = 0.036 mmol NH + 4 per mmol O 2 .

Diagenetic model parameter optimization
Optimization of the diagenetic model parameters lowered the cost function (Eq. 3) significantly compared to the original parameter set (Table 3).NH + 4 profiles and sedimentwater fluxes simulated with the optimized parameters are, in most cases, within 2 standard deviations of the observations (Fig. 2).Simulated O 2 fluxes match the observations at station Z02 but are underestimated somewhat in April and June at station Z03.Observed O 2 fluxes are relatively high in April and June at station Z03 despite low sediment-water nutrient fluxes and NH + 4 concentration in the sediment.Observed O 2 flux had a very large standard deviation in April at station Z03 and therefore did not influence the optimization.NH + 4 and NO − 3 fluxes represent a more difficult problem for the optimization and therefore their cost is larger, especially at station Z03.Overall, sediment-water fluxes are better simulated at station Z02, and therefore station Z03 contributes more to the total cost for the optimized parameter set (Table 3).Temporal variations in NH + 4 and NO − 3 fluxes are in qualitative agreement with observations although the model underestimates their magnitudes (Fig. 2).The model is able to simulate observed NO − 3 flux realistically, in particular the observed NO − 3 flux into the sediment under low bottom O 2 conditions (Fig. 2).Within the sediment, simulated NH + 4 concentrations agree with observations in April and June, but are underestimated in September.High NH + 4 concentrations were observed at station Z02 at this time despite low NH + 4 effluxes from the sediment.Note that the observations have large standard deviations for this case and therefore this NH + 4 sediment profile had only a small influence on the optimization.Some of the observed NH + 4 profiles in April and September display a gradient at depth (Fig. 2) that the diagenetic model might not be able to resolve.There is also a deep negative gradient in the simulated profiles in April, indicating that the model did not reach full steadystate conditions at depth.However, this mismatch at depth has a limited effect on sediment-water fluxes.
Within the optimized parameter set, several parameter values reached the lower or upper edge of their allowed range, which can be informative about the dynamics of the system (Table 2).Except for the bioturbation diffusivity (Dbio 0 ), all other parameters associated with bioturbation reduced the effect of bioturbation on sediment-water fluxes over the course of the optimization: the depth of the bioturbated layer (z bio ) decreased to 1 cm; the optimized Q 10 factor for bioturbation (θ bio ) moved to the lower limit of the Q 10 range (2 < θ < 3); and the non-local mixing coefficient (α 0 ) was reduced to a small value, essentially removing the influence of non-local mixing from the system.In addition to the reduction in bioturbation, permanent burial of ODUs does not occur in the optimized model (PB = 0, Table 2).Conversely, the optimized Q 10 factors for the remineralization of the slow (θ r1 ) and fast (θ r2 ) decaying pools of organic matter are at their upper limits indicating a strong dependence of remineralization on temperature (Table 2).For denitrification, the optimized value for the inhibition effect of NO − 3 (k dnf ) is low compared to the original parameter, whereas the inhibition effect of O 2 Table 3. Cost F (p), calculated using Eq. ( 3), for each variable type at stations Z02 and Z03.Simulations were run with the parameter set from Soetaert et al. (1996a;original) and with the optimized parameter set (baseline).Additional optimizations were carried out for each station independently (site-specific), for each station using sediment-water fluxes only (site-specific, fluxes only), and including POM depositional flux in the optimization (site-specific, fluxes only, +F POM ).(kin dnf ) is high (Table 2).The inhibition effect of O 2 on nitrification (k nit ) and of NO − 3 (kin anox ) and O 2 (kin odu ) on anaerobic remineralization is small in comparison to the original parameters.The maximum rate of nitrification (Nit) is significantly higher than in the original parameter set (Table 2).

Optimization
We examined the sources of model-data discrepancies by sequentially releasing part of the constraints on the parameter optimization (Fig. 2, Table 3).Optimizing stations Z02 and Z03 separately improves the total cost by decreasing the cost associated with NH + 4 and NO − 3 fluxes (Table 3), in particular for NO − 3 at station Z02 (Fig. 3, Table 3).Removing the constraint of sediment NH + 4 profiles from the optimization improves the total cost further (Table 3).This is due, in part, to the absence of NH + 4 profiles from the cost calculation, but also to somewhat improved sediment-water fluxes (Fig. 2).The best agreement between simulated and observed sediment-water fluxes is achieved by including POM depositional fluxes as an additional parameter to optimize (Fig. 3, Table 3).In this case POM deposition is increased in June (×2 and ×1.3 at stations Z02 and Z03, respectively) and reduced in spring (×0.5 and ×0.25 at stations Z02 and Z03, respectively) and fall (×0.5 at station Z03), and the cost as- sociated with NO − 3 and NH + 4 fluxes decreases significantly (Table 3).However, when NH + 4 profiles are not included in the cost calculation, there is a large deviation between observed and modeled sediment NH + 4 concentrations (not included in the cost).The root mean square error for the sediment profiles increases from 87.59 mmol N m −2 d −1 for the baseline case to 174.45 mmol N m −2 d −1 (site-specific, flux only) and 111.86 mmol N m −2 d −1 (site-specific, flux only +F POM ).Since the parameter set with all constraints best represents sediment-water fluxes and NH + 4 sediment concentrations throughout the Louisiana Shelf, it is used subsequently to parameterize sediment-water fluxes and is referred to as baseline.
For most of the parameter set, the optimized model is insensitive to parameter variation (Fig. 5).The most sensitive process in the diagenetic model is the remineralization of the fast decaying organic matter pool, since the optimized model is sensitive to all the associated parameters, namely the remineralization of the fast decaying organic matter pool (R 2 (T )), the base temperature (T b ) and the Q 10 factor for fast decaying organic matter (θ r1 ) in the Q 10 relationship.The optimized model is also sensitive to the variation in POM deposition rates at station Z03 (F POM 3 x ), mainly in June.Variation in deposition rates at station Z02, however, does not influence the overall cost.The sensitivity to parameters or model forcing related to organic matter is not surprising given the high magnitude and large temporal and spatial variations in POM deposition in this region.Nonetheless, it highlights the overall uncertainty in the optimized model due to the lack of observations on depositional flux.The difference in sensitivities to the depositional flux at stations Z02 and Z03 can be explained by the magnitude of the total cost, which is higher at station Z03 (Table 3).The cost at station Z02 is sensitive to the POM deposition rate (e.g., > 300 % increase in April), but since the cost at station Z03 is much higher, the effect on the total cost is small.The uncertainty associated with POM deposition rates is then larger at station Z03.To a lesser extent, the optimized model is sensitive to the bioturbation diffusivity (Dbio 0 ) and to the maximum rate of nitrification (Nit).The cost is largest for NO − 3 flux (Table 3), which indicates that the optimization has more difficulty fitting the observations for this flux.The sensitivity of the optimized value for nitrification rate, which influence NO − 3 flux, is therefore higher.

Meta-modeling parameterization
A meta-model of sediment-water fluxes was derived using simulations with the optimized diagenetic model, as described in Sect.2.2.3.The coefficients of the meta-model parameterizations for O 2 , NH + 4 and NO − 3 sediment-water fluxes and the range of bottom water conditions used for the parameterization are presented in Table 4.Each parameterization is able to reproduce the sediment-water fluxes simulated with the diagenetic model (Fig. 6).The spatially resolved correlation coefficients are above 0.8 for most of the Louisiana Shelf for O 2 and NH + 4 fluxes and above 0.6 for NO − 3 fluxes (Fig. 6).The parameterization fails to retrieve the simulated fluxes in some limited areas near the offshore limit of the shelf.Bottom water conditions for depths greater than 50 m were not included in the meta-modeling parameterization, which explains why the meta-model does not perform well in a few limited areas along the 50 m isobath.
Overall, the main contributors to the meta-model are temperature, salinity and O 2 (Table 4).The average contribution of POM deposition is low (Table 4, Fig. 7).The time dependency between POM deposition and sediment-water fluxes is implicit in the meta-model, and therefore instant POM deposition is not a good predictor of sediment-water fluxes.Temperature is the largest contributor for all fluxes (Table 4) and is associated with the seasonal variation in sediment-water fluxes.Salinity is not included in the diagenetic model, but is a significant contributor in the meta-model because it is associated with the spatial variation in sediment-water fluxes on the Louisiana Shelf.Bottom water O 2 has a growing effect on NH + 4 and NO − 3 flux under hypoxic conditions (Table 4, Fig. 6).When bottom water O 2 is low, NH + 4 flux increases with decreasing O 2 .More deposited particulate organic N is thus returned to the water column as NH + 4 .O 2 concentration controls both the direction and intensity of NO − 3 flux in the meta-model.With oxygenated bottom waters, NO − 3 flux depends on bottom NO − 3 concentration due to NO − 3 diffu- ) and NO − 3 flux (F NO −

3
).The form of the relationship is given in Eq. ( 4).For each flux, the average contribution of each input variable is indicated as well as the dominant direction of its effect.A positive effect promotes a weaker flux into the sediment or a larger flux to the water column (depending on the direction of the flux), whereas a negative effect leads to a larger sink into the sediment or a weaker flux to the water column.± indicates that the effect's direction varies as a function of the variable.The contributions were calculated from standardized coefficients.Bold values indicate variables contributing > 10 % on average.sion across the sediment-water interface.NO − 3 flux is into the sediment when the bottom water NO − 3 concentration is high and out of the sediment when the bottom water NO − 3 concentration is low.When bottom waters are hypoxic, NO − 3 flux is oriented into the sediment, which then becomes a sink for water column NO − 3 (Fig. 7).
By using simulated bottom water conditions from our biogeochemical circulation model as input for the meta-model, we can assess the spatial and temporal variability in parameterized sediment-water fluxes over the Louisiana Shelf (see Figs. 8 and 9).Sediment-water fluxes were computed from the meta-model in mid-August 2009 (Fig. 8) and throughout 2009 at stations Z02 and Z03 (Fig. 9).Bottom water conditions are presented in Fig. 3 of parameterized O 2 and NH + 4 fluxes are somewhat similar (Fig. 8), with large fluxes near Atchafalaya Bay and the Mississippi River delta where POM deposition is high in late spring (> 5 mmol N m −2 d −1 , Fig. 3).Patches of moderate to high NH + 4 flux (1-4 mmol N m −2 d −1 ) occur southwest of Terrebonne Bay and further west on the shelf where bottom waters are hypoxic (Fig. 8).NO − 3 flux follows the distribution of bottom water O 2 on the shelf with flux into the sediment in hypoxic areas and flux out of the sediment elsewhere (Fig. 8).
The time series at stations Z02 and Z03 indicate high temporal variability in parameterized sediment-water fluxes in summer (Fig. 9) that are driven by rapid changes in bottom water conditions (Fig. 3).The difference in the magnitude of O 2 flux is large between the two stations and coincides with the distinct POM deposition rate at the two stations in spring and early summer (Fig. 9).This time-dependent effect is implicit in the meta-model.A similar pattern occurs for NH + 4 flux at station Z02 (Fig. 9).The annual peak in NH + 4 flux occurs under hypoxic conditions.In late summer and fall, transient hypoxic conditions at station Z03 result in enhanced NH + 4 flux to the water column.The direction and magnitude of NO − 3 fluxes closely follows the O 2 concentration in bottom water.Hypoxic conditions starting in early July at station Z02 result in a switch from efflux of NO − 3 from the sediment to influx of NO − 3 into the sediment (Fig. 9).As for NH + 4 , rapid reversal in NO − 3 flux direction in late summer and fall at station Z03 is associated with changes between oxic and hypoxic conditions.

Comparison with other parameterizations
Here we explore the differences between the meta-models and the three sediment-water flux parameterizations we used previously in our ROMS models for the Louisiana Shelf, i.e., IR, which assumes instant remineralization of deposited POM, and H&D and M&L, which are functions of bottom temperature and O 2 concentration only.In contrast to the H&D and M&L parameterizations, O 2 flux has a relatively weak sensitivity to bottom water O 2 concentrations in the meta-model (Fig. 10).O 2 flux decreases at low bottom water O 2 concentration but does not stop in anoxic conditions, as is the case for H&D and M&L.In the model, at low O 2 , ODUs become the dominant O 2 sink (due to ODU oxidation in the water column), and therefore the O 2 sink can be significant despite the lack of O 2 in bottom waters.Similar to the IR parameterization, O 2 flux increases with PON depositional flux, but this effect is much weaker in the meta-model (Fig. 10).
The NH + 4 flux parameterized with the meta-model falls within the range of the H&D and M&L parameterizations when O 2 is available (O 2 > 50 mmol O 2 m −3 , Fig. 11).However, the meta-model differs significantly from H&D and M&L in hypoxic conditions; NH + 4 flux increases with decreasing O 2 , opposite to the H&D and M&L parameterizations.As for O 2 flux, the increase in NH + 4 flux with PON deposition is weaker than in the IR parameterization (Fig. 11).
In the meta-model, the NH + 4 flux is larger than in IR under hypoxic conditions and low PON deposition, and lower than in IR at high deposition.
Sediment-water fluxes were calculated by applying the meta-models to output from the biogeochemical circulation model and are compared to those parameterized with the H&D parameterization (Fig. 12).O 2 fluxes are larger in the meta-model in the areas of hypoxia near the Mississippi and Atchafalaya river mouths and on the mid-shelf (see Fig. 7).O 2 fluxes are smaller in the meta-model in other regions, especially on the western Louisiana Shelf, where bottom water salinity and O 2 concentrations are elevated.NH + 4 flux is also much higher in the meta-model in regions where hypoxia occurs (Fig. 12).In the other areas, NH + 4 flux is slightly lower in the meta-model.

Discussion
The meta-model procedure for parameterizing sedimentwater fluxes requires a diagenetic model that realistically represents sediment processes.In order to obtain such a realistic diagenetic model for the Louisiana Shelf, we optimized a modified version of Soetaert et al.'s model (1996a), which captures the main temporal variations in sediment biogeochemistry, sediment NH + 4 concentration and sediment-water fluxes at the two sampling locations on the eastern and western Louisiana Shelf.An issue with the optimization of large parameter sets in diagenetic models is the poor identifiability of some parameters that results in a large uncertainty in their value (Soetaert et al., 1998).This caveat in our optimization approach would not be alleviated by using a different type of optimization.Several methods have been proposed to estimate parameter identifiability and uncertainty (Soetaert et al., 1998;Soetaert and Petzoldt, 2010;Fennel et al., 2001).However, a more complete set of observations would be necessary.The available observations were also not sufficient to allow running of the diagenetic model in a time-dependent mode, and therefore the optimization was carried out with constant forcing conditions.To evaluate the effect of parameter variations (i.e., uncertainty) on the model results, we carried out a sensitivity analysis on the optimized model.A key driver of diagenetic processes is POM deposition, and the remineralization of the labile deposited POM is the most sensitive parameter in the model.Observations of POM deposition were not available, and using average rates of POM deposition from a biogeochemical model, as we have done here, is an additional source of uncertainty.This is demonstrated by the improved agreement between simulated and observed sediment-water fluxes when including POM deposition in the optimization.
Some of the discrepancies between model and observations can also be attributed to the imposition of a single parameter set.For example, sediment porosity and bioturbation are interdependent (Mulsow et al., 1998) and influence sediment-water fluxes (Aller, 1982).They are known to vary spatially on the Louisiana Shelf (Lehrter et al., 2012;Briggs et al., 2014), which is not represented in the optimized parameter set.This limitation could be resolved by introducing spatially dependent bioturbation and porosity coefficients; however, a much larger spatially resolved data set would be necessary to obtain these dependencies.Another limitation is the observed deep gradient in some of the NH + 4 profiles (e.g., in April), whereas the diagenetic model imposes a no-gradient boundary condition at depth.Some mismatch between model and observations may also be gener- ated by missing processes in the diagenetic model.As in earlier studies of the Louisiana Shelf (Morse and Eldridge, 2007;Eldridge and Morse, 2008), the diagenetic model does not represent DNRA and anammox.Although DNRA can be an important contributor to the N cycle under severe hypoxia (Dale et al., 2013), there is a poor understanding of the importance of DNRA on the Louisiana Shelf due to the lack of observations (Dagg et al., 2007).High porewater sulfide concentrations near the sediment-water interface are not reported for sediments of the Louisiana Shelf (Lin and Morse, 1991;Morse and Eldridge, 2007), which tend to minimize the importance of DNRA.However, the large NH + 4 porewater concentrations observed at station Z02 in September (Fig. 2) could be explained by the occurrence of DNRA.Anammox may also be a sink for bottom water NH + 4 on the Louisiana Shelf (Lin et al., 2011). McCarthy et al. (2015) found that anammox may represent, at times, up to 30 % of denitrification (including anammox) in some locations of the Louisiana Shelf.As a result, NH + 4 flux to the water column may be overestimated by the diagenetic model, and in the parameterization, under low bottom O 2 conditions.
Overall, despite some discrepancies with observations primarily due to uncertainty about POM deposition, diagenetic processes are represented reasonably well in the optimized model.Therefore, we deemed the optimized model to be an appropriate framework for representing the main diagenetic processes on the Louisiana Shelf.Further development of the diagenetic model may include explicit anaerobic reactions, including DNRA and anammox.However, this is beyond the scope of this work.
Comparing optimized parameters to the original parameter set used by Soetaert et al. (1996a) is informative about sediment biogeochemistry on the Louisiana Shelf .The optimization minimized the influence of bioturbation, likely a reflection of the negative impact of hypoxia on sediment biota (Diaz and Rosenberg, 1995;Middelburg and Levin, 2009).This result is also consistent with the dominance of bacteria over invertebrates in the sediment community as observed by Rowe et al. (2002).The small O 2 and NO − 3 inhibition parameters for anaerobic remineralization emphasize the importance of anaerobic processes in the area (Morse and Berner, 1995).This is consistent with observations for Mississippi River plume sediments that suggest a substantial production of reduced substances under low O 2 conditions throughout the Louisiana Shelf (Rowe et al., 2002;Lehrter et al., 2012) and reflects the important role of ODU in the O 2 flux metamodel.The small optimized value for NO − 3 limitation of denitrification indicates that direct denitrification is an important process on the Louisiana Shelf when low O 2 limits coupled nitrification-denitrification (Nunnally et al., 2013).Direct denitrification occurs when NO − 3 is available in bottom waters and tends to increase with increasing NO − 3 concentration (Fennel et al., 2009).The small optimized value of O 2 inhibition on nitrification and the relatively high maximum rate of nitrification compared to the original parameter values are also indications that sediment nitrification is an important process on the Louisiana Shelf, contributing to O 2 consumption in the sediment.This result is also consistent with earlier observations (Lehrter et al., 2012).
We added temperature dependence of remineralization to the original model from Soetaert et al. (1996a).Model results were very sensitive to changes in the remineralization of the fast decaying organic matter pool (R 2 (T )).The optimum temperature of remineralization (T b ), the remineralization at base temperature (R T b 2 ) and the Q 10 parameter for the fast decaying organic matter pool (θ 2 ) all influence R 2 (T ), and therefore model results are very sensitive to variations in these parameter values.
The meta-model reproduced the results from the optimized diagenetic model remarkably well, suggesting that it is possible to use such parameterizations in place of a full, vertically resolved diagenetic model to prescribe sediment-water boundary conditions in biogeochemical circulation models.Previous meta-model parameterizations of diagenetic rates (Middelburg et al., 1996;Soetaert et al., 2000;Gypens et al., 2008) and perturbation response experiments (Rabouille et al., 2001) had similar success.The present method is somewhat different because the goal is to parameterize sedimentwater exchanges directly as a function of bottom water conditions.This simplified parameterization method does not require an additional, vertically integrated sediment layer to track deposited POM as in the method proposed by Soetaert et al. (2000).Although the meta-model is not mass conservative, violation of mass conservation should be minor if the meta-model is used for the system and within the range of conditions that were used for its development.The result-ing meta-model exhibits realistic dynamics such as the increase in sediment-water fluxes in summer due to warmer temperature and the time delay between POM deposition and remineralization, the decrease in coupled nitrificationdenitrification at low bottom O 2 concentrations and the prominent role of reduced substances (represented by the ODU pool) as an O 2 sink in suboxic conditions.
Perhaps a key difference to other sediment-water parameterizations is the importance of ODU at low O 2 , which results in a relatively weak relationship between O 2 flux and bottom O 2 concentration in hypoxic conditions, and the occurrence of O 2 flux in anoxic conditions; in the meta-model, ODU is the dominant source of O 2 consumption in hypoxic conditions and at high temperature (i.e., in summer), independently of bottom O 2 concentration.Previous parameterizations of sediment-water O 2 flux on the Louisiana Shelf considered only SOC, and therefore O 2 flux decreased toward zero with decreasing bottom O 2 in the hypoxic range (with a zero intercept for anoxic conditions).However, Lehrter et al. (2012) found an increase in the DIC / O 2 flux ratio with bottom O 2 depletion that they attributed to anaerobic metabolism, i.e., the production of reduced chemical species that accumulate in the sediment, diffuse back and reoxidize in the water column when O 2 becomes available.Justić and Wang (2014) considered the effect of reduced chemical species on biological oxygen demand in their hypoxia model.It represents a significant O 2 sink in bottom waters and needs to be accounted for in the sediment-water O 2 flux parameterization.The O 2 flux meta-model combines SOC and ODU fluxes and is therefore a more realistic representation of O 2 consumption at the sediment-water interface.This formulation assumes instant ODU oxidation in the water column, even in anoxic conditions, whereas oxidation occurs in oxygenated waters only.The time delay between ODU flux and oxidization is therefore missing in the meta-model, but is accounted for if the coupled biogeochemical-circulation model carries an O 2 debt in anoxic conditions, as is the case in the models of Fennel et al. (2009Fennel et al. ( , 2013) ) and Laurent and Fennel (2014).
The meta-model simulates both the O 2 dependence of coupled nitrification-denitrification and direct denitrification, which are also key differences to simple parameterizations of sediment-water fluxes in biogeochemical models.The inhibition of coupled nitrification-denitrification at low O 2 stimulates eutrophication and therefore represents a positive feedback of hypoxia, as observed in Chesapeake Bay and other eutrophic systems (Kemp et al., 1990) and estimated for the global coastal ocean (Rabouille et al., 2001).It is essential to represent this feedback in high N/low O 2 systems such as the Louisiana Shelf.In the NO − 3 meta-model, the inhibition of coupled nitrification-denitrification in hypoxic conditions is partly compensated for by the increase in direct denitrification in areas where NO − 3 is available in bottom waters, which results in a nitrate flux to the sediment.On the Louisiana Shelf, this is the case in areas near the Mississippi-Atchafalaya River source, especially in the shallow area near Atchafalaya Bay.The parameterized nitrate uptake by the sediment agrees with observations from the Louisiana Shelf (Gardner et al., 1993;Nunnally et al., 2013).Nunnally et al. (2013) suggest a limited coupling between nitrification and denitrification in the Louisiana Shelf hypoxic zone.Nonetheless, the magnitude of this NO − 3 sink remains much smaller than the NH + 4 flux to the water column, and therefore the overall effect of low bottom O 2 is an enrichment of N in the water column, i.e., a positive feedback on eutrophication.
The meta-model method can be easily implemented in biogeochemical circulation models.However, the method should be applied only on regional scales because different types of bacterial, meio-or macro-faunal communities with various levels of bioturbation are associated with distinct types of substrate, porosity and POM quality and quantity affect POM recycling, and thus influence the rates of sediment diagenetic processes locally (Herman et al., 1999).In other words, diagenetic models are region-specific.

Summary and conclusions
Benthic-pelagic coupling in biogeochemical circulation models is usually implemented through simple parameterizations or with a diagenetic model.These methods are either too simplistic or computationally very costly.Soetaert et al. (2000) proposed an intermediate method to improve the efficiency of benthic-pelagic coupling in biogeochemicalcirculation models.Here we presented a simplified version computing a meta-model of sediment-water fluxes for use in a regional biogeochemical model through optimization of a diagenetic model.The method results in a realistic and computationally efficient representation of sediment-water fluxes.Applied to the Louisiana Shelf, the method provides insight into the sediment biogeochemistry of the region, such as the importance of anaerobic processes and reduced substances, the limited level of bioturbation, the occurrence of direct denitrification and the inhibition of coupled nitrification-denitrification under hypoxic conditions.The meta-models represent these Louisiana Shelf processes, resulting in more realistic, nonlinear interactions between POM deposition, bottom water concentrations and sedimentwater fluxes, in particular under hypoxic conditions.A potential limitation of the method is the need for local observations to optimize the diagenetic model.

Figure 1 .
Figure 1.Map of the Louisiana Shelf showing the location of sample collection sites Z02 and Z03.

Figure 2 .
Figure 2. Model-data comparison of sediment-water fluxes (top row) and NH + 4 profiles (bottom row) for sites Z02 and Z03.Simulations use the optimized parameter set (baseline).

Figure 3 .
Figure 3. Spatial (top) and temporal (bottom) POM depositional flux and bottom water O 2 , NH + 4 and NO − 3 concentrations in the biogeochemical circulation model.The upper panels represent a snapshot of bottom water conditions on 15 August 2009 and the lower panels time series at stations Z02and Z03.This data set is used to force the diagenetic model in the meta-modeling procedure (Sect.2.2), to compute spatial fluxes with the meta-model (Fig.8) and to compare the meta-model and H&D parameterizations (Fig.12).

Figure 4 .
Figure 4. Model-data comparison of sediment-water fluxes at stations Z02 and Z03 for several different optimization schemes (baseline includes all constraints).

Figure 5 .
Figure 5. Sensitivity of model results to parameter variation.

Figure 6 .
Figure 6.Correlation coefficients between time-dependent diagenetic model simulations and the parameterized fluxes for each location on the Louisiana Shelf.

Figure 7 .
Figure 7. Influence of selected contributors to O 2 , NH + 4 and NO − 3 fluxes.Negative fluxes (blue shades) are into the sediment and positive fluxes (orange shades) are out of the sediment.

Figure 11 .
Figure 11.NH + 4 flux in the meta-model compared with that from the IR, H&D and M&L parameterizations.NH + 4 flux is represented as a function of (left) bottom O 2 concentration and (right) PON depositional flux.The grey area and the black line in the left panel correspond to the variation in O 2 flux when 1 < F POM < 10 mmol N m −2 d −1 and F POM = 5 mmol N m −2 d −1 , respectively.The black lines on the right indicate the O 2 flux at bottom O 2 concentrations of 0, 50 and 250 mmol O 2 m −3 .

Figure 12 .
Figure 12.Difference between parameterized oxygen (top panel) and ammonium (bottom panel) fluxes and fluxes simulated with the H&D parameterization on 15 August 2009.

Table 1 .
Bottom water conditions at stationsZ02 and Z03 in 2006.These data are used as forcing conditions during the optimization of the diagenetic model.POM deposition flux (F POM ) was not measured; F POM monthly climatologies were calculated for stations Z02 and Z03 from a multiyear simulation with a biogeochemical circulation model (see Sect. 2.3).

Table 2 .
Soetaert et al. (1996a)ters.The 20 parameters that were optimized are indicated with a + sign.The original values are fromSoetaert et al. (1996a); an asterisk indicates values that are identical in the optimized parameter set.

Laurent et al.: Parameterization of sediment-water fluxes on the Louisiana Shelf each
time and location where observations were available.During the optimization, the model was forced with observed bottom water conditions, namely salinity, temperature, NH + 4 , NO − 3 , and O 2 (Table www.biogeosciences.net/13/77/2016/Biogeosciences, 13, 77-94, 2016 A.