Qual2K: a modeling Framework for Simulating River and Stream Water Quality (Version 11)

SOD/Nutrient Flux Model

Sediment nutrient fluxes and sediment oxygen demand (SOD) are based on a model developed by Di Toro (Di Toro et al. 1991, Di Toro and Fitzpatrick. 1993, Di Toro 2001). The present version also benefited from James Martin’s (Mississippi State University, personal communication) efforts to incorporate the Di Toro approach into EPA’s WASP modeling framework.

A schematic of the model is depicted in Figure . As can be seen, the approach allows oxygen and nutrient sediment-water fluxes to be computed based on the downward flux of particulate organic matter from the overlying water. The sediments are divided into 2 layers: a thin (? 1 mm) surface aerobic layer underlain by a thicker (10 cm) lower anaerobic layer. Organic carbon, nitrogen and phosphorus are delivered to the anaerobic sediments via the settling of particulate organic matter (i.e., phytoplankton and detritus). There they are transformed by mineralization reactions into dissolved methane, ammonium and inorganic phosphorus. These constituents are then transported to the aerobic layer where some of the methane and ammonium are oxidized. The flux of oxygen from the water required for these oxidations is the sediment oxygen demand. The following sections provide details on how the model computes this SOD along with the sediment-water fluxes of carbon, nitrogen and phosphorus that are also generated in the process.

Figure Schematic of SOD-nutrient flux model of the sediments.

1. Diagenesis

As summarized in Figure , the first step in the computation involves determining how much of the downward flux of particulate organic matter (POM) is converted into soluble reactive forms in the anaerobic sediments. This process is referred to as diagenesis. First the total downward flux of carbon, nitrogen and prosphorus are computed as the sums of the fluxes of settling phytoplankton and organic matter delivered to the sediments from the water column

()

where J_POC = the downward flux of POC [gC m^?2 d^?1], r_ca = the ratio of carbon to chlorophyll a [gC/mgA], v_a = phytoplankton settling velocity [m/d], a_p = phytoplankton concentration [mgA/m³], r_cd = the ratio of carbon to dry weight [gC/gD], v_dt = detritus settling velocity [m/d], m_o = detritus concentration [gD/m³], J_PON = the downward flux of PON [mgN m^?2 d^?1], q_Np = the phytoplankton nitrogen cell quota [mgN mgA^?1], v_on = organic nitrogen settling velocity [m/d], n_o = organic nitrogen concentration [mgN m^–3], J_POP = the downward flux of POP [mgP m^?2 d^?1], q_Pp = the phytoplankton phosphorus cell quota [mgP mgA^?1], v_op = organic phosphorus settling velocity [m/d], and p_o = organic phosphorus concentration [mgP m^–3] .



Figure Representation of how settling particulate matter is transformed into fluxes of dissolved carbon (J_C) , nitrogen (J_N) and phosphorus (J_P) in the anaerobic sediments.
Note that for convenience, we express the particulate organic carbon (POC) as oxygen equivalents using the stoichiometric coefficient r_oc. Each of the nutrient fluxes is further broken down into three reactive fractions: labile (G1), slowly reacting (G2) and non-reacting (G3).
These fluxes are then entered into mass balances to compute the concentration of each fraction in the anaerobic layer. For example, for labile POC, a mass balance is written as
()
where H₂ = the thickness of the anaerobic layer [m], POC_2,G1 = the concentration of the labile fraction of POC in the anaerobic layer [gO₂/m³], J_POC,G1 = the flux of labile POC delivered to the anaerobic layer [gO₂/m²/d], k_POC,G1 = the mineralization rate of labile POC [d^?1], ?_POC,G1 = temperature correction factor for labile POC mineralization [dimensionless], and w₂ = the burial velocity [m/d]. At steady state, Eq. 197 can be solved for
()
The flux of labile dissolved carbon, J_C,G1 [gO₂/m²/d], can then be computed as
()
In a similar fashion, a mass balance can be written and solved for the slowly reacting dissolved organic carbon. This result is then added to Eq. 199 to arrive at the total flux of dissolved carbon generated in the anaerobic sediments.
()
Similar equations are developed to compute the diagenesis fluxes of nitrogen, J_N [gN/m²/d], and phosphorus J_P [gP/m²/d].

2. Ammonium

Based on the mechanisms depicted in Figure , mass balances can be written for total ammonium in the aerobic layer and the anaerobic layers,

()
()
where H₁ = the thickness of the aerobic layer [m], NH_4,1 and NH_4,2 = the concentration of total ammonium in the aerobic layer and the anaerobic layers, respectively [gN/m³], n_a = the ammonium concentration in the overlying water [mgN/m³], ?_NH4,1 = the reaction velocity for nitrification in the aerobic sediments [m/d], ?_NH4 = temperature correction factor for nitrification [dimensionless], K_NH4 = ammonium half-saturation constant [gN/m³], o = the dissolved oxygen concentration in the overlying water [gO₂/m³], and K_NH4,O2 = oxygen half-saturation constant [mgO₂/L], and J_N = the diagenesis flux of ammonium [gN/m²/d].
The fraction of ammonium in dissolved (f_dai) and particulate (f_pai) form are computed as
()
()
where m_i = the solids concentration in layer i [gD/m³], and ?_ai = the partition coefficient for ammonium in layer i [m³/gD].
The mass transfer coefficient for particle mixing due to bioturbation between the layers, ?₁₂ [m/d], is computed as
()
where D_p = diffusion coefficient for bioturbation [m²/d], ?_Dp = temperature coefficient [dimensionless], POC_R = reference G1 concentration for bioturbation [gC/m³] and K_M,Dp = oxygen half-saturation constant for bioturbation [gO₂/m³].
The mass transfer coefficient for pore water diffusion between the layers, K_L12 [m/d], is computed as,
()
where D_d = pore water diffusion coefficient [m²/d], and ?_Dd = temperature coefficient [dimensionless].
The mass transfer coefficient between the water and the aerobic sediments, s [m/d], is computed as
()
where SOD = the sediment oxygen demand [gO₂/m²/d].
At steady state, Eqs. 201 and 202 are two simultaneous nonlinear algebraic equations. The equations can be linearized by assuming that the NH_4,1 term in the Monod term for nitrification is constant. The simultaneous linear equations can then be solved for NH_4,1 and NH_4,2. The flux of ammonium to the overlying water can then be computed as
()

3. Nitrate

Mass balances for nitrate can be written for the aerobic and anaerobic layers as

()
()
where NO_3,1 and NO_3,2 = the concentration of nitrate in the aerobic layer and the anaerobic layers, respectively [gN/m³], n_n = the nitrate concentration in the overlying water [mgN/m³], ?_NO3,1 and ?_NO3,2 = the reaction velocities for denitrification in the aerobic and anaerobic sediments, respectively [m/d], and ?_NO3 = temperature correction factor for denitrification [dimensionless].
In the same fashion as for Eqs. 201 and 202, Eqs. 209 and 210 can be linearized and solved for NO_3,1 and NO_3,2. The flux of nitrate to the overlying water can then be computed as
()
Denitrification requires a carbon source as represented by the following chemical equation,
()
The carbon requirement (expressed in oxygen equivalents per nitrogen) can therefore be computed as
()
Therefore, the oxygen equivalents consumed during denitrification, J_O2,dn [gO₂/m²/d], can be computed as
()

4. Methane

The dissolved carbon generated by diagenesis is converted to methane in the anaerobic sediments. Because methane is relatively insoluble, its saturation can be exceeded and methane gas produced. As a consequence, rather than write a mass balance for methane in the anaerobic layer, an analytical model developed by Di Toro et al. (1990) is used to determine the steady-state flux of dissolved methane corrected for gas loss delivered to the aerobic sediments.

First, the carbon diagenesis flux is corrected for the oxygen equivalents consumed during denitrification,
()
where J_CH4,T = the carbon diagenesis flux corrected for denitrification [gO₂/m²/d]. In other words, this is the total anaerobic methane production flux expressed in oxygen equivalents.
If J_CH4,T is sufficiently large (? 2K_L12C_s), methane gas will form. In such cases, the flux can be corrected for the gas loss,
()
where J_CH4,d = the flux of dissolved methane (expressed in oxygen equivalents) that is generated in the anaerobic sediments and delivered to the aerobic sediments [gO₂/m²/d], C_s = the saturation concentration of methane expressed in oxygen equivalents [mgO₂/L]. If J_CH4,T < 2K_L12C_s, then no gas forms and
()
The methane saturation concentration is computed as
()
where H = water depth [m] and T = water temperature [^oC].
A methane mass balance can then be written for the aerobic layer as
()
where CH_4,1 = methane concentration in the aerobic layer [gO₂/m³], c_f = fast CBOD in the overlying water [gO₂/m³], ?_CH4,1 = the reaction velocity for methane oxidation in the aerobic sediments [m/d], and ?_CH4 = temperature correction factor [dimensionless]. At steady, state, this balance can be solved for
()
The flux of methane to the overlying water, J_CH4 [gO₂/m²/d], can then be computed as
()

5. SOD

The SOD [gO₂/m²/d] is equal to the sum of the oxygen consumed in methane oxidation and nitrification,

()
where CSOD = the amount of oxygen demand generated by methane oxidation [gO₂/m²/d] and NSOD = the amount of oxygen demand generated by nitrification [gO₂/m²/d]. These are computed as
()
()
where r_on = the ratio of oxygen to nitrogen consumed during nitrification [= 4.57 gO₂/gN].

6. Inorganic Phosphorus

Mass balances can be written total inorganic phosphorus in the aerobic layer and the anaerobic layers as

()
()
where PO_4,1 and PO_4,2 = the concentration of total inorganic phosphorus in the aerobic layer and the anaerobic layers, respectively [gP/m³], p_i = the inorganic phosphorus in the overlying water [mgP/m³], and J_P = the diagenesis flux of phosphorus [gP/m²/d].
The fraction of phosphorus in dissolved (f_dpi) and particulate (f_ppi) form are computed as
()
()
where ?_pi = the partition coefficient for inorganic phosphorus in layer i [m³/gD].
The partition coefficient in the anaerobic layer is set to an input value. For the aerobic layer, if the oxygen concentration in the overlying water column exceeds a critical concentration, o_crit [gO₂/m³], then the partition coefficient is increased to represent the sorption of phosphorus onto iron oxyhydroxides as in
()
where ??_PO4,1 is a factor that increases the aerobic layer partition coefficient relative to the anaerobic coefficient.
If the oxygen concentration falls below o_crit then the partition coefficient is decreased smoothly until it reaches the anaerobic value at zero oxygen,
()
Equations 225 and 226 can be solved for PO_4,1 and PO_4,2. The flux of phosphorus to the overlying water can then be computed as
()

7. Solution Scheme

Although the foregoing sequence of equations can be solved, a single computation will not yield a correct result because of the interdependence of the equations. For example, the surface mass transfer coefficient s depends on SOD. The SOD in turn depends on the ammonium and methane concentrations which themselves are computed via mass balances that depend on s. Hence, an iterative technique must be used. The procedure used in QUAL2K is

1. 
   Determine the diagenesis fluxes: J_C, J_N and J_P.
   

2. 
   Start with an initial estimate of SOD,
   

()
where r_on’ = the ratio of oxygen to nitrogen consumed for total conversion of ammonium to nitrogen gas via nitrification/denitrification [= 1.714 gO₂/gN]. This ratio accounts for the carbon utilized for denitrification.

3. 
   Compute s using
   

()

4. 
   Solve for ammonium, nitrate and methane, and compute the CSOD and NSOD.
   

5. 
   Make a revised estimate of SOD using the following weighted average
   

()

6. 
   Check convergence by calculating an approximate relative error
   

()

7. 
   If ?_a is greater than a prespecified stopping criterion ?_s then set SOD_init = SOD and return to step 2.
   

8. 
   If convergence is adequate (?_a ? ?_s), then compute the inorganic phosphorus concentrations.
   

9. 
   Compute the ammonium, nitrate, methane and phosphate fluxes.
   

8. Supplementary Fluxes

Because of the presence of organic matter deposited prior to the summer steady-state period (e.g., during spring runoff), it is possible that the downward flux of particulate organic matter is insufficient to generate the observed SOD. In such cases, a supplementary SOD can be prescribed,

()
where SOD_t = the total sediment oxygen demand [gO₂/m²/d], and SOD_s = the supplementary SOD [gO₂/m²/d]. In addition, prescribed ammonia and methane fluxes can be used to supplement the computed fluxes.