Bottom RedOx Model ( BROM , v . 1 . 1 ) : a coupled benthic-1 pelagic model for simulation of water and sediment 2 biogeochemistry

Interactions between seawater and benthic systems play an important role in global biogeochemical 15 cycling. Benthic fluxes of some chemical elements (e.g. C, N, P, O, Si, Fe, Mn, S) alter the redox state and 16 marine carbonate system (i.e. pH and carbonate saturation state), which in turn modulate the functioning of 17 benthic and pelagic ecosystems. The redox state of the near bottom layer in many regions can change with time, 18 responding to the supply of organic matter, physical regime and coastal discharge. We developed a model 19 (BROM) to represent key biogeochemical processes in the water and sediments and to simulate changes 20 occurring in the bottom boundary layer. BROM consists of a transport module (BROM-transport) and several 21 biogeochemical modules that are fully compatible with the Framework for the Aquatic Biogeochemical Models, 22 allowing independent coupling to hydrophysical models in 1D, 2D or 3D. We demonstrate that BROM is 23 capable of simulating the seasonality in production and mineralization of organic matter as well as the mixing 24 that leads to variations in redox conditions. BROM can be used for analyzing and interpreting data on sediment25 water exchange, and for simulating the consequences of forcings such as climate change, external nutrient 26 loading, ocean acidification, carbon storage leakage, and point-source metal pollution. 27


454
E. V. Yakushev et al.: Bottom RedOx Model (BROM v.1.1) The redox state and oxygenation of near-bottom water varies due to the transport of oxidized and reduced species across the SWI and biogeochemical processes occurring in the sediments (Cooper and Morse, 1996;Jorgensen et al., 1990;Roden and Tuttle, 1992;Sell and Morse, 2006). The sediments generally consume oxygen due to the deposition of labile OM and the presence of reduced forms of chemical elements. Their capacity to exchange oxygen with the pelagic layer is limited, as near-bottom water is usually characterized by low water velocity and reduced mixing in the vicinity of the SWI (Glud, 2008). In some cases, a high benthic oxygen demand (BOD) associated with local OM mineralization and low mixing rates can cause anoxia in the bottom water. This may lead to death, migration, or changed behavior of the benthic macro-and meiofaunal organisms responsible for bioturbation and bioirrigation (Blackwelder et al., 1996;Sen Gupta et al., 1996;Morse and Eldridge, 2007), which in turn can greatly slow down the transport of solid and dissolved species inside the sediments and therefore the rates of oxidative reactions. Under such conditions, sedimentary sulfides can build up, and dissolution of carbonate minerals may come to a halt (Morse and Eldridge, 2007). When oxic conditions return, there can be an "oxygen debt" of reduced species in the water column  which may buffer and delay reoxygenation of the sediments (Morse and Eldridge, 2007).
In areas experiencing seasonal hypoxia/anoxia, the processes taking place in the water column and in the sediments are tightly coupled to each other, as well as to the fluxes and exchanges of organic matter over a range of timescales. An accurate understanding of physical, chemical, and biological processes driving changes in redox conditions is needed to predict the distribution of hypoxia/anoxia in a given environment. This "benthic-pelagic coupling" broadly encompasses the fluxes of OM to the sediments and the return fluxes of inorganic nutrients to the water column. Variations in supply, dynamics, and reactivity of OM affect benthic communities (Pearson and Rosenberg, 1978), sediment and porewater geochemistry (Berner, 1980), and nutrient and oxygen fluxes at the SWI (Boudreau, 1997).
Many previous studies have demonstrated the capability of sophisticated reactive transport codes for integrated modeling of biogeochemical cycles in sediments (Boudreau, 1996;Van Cappellen and Wang, 1996;Couture et al., 2010;Jourabchi et al., 2005;Katsev et al., 2006Katsev et al., , 2007Paraska et al., 2014;Soetaert et al., 1996). The water column redox interface was also specifically targeted in the models of Konovalov et al. (2006) and Yakushev et al. (2006Yakushev et al. ( , 2007. However, the process of integrating such models with pelagic biogeochemical models to produce benthic-pelagic coupled models has only begun in recent years. As of the year 2000, benthic-pelagic coupling was either neglected or crudely approximated in many pelagic biogeochemical and early diagenetic models (Soetaert et al., 2000). One of the first fully coupled physical-pelagic-benthic bio-geochemical modes was developed for the Goban Spur shelf break area to examine the impact of in situ atmospheric conditions on ecosystem dynamics, to understand biogeochemical distributions in the water column and the sediments, and to derive a nitrogen budget for the area. This model was most suited for testing the impact of short-term physical forcing on the ecosystem (Soetaert et al., 2001).
Later, several coupled benthic-pelagic models were produced with an emphasis on studying eutrophication (Cerco et al., 2006;Fennel et al., 2011;Soetaert and Middelburg, 2009) or hypoxia in various locations including Tokyo Bay (Sohma et al., 2008), the Baltic Sea (Reed et al., 2011), the North Sea oyster grounds (Meire et al., 2013), and the Southern Bight (Lancelot et al., 2005). Another model was created to investigate early diagenesis of silica in the Scheldt estuary, with benthic-pelagic coupling only of silica (Arndt and Regnier, 2007).
By coupling two quite sophisticated models ECOHAM1 and C.CANDI, a 3-D model for the North Sea was created where pelagic model output was used to force a benthic biogeochemical module (Luff and Moll, 2004). Another physical-biological model for the North Sea, PROWQM, is more complex than ECOHAM1 and has been coupled to a benthic module to simulate seasonal changes of chlorophyll, nutrients, and oxygen at the PROVESS north site, southeast of the Shetland Islands (Lee et al., 2002). Brigolin et al. (2011) developed a spatially explicit model for the northwestern Adriatic coastal zone by coupling a 1-D transient early diagenesis model with a 2-D reaction-transport pelagic biogeochemical model. Currently, the most known and established coupled model is ERSEM -the European Regional Seas Ecosystem Model, which was initially developed as a coastal ecosystem model for the North Sea and which has evolved into a generic tool for ecosystem simulations from shelf seas to the global ocean (Butenschön et al., 2016).
The BROM model described herein is a fully coupled benthic-pelagic model with a special focus on deoxygenation and redox biogeochemistry in the sediments and benthic boundary layer (BBL). The BBL is "the part of the marine environment that is directly influenced by the presence of the interface between the bed and its overlying water" (Dade et al., 2001). Physical scientists tend to prefer the term "bottom boundary layer", but this is largely synonymous with the BBL (Thorpe, 2005). Within BROM, the term BBL refers to the lower parts of the fluid bottom boundary layer where bottom friction strongly inhibits current speed and vertical mixing, hence including the viscous and logarithmic sublayers up to at most a few meters above the sediment. This calm-water layer plays a critical role in mediating the interaction of the water column and sediment biogeochemistry and in determining, e.g., near-bottom oxygen levels, yet it remains poorly resolved in most physical circulation models. For BROM, we have developed an accompanying offline transport module (BROM-transport) that uses output from hydrodynamic water column models but solves the transport-2 BROM description BROM consists of two modules, BROM-biogeochemistry and BROM-transport. BROM-biogeochemistry builds on ROLM (RedOx Layer Model), a model constructed to simulate basic biogeochemical structure of the water column oxic/anoxic interface in the Black Sea, Baltic Sea, and Norwegian fjords (He et al., 2012;Stanev et al., 2014;Yakushev et al., 2009Yakushev et al., , 2006Yakushev et al., , 2007. In BROM-biogeochemistry, we extended the list of modeled compounds and processes (Fig. 1). BROM considers interconnected transformations of species of N, P, Si, C, O, S, Mn, and Fe, and resolves OM in nitrogen currency. OM dynamics include parameterizations of OM production (via photosynthesis and chemosynthesis) and OM decay via oxic mineralization, denitrification, metal reduction, sulfate reduction, and methanogenesis. To provide a detailed representation of changing redox conditions, OM in BROM is mineralized by several different electron accep-tors and dissolved oxygen is consumed during both mineralization of OM and oxidation of various reduced compounds. Process inhibition in accordance with redox potential is parameterized by various redox-dependent switches. BROM also includes a module describing the carbonate equilibria; this allows BROM to be used to investigate acidification and impacts of changing pH and saturation states on water and sediment biogeochemistry.
The physical domain of BROM-transport spans the water column, BBL, and upper layers of the sediments in a continuous fashion. This allows for an explicit, high-resolution representation of the BBL and upper sediments, while also allowing the boundary conditions to be moved as far as possible from these foci of interest, i.e., to the air-sea interface and to deep in the sediment.
BROM is integrated into an existing modular platform (FABM) and is therefore coded as a set of reusable "LEGO brick" components, including the offline transport driver BROM-transport and modules for ecology, redox chemistry, and carbonate chemistry. This means that BROM-transport can be used with all biogeochemical modules available in FABM, including, e.g., the modules comprising ERSEM, and that BROM biogeochemical modules can be used in all other 1-D and 3-D hydrodynamic models supported by FABM (e.g., GOTM, GETM, MOM5, NEMO, FVCOM). Individual BROM modules can also be coupled to existing ecological models to expand their scope, e.g., by providing descriptions of redox and carbonate chemistry. Using the FABM framework thus facilitates the transparent and consistent setup of complex biogeochemical reaction networks for the prediction of hypoxia/anoxia while harnessing the capabilities of various hydrophysical drivers.
2.1 Biogeochemical module 2.1.1 General description BROM-biogeochemistry consists of three biogeochemical submodules: BROM_bio (ecological model), BROM_redox (redox processes), and BROM_carb (carbonate system). Interactions between modeled variables are either kinetic (e.g., OM degradation) or equilibrium processes (e.g., carbonate system equilibration) (Boudreau, 1996;Jourabchi et al., 2008;Luff et al., 2001). In general, the redox reactions are fast in comparison with the other processes and a typical model time step. Species involved in such reactions are therefore set to equilibrium concentrations using mass action laws and equilibrium constants for seawater (Millero, 1995). Total scale pH is also diagnosed at every time step, mainly as a function of dissolved inorganic carbon (DIC) and total alkalinity (Alk) which are both prognostic (state) variables.
The model has 33 state variables (Table 1), including frequently measured components such as hydrogen sulfide (H 2 S) and phosphate (PO 4 ), as well as rarely measured variables such as elemental sulfur (S 0 ), thiosulfate (S 2 O 3 ), triva-456 E. V. Yakushev et al.: Bottom RedOx Model (BROM v.  lent manganese species Mn(III), and bacteria. We acknowledge that for many of these, site-specific estimates of associated model parameters and initial/boundary conditions may be difficult or impossible to obtain, and may in practice require some crude assumptions and approximations (e.g., universal default parameter values, no-flux boundary conditions, and initial conditions from a steady annual cycle). Nevertheless, we believe that for many applications this caveat will be acceptable given the additional process resolution and realism provided by BROM for important biogeochemical processes in the BBL and sediments. The equations and parameters employed in BROM are given in Tables 2 and 3, and a flow chart is shown in Fig. 1.

Ecosystem and redox models
The overall goal of the ecosystem representation is to parameterize the key features of OM production and decomposition, following Redfield and Richards stoichiometry (Richards, 1965). We divide all the living OM (biota) into Phy (photosynthetic biota), Het (non-microbial heterotrophic biota), and four groups of "bacteria" which may be considered to include microbial fungi. These latter are Baae (aerobic chemoautotrophic bacteria), Baan (anaerobic chemoautotrophic bacteria), Bhae (aerobic heterotrophic bacteria), and Bhan (anaerobic heterotrophic bacteria). OM is produced photosynthetically by Phy and chemosynthetically by bacteria, specifically by Baae in oxic conditions and by Baan in anoxic conditions. Growth of heterotrophic bacteria is tied to mineralization of OM, favoring Bhae in oxic conditions and Bhan in anoxic conditions. Secondary production is represented by Het, which consumes phytoplankton as well as all types of bacteria and dead particulate organic matter (detritus, which is also explicitly modeled). The effect of suboxia and anoxia is parameterized by letting the mortality of aerobic organisms depend on the oxygen availability. A detailed account of processes representing the inorganic cycling of N, S, Mn, Fe, and P is given in the description of ROLM (Yakushev et al., 2007(Yakushev et al., , 2013a, while the process parameterization, chemical reactions, rates, and stoichiometric constants values are summarized in Tables 2-4.  Table 2 also describes the redox-dependent switches, nutrient limitation, and substrate consumption rates for heterotrophs. The redox-dependent switches are mostly based on hyperbolic tangent functions which improve system stability compared with discrete switches. The nutrient limitation and heterotrophic transfer functions are based on squared Monod laws for nutrient-biomass ratio, which also stabilizes the system compared with Michaelis-Menten and Ivlev formulations. Here, we describe the parameterization of carbon that was not considered in ROLM and was not described in Yakushev (2013).

Total alkalinity
Total alkalinity, A T , is a model state variable. Following the formal definition of A T (Dickson, 1992;Wolf-Gladrow et al., 2007;Zeebe and Wolf-Gladrow, 2001), the following alkalinity components were considered:  Luff et al. (2001) and Volkov (1984). The boric alkalinity A B = [B(OH) − 4 ] was estimated from total dissolved boron, which in turn was calculated from salinity.
Biogeochemical processes can lead to either increase or decrease of alkalinity, and alkalinity can be used as an indicator of specific biogeochemical processes (Soetaert et al., 2007). Organic matter production can affect alkalinity via the "nutrient-H + compensating principle" formulated by Wolf-Gladrow et al. (2007): during uptake or release of charged nutrient species, electroneutrality is maintained by consumption or production of a proton (i.e., during uptake of nitrate for photosynthesis or denitrification, or production of nitrate by nitrification).

Carbonate system
Equilibration of the carbonate system was considered as a fast process occurring within seconds (Zeebe and Wolf-Gladrow, 2001). Accordingly, the equilibrium solution was calculated at every time step using an iterative procedure. The carbonate system was described using standard approaches (Lewis and Wallace, 1998;Munhoven, 2013;Roy et al., 1993;Wanninkhof, 2014;Wolf-Gladrow et al., 2007;Zeebe and Wolf-Gladrow, 2001). The set of constants of Roy et al. (1993) was used for carbonic acid. Constants for boric, hydrofluoric, and hydrogen sulfate alkalinity were calculated according to Dickson (1992), for silicic alkalinity according to Millero (1995), for ammonia alkalinity according to Luff et al. (2001), and for hydrogen sulfide alkalinity according to Luff et al. (2001) and Volkov (1984). The ion product of Dependence of photosynthesis on P LimP = Excretion rate of phytoplankton ExcrPhy = K_phy_exc × Phy Phytoplankton mortality rate

Bacteria
Growth rate of bacteria aerobic autotrophic Growth rate of bacteria aerobic heterotrophic Rate of mortality of bacteria aerobic heterotrophic Rate of mortality of bacteria anaerobic autotrophic MortBaan = K_Baan_mrt × Baan Growth rate of bacteria anaerobic heterotrophic HetBhan = (DcPM_NOX + DcDM_NOX + DcDM_Mn + DcPM_Mn + DcDM_Fe + DcPM_Fe + DcDM_SO4 Rate of mortality of bacteria anaerobic heterotrophic Yakushev et al.: Bottom RedOx Model (BROM v.1.1) water was calculated according to Millero (1995). Total scale pH was calculated using the Newton-Raphson method with the modifications proposed in Munhoven (2013). Precipitation and dissolution of calcium carbonate were modeled following the approach of Luff et al. (2001) (Table 2).

Physical environment
BROM-biogeochemistry can be coupled online with various hydrodynamic models using FABM, but this may require extensive adaptation of the hydrodynamic model to resolve the BBL and upper sediments. We have therefore developed a simple 1-D offline transport-reaction model, BROMtransport, whose model domain spans the water column, BBL, and upper layers of the sediments, with enhanced spatial resolution in the BBL and sediments. All options and parameter values for BROM-transport are specified in a runtime input file brom.yaml. A step-by-step guide to running BROM-transport is provided in Appendix A.

BROM-transport model formulation
The time-space evolution of state variables in BROMtransport is described by a system of 1-D transport-reaction equations in Cartesian coordinates. In the water column, the dynamics are whereĈ i is the concentration in units [mmol m −3 total volume] of the ith state variable, D(z, t) is the vertical diffusivity, v i is the settling or sinking velocity, ε h (z, t) is a rate of horizontal mixing with an external concentrationĈ 0i (z, t) (or alternatively, a restoring rate to a climatological concentration), T birr(i) is a tendency due to bioirrigation (only nonzero for dissolved substances in the bottom layer of the water column; see below), and R i is the combined sources minus sinks (in this study provided by BROM-biogeochemistry, but in principle any biogeochemical model in FABM could be used). Values for D, ε h ,Ĉ 0i , and other forcings used by R i are configured at runtime through input files (see Sect. 2.2.7). Sinking velocities v i are non-zero only for particulate (non-dissolved) variables and are determined at each time step by the biogeochemical module (through FABM). BROM-biogeochemistry assumes constant sinking velocities for phytoplankton, zooplankton, bacteria, detritus, and inorganic particles (Table 3e).
In the sediments, dissolved substances or solutes obey the dynamics where ϕ is the porosity (assumed constant in time), D C is the total solute diffusivity, u is the solute burial velocity, and C i is the porewater concentration in units [mmol m −3 porewater]. Particulate substances become part of the solid matrix in the sediments. These obey where D B is the particulate (bioturbation) diffusivity, w is the particulate burial velocity, and B i is the particulate concentration in units [mmol m −3 total solids]. The porosity ϕ(z) in Eqs. (2) and (3) is prescribed as an exponential decay, following Soetaert et al. (1996): where ϕ ∞ is the deep (compacted) porosity, ϕ 0 is the sediment surface porosity, z SWI is the depth of the SWI, and δ is a decay scale defining the rate of compaction. Diffusion within the sediments is assumed to be strictly "intraphase" (Boudreau, 1997), hence the Fickian gradients in Eqs. (2) and (3) are formed using the concentration per unit volume porewater for solutes and per unit volume total solids for particulates. The total solute diffusivity D C = D m + D B , where D m is the apparent molecular/ionic diffusivity and D B is the bioturbation diffusivity due to animal movement and ingestion/excretion. The apparent molecular diffusivity D m (z) = θ −2 D 0 µ 0 µ sw is derived from the infinite-dilution molecular diffusivity D 0 (an input parameter) assuming a constant relative dynamic viscosity µ 0 µ sw (default value 0.94, cf. Boudreau, 1997, Table 4.10) and a tortuosity parameterized as θ 2 = 1−2 ln ϕ from Boudreau (1997, Eq. 4.120). The bioturbation diffusivity D B (z, t) is modeled as a Michaelis-Menten function of the dissolved oxygen concentration in the bottom layer of the water column: where D Bmax (z) is a constant over a fixed mixed layer depth in the surface sediments, then decays to zero with increasing depth, and K O2s is a half-saturation constant. The rationale for Eq. (5) is that the benthic animals that cause bioturbation require a source of oxygen at the sediment surface for respiration. Diffusion between the sediments and water column, i.e., across the SWI, raises a subtle issue in regard to particulates. Here, any diffusive flux cannot be strictly intraphase, because particulates are modeled as [mmol m −3 total solids] in the sediments but as [mmol m −3 total volume] in the water column. In BROM-transport, the bottom layer of the water column is considered a "fluff layer"; particles enter through the upper interface at their sinking velocity and leave through the SWI at the particulate burial velocity. It follows that a portion of the particulate matter in the fluff layer must be considered as settled fluff, but that portion is not predicted by the  (Yakushev et al., 2007) model. BROM-transport therefore offers two approaches. In the first approach, the bioturbation diffusivity is set to zero on the SWI, so that only solutes can diffuse across the SWI by molecular diffusion. Since the present version of BROMtransport does not parameterize resuspension through the SWI due to fluid turbulence, the SWI thus becomes a oneway street for particulate matter, whose components can only reenter the water column after dissolution. In the second approach, the bioturbation diffusivity is given by Eq. (5) on the SWI, but the bioturbation flux is interphase, mixing concen-trations in units [mmol m −3 total volume] for both solutes and particulates. This approach is appropriate if bioturbation can be assumed to exchange fluff and sediment, or if it contributes significantly to particulate resuspension. The burial velocities u and w in Eqs.
(2) and (3) can be inferred from the porosity profile under the assumptions of steady-state compaction (ϕ constant in time) and no externally impressed porewater flow (Berner, 1971(Berner, , 1980Boudreau, 1997;Meysman et al., 2005). Here, BROMtransport again offers two approaches. In the first approach,  3 × 10 −4 10 −4 -10 −2 mol g −1 yr −1 ; (Wersin, 1990;Wollast, 1990) Specific rate of MnCO 3 dissolution K_mnco3_ diss d −1 7 × 10 −4 10 −2 -10 3 yr −1 ; (Wersin, 1990;Wollast, 1990)  K_feco3_ diss d −1 7 × 10 −4 2.5 × 10 −1 -10 −2 yr −1 ; (Wersin, 1990;Wollast, 1990) Specific rate of FeCO 3 formation K_feco3_form d −1 3.4 × 10 −4 10 −6 -10 −2 mol/g yr; (Boudreau, 1996;Wersin, 1990;Wollast, 1990     the reactions of particles in the sediments are assumed to have negligible impact on the volume fraction of total solids, and the deep particulate burial velocity w ∞ in compacted sediments (where ϕ = ϕ ∞ ) is assumed to be a known constant w b∞ (an input parameter). Since compaction ceases at this (possibly infinite) depth, the solute burial velocity must here equal the particulate burial velocity (u ∞ = w b∞ ). Steady state then implies the following burial velocities (Appendix B): where D inter B is the interphase bioturbation diffusivity, nonzero only at the SWI and only if bioturbation across the SWI is enabled. In the second approach, the reactions of the modeled particulate substances in the sediments modify the total solid volume fraction, and the modeled sinking fluxes from the water column modify the flux of solid volume at the SWI. The velocities in Eqs. (6) and (7) then define background velocities (w b , u b ) due to non-modeled particulates. Assuming steady-state compaction leads to the following corrections to the background burial velocities (see Appendix B): where w = w − w b , u = u − u b , N p is the number of particulate variables, ρ i is the density of the ith particle type, v f(i) is the sinking velocity in the fluff layer,Ĉ sf(i) is the suspended particulate concentration in the fluff layer, R i is the particulate reaction term, and w ∞ is the correction to the deep particulate burial velocity, in practice approximated by the deepest value of w . Since the suspended portionĈ sf(i) is not explicitly modeled, it is approximated as the minimum of the particulate concentrations in the fluff layer and the layer immediately above. In our applications, we have found that Eqs. (8) and (9) can improve the realism of sediment organic matter distributions, mainly by increasing the burial rate following pelagic production and export events such as the spring bloom. Finally, the process of bioirrigation, whereby benthic organisms flush out their burrows with water from the sediment surface, is modeled as a non-local solute exchange (following Aller, 2001;Meile et al., 2001;Rutgers Van Der Loeff and Boudreau, 1997;Schlüter et al., 2000): where α (z) is the bioirrigation rate in oxic conditions,Ĉ f(i) is the flushing concentration of solute in the fluff layer, and the Michaelis-Menten function again accounts for the suppression of worm activity in anoxic conditions. The oxic bioirrigation rate α (z) is parameterized as an exponential decay from the sediment surface as in Schlüter et al. (2000). The total mass transfer to/from the sediment column must be balanced by a flux into/out of the fluff layer (see Eq. 1): where h f is the thickness of the fluff layer and z max is the depth of the bottom of the modeled sediment column. T birrC(i) , T birr(i) = 0 for all particulate variables.

BROM-transport numerical integration
Equations (1)-(3) are integrated numerically over a single combined grid (water column plus sediments) and using the same model time step in both water column and sediments. All concentrations are stored internally and input/output in units [mmol m −3 total volume]. Time stepping follows an operator splitting approach (Butenschön et al., 2012): concentrations are successively updated by contributions over one time step of diffusion, bioirrigation, reaction, and sedimentation, in that order. If any state variable has any "not-anumber" values at the end of the time step then the program is terminated.
Diffusive updates are calculated either by a simple forward-time central-space (FTCS) algorithm or by a semiimplicit central-space algorithm adapted from a routine in the General Ocean Turbulence Model, GOTM (Umlauf et al., 2005). Bioirrigation and reaction updates are calculated from forward Euler time steps, using FABM to compute R i , and sedimentation updates are calculated using a simple first-order upwind differencing scheme. After each update, Dirichlet boundary conditions (see below) are reimposed and all concentrations are low bounded by a minimum value (default = 10 −11 µM) to avoid negative values. Maximum diffusive and advective Courant numbers can optionally be output after every time step or when/if a not-a-number value is detected. Before starting the integration, the program calculates Courant numbers due to eddy/molecular diffusion and returns a warning message if maximum values on any given day exceed 0.5 and the FTCS option is selected.
BROM-transport also provides an option to divide the diffusion and sedimentation updates into smaller time steps related to the sources-minus-sinks time step by fixed factors, since the physical transport processes are often numerically limiting (Butenschön et al., 2012). The default time step is 0.0025 days or 216 s, which is much longer than the characteristic equilibration timescale of the CO 2 kinetics (Zeebe and Wolf-Gladrow, 2001).

BROM-transport vertical grid
The vertical grid in BROM-transport is divided into the pelagic water column, the BBL, and the sediments. The pelagic water column grid is either set as uniform with height/spacing set by the brom.yaml file (see Sect. S1 in the Supplement), or it is read from the NetCDF forcing input file (see below), with an option to decrease resolution by subsampling. In principle, the NetCDF input from the hydrodynamic model may already include a fully resolved BBL, but in practice we find this is rarely the case. BROM-transport therefore allows the user to "insert" a high-resolution BBL into the bottom of the input water column. This BBL has non-uniform grid spacing with layer thickness decreasing geometrically towards the SWI, reaching O (cm) thickness for the fluff layer, based on parameters from the brom.yaml file. For the upper sediments, the layer thickness is increased geometrically moving down from the SWI, from O (mm) thickness in the surface layer to O (cm) thickness deeper in the sediments, again based on brom.yaml parameters. The result is a full grid with non-uniform spacing and maximum resolution near the SWI. As in many ocean models (e.g., ROMS, GOTM) the vertical grid in BROM-transport is staggered: temperature, salinity, and biogeochemical concentrations are defined at layer midpoints, while diffusivities, sinking/burial velocities, and resulting transport fluxes are all defined on layer interfaces.

BROM-transport initial conditions
Initial conditions for all concentrations in Eqs. (1)-(3) can be provided by either using the initialization values defined in the fabm.yaml file (see Sect. S2 in the Supplement) as uniform initial conditions for each variable, or by providing the initial conditions for all variables at every depth in a text file with a specific format. Typically, these initial-condition text files are generated by running the model to a steady state annual cycle and saving the final values as the desired start date. Alternatively, they could be generated by interpolating/smoothing data, in which case the user should note that the input concentrations must be in units [mmol m −3 total volume].

BROM-transport boundary conditions
BROM-transport presently allows the user to choose between four different types of boundary conditions for each variable and for upper and lower boundaries: (1) no gradient at the bottom boundary (no diffusive flux) or no flux at the surface boundary, except where parameterized by the FABM biogeochemical model (i.e., for O 2 and DIC in the case of BROMbiogeochemistry); (2) a fixed constant value; (3) a fixed sinusoidal variation in time defined by amplitude, mean value, and phase parameters; or (4) an arbitrary fixed variation in time read from the input NetCDF file. All boundary condition options and parameters are set in the brom.yaml file (see Sect. S1). Note that options 2-4 are Dirichlet boundary conditions which define implicit fluxes of matter into and out of the model domain, and that all boundary concentrations should be in units [mmol m −3 total volume (water+solids)]. The default option 1 is generally the preferred choice, but the Dirichlet options can also be useful to allow a simple representation of, e.g., fluxes of nutrients into and out of the surface layer due to lateral riverine input. A possible alternative is to use the forcings' parameters for horizontal mixing (see Eq. 1) to specify horizontal exchanges or restoring terms to observed climatology (see Sect. 2.2.7).
Under option 1, and using BROM-biogeochemistry, a surface O 2 flux representing exchange with the atmosphere is parameterized as where O 2sat is the oxygen saturation as a function of temperature and salinity, according to UNESCO (1986), Sc is the Schmidt number for oxygen (Raymond et al., 2012), and k 660 is the reference gas-exchange transfer velocity, parameterized as k 660 = 0.365 u 2 + 0.46 u (Schneider et al., 2002) where u is the wind speed 10 m above the sea surface (m s −1 ). Air-sea exchange of CO 2 in BROMbiogeochemistry is parameterized using the partial pressures in water (pCO water 2 ) and air (pCO air 2 ) following the formulation and coefficients in Butenschön et al. (2016): where F wind = (0.222 u 2 + 0.333 u)(Sc/660) −0.5 is a wind parameter (Nightingale et al., 2000), u is the wind speed, and Sc is the Schmidt number for CO 2 (Raymond et al., 2012).

BROM-transport irradiance model
BROM-transport includes two simple Beer-Lambert attenuation models to calculate in situ 24 h average photosynthetically active radiation (PAR) as needed by BROMbiogeochemistry and many other biogeochemical models. The first is derived from the current ERSEM default model (Blackford et al., 2004;Butenschön et al., 2016) and models the total attenuation as where k 0 is the background attenuation of seawater, k Phy and k PON are the specific attenuations due to phytoplankton and detritus, respectively, and k s is the specific attenuation due to "other" optically active substances with concentration S (currently a constant input parameter). The second model includes attenuation due to other optically active concentrations that are modeled by BROM-biogeochemistry: where B is the total bacterial concentration (Baae + Baan + Bhae + Bhan) and PIV is the total volume fraction of modeled inorganic particles, calculated from the concentrations using input densities of each inorganic solid. The final irradiance is scaled by a constant parameter representing either the photosynthetically active fraction of the in situ irradiance or the relationship between surface PAR in water and the forcing surface irradiance (Mobley and Boss, 2012). The forcing surface irradiance Eair(t) can be read from NetCDF input or otherwise calculated using a sinusoidal function (Yakushev et al., 2013b). In addition, the surface attenuation due to ice cover can be accounted for as a simple linear function of a NetCDF input ice thickness variable hice(t).

BROM-transport input forcings
BROM-transport requires forcing inputs at least for temperature, salinity, and vertical diffusivity at all depths in the pelagic water column and for each day of the simulation. These may be provided from an input subroutine that creates simple, hypothetical profiles, or from text/NetCDF files containing data from interpolations of measurements or hydrodynamic model output. Forcing time series of surface irradiance and ice thickness may also be read as NetCDF input. BROM-transport then uses these inputs in combination with parameters set in the runtime input file brom.yaml (see Sect. S1) to solve the transport-reaction equations on a "full" vertical grid including pelagic water column, BBL, and sediment subgrids. In order to run, BROM-transport must extend the input pelagic (temperature, salinity, diffusivity) forcings over the full grid. Temperature and salinity in the BBL and sediments are set as uniform and equal to the values at the bottom of the input pelagic water column for each day. The vertical diffusivity needs a more careful treatment, as it is the main defining characteristic of the pelagic vs. BBL vs. sediment environments. Within the water column, the total vertical diffusivity D = D m + D e for solutes and D = D e for particulates, where D m is a constant molecular diffusivity at infinite dilution, and D e is the eddy diffusivity read from the input file for the pelagic water column. For the BBL, D e can be defined as "dynamic", in which case it is linearly interpolated for each day between the deepest input forcing value above the SWI and zero at a depth h DBL above the SWI, where h DBL is the diffusive boundary layer (DBL) thickness (default value 0.5 mm). This option is likely appropriate for shallow-water applications where D e may be strongly time dependent within the user-defined BBL (default thickness 0.5 m). Alternatively, a static, fixed profile D eBBL (z) may be more appropriate for deep-water BBLs, where time dependence may be weak and deepest values from hydrodynamic models may be relatively far above the SWI. In this case, BROM-transport offers two options for D eBBL (z): (1) a constant value, dropping to zero in the DBL, or (2) a linear variation between a fixed value at the top of the BBL and zero at the top of the DBL. Option 1 defines a simplest-possible assumption, while option 2 corresponds to the assumption of a log layer for the current speed (e.g., Boudreau and Jorgensen, 2001;Holtappels and Lorke, 2011). Eddy diffusivity is strictly zero in the DBL, on the SWI, and within the sediments. Diffusivity in the sediments is due to molecular diffusion and bioturbation and is parameterized as described in Sect. 2.2.1.
Optional forcings for BROM-transport include 24 h average surface irradiance Eair(t), which is often supplied by hydrodynamic models (e.g., ROMS), a surface ice thickness forcing hice(t), and depth-time arrays of horizontal mixing rates ε h (z,t) and horizontal mixing concentrationsĈ 0i (z,t) (see Eq. 1). Horizontal mixing rates within the inserted BBL and sediments are set to zero. Note that these horizontal mixing forcings can also be used to define relaxation or restoring fluxes to climatological values within the pelagic water column, which may in some cases provide a valid means of accounting for horizontal flux divergence effects that are missing in the 1-D model.

Model setup
A North Sea hydrodynamic scenario was used to demonstrate the ability of BROM to reproduce the biogeochemical mechanisms of oxic/anoxic transformations. Complete lists of the model options and parameter values used are given in Sect. S1 (brom.yaml input file for BROM-transport) and Sect. S2 (fabm.yaml input file for BROM-biogeochemistry).
The BROM-transport water column extended from 0 to 110 m, with a pelagic spatial resolution of 1 m inherited from the GOTM hydrodynamic model used to provide forcings. A high-resolution BBL was inserted from 109.5 to 110 m, with layer thickness decreasing from approximately 25 to 3 cm in the fluff layer. Sediment grid points were added to cover the upper 10 cm of sediments with layer thickness increasing from 0.5 mm in the surface layer to 1 cm at depth. This choice of grid does not explicitly resolve the DBL (default thickness 0.5 mm) but the main DBL function of limiting solute exchange between the BBL and sediments is largely fulfilled by the fluff layer (thickness 3 cm) and upper sediment layer (thickness 0.5 mm). The model time step for BROMtransport was set to 0.0025 days (216 s).
Upper boundary conditions included sinusoidal, timevarying Dirichlet boundary conditions for nitrate, phosphate, and silicate, implying net influxes and outfluxes of surface nutrients, as well as the default parameterized air-sea fluxes of O 2 and DIC (see Sect. S1). Lower boundary conditions assumed (by default) zero diffusive flux for all reduced components (i.e., hydrogen sulfide, solid-phase concentrations of metal sulfides and carbonates, silicon, and OM). The simulation therefore focuses on the consequences of the supply of fresh OM as a main reducer in both water column and sediments.
The pelagic water column was forced by output from a GOTM hydrodynamical simulation for temperature, salinity, and vertical diffusivity (taken from the salinity diffusivity) and surface irradiance calculated using the sinusoidal option. We aimed for a solution representative of "present day" and therefore treated the GOTM forcing as representative of a "normal year". BROM-transport was spun up from vertically homogeneous initial conditions for 100 model years with repeated-year forcings and boundary conditions. After this time, a quasi-stationary solution with seasonally forced oscillations of the biogeochemical variables had been reached.
The results of these calculations were written to an output file in NetCDF format, including the daily vertical distributions of model state variables, diagnostic rates of biogeochemical transformations, and fluxes associated with diffusion and sedimentation. This output can be visualized by any NetCDF-compatible software.

Results
The model simulated the periodic replacement of oxic with anoxic conditions in the BBL following seasonal mixing and OM production. The simulation demonstrates the character-istic features of biogeochemical profiles in the water column, BBL, and upper sediments, as well as their variability under changing redox conditions (Figs. 2-4).
During intensive mixing conditions in winter, the water column is well oxygenated and the oxic/anoxic interface is located at a depth of several centimeters in the sediments (Figs. 2, 3). In summer, just after the spring bloom, an enrichment of the sediment surface with fresh OM and a restricted oxygen supply leads to the consumption of O 2 by OM mineralization and close to suboxic conditions (Fig. 2). The second bloom in autumn leads to a further decrease of oxygen concentrations to complete depletion. There is a concomitant increase in reduced forms of N, Mn, and Fe, and finally of hydrogen sulfide in the bottom water (Figs. 2, 4). The redox interface thus moves from the sediment to the BBL. Figure 5 shows the rate of OM mineralization with a variety of electron acceptors. Oxygen is consumed during OM mineralization in summer and autumn and, after its complete depletion, denitrification dominates, with both nitrate reduction and nitrite reduction playing significant roles. The rate 474 E. V. Yakushev et al.: Bottom RedOx Model (BROM v.  of mineralization of OM with Mn and Fe oxides is small, but as these processes prevent mineralization with sulfate, they cause a lag of a few days between the depletion of oxygen and the appearance of hydrogen sulfide in the water column (Figs. 2, 5). The amount of labile degradable OM is relatively small and mineralization with sulfate completely removes the remaining OM, thus preventing methanogenesis (Fig. 5).
The seasonal variability of the sediment-water fluxes clearly demonstrates the appearance in the bottom water of reduced forms of N, Mn, Fe, and phosphate (Fig. 6).
Generally, the concentrations, vertical distributions, and benthic-pelagic fluxes of the parameters considered in the model are reasonable and are within observed ranges for the North Sea (Queirós et al., 2014) and some other regions with temporary bottom anoxia (Almroth et al., 2009;McCarthy et al., 2008;Morse and Eldridge, 2007;Pakhomova et al., 2007;Queirós et al., 2014;Yu et al., 2015).

Conclusion and future work
This paper presents a description of BROM, a fully coupled pelagic-benthic model that provides an integrated framework to study the biogeochemistry of a water column and upper sediments. BROM simulates changes in redox conditions and their impact on the distributions of a wide range of biogeochemical variables. In particular, BROM provides a detailed description of the fate and availability of dissolved oxygen and hydrogen sulfide: the former essential for macroscopic marine life, the latter highly toxic to it. BROM can therefore provide valuable information to ecological studies, particularly in the context of multistressor impacts. The model suggests that the timing of hydrogen sulfide release into the pelagic is linked to the dynamics of several electron acceptors that are themselves of limited interest for biogeochemical and ecological purposes, and that are therefore rarely included in models. The ability of BROM to simulate and forecast H 2 S toxicity is in fact the direct result of its inclusion of several of these rarely modeled chemical compounds (e.g., Mn(IV), Fe(III)).
This paper was not devoted to a detailed validation of BROM with in situ data; we plan to explore this in future work. A qualitative analysis of the model results (Sect. 3) suggests that the model can produce realistic distributions and fluxes of key biogeochemical variables during periodic changes in redox conditions.
In summary, we present a new benthic-pelagic biogeochemical model (BROM) that combines a relatively simple pelagic ecosystem model with a detailed biogeochemical model of the coupled cycles of N, P, Si, C, O, S, Mn, and Fe in the water column, benthic boundary layer, and sediments, with a focus on oxygen and redox state. BROM should be of interest for the study a range of environmental applications in addition to hypoxia, such as benthic nutrient recycling, redox biogeochemistry, eutrophication, industrial pollution from trace elements, organic loading, and ocean acidification.

Code availability
The model as presented consists of two components. The first is a set of biogeochemical modules (brom/redox, brom/bio, brom/carb, brom/eqconst), available as part of the official FABM distribution (http://fabm.net); BROM-specific files are located in subdirectory src/models/niva/brom). The second is a hydrophysical driver (BROM-transport) that provides the 1-D vertical context and resolves transport; this is available separately from https://github.com/e-yakushev/ brom-git.git. When combined, the 1-D BROM model as presented is obtained.
Both FABM and BROM-transport are coded in objectoriented Fortran 2003, have a build system based on CMake (https://cmake.org), and use YAML files (http://yaml.org) for runtime configuration. The code is platform independent and only requires a Fortran 2003-capable compiler, e.g., gfortran 4.7 or higher, or the Intel Fortran compiler version 12.1 or higher. BROM-transport includes facilities for producing results as NetCDF files, which can be read by a variety of software on different platforms.
As BROM's biogeochemical modules are built on FABM, they can be used from a wide range of 1-D and 3-D hydrodynamic models, including GOTM, GETM, ROMS, MOM, NEMO, and FVCOM (a ROMS-FABM coupler has been developed by P. Wallhead; NEMO-FABM and FVCOM-FABM couplers have been developed by the Plymouth Marine Laboratory; contact J. Bruggeman for information).
Results shown in this paper were produced with BROMtransport tag v1.1 and the BROM-biogeochemistry code in FABM tag v0.95.3, available from the above repositories. The simulation was run using the netCDF/.yaml input files found in the data folder of the BROM-transport repository. However, we envisage BROM to be further developed in a backward compatible manner, and encourage users to adopt the latest version of the code.
Step-by-step instructions for running BROM are found in Appendix A. Both FABM and BROM-transport are distributed under the GNU General Public License (http://www.gnu.org/licenses/). As a component of FABM, BROM-biogeochemistry is licensed under the same conditions as FABM. 2. Preparation of input files consists of the model reading two .yaml files with the model parameters (fabm.yaml and brom.yaml), as well as a NetCDF or text file with the hydrophysical forcing data. Optionally, the biogeochemical initial conditions can be read from a text file start.dat; this may be a file written by a previous simulation (the final model state is written to a file named finish.dat at the end of every simulation).
i. The brom.yaml (see Sect. S1) file specifies the values of transport model parameters as well as various option switches and input/output file and variable names. Text comments provide guidance and references for setting parameter values. If using NetCDF input, the user should pay careful attention to the NetCDF input parameters and names, ensuring that this information is consistent with the input NetCDF file. The selected-year parameter year must refer to a year that is covered by the input forcing data.
ii. The fabm.yaml (see Sect. S2) file specifies the values of biogeochemical model parameters, default initial values for state variables, and the coupling of FABM modules. Text comments provide annotation and references.
iii. The nns_annual.nc (in the example) file contains input forcing data that may be derived from observations or hydrodynamical model output (GOTM in our demonstration). It can be replaced by a text (.dat) file if this is the format of the hydrodynamical model output.
iv. The start.dat is the text file with initial values for model state variables at every depth. This file may be created by renaming the output of a previous simulation (finish.dat is the state on 1 January of the last modeled year).
3. Output files are NetCDF and headed text files generated automatically by the model during the simulation. Output files can be readily imported into various software packages for visualization and further analysis. Certain output files (Vertical_grid.dat and Hydrophysics.dat) are generated early in the simulation and should be checked by the user to ensure that the model grid and hydrophysical forcings are set up as intended.
i. Vertical_grid.dat is the text file with model layer indices, midpoint depths, increments between midpoint depths, and thicknesses.
iii. The finish.dat is the text file with the state variables for the 1 January of the last modeled year. It can be used for visualization or as initial conditions for further calculations.
iv. The output_NNday.dat is the optional text file with the state variables and diagnostic variables for day NN to make plots of vertical distributions (e.g., Fig. 3) v. BROM_out.nc is the NetCDF file with daily profiles of state variables, rates of biogeochemical transformations, and vertical fluxes.
4. For visualization of NetCDF output files, any software with NetCDF input can be used. In the example, we used PyNcView for 2-D and BROM_pictures for 1-D (available at https://github.com/BottomRedoxModel/ brom_pictures).

Appendix B: Derivation of burial velocities
The conservation equations for liquid and total solid volume fractions in the sediments can be written as where D inter B is the interphase bioturbation diffusivity (possibly non-zero only at the SWI), ρ i is the density of the ith particulate substance, and R i is the corresponding reaction term. Equations (B1) and (B2) assume that the densities of liquid and total solid are both constant, and they retain the net contributions of reactive terms although these are often considered negligible, e.g., (Boudreau, 1997;Meysman et al., 2005). Summing Eqs. (B1) and (B2) and integrating over depth gives a useful and quite general relationship: where U (t) is only a function of time. If we now assume no externally impressed porewater flow, it follows that at some (possibly infinite) depth where compaction ceases ( ∂ϕ ∂z = 0, ϕ = ϕ ∞ ), the solute burial velocity u must here equal the particulate burial velocity w, hence u ∞ = w ∞ . Equation (B3) becomes Now assuming steady state compaction ( ∂ϕ ∂t = 0), Eq. (B2) can be integrated from the SWI to a depth z within the sediments: To determine the first term on the right-hand side of Eq. (B5), we assume that the total solid volume flux across the SWI is equal to the total solid volume flux from the sinking of suspended particulate matter in the fluff layer: (1 − ϕ SWI )w SWI + D inter where F b defines a constant background solid volume flux due to non-modeled particles, v f(i) is the sinking velocity in the fluff layer, andĈ sf(i) is the suspended particulate concentration in the fluff layer. Substituting into Eq. (B5), we have Since D inter B ∂ϕ ∂z is zero at depth, the constant surface flux term is given by F b = (1 − ϕ ∞ ) w b∞ , where both ϕ ∞ and w b∞ are input parameters. Hence, we have Equation (6) directly follows from Eq. (B8) by neglecting the modeled settling flux and reaction terms, then Eq. (7) follows by application of Eq. (B4). Equations (8) and (9) follow by considering the additional particulate burial velocity due to fluxes and reactions (from the last term in Eq. B8) and applying Eq. (B4) to obtain the additional solute burial velocity.