## ABSTRACT

We developed a user-friendly numerical model, Shoreline Bioremediation Model (SBM), to simulate the biodegradation and bioremediation of oil entrapped within shorelines. The model takes oil properties and environmental conditions as input and produces variation of oil concentration with time, up to several years from the time of the spill. SBM is equipped with a user-friendly graphical user interface (GUI). The accessibility and easy-to-use interface allow the user to quickly produce several biodegradation and bioremediation scenarios before they are implemented at the contaminated shoreline. The model has been calibrated to predict the biodegradation rate of saturates and aromatics, but it can be also used to predict the biodegradation rates of individual oil components and to decide on bioremediation studies in shorelines based on the enhancement due to biostimulation (addition of nutrients) or bioaugmentation (addition of hydrocarbon degrading organisms). The GUI provides the oil concentration with time along with best case and worst case scenarios, which is commonly needed for decision making.

## INTRODUCTION:

Accidental oil spills such as the Lakeview Gusher spill in California (Harvey, 2010), the Exxon Valdez oil spill in Prince William Sound (PWS) of Alaska, the Deepwater Horizon oil spill in the Gulf of Mexico and others have resulted in the total spillage of 37 billion barrels of crude oil which is more than the annual consumption of oil in 2006 (30 billion barrels) and have resulted in contamination of several miles of shorelines (Wang *et al.*, 2011). Over the last several years, bioremediation has been used by researchers in order to improve the natural degradation process of residual oil in order to diminish the ecological impact, specifically on intertidal communities, which have been affected by major oil spills.

Several studies have been conducted on bioremediation and the factors that affect the rate of biodegradation (Atlas, 1981; Atlas and Hazen, 2011; Jones *et al.*, 2008; Karamalidis *et al.*, 2010; Kostka *et al.*, 2011; Leahy and Colwell, 1990; Lee *et al.*, 2011). The surface area of the oil is important because the growth of microorganisms occurs at the interface of the oil and water (Atlas and Bartha, 1992; Vilcáez *et al.*, 2013). The formation of oil-in-water emulsions is important in the resulting uptake of hydrocarbons by the oil degrading microorganisms (Mohanty *et al.*, 2013; Singer and Finnerty, 1984). Venosa *et al.* (2010) conducted laboratory experiments to study the biodegradation of 19 year old lingering weathered oil in Prince William Sound (PWS) from the Exxon Valdez Oil Spill at three different islands. They observed up to 80% biodegradation of the polycyclic aromatic hydrocarbons (PAHs), the compounds typically of concern due to their high toxicity. The results obtained from these experiments can be used to gain a better understanding of the kinetic processes involved in biodegradation when compared with numerical models.

Numerical models are routinely used to estimate the biodegradation rates of hydrocarbons and the biodegradation potential at a site. Several researchers have used first order decay models to simulate the biodegradation of oil (Chen *et al.*, 2008; Venosa *et al.*, 2010; Venosa *et al.*, 1996; Xu and Obbard, 2004). However, these models cannot be used to predict as the fitted rates cannot be related to the biochemistry of the system. In terms of microbial growth, it is widely accepted that increase in biomass follows Monod kinetics (Desai *et al.*, 2008; Essaid *et al.*, 2003; Geng *et al.*, 2014; Geng *et al.*, 2013). Desai *et al.* (2008) investigated the biodegradation kinetics of PAHs and their studies indicated that the prediction of natural or enhanced biodegradation of PAHs cannot be based on single compound kinetics as this assumption would overestimate the rate of disappearance. Some researchers have also developed graphical user interfaces (GUI) to simulate the biodegradation process eg; (BIOMOC; Essaid and Bekins, 1997). In addition, EPISUITE (EPA, 2012) contains tools such as BIOWIN and BioHCwin that estimate the biodegradation rates of different hydrocarbons in their NAPL phase and their predicted biodegradability. These tools do not visually show the biodegradation of hydrocarbons at a specific site and require other modeling tools to use the parameters presented in these databases.

In this study, we present the details of Shoreline Bioremediation Model (SBM) that is capable of predicting the biodegradation of hydrocarbons at contaminated shorelines and other sites. We use the experimental dataset presented in Venosa *et al.* (2010) as a case study to show the robustness of the SBM model.

## METHODS:

### Numerical model

The numerical model SBM provides a front-end to the Monod kinetic model, BIOB (Geng *et al.*, 2014; Geng *et al.*, 2013) which simulates the growth of bacteria in the presence of nutrients and the decay of the hydrocarbon attached to the sediment. The decay of hydrocarbons can be mathematically expressed as follows (Geng *et al.*, 2014):

Where S is the concentration of the PAHs (mg S/kg sediment), X is the concentration of the biomass (mg X/kg sediment), Y_{x} is the biomass yield coefficient for growth on the hydrocarbon (mg X/mg S) and μ is the growth rate of the biomass (day^{−1}) given by (Geng *et al.*, 2014):

Where μ_{max} is the maximum growth rate (day^{−1}), X_{max} is the maximum allowable microbial concentration (mg X/kg sediment), K_{s} is the half saturation concentration of hydrocarbon (mg S/kg sediment), N is the nitrogen-based nutrient concentration (nitrate+nitrite+ammonia) (mg-N/L of pore water), and K_{N} is the half-saturation concentration for nitrogen consumption (mg-N/L of pore water). The term was introduced by (Geng *et al.*, 2014) to account for the decrease in biomass accumulation when the microbial concentration approaches its maximum value. Geng *et al.* (2014) also included an expression for hopane removal due to physical processes in Eq. 1. However, it was noted in Venosa *et al.* (2010) study that hopane concentration remained constant throughout the experiment and hence this term was not included in this study. BIOB also includes equations to solve for the consumption of nutrients and O_{2}, and the production of CO_{2}. In the Venosa *et al.* (2010) study, the nutrient concentration, O_{2} and CO_{2} levels were not measured as they were always present in excess, and hence these equations were not considered in the current model.

The biomass growth can be expressed mathematically as follows (Geng *et al.*, 2014):

Where k_{d} is the endogenous biomass decay rate (day^{−1}).

The model can be solved for the concentration of biomass and PAH at different times by numerically solving Equations (1) through (3) using an ordinary differential equation (ODE) solver. BIOB uses a Runge-Kutta Felhberg adaptive time-stepping scheme to solve the equations (Chapra and Canale, 1998).

Figure 1 shows the layout of the SBM model. The GUI uses textboxes for input where the concentration of oil, nutrients, temperature, number of days of simulation and biomass concentration can be input by the user. The GUI also has a help button which shows the different constants used by the model and the references from which they have been obtained from. The model can be run by pressing the calculate button; the results are populated in the charts presented below. The oil biodegradation with time is presented in the bottom left hand graph, whereas the biomass concentration (growth) with time is presented in the bottom right hand graph. The graphs also show variation of oil biodegradation and biomass (growth) for a 50% increase and decrease in the maximum growth rate (μ_{max}) value.

### Parameter sensitivity

Sensitivity analyses were performed to test the sensitivity of the model to the parameter estimates and to identify the kinetic parameters that affect the hydrocarbon biodegradation process. These analyses also ensured that the best fit parameters estimated for each dataset were indeed the global minimum for that dataset. The covariance matrix, V_{x}, of the parameters was evaluated according to the equation (Bard, 1974):

where H is the Hessian matrix where terms are the second derivative of the objective function with respect to the parameters. The variance of errors, σ^{2}, can be expressed by the following formula (Bard, 1974):

where F is the value of the objective function at the optimum, n is the number of observations and p is number of estimated parameters.

### Inverse problem – Parameter estimation

The Monod kinetic model requires several parameters that define the system as shown in Equations (1) – (3). Therefore, we need to estimate these model parameters before using SBM. A parameter estimation tool based on Genetic Algorithms (GA) was used to estimate the parameters required for simulating the experimental dataset (Torlapati, 2013). Since no microbial concentration was included in Venosa *et al.* (2010), the maximum microbial concentration (X_{max}) and the initial microbial concentration (X) were used as adjustable parameters. Therefore, a total of 7 parameters were estimated using a parameter estimation tool. Table 1 shows the parameters that were estimated along with their higher and lower bounds. The general procedure of a parameter estimation process using GA is described below.

GAs are a branch of evolutionary algorithms which are based on the “survival of the fittest” paradigm. The six key steps involved in a traditional GA are: encoding, population generation, selection, crossover, mutation, and termination (Holland, 1975). The GA starts with a randomly-generated initial set of solutions between the bounds provided by the user (also known as chromosomes); this is called the initial population. The objective function is calculated for each chromosome and is assigned as the fitness. A lower value of fitness is desirable for a minimization problem. Based on the fitness values, two parents are selected using a selection process and the selected parents undergo a crossover, where the genetic information is exchanged between the parents using a crossover function. It is also possible that an offspring generated from the crossover of the parents could undergo a mutation operation governed by a mutation probability. The fitness of the offspring is calculated and is combined with the entire population. The individuals with poor fitness are removed from the population (death) at the end of the generation. Specific details and methods used for the GA used in the parameter estimation process for reactive contaminant models are available in Torlapati (2013).

## BIODEGRADATION OF EXXON VALDEZ OIL:

Venosa *et al.* (2010) conducted laboratory experiments to test the biodegradability of 19 year old lingering oil due to the 1989 Exxon Valdez Oil Spill. Samples were collected from three different gravel beaches: Eleanor Island (EL107), Knight Island (KN114A), and Smith Island (SM006B) whose respective mass weathering indices (MWI) were 30%, 76%, and 60%. Samples of these oil-contaminated sediment were collected from each site by digging up to 10–40 cm below the surface. The volume of each sample was about 50 L and the collected samples were used for aerobic microcosm experiments. Seawater samples from each island were also collected for microcosm experiments.

The collected samples were then used for experiments to study: a) natural attenuation, the biodegradation due to the nutrients naturally present in seawater/sediment; and b) nutrient-amended treatment, where nutrients were added to the samples and the concentrations of the nutrients were maintained by adding 10 mg/L of KNO^{3} to ensure no nutrient limitation occurs (Boufadel *et al.*, 1999; Du *et al.*, 1999; Venosa *et al.*, 1996). The samples were collected at 0, 14, 28, 56, 112 and 168 days and the temperature during the experimental study was similar to the conditions *in situ* (15 °C).

The PAH concentrations were normalized using the biomarker concentration (hopane). Biomarkers are compounds that biodegrade very slowly within the time frame of interest and have been used to evaluate the effectiveness of biodegradation (Bragg *et al.*, 1994; Venosa *et al.*, 1996). The results obtained from these experimental analyses indicated that there was significant biodegradation in the nutrient-amended samples, even when the mass weathering index (MWI) was high. The experimental results were fitted by Venosa *et al.* (2010) using a first-order decay model. The results for natural attenuation also showed significant amount of biodegradation, although the rates were lower than those of the nutrient-amended experiments.

## RESULTS AND DISCUSSION:

The initial hopane-normalized PAH concentration used was 58.68 mg/mg hopane and we set the concentration of nitrogen (N) in Eq. 2 at 10 mg/L for the complete duration of 168 days since Venosa *et al.* (2010) study ensured that there was sufficient nutrient concentration available for nutrient-amended experiments. The goal was to fit the model to the observed data to estimate parameter values. Since the experiments were performed in triplicates, the average value was obtained at each sample measurement and was used as an input to the GA. The fitness of each solution was then evaluated by calculating the error (difference) between the model concentration and the input experimental concentration that was provided. This error was normalized by the input concentration at that time. This normalized error was squared and added for all the experimental data and assigned as fitness for that population. The objective of the GA was to minimize the weighted least squares (WLS). All the simulations were performed on an Intel Core i7 desktop computer with 8 cores and a total clock speed of 3.40 GHz. The time step used for all the simulations was 3 hours. The estimated kinetic parameters for the experimental dataset from the Smith Island are shown in Table 1. Based on Eq. 5, the variance of errors (σ^{2}) at the optimum was 0.0017 for the nutrient amended experiments. To evaluate the standard deviation using Eq. 4, the objective function for each parameter was evaluated by varying the best fit parameter between –20% and +50% while the rest of the parameters were kept constant; and the results were plotted (not shown here). A quadratic function was fitted to these “objective function vs deviation from the best fit parameter” plots and the second derivative of this fitted quadratic function was evaluated (Hessian matrix; H in Eq. 5). For the sake of simplicity, only diagonal terms were evaluated. The standard deviation was small for most of the parameters except X_{0} which indicates the low sensitivity of the system to changes in initial biomass concentration. The estimated parameters are well within the range of the literature values.

Subsequently, the best-fit parameters obtained for the nutrient amended experiments (Table 1) were used to estimate the nutrient concentration in the natural attenuation experiment. The microcosms were gently mixed and supplied with seawater collected from each site every 4 hours to provide reaeration. The concentration of nitrogen for natural attenuation estimated by the parameter estimation tool was 2.21 mg/L.

The time-varying hopane-normalized PAH concentrations for the nutrient amended experiments and natural attenuation obtained from the parameter estimation tool are shown in Figures 2a and 2b respectively. Comparisons were made against the experimental data collected on different days using an exponential fit (standard deviations are included). The model simulations predicted the experimental results well. In order to compare the goodness of fit, the objective function (WLS) was calculated (0.014 and 0.029 for nutrient amended and natural attenuation experiments, respectively, vs 0.06 and 0.012 for the exponential fit). The lower value of the objective function indicated that the SBM provided a better fit for the nutrient amended experiments than the exponential fit for the natural attenuation data. The exponential fit was marginally better for Smith Island, however. The SBM's prediction for the Smith Island's nutrient amended experiments is more intuitive as there is a slight lag before the actual biodegradation begins, whereas the exponential fit indicates an immediate biodegradation. This is not accurate because the biomass requires some time to grow before they are able to biodegrade the PAHs and this observation is consistent with the experimental data.

## SENSITIVITY ANALYSIS:

### Sensitivity to nutrient concentration

Nutrient availability is an important parameter in biodegradation, as it has been observed by several researchers that the nitrate concentration should be about 2–10 mg N/L to achieve maximum biodegradation rate (Atlas, 1981; Atlas and Bartha, 1972; Boufadel *et al.*, 1999; Du *et al.*, 1999; Wrenn *et al.*, 2006). Simulations were performed with different nutrient concentrations between 1 and 10 mg/L (Figure 3). The nutrient concentration for the base case scenario was 10 mg/L (Figure 3). The results are presented in Figure 3 and it can be observed from the figures that for all datasets that, as the nutrient concentration decreases, the biodegradation rate decreases. Compared to the base case scenario, the decrease in nutrient concentration resulted in larger oil values at the end of the experiments; there was an decrease in PAH biodegradation at the end of 168 days experiment by 30%, 81%, 161%, 201%, 258%, 346% and 429% for nutrient concentrations of 5, 3, 2, 1.75, 1.5, 1.25 and 1 mg/L respectively. Therefore, nutrient concentration plays a significant role in biodegradation of PAHs on the Smith Island.

### Sensitivity to initial biomass concentration

Different scenarios were chosen in addition to the base case scenario to study the sensitivity of the system to the initial biomass concentration. The initial biomass concentration for the base case scenario was 2.46E-03 mg X/kg of sediments. The initial biomass concentrations considered were 2.46E-02, 2.46E03, 2.46E-05, and 2.46E-07 mg X/kg of sediments (Figure 4). The initial lag before biodegradation started was larger when the initial biomass concentration was lower. Compared to the base case scenario, there was a decrease in the PAH biodegradation by 14% and 20% when the biomass concentration was decreased by 100 and 10000 fold respectively. There was no difference in the final PAH concentration when the biomass concentration was increased 10 fold. The initial biomass concentration seems to play a significant role.

### Sensitivity to temperature

In general, the rate of oil biodegradation increases with temperature. However, there are exceptions due to the fact that some microorganisms thrive at low temperature. In this study, we use the same approach used for wastewater treatment and the effect of temperature on the reaction rate of a biological process is usually expressed as (Metcalf and Eddy, 1991, p 372):

Where r_{t} is the reaction rate at t °C, r_{20} is the reaction rate at 20 °C, , t is the temperature in °C, θ is the temperature-activity coefficient given by 1.035 (Metcalf and Eddy, 1991) for trickling filters which is a fair approximation to the gravel microcosms used in the Venosa *et al.* (2010) experiments.

The ambient temperature in the experiments conducted by Venosa *et al.* (2010) was 15 °C and the maximum biomass growth rate (μ_{max}) obtained from parameter estimation tool at this temperature was 1.35 day^{−1}. Using Eq. 6, the biomass growth rate at temperatures of 10, 20 and 25 °C were calculated to be 1.13, 1.60, 1.90 day^{−1}(Figure 5). As the temperature increased, the biodegradation rate increased. Compared to the base case scenario, the hopane-normalized PAH concentration after 168 days there was an increase in PAH biodegradation by 26% and 43% for temperatures of 20 and 25 °C respectively whereas for a temperature of 10 °C, there was decrease in PAH biodegradation by 39%.

## CONCLUSION:

In this study, the biodegradation model SBM was used to simulate the experimental data from nutrient amended and the natural attenuation treatments from Venosa *et al.* (2010) study of the Smith Island. The kinetic parameters required for the SBM were estimated using a parameter estimation tool based on GA using the experimental data from the nutrient amended treatments. Subsequently, the best-fit parameters obtained from the nutrient amended experiments were used to estimate the nutrient concentration for the natural attenuation experiments. The estimated parameters were able to predict the experimental results well for both datasets and comparison with the exponential decay model showed that the results from SBM had better goodness of fit for nutrient amended experiments whereas for natural attenuation treatment, the exponential fit provided a marginally better fit. In particular, SBM was able to capture the initial slow biodegradation due to the lag phase in microbial growth. Sensitivity analyses studies were performed to study the system response to nutrient concentration, initial biomass concentration and ambient temperature. The system was sensitive to small changes in the nutrient concentration. When the initial biomass concentration was varied, it affected the oil biodegradation significantly during the initial phase of the biodegradation process. As the temperature increased, the rate of biodegradation also increased. The effect of temperature on biodegradation rate needs to be validated with published data, but that is beyond the scope of this manuscript. We determine that the Monod kinetic model, SBM is a useful tool that can be used to predict the extent of biodegradation and the biodegradation potential in marine environments (e.g., beaches) where the oil is trapped in sediments with no significant transport or dissolution.