Process-Based Model to Describe Treatment of Nitrate-Rich Groundwater Using Emulsified Oil

Abstract

Permeable reactive barriers (PRBs) are being considered for treatment where the discharge of nitrate plumes contributes to eutrophication (e.g., Cape Cod, MA). PRBs enhance denitrification through the addition of carbon-based amendments such as the injection of emulsified vegetable oil (EVO). The use of EVO to stimulate denitrification foregrounds aspects of carbon utilization, dosing, longevity, and secondary effects in ways that differ from the application of EVO at hazardous waste sites. The overall objective of this study was to develop and evaluate a process-based modeling approach for simulating denitrification stimulated and supported by EVO. A series of one-dimensional column experiments assessed emulsion retention, production of soluble substrate, and utilization of carbon for nitrate reduction. Retention of 5.5 g dispersed phase emulsion resulted in sustained reduction of nitrate (∼43 mg/day) at ∼2 m/day porewater velocity. Biokinetic processes underlying the model are based on the two-step denitrification description of the Activated Sludge Model (ASM) No. 3. Biokinetic processes were integrated within the flow and transport simulator COMSOL to simulate the column experiment. The model capitalizes on the biokinetic parameters available in the ASM literature to limit the number of site-specific fits of model parameters. Simulation results demonstrate how this approach can result in reasonable predictions, although model performance was enhanced by fitting two parameters—yield coefficients for nitrate and nitrite. Comparisons with existing biokinetic transport models that were similarly fit to the column data suggest that the better overall descriptions of the column data using the process-based model stem from a more robust handling of spatial and temporal distribution of biomass. Sensitivity analyses highlight the importance of accurately describing the transformation of complex carbon into soluble substrate, and the subsequent utilization of that substrate. This research establishes a foundation for exploring implications of carbon processing on dosing, longevity, and effectiveness in denitrifying PRBs.

Introduction

Nitrate in groundwater is a substantial water quality concern in areas where septic systems and fertilizer applications are common. Recently, communities have begun exploring the enhancement of denitrification within the subsurface near locations where groundwater discharges contribute to eutrophication. To enhance denitrification within the subsurface, carbon sources are added to stimulate indigenous bacteria that are capable to perform each step in the sequential reduction NO₃⁻ → NO₂⁻ → NO → N₂O → N₂ (Supplementary Fig. S1). When treating nitrate-rich groundwater, carbon amendments can be added via injection, mixing, or trenching (Huno et al., 2018). Permeable reactive barriers (PRBs) established with injectable carbon in the form of oil-in-water emulsions offer potential to address nonpoint sources of nitrate (Parece et al., 2019). In fact, this type of treatment is being tested on Cape Cod, MA to reduce nitrate discharges where groundwater flows into embayments.

The oil-in-water emulsions used in these applications comprise vegetable oil, typically soybean oil, emulsifiers, and sometimes also contain lactate (Dombrowski et al., 2017). Demonstration and testing of emulsified vegetable oil (EVO) to establish PRBs on Cape Cod have identified the need for better capabilities in determining the dosing and longevity of the emulsified oil.

Denitrification in porous media has been examined and modeled at the batch (Calderer et al., 2010; Mekala and Nambi, 2016), column (Clement et al., 1997; Glass and Silverstein, 1998; Lu et al., 2020; Rodríguez-Escales et al., 2016; Sun et al., 1998; Wang et al., 2020), and field (Killingstad et al., 2002; Smith et al., 2001) scales. Monod-type expressions are commonly employed to describe biological processes related to denitrification (Li et al., 2017). Models describing denitrification are often characterized as one-step (NO₃⁻ to N₂) or two-step (NO₃⁻ to NO₂⁻ to N₂), although explicit consideration of each step in the reduction pathway (Supplementary Fig. S1) may provide benefit where accumulation or release of NO and N₂O is of concern (Kaelin et al., 2009; Pan et al., 2015).

Clement et al. (1997) describe a one-step denitrification process supported by the utilization of a simple carbon (acetate) soluble substrate. The model includes biomass growth but neglects inhibition, endogenous respiration, O₂, and N₂. Killingstad et al. (2002) used the two-step denitrification model to simulate nitrate reduction using formate. The increased model complexity of the two-step model affords consideration of the role nitrite has in inhibiting other processes (Du et al., 2016). Killingstad et al. (2002) include nitrite inhibition but exclude endogenous respiration and O₂. However, both of these models, along with the others noted above, have less applicability to PRBs established using EVO as each lacks the capability for modeling the biotransformation of complex carbon (e.g., soybean oil) into soluble substrate (e.g., acetate), and each relies on fitting model coefficients to cite specific data.

Our overarching objective is to develop a process-based modeling approach that can meaningfully contribute to understanding amendment longevity when EVO is used to create a PRB for the purposes of denitrification. PRB performance may also be influenced by nitrite and dissolved oxygen, and the contribution of endogenous respiration to nitrate reduction may be important over longer periods of biomass growth and decay. Our previous work has demonstrated the ability of the Activated Sludge Model (ASM) No. 3 framework (Kaelin et al., 2009) to describe microcosm experiments conducted to explore denitrification stimulated with EVO (Gonsalez et al., 2023). The ASM framework was developed and refined in the context of wastewater treatment to standardize descriptions of biological processes (Brun et al., 2002; Gujer et al., 1999; Henze et al., 1999; Hiatt and Grady, 2008; Kappeler and Gujer, 1992).

ASMs use Monod kinetics to describe a single biomass population growing under ideal mixing conditions in aerobic, anoxic, and autotrophic environments. Saccani et al. (2020) used the ASM framework to simulate denitrification in a biofilter. Our work aims to couple the process-based description of the biokinetics offered by the ASM framework with transport in porous media to describe the stimulation and sustention of denitrification via the introduction of emulsified oil in the subsurface environment. Such a model may be subsequently combined with emulsion transport models to provide a fuller understanding of the interplay between emulsion retention, geometry of the treatment zone, amendment longevity, and fate of the added carbon.

Materials and Methods

Materials

Saturated aquifer material was collected from Orleans, MA. The material was refrigerated at 4°C until use and homogenized before column packing. The Orleans material is a medium sand having 10% wt. >4,250 μm, 50% wt. >1,500 μm, and 90% wt. >287.5 μm. The sand has 3,600 mg/kg organic matter (loss on ignition to 550°C) and contains 3,100 mg/kg total iron and 28.8 mg/kg total manganese. Stream water from Naaman Creek, Claymont, Delaware, was used to represent a natural water containing bacteria capable of denitrification (water quality provided in Supplementary Table S1). SRS^®-NR (60% soybean oil, 10.9% surfactants, and 4.0% sodium lactate) was provided by Terra Systems, Inc.

An emulsion comprising 25 g/L dispersed phase (DP) content was created by diluting SRS-NR with stream water, which represents ∼10% of the dosage typically used in Cape Cod demonstration tests (Dombrowski et al., 2017). The diluted EVO results in a measured total chemical oxygen demand (COD) concentration of 15,550 mg/L.

Column methods

Two column experiments (15 cm length × 5.1 cm diameter, clear polyvinyl chloride) were conducted at 22°C ± 2°C using Orleans aquifer material (bulk density: 1.83 g/cm³; porosity: 0.31). Column A served as the treatment column that received both emulsion and nitrate as described below. Column B served as a control column for the delivery of the emulsion. Columns were treated identically during simultaneous emulsion injections until such time as Column B was sectioned to obtain the retention profile of the oil and biomass.

Two hundred twenty milliliters of emulsion was introduced to Columns A and B for a duration of 23 min. At the conclusion of the emulsion pulse, the influent solution was switched to the background solution of 500 mg/L calcium chloride for 1.6 pore volumes. Effluent samples collected from both columns during this phase of the experiment were analyzed for turbidity (to estimate DP emulsified oil concentration), specific conductivity, COD, and total organic carbon.

After emulsion delivery, Column B was extruded and sampled at 1-cm intervals to quantify DP emulsified oil retention using the loss on ignition approach described in Muller et al. (2018). Column A did not receive flow for a period of 1.7 days before the introduction of ∼30 mg-N/L nitrate solution at ∼1 mL/min (Cole Parmer model 7553-80 pump; Darcy velocity: 0.7 m/day) using water purged with nitrogen gas. Influent concentrations were found to vary over the course of the 33-day experiment (20 ≤ nitrate mg-N/L ≤ 30; 14 ≤ mg-COD/L ≤ 41). Measured flow rates became less stable around day 16 due to technical issues with the pump (Supplementary Fig. S2). At the conclusion of the experiment, the column was extruded and sampled in 1-cm intervals for emulsion and biomass.

Analytical methods

While the use of COD precludes determination of the breakdown pathways of the complex carbon substrate, it is well aligned with the field and modeling approaches used to monitor denitrification (Henze et al., 1987). COD was quantified in the effluent samples using HACH method TNT 822. Effluent samples were submitted to Eurofins Lancaster Laboratories for analysis of nitrate and nitrite (ion chromatography method 300). Turbidity was quantified using Hach 890 colorimeter method 8277. Specific conductivity and pH were quantified using Oakton Con6 and Thermo 290A meters. Dissolved oxygen (DO) was quantified using YSI ProODO. Biomass quantification was performed by Gene-Trac^® testing using quantitative polymerase chain reaction (qPCR) at SiREM Laboratory (Guelph, Canada).

Total bacteria were quantified by the 16S rRNA gene count, present in all bacteria (copies per gram of wet soil), and denitrification target genes (i.e., nirS, nirK, and narG) quantified the denitrifying community. Biomass concentrations were calculated by assuming that each bacterium possesses only one 16S rRNA gene (Jenkins et al., 2012).

Process-based modeling

Columns were conceptualized as a one-dimensional, axisymmetric transport [Eq. (1)] of the following components (i): biomass, emulsion, oxygen, soluble substrate, nitrate, nitrite, and nitrogen gas (Table 1).

Table 1.

Transport Model Components

Component	Units	Initial condition	Influent condition	Transport condition	Reaction ID
Nitrate (NO₃⁻)	mg-N/L	0.8	29.9	Mobile	R_NO3−
Nitrite (NO₂⁻)	mg-N/L	0	0	Mobile	R_NO2−
Soluble substrate (S)	mg-COD/L	39.2	39.2	Mobile	R_S
Oxygen (DO)	mg/L	0	0	Mobile	R_DO
Nitrogen gas (N₂)	mg-N/L	0	0	Mobile	R_N2
Dispersed phase emulsion (DP)	g-COD/kg	Spatially variable^a	0	Immobile	R_DP
Biomass^b (X)	mg-COD/kg	0.33	0	Immobile	R_X

Emulsified vegetable oil concentration as a function of the fractional position in the column (x/L) is $530 + 319987.5 \cdot e^{(- 11.5 \frac{x}{L})}$ .

Biomass is assumed to be the fraction of heterotrophic bacteria contributing to denitrification and was established using quantitative polymerase chain reaction results for the aquifer material used in the column experiment.

COD, chemical oxygen demand.

\frac{\partial C_{i}}{\partial t} = D_{H} \frac{\partial^{2} C_{i}}{\partial x^{2}} - v_{x} \frac{\partial C_{i}}{\partial x} + R_{i}

(1)

where C_i in Equation (1) represents the aqueous concentration of component i, [M∙L⁻³], t is time [T], D_H is hydrodynamic dispersion [L²∙T], x is axial position [L], v_x is pore water velocity [L∙T], and R_i is the reaction term for component i. Modeling transport using Equation (1) relies on several simplifying assumptions: porosity and aqueous-phase saturation do not vary in space in time; solid-phase interactions are negligible; and reactions occur in the aqueous phase. Once retained, the DP of well-designed oil-in-water emulsions is immobile and has limited influence on hydraulic conductivity (Coulibaly and Borden, 2004; Crocker et al., 2007; Muller et al., 2018).

Thus, we elected to model emulsion as an immobile aqueous-phase component rather than increase model complexity with a description of the exchange between the solid and aqueous phases. The initial condition for the emulsion within the column was based on the measured retention profile from Column B (see distribution in the notes of Table 1). The retained component of the emulsion was the DP, which constitutes the droplets of vegetable oil that make up the oil-in-water emulsion. The retained DP emulsion represents a source of slow-release soluble substrate (S_S) produced via hydrolysis (P1 in Table 2, and reaction terms R_DP and R_s in Supplementary Table S4). Biomass was similarly conceptualized as an immobile aqueous-phase component (Killingstad et al., 2002) and growth was limited by a threshold term, 1-X/X_max (Kildsgaard and Engesgaard, 2001).

Table 2.

Model Process Rates

ID	Process	Rate equation
Process-based model in this work
P1	Hydrolysis	$k_{D P} \cdot \frac{D P ∕ X}{K_{D P} + D P ∕ X} \cdot (1 - \frac{D P_{m i n}}{D P}) \cdot X$
P2	Aerobic growth on S_S	$μ_{H} \cdot \frac{D O}{K_{D O} + D O} \cdot \frac{S_{S}}{K_{S} + S_{S}} \cdot (1 - \frac{X}{X_{m a x}}) \cdot X$
P3	Anoxic growth on S_S using nitrate	$μ_{H} \cdot η \cdot \frac{K_{D O}}{K_{D O} + D O} \cdot \frac{S_{S}}{K_{S} + S_{S}} \cdot \frac{N O_{3}^{-}}{K_{N O_{3}^{-}} + N O_{3}^{-}} \cdot (1 - \frac{X}{X_{m a x}}) \cdot X$
P4	Anoxic growth on S_S using nitrite	$μ_{H} \cdot η \cdot \frac{K_{D O}}{K_{D O} + D O} \cdot \frac{S_{S}}{K_{S} + S_{S}} \cdot \frac{N O_{2}^{-}}{K_{N O_{2}^{-}} + N O_{2}^{-}} \cdot (1 - \frac{X}{X_{m a x}}) \cdot X$
P5	Aerobic endogenous respiration	$b_{H} \cdot \frac{D O}{K_{D O} + D O} \cdot X$
P6	Anoxic endogenous respiration using nitrate	$b_{H} \cdot η_{} \cdot \frac{K_{D O}}{K_{D O} + D O} \cdot \frac{N O_{3}^{-}}{K_{N O_{3}^{-}} + N O_{3}^{-}} \cdot X$
P7	Anoxic endogenous respiration using nitrite	$b_{H} \cdot η_{} \cdot \frac{K_{D O}}{K_{D O} + D O} \cdot \frac{N O_{2}^{-}}{K_{N O_{2}^{-}} + N O_{2}^{-}} \cdot X$
Inferred from the model of Clement et al. (1997)
C1	Hydrolysis	$k_{D P} \cdot \frac{D P ∕ X}{K_{D P} + D P ∕ X} \cdot (1 - \frac{D P_{m i n}}{D P}) \cdot X$
C2	Anoxic growth on S_S using nitrate	$\hat{q} m a x [\frac{C_{N}}{k_{N} + C_{N}}] [\frac{C_{A}}{k_{A} + C_{A}}] X (1 - \frac{X}{c x M A X})$
C3	General expression for biomass decay	$K_{d} X$

A maximum concentration of biomass, X_max, of 205 mg-COD/kg was determined from qPCR measurements obtained from column sections collected at the conclusion of Column A, and is consistent with that reported by Kildsgaard and Engesgaard (2001). As the biomass concentration approaches X_max, the threshold term serves to slow the rate of biomass growth to a pseudosteady limit where growth equals decay. EVO retention and biomass growth increase dispersivity (Rodríguez-Escales et al., 2016; Muller and Ramsburg, 2018; Muller et al., 2018). In the work described herein, the decrease in dispersivity due to consumption of the retained emulsion and the increase in dispersivity due to biomass growth were not determined experimentally. Thus, we assume that decreases in dispersivity due to consumption of EVO are offset by increases to dispersivity due to biomass growth such that dispersivity remains constant at 0.025 m (Muller et al., 2018).

Parameters related to the simulation of the column are shown in Supplementary Table S2. Sorption of nitrogen and soluble substrate (here taken to be acetate produced from degradation of the soybean oil) to the sandy medium was assumed negligible (Mekala and Nambi, 2016).

The biochemical processes described in this model are a simplification of the ASM proposed by Kaelin et al. (2009) used for describing the stimulation of denitrification in the subsurface environment (Gonsalez et al., 2023). The biochemical reactions conceptualize denitrification as a two-step process (Kaelin et al., 2009). The modeling approach uses seven process rates to describe hydrolysis of the retained DP to produce soluble substrate, aerobic growth, aerobic endogenous respiration, anoxic growth using nitrate, endogenous respiration using nitrate, anoxic growth using nitrite, and endogenous respiration using nitrite (Table 2). The influence of these processes on the components is shown via the model stoichiometry matrix (Supplementary Table S3).

Following the ASM matrix notation, reaction terms are produced for each component (columns in Supplementary Table S3) by summing the products of each stoichiometric term in the matrix by the corresponding process rate (rows in Supplementary Table S3). This notation focuses on the underlying biochemical processes and highlights differences between models. The resulting reaction terms [R_i in Eq. (1)] are shown in Supplementary Table S4.

Parameterization of the model (Table 3) was based on Kaelin et al. (2009), although with some modifications to account for conditions in the subsurface environment, particularly the enrichment of denitrifiers when stimulated with EVO (Gihiring et al., 2011). Enrichment would result in a reduction factor, $η$ , that varies in time and space. In this study, we assume a constant and uniform value of $η$ = 1, based on the idea that anoxic growth of denitrifiers is what controls temporal variation of the biomass concentration. This simplification greatly diminishes the role of DO in the model, which is appropriate when simulating the column experiments described herein but may be inappropriate when simulating bioactivity within PRBs receiving higher mass loading of DO. The hydrolysis rate (P1) is modulated using a minimum concentration of DP emulsion, DP_min, of 158 mg-COD/kg based on the minimum measured values in the column experiment.

Table 3.

Biokinetic Parameters Used in the Process-Based Model

Symbol	Value	Units	Description
μ_H	3.0	day⁻¹	Maximum specific growth rate on substrate
b_H	0.30	day⁻¹	Endogenous respiration rate for biomass
η	1.0^a	—	Reduction factor for denitrification
k_DP	9.0	day⁻¹	Hydrolysis rate constant
K_DP	1.0	g/g	Hydrolysis saturation constant
K_s	10	mg-COD/L	Saturation coefficient for readily soluble substrate
K_DO	0.20	mg-O₂/L	Saturation/inhibition coefficient for oxygen
$K_{N O_{3}^{-}}$	0.50	mg-N/L	Saturation/inhibition coefficient for nitrate
$K_{N O_{2}^{-}}$	0.50	mg-N/L	Saturation/inhibition coefficient for nitrite
$Y_{D O}$	0.80	—	Yield coefficient for aerobic growth
$Y_{N O_{3}^{-}}$	0.30^b	—	Yield coefficient for anoxic growth using nitrate
$Y_{N O_{2}^{-}}$	0.06^b	—	Yield coefficient for anoxic growth using nitrite
X_max	205^c	mg-COD/kg	Maximum concentration denitrifying biomass
DP_min	0.14^c	g-COD/kg	Minimum concentration of emulsion

Parameter values from Kaelin et al. (2009) except: ^aassumed; ^badjusted; and ^cobserved.

The yield coefficient for heterotrophic bacteria in wastewater typically ranges from 0.46 to 0.69 (Henze et al., 1987). Values reported from denitrification in the subsurface can be as low as 0.2 (Killingstad et al., 2002). Thus, we elected to fit the yield coefficients related to nitrate and nitrite reduction (Koch et al., 2000) and rely on literature values for the remaining biokinetic parameters (Table 3). It should be noted that groundwater flow and transport simulators may require their own calibration via site-specific adjustment of flow and transport parameters.

Model simulations were produced using COMSOL Multiphysics 6.0 with the Transport of Diluted Species in Porous Media module (COMSOL, Inc., Burlington, MA) using the Newton nonlinear solver. Relevant properties of the simulations are provided in Supplementary Table S2. Initial conditions shown in Table 1 are based on Column A at the time nitrate was introduced. All components except emulsion were assumed to be uniformly distributed at time zero. These reflect the characterization of the stream water (nitrate, nitrite, COD), aquifer material (biomass), and Column B (emulsion). The initial DO was assumed to be zero because of the N₂ sparging of the water and availability of excess carbon to support rapid consumption of DO within column.

Initial biomass was taken as 0.33 mg-COD/kg indicative of the median qPCR results in samples obtained from the emulsion control column. Emulsion was distributed based on the results of the emulsion control column (Table 1). Fits were produced by minimizing the summed squared error for nitrate and nitrite using the Nelder–Mead solver (tolerance of 0.001).

Results and Discussion

Column experiments

Effluent concentrations and retention profiles are shown for the experiment in Fig. 1. Concentrations of COD and nitrate in the effluent decreased rapidly over the initial 20 pore volumes (PV) (∼3 days) indicative of the growth of denitrifying biomass. Effluent COD appeared to slowly move through a secondary maximum of 130 mg-COD/L around 70 PV. Nitrite accumulation was limited (<2.5 mg-N/L) throughout the experiment. The lack of nitrite accumulation is attributable to the availability of electron donor and enrichment of the nirS fraction of the total denitrifying gene copies (i.e., nirS, nirK, and narG) from ∼3% to ∼20% over the course of the experiment (Supplementary Table S5). While nitrate is fully transformed for 260 PV (∼20 days) when introduced at a flux of 22 g/(m²·day), this duration should not be misconstrued as an indicator of PRB longevity. In practice, PRBs operate with greater thickness, greater amounts of retained emulsion, and slower flow.

FIG. 1.

Simulations using the process-based model. Left: prediction (i.e., no adjusted parameter); right: fit by adjusting $Y_{N O_{3}^{-}}$ and $Y_{N O_{2}^{-}}$ ; top: effluent concentrations of nitrate, nitrite, soluble substrate; middle: profile of retained DP emulsion; bottom: biomass profile. DP, dispersed phase.

The experimental design used a fast flow (pore water velocity ∼2 m/day) to accelerate the depletion of carbon and breakthrough of nitrate. This approach produced a data set that required the subsequent modeling to capture both the increase in the rate of denitrification in the initial phase of the column and the decrease in the rate of denitrification around 260 PV. The accurate interrelation of the processes responsible for denitrification is necessary to inform decisions related to treatment dosing and longevity. Consumption of the retained oil and biomass growth are apparent in the retention profile data (Fig. 1). The increase in total bacteria and denitrifying gene copies (i.e., nirS, nirK, and narG) is an evidence of bacterial growth, which led to an increase in the relative abundance of denitrifiers. Similar enrichment for denitrifiers was observed by Gihring et al. (2011) and supports the model assumption of $η$ = 1.

Model comparison

The process-based model used herein engages multiple nonlinear biological processes, which foregrounds a question related to the value of increased model complexity. Killingstad et al. (2002) offer similar complexity to the process-based model examined herein but lack the capability of transforming EVO to soluble substrate. While the addition of a reaction to process EVO to soluble substrate is relatively straightforward, our preliminary investigation of the Killingstad et al. (2002) approach suggested that many of the 15 parameters fit by these authors would need to be adjusted to describe our data (Supplementary Fig. S3). While this likely could be accomplished, we chose to emphasize using the ASM framework that includes well-established values for the biokinetic parameters to limit the number of adjustable parameters.

Simulation results using the process-based model in a predictive manner (without adjustable parameters) are shown in the left side of Fig. 1. This simulation uses the parameters shown in Table 3 but with yield coefficients for nitrate (Y_NO3-) and nitrite (Y_NO2-) of 0.65 g/g as described in Kaelin et al. (2009). This simulation represents a reasonable visual prediction of effluent concentrations of nitrate and COD, although the model produces much higher levels of nitrite (∼10 mg/L). Biomass profiles shown for the simulations include initial, final, and 260 PV. Two hundred sixty PV is included here as the biomass data represent presence, not activity. Two hundred sixty PV represents nitrate breakthrough—the point in time when the influence of decay ascends.

The performance of the process-based model executed in a predictive manner is summarized by a sum square normalized error (SSNE) of 825, and small-sample Akaike information criterion (AICc, Burnham and Anderson, 2002) of 145.

The process-based model was fit to the nitrate and nitrite data by adjusting the yield coefficients for nitrate and nitrite from 0.65 g/g for both to 0.30 and 0.06 g/g, respectively (Table 3, Fig. 1—right side). Adjustment of the model enhances the conversion of nitrite such that the product of the leading coefficients in R_NO2- increases from 0.94 day⁻¹ in the prediction to 27 day⁻¹ in the fit (Table 4). Nitrate concentration, as well as the time of nitrate breakthrough, is improved in the fit due to a reduction in the product of leading coefficient for R_NO3- from 1.4 day⁻¹ in the prediction to 27 day⁻¹ in the fit (Table 4).

Table 4.

Comparison of Key Biokinetic Coefficients Between Models

Description	Clement et al. (1997)^a			This work^b
Description	Reaction	Parameter	Value, day⁻¹	Reaction	Parameter	Value, day⁻¹
Biomass growth	${\hat{r}}_{X}$	$Ŷ_{X ∕ N} \cdot q_{m a x}$	2.1	R_X	$μ_{H}$	3.0
Nitrate reduction	${\hat{r}}_{N}$	$q_{m a x}$	7.0	R_NO3-	$\frac{1 - Y_{N O_{3}^{-}}}{1.14 \cdot Y_{N O_{3}^{-}}} μ_{H}$	6.1
Nitrite reduction	Model excludes nitrite			R_NO2-	$\frac{1 - Y_{N O_{2}^{-}}}{1.72 \cdot Y_{N O_{2}^{-}}} μ_{H}$	27
Production of soluble substrate from emulsion	Emulsion reaction added^a			R_DP	k_DP	9.0
Substrate utilization with nitrate	${\hat{r}}_{A}$	$Ŷ_{A ∕ N} \cdot q_{m a x}$	15	R_S	$\frac{1}{Y_{N O_{3}^{-}}} μ_{H}$	10
Substrate utilization with nitrite	Model excludes nitrite			R_S	$\frac{1}{Y_{N O_{2}^{-}}} μ_{H}$	52

The model of Clement et al. (1997) assumes soluble substrate. The hydrolysis reaction, R_DP, was added to allow the model to be fit to the data in this study. Fits were obtained by adjusting $Ŷ_{X ∕ N},$ $Ŷ_{A ∕ N}$ , and $q_{m a x}$ .

The model developed in this study was fit to the data by adjusting $Y_{N O_{3}^{-}}$ and $Y_{N O_{2}^{-}}$ .

In addition, the biomass profile at 260 PV is improved with the fit (biomass SSNE decreased from 1.6 for the prediction to 0.89 for the fit). The decreases in the composite SSNE from 825 to 7.64 and AICc from 145 to 118 indicate that two adjustable parameters improved the overall performance of the process-based model.

The process-based approach is compared with a lower complexity, one-step nitrate reduction formulation established by Clement et al. (1997) (see Table 3, Supplementary Tables S6–S8 and Model Differences in the Supplementary Data). The limited nitrite observed in the effluent suggests that the Clement biokinetic model may have utility in describing the column data. Because the Clement model lacks the ability to process complex carbon into soluble substrate, we added R_DP, which included k_DP and K_DP, DP_min, and X_max (Supplementary Table S4), to establish a reasonable point of comparison. Note that the implementation of the Clement model included immobile DP emulsified oil and biomass to mimic the approach used with the process-based model. The values for all the added parameters are the same as those used in the process-based model (Table 3).

Predictions using the Clement model were obtained using values for the $Ŷ_{X ∕ N},$ $Ŷ_{A ∕ N}$ , and $q_{m a x}$ parameters reported by those authors [0.13 g/g, 0.84 g/g, and 7.2 mg/(mg·day), respectively]. We also fit the model to the data by adjusting three parameters, $Ŷ_{X ∕ N},$ $Ŷ_{A ∕ N}$ , and $q_{m a x}$ , to the values of 0.3 g/g, 2.2 g/g, and 7.0 g/(g·day), respectively. The Clement model prediction overestimates carbon utilization and biomass growth resulting in much earlier than observed breakthrough of the nitrate (Fig. 2—left side). The prediction has an SSNE of 13.5 and AICc of 148. The fitted model has a good visual representation of the effluent data, but overpredicts biomass growth (Fig. 2, right). The reduction in the SSNE and AICc to values of 8.84 and 111, respectively, supports the use of the three adjustable parameters.

FIG. 2.

Simulations using the kinetic model of Clement et al. (1997). Left: prediction (i.e., no adjusted parameter); right: fit by adjusting $Ŷ_{X ∕ N},$ $Ŷ_{A ∕ N}$ , and $q_{m a x}$ ; top: effluent concentrations of nitrate, nitrite, soluble substrate; middle: profile of retained DP emulsion; bottom: biomass profile.

Sensitivity analysis

The sensitivity of the fitted models to effluent nitrate was assessed using the adjoint sensitivity method available in COMSOL. Both the process-based and Clement models were found to be most sensitive to the maximum growth rate parameter [maximum specific growth rate (μ_max) and maximum specific nitrate utilization rate (q_max), respectively] (Table 5). Both models were also sensitive to the hydrolysis rate coefficient (k_DP), with all other parameters having sensitivity values that are several orders of magnitude lower (Table 5). While the highest sensitivity parameters in both models are consistent, the process-based model showed significantly higher sensitivity values suggesting that denitrification in the process-rate model may be more sensitive to the creation and consumption of soluble substrate. Thus, we examined the models further using hypothetical simulations conducted using an order of magnitude larger EVO dosage and/or an order of magnitude slower Darcy velocity. All other parameters and inputs were unchanged.

Table 5.

Model Sensitivity to Nitrate

Clement et al. (1997)^a			This work
Parameter	Description	Sensitivity	Parameter	Description	Sensitivity
k_DP	Hydrolysis rate coefficient	8.3 × 10³	k_DP	Hydrolysis rate coefficient	3.2 × 10⁷
q_max	Maximum specific nitrate utilization rate	−2.7 × 10⁴	μ_H	Maximum specific growth rate on soluble substrate	−1.8 × 10⁸
$Ŷ_{X ∕ N}$	Yield coefficient for cells	4.4 × 10¹	$Y_{N O_{3}^{-}}$	Yield coefficient for anoxic growth using nitrate	1.3 × 10⁴
$Ŷ_{A ∕ N}$	Yield coefficient for acetate (soluble substrate)	1.8 × 10⁰	$Y_{N O_{2}^{-}}$	Yield coefficient for anoxic growth using nitrite	1.9 × 10⁵
K_A	Saturation coefficient for acetate (soluble substrate)	2.1 × 10^–2	K_s	Saturation coefficient for soluble substrate	1.2 × 10³
X_max	Maximum concentration denitrifying biomass	−9.7 × 10^–3	X_max	Maximum concentration denitrifying biomass	−1.2 × 10⁰

The kinetic model of Clement et al. (1997) was modified to include k_DP and X_max for the purposes of comparison herein (see the Model Comparison section).

Given that the role of carbon processing is critical to understanding longer term performance and dosing of denitrifying PRBs established with emulsified oil, it is reasonable to explore scenarios where the amount of emulsion in the system is much greater (i.e., a longer lasting barrier). Use of emulsions at 25% vol. is common in practice. Muller et al. (2018) show that these higher concentrations can lead to DP retention that is about an order of magnitude greater than that observed in the column experiments where we accelerated the diminishment of the retained emulsion by using a relative high velocity. The order of magnitude slower flow is in the range of Darcy velocities observed on Cape Cod. Slower flows represent lower nitrate loadings and thus shift the dynamics related to carbon processing.

Both the hypothetical conditions create scenarios where carbon may be in even greater excess, with the greater emulsion mass representing capacity and the slower flow representing throughput.

Results of these additional simulations are summarized in Table 6. The base case shown in Table 6 comprises the models as calibrated (see the right-hand columns in Figs. 1 [process-based model] and 2 [Clement model]). The metrics shown in Table 6 highlight the similar performance between these calibrated models, with the process-based model exhibiting marginally better carbon utilization (e.g., less substrate discharged, more nitrate reduced, and longer treatment duration—time for which nitrate remained lower than an arbitrary limit of 1 mg-N/L). For simulations conducted with 10 × more emulsion mass, the process-based model shows an increase in the duration of treatment from 274 to 2,142 PV, although with an accompanying increase in the discharge of soluble substrate.

Table 6.

Summary Metrics from Hypothetical Simulations Containing More Emulsion Mass and/or Slower Flow

Metric	Base case		0.1 × Velocity		10 × Emulsion mass		0.1 × Velocity and 10 × emulsion mass
Metric	This work	Clement	This work	Clement	This work	Clement	This work	Clement
Treatment duration defined as time for which effluent nitrate remained <1 mg-N/L (PV)	274	213	461	687	2,142	245	3,455	3,357
Time at which the emulsion mass is exhausted (PV)	408	393	663	890	3,450	3,184	3,448	4,474
Maximum effluent concentration of soluble substrate (t > 5 PV) (mg-COD/L)	139	162	1,063	1,076	251	309	1,248	1,116
Nitrate reduced (g-N) before the appearance of nitrate in effluent	0.724	0.533	1.30	1.94	6.05	0.614	9.82	9.50
Soluble substrate discharged (g-COD) before the appearance of nitrate in effluent	0.937	1.40	20.1	20.4	11.1	4.36	216	156

PV, pore volumes.

This suggests that the use of carbon is less efficient or effective as the duration for which nitrate is treated is about a factor of 8 greater (i.e., not the full 10 × of the additional mass). In contrast, when the Clement model is used with the 10 × amount of emulsion mass, the duration of treatment is largely unchanged from the base case. Soluble substrate in the effluent was higher with greater emulsion mass, although the cumulative total of the substrate discharged downgradient was held lower by the shorter treatment duration. Metrics of the two models are more similar in cases where the flow was slower (Table 6).

To understand these differences in model performance, we explored the temporal evolution of the biokinetic rates in the base case simulations (Supplementary Fig. S4). Both models generally show production of soluble substrate near the influent of the column, where the rates of hydrolysis are the greatest. However, in the case of the Clement model, the larger rates in biomass production (and nitrate and substrate consumption) move the downgradient through the course of the simulation with the biomass concentration ultimately approaching X_max over most of the column length. In contrast, the carbon processing and growth dynamics in the process rate model are more moderate resulting in biomass growing mostly in the first half of the column length before the rate turns negative around the time of nitrate breakthrough in the effluent (∼22 days).

These dynamics in the process-based model work to produce a superior description of the measured biomass profile (Figs. 1 and 2). The uniform biomass profile in the Clement model is thought to result from the lower sensitivity to the hydrolysis rate (Table 5). Because the biomass more quickly and uniformly reaches X_max, the amount of nitrate subsequently reduced is that for maintenance of the biomass concentration (i.e., that needed for the growth rate to match the rate of decay). This maintenance demand on nitrate is limited in the one-step model due to a relatively low decay rate coefficient (Supplementary Table S8) that is not dependent on an electron acceptor (see ${\hat{R}}_{X}$ —Supplementary Table S7). Thus, the limited demand for nitrate means that adding more emulsion mass (i.e., longer term source of soluble substrate) does not influence the treatment duration in these simulations.

However, a different maintenance demand (i.e., different parameters associated with growth and decay), or longer residence times, either through slower flow or longer flow paths can change this dynamic.

Shown in Supplementary Fig. S5 are simulated biomass profiles for the scenario using a slower flow and greater mass of emulsion. In the case of the Clement model, the biomass forms a sharp profile in the first third of the column length. In contrast, the biomass profile from the process-based model is distributed throughout the porous medium. It is important to recognize that the nitrate flux in this simulation is 10% of that in the base case. The carbon processing in the Clement model quickly asserts maximum biomass conditions near the inlet of the column akin to the swift growth to X_max in the base case. Note that the dotted line in Supplementary Fig. S5 represents the approximate time the effluent nitrate concentration exceeds 1 mg-N/L. Both models show a small pulse of biomass growth moving through the column, although the dynamics in the Clement model are more pronounced. In both models, this pulse of biomass growth is a transient condition that arises where soluble substrate concentrations begin decreasing toward zero and nitrate concentrations start increasing from zero. Similar dynamics are observed in the rate plots for the base case (Supplementary Fig. S4).

Implications

This research adapted the ASM3 framework using a two-step denitrification description to produce a process-based model for nitrate reduction within the subsurface environment. While subsurface denitrification has been modeled before, previous model formulations do not consider the transformation of a complex carbon, such as when EVO is used to form a PRB for nitrate reduction. Rather, most prior approaches aim to simulate the utilization of a simple carbon substrate (e.g., acetate, ethanol, formate) in support of a one-step reduction of nitrate to dinitrogen. One benefit of the process-based framework is that the model is readily parameterized based on the rich set of literature on ASM.

In fact, when used in a predictive manner (i.e., without calibration by using literature values for all biokinetic parameters), the model produced herein was able to capture key features of the column experiment such as the release of soluble substrate and rate at which denitrification concludes upon the exhaustion of the carbon supply. Decreasing the yield coefficients for anoxic growth on nitrate and nitrite increased model performance and allowed capture of the time of nitrate breakthrough and produced better description of limited amounts of nitrite accumulation observed in the column effluent.

Biokinetic coefficients are structured differently between the process-based model and the Clement model, but a comparison of the key model coefficients is possible by grouping model parameters (Table 4). Given the similarity in model performance, it is not surprising that the coefficients shown in Table 4 have similar values. The fits using the process-based model offer a lower overall SSNE (7.64) than did the fits using the Clement model (overall SSNE 8.84), with the largest discrepancy between the models found in the biomass residuals (SSNE of 0.89 and 2.4 for the process-rate model and Clement model, respectively). Differences in the distribution of biomass between the two models may explain the difference in performance.

The Clement model biomass profiles are generally sharper with biomass existing in concentrations near the maximum threshold. Understanding the spatial distribution of the amendment and biomass is crucial for stimulating and sustaining a treatment zone within the subsurface. The variations observed in the biomass dynamics between the models indicate differences in the locations and rates of the consumption of soluble substrate and nitrate. These variations suggest that more emphasis should be placed on quantifying and utilizing biomass profiles when evaluating biokinetic models. The process-based model may hold promise in this regard as the simulations herein demonstrate good description of measured biomass concentrations in the porous medium.

Thus, the process-based model offers potential to enhance engineering capability when designing denitrifying PRBs with EVO through simulation of the interplay between transport, biomass growth, and carbon processing when making decisions about amendment dosing. Since EVO degrades into long-chain fatty acids, acetate, and hydrogen, future research may explore further refinement of the carbon and nitrogen processing to better characterize amendment utilization and nitrogen speciation.

Footnotes

Acknowledgments

In kind contributions of materials and supplies were provided by Terra Systems, Inc. The authors express appreciation to Amy Hunter at AECOM for the aquifer materials and Paul Dombrowski at ISOTEC for insights on denitrifying PRBs on Cape Cod. In addition, the authors express appreciation to the anonymous reviewers for their contributions to improving the article.

Authors' Contributions

V.L.G.: investigation, methodology, software, and writing—original draft; M.D.L.: data collection, methodology, conceptualization, resources, and writing—review and editing; K.A.M.: conceptualization, resources, and writing—review and editing; C.A.R.: conceptualization, resources, supervision, and writing—review and editing.

Author Disclosure Statement

V.L.G., Katherine K.A.M., and C.A.R. declare no conflicts of interest. M.D.L. is Vice President of Research and Development at Terra Systems, Inc.

Funding Statement

Funding for this work was provided through discretionary resources managed by C.A.R. at Tufts University.

Supplementary Material

References

Brun

, Kühni

, Siegrist

, et al. Practical identifiability of ASM2d parameters—Systematic selection and tuning of parameter subsets. Water Res, 2002; 36(16):4113–4127; doi: 10.1016/S0043-1354(02)00104-5

Burnham

, Anderson

. Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach. Springer-Verlag: New York, NY, USA; 2002.

Calderer

, Gibert

, Martí

, et al. Denitrification in presence of acetate and glucose for bioremediation of nitrate-contaminated groundwater. Environ Technol, 2010; 31(7):799–814; doi: 10.1080/09593331003667741

Clement

, Peyton

, Skeen

, et al. Microbial growth and transport in porous media under denitrification conditions: Experiments and simulations. J Contam Hydrol, 1997; 24(3):269–285; doi: 10.1016/S0169-7722(96)00014-9

Coulibaly

, Borden

. Impact of edible oil injection on the permeability of aquifer sands. J Contam Hydrol, 2004; 71(1–4):219–237; doi.org/10.1016/j.jconhyd.2003.10.002

Crocker

, Berge

, Ramsburg

. Encapsulated delivery of reactive iron particles using oil-in-water emulsions. In: GQ07: Securing Groundwater Quality in Urban and Industrial Environments. Proceedings of 6th International Groundwater Quality Conference held in Fremantle, Western Australia, December 2–7, 2007. IAHS Publ. no. 324.2008.

Dombrowski

, Parece

, Begley, et al. Denitrifying permeable reactive barriers on Cape Cod: Bench scale studies and implementation of the first in-situ EVO PRB. In: Bioremediation and Sustainable Environmental Technologies—2017. Fourth International Symposium on Bioremediation and Sustainable Environmental Technologies, Miami, FL, USA; 2017.

, Peng

, Cao

, et al. Mechanisms and microbial structure of partial denitrification with high nitrite accumulation. Appl Microbiol Biotech, 2016; 100:2011–2021; doi: 10.1007/s00253-015-7052-9

Gihring

, Zhang

, Brandt

, et al. A limited microbial consortium is responsible for extended bioreduction of uranium in a contaminated aquifer. Appl Environ Microbiol, 2011; 77(17):5955–5965; doi: 10.1128/AEM.00220-11

10.

Glass

, Silverstein

. Denitrification kinetics of high nitrate concentration water: PH effect on inhibition and nitrite accumulation. Water Res, 1998; 32(3):831–839; doi: 10.1016/S0043-1354(97)00260-1

11.

Gonsalez

, Dombrowski

, Lee

, et al. The influence of pH on subsurface denitrification stimulated with emulsified vegetable oil. Water, 2023; 15(5):883; doi: 10.3390/w15050883

12.

Gujer

, Henze

, Mino

, et al. Activated sludge model no. 3. Water Sci Technol, 1999; 39(1):183–193; doi: 10.1016/S0273-1223(98)00785-9

13.

Henze

, Grady

CPL

, Gujer

, et al. A general model for single-sludge wastewater treatment systems. Water Res, 1987; 21(5):505–515; doi: 10.1016/0043-1354(87)90058-3

14.

Henze

, Gujer

, Mino

, et al. Activated sludge model no.2D, ASM2D. Water Sci Technol, 1999; 39(1):165–182; doi: 10.1016/S0273-1223(98)00829-4

15.

Hiatt

, Grady

CPL

. An updated process model for carbon oxidation, nitrification, and denitrification. Water Environ Res, 2008; 80(11):2145–2156; doi: 10.2175/106143008X304776

16.

Huno

SKM

, Rene

, van Hullebusch

, et al. Nitrate removal from groundwater: A review of natural and engineered processes. J Water Supply Res Techonol Aqua, 2018, 67(8):885–902; doi: 10.2166/aqua.2018.194

17.

Jenkins

, Ling

, Ciesielczuk

, et al. Detection and identification of bacteria in clinical samples by 16S RRNA gene sequencing: Comparison of two different approaches in clinical practice. J Med Microbiol, 2012; 61(Pt_4):483–488; doi: 10.1099/jmm.0.030387-0

18.

Kaelin

, Manser

, Rieger

, et al. Extension of ASM3 for two-step nitrification and denitrification and its calibration and validation with batch tests and pilot scale data. Water Res, 2009; 43(6):1680–1692; doi: 10.1016/j.watres.2008.12.039

19.

Kappeler

, Gujer

. Estimation of kinetic parameters of heterotrophic biomass under aerobic conditions and characterization of wastewater for activated sludge modelling. Water Sci Technol, 1992; 25(6):125–139; doi: 10.2166/wst.1992.0118

20.

Kildsgaard

, Engesgaard

. Numerical analysis of biological clogging in two-dimensional sand box experiments. J Contam Hydrol, 2001; 50(3–4):261–285; doi: 10.1016/S0169-7722(01)00109-7

21.

Killingstad

, Widdowson

, Smith

. Modeling enhanced in situ denitrification in groundwater. J Environ Eng, 2002; 128(6):491–504; doi: 10.1061/(ASCE)0733-9372(2002)128:6(491)

22.

Koch

, Kühni

, Gujer

, et al. Calibration and validation of activated sludge model no. 3 for Swiss municipal wastewater. Water Res, 2000; 34(14):3580–3590; doi: 10.1016/S0043-1354(00)00105-6

23.

, Qian

, Gao

, et al. Functional enzyme-based approach for linking microbial community functions with biogeochemical process kinetics. Environ Sci Technol, 2017; 51(20):11848–11857; doi: 10.1021/acs.est.7b03158

24.

, Liu

, Dong

, et al. Quantifying fate and transport of nitrate in saturated soil systems using fractional derivative model. Appl Math Model, 2020; 81:279–295; doi: 10.1016/j.apm.2019.12.005

25.

Mekala

, Nambi

. Transport of ammonium and nitrate in saturated porous media incorporating physiobiotransformations and bioclogging. Bioremediation J, 2016; 20(2):117–132; doi: 10.1080/10889868.2015.1113925

26.

Muller

, Ramsburg

. Influence of nonwetting phase saturation on dispersivity in laboratory-scale sandy porous media. Environ Eng Sci, 2018; 35(10):1062–1074; doi: org/10.1089/ees.2017.0444

27.

Muller

, Esfahani

, Chapra

, et al. Transport and retention of concentrated oil-in-water emulsions in porous media. Environ Sci Technol, 2018; 52:4256–4264; doi: 10.1021/acs.est.7b06012

28.

Pan

, Ni

, Lu

, et al. Evaluating two concepts for the modelling of intermediates accumulation during biological denitrification in wastewater treatment. Water Res, 2015; 71: 21–31; doi: 10.1016/j.watres.2014.12.029

29.

Parece

, Begley

, Shreve-Gibb

, et al. Orleans Massachusetts permeable reactive barriers demonstration project: A non-traditional technology to meet estuary nitrogen TMDLs. In: Proceedings of WEFTEC 2019—92nd Annual Water Environment Federation's Technical Exhibition and Conference, 2019; pp. 3608–3623.

30.

Rodríguez-Escales

, Folch

, van Breukelen

, et al. Modeling long term enhanced in situ biodenitrification and induced heterogeneity in column experiments under different feeding strategies. J Hydrol, 2016; 538:127–137; doi: 10.1016/j.jhydrol.2016.04.012

31.

Saccani

, Hakanen

, Sindhya

, et al. Potential of interactive multiobjective optimization in supporting the design of a groundwater biodenitrification process. J Environ Manage, 2020; 254:109770; doi: 10.1016/j.jenvman.2019.109770

32.

Smith

, Miller

, Brooks

, et al. In situ stimulation of groundwater denitrification with formate to remediate nitrate contamination. Environ Sci Technol, 2001; 35(1):196–203; doi: 10.1021/es001360p

33.

Sun

, Petersen

, Clement

, et al. Effect of reaction kinetics on predicted concentration profiles during subsurface bioremediation. J Contam Hydrol, 1998; 31(3):359–372; doi: 10.1016/S0169-7722(97)00072-7

34.

Wang

, Ma

, Guo

, et al. Experiment and multicomponent model based analysis on the effect of flow rate and nitrate concentration on denitrification in low-permeability media. J Contam Hydrol, 2020; 235:103727; doi: 10.1016/j.jconhyd.2020.103727

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.42 MB