Field site example: Keesler Air-Force Base¶

This notebooks outlines analysis of a field site example of contaminant transport at the Keesler Air Force Base in Mississippi, USA$^1$. This example data has been introduced along the distribution of BIOSCREEN$^2$ v1.4. We use this example here specifically to showcase the comparison of mibitrans with BIOSCREEN. Note that the data is sourced from the USA, as a results units are mixed imperial and metric. For the modelling below, all parameters are converted to metric.

Authors: Jorrit Bakker, Alraune Zech

Site description¶

At the site, underground storage tanks (UST) with a mixture of BTEX were stored. During their removal in 1987, multiple USTs leaked into the groundwater. Remediation efforts included a 'bioventing' system and an in well aeration system. Monitoring wells were used to evaluate natural attenuation.

Flow and contaminant transport is within an unconfined aquifer with groundwater levels varying between 5 and 9 ft below ground level. Groundwater flow is to the north-east, which eventually discharges into the Back Bay of Biloxi, 2100ft downgradient.

Local geology is a fine- to medium-grained sand, underlain by a clay layer at 20ft depth. At the bottom of the sandy layer peat is locally present. Thickness and continuity of the clay layer are unknown, since not every boring showed the clay layer at 20ft. During site modelling, a continuous clay layer was assumed.

Transport Field Observations

Field observations are available at different distances from the spill source after 6 years with contaminant concentrations of

distance x in ft	0	32	64	192	288
concentration	12	5	1	0.5	0.001

Modeling parameters¶

Parameters used in the BIOSCREEN$^2$ modelling differ at several occasions from values reported in $^{1,3}$. For consistency, we use the values reported in BIOSCREEN$^2$.

Hydrogeological parameters¶

Values of hydraulic gradient $i$ and conductivity $K$ reported by $^{1,3}$ vary, with values of $i$ between 0.003 and 0.0083 ft/ft and of $K$ in the range od 40, 61, 32 ft/day$^3$; the same authors estimate effective porosity with $0.25$, resulting in groundwater velocity of 0.8 ft/day. BIOSCREEN$^2$ uses the values of $i=0.003$ ft/ft and $k=31.2$ ft/day ($0.011$ cm/s), resulting in groundwater velocity of 0.31 ft/day.

Dispersivity values have not been identified from field data for the site. BIOSCREEN uses initial dispersivity values based on estimated plume length of 280ft, with 13.3ft, 1.3ft and 0ft for longitudinal, transverse horizontal and transverse vertical dispersivity, respectively. Values have been adapted during calibration to field data to $\alpha_L = 32.5$ft, $\alpha_T= 3.25$ft and $\alpha_V=0$ft. These values are at the higher end for observed reliable longitudinal macrodispersivities and beyond any reliable transverse dispersivity values$^{4,5}$.

We provide a sensitivity analysis showing the impact of dispersivity values below. For consistancy, we stick to the BIOSCREEN - calibrated value, although they are too high.

Attenuation parameters¶

The fraction of organic carbon is $0.000057$ based of lab analysis, soil bulk density is estimated to be $1.7$ (kg/L)$^2$. Partition coefficient of BTEX differs by order of magnitude (38, 135, 95, 240 L/kg for B, T, E, X). BIOSCREEN uses the value of benzene; $38$ L/kg.

Instant reaction model parameters¶

Electron acceptor concentrations are based on groundwater sampling performed in 1995. For oxygen, nitrate and sulfate, the difference between background concentration in the aquifer and minimum concentrations in the plume are used. For ferrous iron and methane, average concentrations in the plume area are used.

Source conditions¶

Conditions and parameters of the contaminant source were inferred from BTEX monitoring data and geologic logs. The source thickness is estimated to be 10ft, with a total 2000kg of pure BTEX. Total source width of 130ft, with 28ft zone of 0.057 g/m3, 30ft zone of 2.508 g/m3, 14ft zone of 13.68, 30ft zone of 2.508 g/m3, 28ft zone of 0.057 g/m3.

Note that BIOSCREEN adjusted source concentrations$^2$ to fit observations and results of the instant reaction model: The input source concentration (representing concentrations 6 yrs ago) were adjusted so the source concentration at $x=0$ for the instant reaction model (representing concentrations now) were equal to 12 mg/L. The final result was a source concentration of 13.68 mg/L in the center of the source zone resulting in an source concentration (at $x=0$) after $6$ years of 12.021 mg/L for the instant reaction model. Note that the source decay term is different for the first order decay and instantaneous reaction models.

Model parameters¶

Model dimensions used for the BIOSCREEN model are:

• length = 500ft = 150m,
• width = 200ft = 60m,
• duration = 6y.

In BIOSCREEN resolution is always a fraction of model dimensions, 1/10 for length, 1/4 for width and 1/10 for time. For the mibitrans package, model resolution is adaptable. We chose a finer resolution as in BIOSCREEN.

References¶

$^1$ Cost and performance case study report - Keesler Air Force Base Base Exchange Service Station

$^2$ Newell, C. J., McLeod, R. K., & Gonzales, J. R. (1997). BIOSCREEN natural attenuation decision support system version 1.4 revisions, Tech. rep., U.S. EPA.

$^3$ Corrective action plan for the risk-based closure of the base exchange service station, area of concern - A (ST-06) Keesler Air Froce Base, Mississippi (make actual citation)

$^4$ Zech, A., Attinger, S. ... (2015). Is unique scaling of aquifer macrodispersivity supported by field data?, Water Resources Research 51(9)

$^5$ Zech, A., Attinger, S. ... (2019). A Critical Analysis of Transverse Dispersivity Field Data, Ground Water, 57(4)

In [1]:

import matplotlib.pyplot as plt
import numpy as np
import mibitrans as mbt

ft = 3.281 #Conversion factor ft/m

import matplotlib.pyplot as plt import numpy as np import mibitrans as mbt ft = 3.281 #Conversion factor ft/m

Setup Data Input¶

In [2]:

field_data_x = np.array([0, 32/ft, 64/ft, 192/ft, 288/ft])
field_data_c = np.array([12,5,1,0.5,0.001])

field_data_x = np.array([0, 32/ft, 64/ft, 192/ft, 288/ft]) field_data_c = np.array([12,5,1,0.5,0.001])

In [3]:

hydro = mbt.HydrologicalParameters(
    h_conductivity = 0.011 / 100 * 3600 * 24, #m/d
    h_gradient = 0.003, # m/m
    porosity = 0.3,
    alpha_x = 32.5 / ft, #m dispersivity BIOSCREEN after calibration (very high)
    alpha_y = 3.25 / ft, #m dispersivity BIOSCREEN after calibration (very high)
    alpha_z = 0
)
print("Hydraulic conductivity K = {:.2f} m/d".format(hydro.h_conductivity))
print("Flow velocity v = {:.4f} m/d".format(hydro.velocity))

hydro = mbt.HydrologicalParameters( h_conductivity = 0.011 / 100 * 3600 * 24, #m/d h_gradient = 0.003, # m/m porosity = 0.3, alpha_x = 32.5 / ft, #m dispersivity BIOSCREEN after calibration (very high) alpha_y = 3.25 / ft, #m dispersivity BIOSCREEN after calibration (very high) alpha_z = 0 ) print("Hydraulic conductivity K = {:.2f} m/d".format(hydro.h_conductivity)) print("Flow velocity v = {:.4f} m/d".format(hydro.velocity))

Hydraulic conductivity K = 9.50 m/d
Flow velocity v = 0.0950 m/d

Attenuation parameters for no decay model setup.

In [4]:

att = mbt.AttenuationParameters(
    bulk_density=1700,             # kg/m^3
    partition_coefficient=0.038,   # m^3/kg (units cancel out with bulk density)
    fraction_organic_carbon=0.000057,
    decay_rate = 0,                # 1/day, no decay
)
att.calculate_retardation(hydro.porosity)
print("Retardation based on input values =", att.retardation)

att = mbt.AttenuationParameters( bulk_density=1700, # kg/m^3 partition_coefficient=0.038, # m^3/kg (units cancel out with bulk density) fraction_organic_carbon=0.000057, decay_rate = 0, # 1/day, no decay ) att.calculate_retardation(hydro.porosity) print("Retardation based on input values =", att.retardation)

Retardation based on input values = 1.012274

We will modify the attenuation parameter for linear decay model setup later by defining a decay rate of $0.15$ years.

Additional parameters for instant reaction model setup.

In [5]:

### input dataclass for instant reaction model
ea_concentrations = mbt.ElectronAcceptors(
    delta_oxygen=2.05 - 0.4, #background conc - minimum conc, [g/m3]
    delta_nitrate=0.07 - 0, #background conc - minimum conc, [g/m3]
    ferrous_iron=16.6, #average conc, [g/m3]
    delta_sulfate=26.2 - 3.8, #background conc - minimum conc, [g/m3]
    methane=6.6, #average conc, [g/m3]
)

### input dataclass for instant reaction model ea_concentrations = mbt.ElectronAcceptors( delta_oxygen=2.05 - 0.4, #background conc - minimum conc, [g/m3] delta_nitrate=0.07 - 0, #background conc - minimum conc, [g/m3] ferrous_iron=16.6, #average conc, [g/m3] delta_sulfate=26.2 - 3.8, #background conc - minimum conc, [g/m3] methane=6.6, #average conc, [g/m3] )

In [6]:

source = mbt.SourceParameters(# Note; as plume is symmetric, source is described from center outwards in y-direction.
    source_zone_boundary=np.array([7/ft, 37/ft, 65/ft]), #[m]
    source_zone_concentration=np.array([13.68, 2.508, 0.057]), #[g/m3]
    depth=10/ft, #[m]
    total_mass=2000000, #[g]
)

source = mbt.SourceParameters(# Note; as plume is symmetric, source is described from center outwards in y-direction. source_zone_boundary=np.array([7/ft, 37/ft, 65/ft]), #[m] source_zone_concentration=np.array([13.68, 2.508, 0.057]), #[g/m3] depth=10/ft, #[m] total_mass=2000000, #[g] )

In [7]:

# Visualize source concentration distribution (e.g. for input check):
source.visualize_source_zone()
plt.show()

# Visualize source concentration distribution (e.g. for input check): source.visualize_source_zone() plt.show()

In [8]:

model = mbt.ModelParameters(model_length=500/ft,
                                 model_width=200/ft,
                                 model_time=6*365,
                                 dx=2/ft,
                                 dy=1/ft,
                                 dt=365/5)

model = mbt.ModelParameters(model_length=500/ft, model_width=200/ft, model_time=6*365, dx=2/ft, dy=1/ft, dt=365/5)

Field data¶

Let's plot the field data from observations, 6 years after the spill. Later, we will compare the measurements to the different models.

In [9]:

plt.figure(figsize = [6,3])
plt.plot(field_data_x, field_data_c, color="black",marker = 'o',ms=10, ls = ':')
plt.title('Keesler field observation on contaminant concentration')
plt.show()

plt.figure(figsize = [6,3]) plt.plot(field_data_x, field_data_c, color="black",marker = 'o',ms=10, ls = ':') plt.title('Keesler field observation on contaminant concentration') plt.show()

Models¶

Three models are run with all three options implemented in mibitrans:

• the Bioscreen model class: reproduced the calculations of BIOSCREEN
• the Anatrans model class: uses the untruncated form of the Bioscreen model and include source decay fully analytically, the additional term increases accuracy for small times and distances
• the Mibitrans model class: uses the fully exact solution to the transport situation

All models come with three options to include degredation during plume movement:

• no decay (i.e. no degradation at all)
• instantaneous reaction
• linear decay

No Decay¶

In [10]:

bio_model = mbt.Bioscreen(hydrological_parameters=hydro,
                          attenuation_parameters=att,
                          source_parameters=source,
                          model_parameters=model)
bio_keesler_nodecay = bio_model.run()

bio_model = mbt.Bioscreen(hydrological_parameters=hydro, attenuation_parameters=att, source_parameters=source, model_parameters=model) bio_keesler_nodecay = bio_model.run()

In [11]:

ana_model = mbt.Anatrans(hydro,att,source,model)
ana_keesler_nodecay = ana_model.run()

ana_model = mbt.Anatrans(hydro,att,source,model) ana_keesler_nodecay = ana_model.run()

In [12]:

mbt_model = mbt.Mibitrans(hydro,att,source,model)
mbt_keesler_nodecay =mbt_model.run()

mbt_model = mbt.Mibitrans(hydro,att,source,model) mbt_keesler_nodecay =mbt_model.run()

In [13]:

gamma = mbt_keesler_nodecay.k_source
print("Source depletion rate gamma (1/d): ",gamma*365)

gamma = mbt_keesler_nodecay.k_source print("Source depletion rate gamma (1/d): ",gamma*365)

Source depletion rate gamma (1/d):  0.0016685484141373108

In [14]:

time_point_1 = 365
time_point_2 = 365 * 6
cmap = plt.get_cmap('tab20b')
colors_1 = cmap.colors   # list of RGB tuples

time_point_1 = 365 time_point_2 = 365 * 6 cmap = plt.get_cmap('tab20b') colors_1 = cmap.colors # list of RGB tuples

In [15]:

bio_keesler_nodecay.centerline(time=time_point_1, color=colors_1[0+12],lw = 2,ls = '--',label="Bioscreen $t_1$")
ana_keesler_nodecay.centerline(time=time_point_1, color=colors_1[0+8],lw = 2,ls = '--', label="Anatrans $t_1$")
mbt_keesler_nodecay.centerline(time=time_point_1, color=colors_1[0],lw = 2,ls = '--', label="Mibitrans $t_1$")

bio_keesler_nodecay.centerline(color=colors_1[2+12],lw = 3,label="Bioscreen  $t_2$")
ana_keesler_nodecay.centerline(color=colors_1[2+8],lw = 2,label="Anatrans $t_2$")
mbt_keesler_nodecay.centerline(color=colors_1[2],lw = 2,label="Mibitrans  $t_2$")

plt.xlim((0,100)) # Reduce the x-axis limit to better observe differences
plt.legend(ncols = 2)
plt.title('No decay, $t_1={}y$, $t_2={}$y'.format(time_point_1/365,time_point_2/365))
plt.show()

bio_keesler_nodecay.centerline(time=time_point_1, color=colors_1[0+12],lw = 2,ls = '--',label="Bioscreen $t_1$") ana_keesler_nodecay.centerline(time=time_point_1, color=colors_1[0+8],lw = 2,ls = '--', label="Anatrans $t_1$") mbt_keesler_nodecay.centerline(time=time_point_1, color=colors_1[0],lw = 2,ls = '--', label="Mibitrans $t_1$") bio_keesler_nodecay.centerline(color=colors_1[2+12],lw = 3,label="Bioscreen $t_2$") ana_keesler_nodecay.centerline(color=colors_1[2+8],lw = 2,label="Anatrans $t_2$") mbt_keesler_nodecay.centerline(color=colors_1[2],lw = 2,label="Mibitrans $t_2$") plt.xlim((0,100)) # Reduce the x-axis limit to better observe differences plt.legend(ncols = 2) plt.title('No decay, $t_1={}y$, $t_2={}$y'.format(time_point_1/365,time_point_2/365)) plt.show()

In [16]:

ax = mbt_keesler_nodecay.plume_3d(cmap="vanimo_r")
ax.view_init(elev=15, azim=340)

ax = mbt_keesler_nodecay.plume_3d(cmap="vanimo_r") ax.view_init(elev=15, azim=340)

Interpretation

Comparing the three models at the first time step shows that the Bioscreen model underestimates the source concentrations. The Anatrans and Mibitrans models start at the correct source zone concentration, but show different behaviour along the domain.

Animation

You can easily create animated figures showing the temporal evolution of the plume center location for all three models by uncommenting the following code in the code-cell. Note that animated plots sometimes causes issues with Jupyter Notebook. In such a case, comment the code and restart the kernel of the notebook.

You will see that over time, the source concentration for the Bioscreen model increases to the expected value, and the results resolve to that of the Anatrans model. Differences between the Anatrans and Mibitrans models decrease when plume front is transported away from source.

In [17]:

### Needed to show animations in Jupyter Notebook
# %matplotlib ipympl

# anim = mbt.centerline(
#     [mbt_results, ana_results, bio_results],
#     legend_names=["Mibitrans", "Anatrans", "Bioscreen"],
#     animate=True)
# plt.show()

### Needed to show animations in Jupyter Notebook # %matplotlib ipympl # anim = mbt.centerline( # [mbt_results, ana_results, bio_results], # legend_names=["Mibitrans", "Anatrans", "Bioscreen"], # animate=True) # plt.show()

Instant Reaction Model¶

The no-decay models above do not consider biodegradation. Using the class method instant_reaction allows modelling biodegradation as instantaneous reaction based on electron acceptor concentrations that are available for the field site. Electron acceptor concentrations are translated into a biodegradation capacity (BC) and applied to the model. For more theoretical background on the instant reaction model, see the package documentation.

In [18]:

bio_model.instant_reaction(electron_acceptors=ea_concentrations)
ana_model.instant_reaction(ea_concentrations)
mbt_model.instant_reaction(ea_concentrations)
bio_keesler_instant = bio_model.run()
ana_keesler_instant = ana_model.run()
mbt_keesler_instant = mbt_model.run()

bio_model.instant_reaction(electron_acceptors=ea_concentrations) ana_model.instant_reaction(ea_concentrations) mbt_model.instant_reaction(ea_concentrations) bio_keesler_instant = bio_model.run() ana_keesler_instant = ana_model.run() mbt_keesler_instant = mbt_model.run()

In [19]:

bio_keesler_instant.centerline(time=time_point_1,color=colors_1[0+12],lw = 2,ls = '--', label="Bioscreen $t_1$")
ana_keesler_instant.centerline(time=time_point_1,color=colors_1[0+8],lw = 2,ls = '--', label="Anatrans $t_1$")
mbt_keesler_instant.centerline(time=time_point_1,color=colors_1[0],lw = 2,ls = '--', label="Mibitrans $t_1$")

bio_keesler_instant.centerline(time=time_point_2,color=colors_1[2+12],lw = 3, label="Bioscreen $t_2$")
ana_keesler_instant.centerline(time=time_point_2,color=colors_1[2+8],lw = 2, label="Anatrans $t_2$")
mbt_keesler_instant.centerline(time=time_point_2,color=colors_1[2],lw = 2, label="Mibitrans $t_2$")
plt.xlim((-1,100)) # Reduce the x-axis limit to better observe differences
plt.title('Instant reaction, $t_1={}y$, $t_2={}$y'.format(time_point_1/365,time_point_2/365))
plt.legend(ncol=2)
plt.show()

bio_keesler_instant.centerline(time=time_point_1,color=colors_1[0+12],lw = 2,ls = '--', label="Bioscreen $t_1$") ana_keesler_instant.centerline(time=time_point_1,color=colors_1[0+8],lw = 2,ls = '--', label="Anatrans $t_1$") mbt_keesler_instant.centerline(time=time_point_1,color=colors_1[0],lw = 2,ls = '--', label="Mibitrans $t_1$") bio_keesler_instant.centerline(time=time_point_2,color=colors_1[2+12],lw = 3, label="Bioscreen $t_2$") ana_keesler_instant.centerline(time=time_point_2,color=colors_1[2+8],lw = 2, label="Anatrans $t_2$") mbt_keesler_instant.centerline(time=time_point_2,color=colors_1[2],lw = 2, label="Mibitrans $t_2$") plt.xlim((-1,100)) # Reduce the x-axis limit to better observe differences plt.title('Instant reaction, $t_1={}y$, $t_2={}$y'.format(time_point_1/365,time_point_2/365)) plt.legend(ncol=2) plt.show()

In [20]:

ax = mbt_keesler_instant.plume_3d(cmap="vanimo_r")
ax.view_init(elev=15, azim=340)

ax = mbt_keesler_instant.plume_3d(cmap="vanimo_r") ax.view_init(elev=15, azim=340)

Interpretation

At the plume fringes, electron acceptor concentrations remain high enough to allow for full biodegradation of the contaminant plume.

Animation

Uncomment the following cell, to see results as temporal animation.

In [21]:

# anim = mbt.centerline(
#     [mbt_keesler_instant, ana_keesler_instant, bio_keesler_instant],
#     legend_names=["Mibitrans", "Anatrans", "Bioscreen"],
#     animate=True
# )
# plt.show()

# anim = mbt.centerline( # [mbt_keesler_instant, ana_keesler_instant, bio_keesler_instant], # legend_names=["Mibitrans", "Anatrans", "Bioscreen"], # animate=True # ) # plt.show()

Linear decay¶

The other option for modelling biodegradation is using a linear decay of contaminants during transport. However, direct measurements of a decay rate/half-life time are not available. Consequently, the parameter of solute half-life (i.e. half life time of solute decay during transport) was estimated in BIOSCREEN$^2$ with $0.15$ years based on calibration of the decay model to the field data. We reproduce these results here.

When updating the model input parameters and the model mode, the setup can be rerun using the three implemented models now allowing for linear decay of contaminants during transport.

In [22]:

half_life = 0.15*365
mbt_model.attenuation_parameters.half_life = half_life
ana_model.attenuation_parameters.half_life = half_life
bio_model.attenuation_parameters.half_life = half_life
mbt_model.mode = "linear"
ana_model.mode = "linear"
bio_model.mode = "linear"
bio_keesler_decay = bio_model.run()
ana_keesler_decay = ana_model.run()
mbt_keesler_decay = mbt_model.run()

half_life = 0.15*365 mbt_model.attenuation_parameters.half_life = half_life ana_model.attenuation_parameters.half_life = half_life bio_model.attenuation_parameters.half_life = half_life mbt_model.mode = "linear" ana_model.mode = "linear" bio_model.mode = "linear" bio_keesler_decay = bio_model.run() ana_keesler_decay = ana_model.run() mbt_keesler_decay = mbt_model.run()

In [23]:

bio_keesler_decay.centerline(time=time_point_1,color=colors_1[0+12],lw = 3,ls = '--', label="Bioscreen, $t_1$")
ana_keesler_decay.centerline(time=time_point_1,color=colors_1[0+8],lw = 2,ls = '--', label="Anatrans, $t_1$")
mbt_keesler_decay.centerline(time=time_point_1,color=colors_1[0],lw = 2,ls = '--', label="Mibitrans, $t_1$")

bio_keesler_decay.centerline(time=time_point_2,color=colors_1[2+12],lw = 3, label="Bioscreen, $t_2$")
ana_keesler_decay.centerline(time=time_point_2,color=colors_1[2+8],lw = 2, label="Anatrans, $t_2$")
mbt_keesler_decay.centerline(time=time_point_2,color=colors_1[2],lw = 2, label="Mibitrans, $t_2$")
plt.xlim(0,50)
plt.title('Linear decay, $t_1={}y$, $t_2={}$y'.format(time_point_1/365,time_point_2/365))
plt.legend(ncols = 2)
plt.show()

bio_keesler_decay.centerline(time=time_point_1,color=colors_1[0+12],lw = 3,ls = '--', label="Bioscreen, $t_1$") ana_keesler_decay.centerline(time=time_point_1,color=colors_1[0+8],lw = 2,ls = '--', label="Anatrans, $t_1$") mbt_keesler_decay.centerline(time=time_point_1,color=colors_1[0],lw = 2,ls = '--', label="Mibitrans, $t_1$") bio_keesler_decay.centerline(time=time_point_2,color=colors_1[2+12],lw = 3, label="Bioscreen, $t_2$") ana_keesler_decay.centerline(time=time_point_2,color=colors_1[2+8],lw = 2, label="Anatrans, $t_2$") mbt_keesler_decay.centerline(time=time_point_2,color=colors_1[2],lw = 2, label="Mibitrans, $t_2$") plt.xlim(0,50) plt.title('Linear decay, $t_1={}y$, $t_2={}$y'.format(time_point_1/365,time_point_2/365)) plt.legend(ncols = 2) plt.show()

In [24]:

ax = mbt_keesler_decay.plume_3d(cmap="vanimo_r")
ax.view_init(elev=15, azim=340)

ax = mbt_keesler_decay.plume_3d(cmap="vanimo_r") ax.view_init(elev=15, azim=340)

Comparison with field data¶

BIOSCREEN$^2$ was used to try to reproduce the movement of the plume from 1989 to 1995. The instantaneous reaction model was used as the primary model to try to reproduce the plume length of about $280$ ft. BIOSCREEN$^2$ reports that dispersivity has been calibrated for a better fit of the instant reaction model: "The initial run of the instantaneous reaction model indicated that the plume was too long. This indicates that there is more mixing of hydrocarbon and electron acceptors at the site than is predicted by the model. Therefore the longitudinal dispersivity was adjusted upwards (more mixing) until BIOSCREEN matched the observed plume length. The final longitudinal dispersivity was 32.5 ft."

Prediction of field behaviour for no-decay and instant reaction model is off from the field data. The initial run of the instantaneous reaction model (with smaller dispersivity values) indicated that the plume was too long and mixing of hydrocarbon and electron acceptors at the site is higher than assumed. Dispersivity values have been adapted to match plume length. But still, the instantaneous reaction model does not fit well the field observations, especially at $20$ and $60$m distance. Consequently, a different model for degradation behaviour is indicated.

In [25]:

bio_keesler_nodecay.centerline(time=time_point_2,color=colors_1[2+12],lw = 3,ls = '-', label="Bioscreen, no decay")
ana_keesler_nodecay.centerline(time=time_point_2,color=colors_1[2+8],lw = 2,ls = '-', label="Anatrans, no decay")
mbt_keesler_nodecay.centerline(time=time_point_2,color=colors_1[2],lw = 2,ls = '-', label="Mibitrans, no decay")

bio_keesler_instant.centerline(time=time_point_2,color=colors_1[2+12],lw = 3,ls = '--', label="Bioscreen, Inst react")
ana_keesler_instant.centerline(time=time_point_2,color=colors_1[2+8],lw = 2,ls = '--', label="Anatrans, Inst react")
mbt_keesler_instant.centerline(time=time_point_2,color=colors_1[2],lw = 2,ls = '--',  label="Mibitrans, Inst react")

bio_keesler_decay.centerline(time=time_point_2,color=colors_1[2+12],lw = 3,ls = ':', label="Bioscreen, lin decay")
ana_keesler_decay.centerline(time=time_point_2,color=colors_1[2+8],lw = 2,ls = ':', label="Anatrans, lin decay")
mbt_keesler_decay.centerline(time=time_point_2,color=colors_1[2],lw = 2,ls = ':',  label="Mibitrans, lin decay")
plt.scatter(field_data_x, field_data_c, label="field data", color="black")

plt.title('Model comparison at $t={}$y'.format(time_point_2/365))
plt.xlim(-1,100)
plt.ylim(-0.2,13)
plt.legend(ncols = 2)
plt.show()

bio_keesler_nodecay.centerline(time=time_point_2,color=colors_1[2+12],lw = 3,ls = '-', label="Bioscreen, no decay") ana_keesler_nodecay.centerline(time=time_point_2,color=colors_1[2+8],lw = 2,ls = '-', label="Anatrans, no decay") mbt_keesler_nodecay.centerline(time=time_point_2,color=colors_1[2],lw = 2,ls = '-', label="Mibitrans, no decay") bio_keesler_instant.centerline(time=time_point_2,color=colors_1[2+12],lw = 3,ls = '--', label="Bioscreen, Inst react") ana_keesler_instant.centerline(time=time_point_2,color=colors_1[2+8],lw = 2,ls = '--', label="Anatrans, Inst react") mbt_keesler_instant.centerline(time=time_point_2,color=colors_1[2],lw = 2,ls = '--', label="Mibitrans, Inst react") bio_keesler_decay.centerline(time=time_point_2,color=colors_1[2+12],lw = 3,ls = ':', label="Bioscreen, lin decay") ana_keesler_decay.centerline(time=time_point_2,color=colors_1[2+8],lw = 2,ls = ':', label="Anatrans, lin decay") mbt_keesler_decay.centerline(time=time_point_2,color=colors_1[2],lw = 2,ls = ':', label="Mibitrans, lin decay") plt.scatter(field_data_x, field_data_c, label="field data", color="black") plt.title('Model comparison at $t={}$y'.format(time_point_2/365)) plt.xlim(-1,100) plt.ylim(-0.2,13) plt.legend(ncols = 2) plt.show()

Interpretation

Prediction of field behaviour for no-decay and instant reaction model is off from the field data. The initial run of the instantaneous reaction model (with smaller dispersivity values) indicated that the plume was too long and mixing of hydrocarbon and electron acceptors at the site is higher than assumed. Dispersivity values have been adapted to match plume length. But still, the instantaneous reaction model does not fit well the field observations, especially at $20$ and $60$m distance. Consequently, a different model for degradation behaviour is indicated.

Note the overlap of the Bioscreen and the Anatrans model which is expected at the large Peclet number given this field sitation.

Summary of Results

Condensing to the crucial outcome of the figure (note the changed color scheme to highlight the differences in degradation model choice):

In [26]:

mbt_keesler_decay.centerline(time=time_point_2, label="Mbt, linear decay", color="C0", linestyle="-",lw = 2)
mbt_keesler_instant.centerline(time=time_point_2, label="Mbt, inst reaction", color="C2", linestyle="-",lw = 2)
mbt_keesler_nodecay.centerline(time=time_point_2, linestyle="-", label="Mbt, no decay", color="C1",lw = 2)
bio_keesler_decay.centerline(time=time_point_2, label="Bioscreen", color="C0", linestyle="--",lw = 2)
bio_keesler_instant.centerline(time=time_point_2, label="Bioscreen", color="C2", linestyle="--",lw = 2)
plt.scatter(field_data_x, field_data_c, linestyle="--", label="field data", color="black",s = 50,zorder = 5)

plt.title('Model comparison at $t={}$y'.format(time_point_2/365))
plt.xlim(-2,100)
plt.ylim(-0.2,13)
plt.legend(ncol=2)

mbt_keesler_decay.centerline(time=time_point_2, label="Mbt, linear decay", color="C0", linestyle="-",lw = 2) mbt_keesler_instant.centerline(time=time_point_2, label="Mbt, inst reaction", color="C2", linestyle="-",lw = 2) mbt_keesler_nodecay.centerline(time=time_point_2, linestyle="-", label="Mbt, no decay", color="C1",lw = 2) bio_keesler_decay.centerline(time=time_point_2, label="Bioscreen", color="C0", linestyle="--",lw = 2) bio_keesler_instant.centerline(time=time_point_2, label="Bioscreen", color="C2", linestyle="--",lw = 2) plt.scatter(field_data_x, field_data_c, linestyle="--", label="field data", color="black",s = 50,zorder = 5) plt.title('Model comparison at $t={}$y'.format(time_point_2/365)) plt.xlim(-2,100) plt.ylim(-0.2,13) plt.legend(ncol=2)

Out[26]:

<matplotlib.legend.Legend at 0x7fc97482d450>

Impact of dispersivity¶

BIOSCREEN$^2$ reports that dispersivity has been calibrated for a better fit of the instant reaction model.: "The initial run of the instantaneous reaction model indicated that the plume was too long. This indicates that there is more mixing of hydrocarbon and electron acceptors at the site than is predicted by the model. Therefore the longitudinal dispersivity was adjusted upwards (more mixing) until BIOSCREEN matched the observed plume length. The final longitudinal dispersivity was 32.5 ft."

However, we saw that the choice of degradation model has a significant inpact of model fit to observation data. The linear decay model matches observation much better. Furthermore, dispersivities of 32.5ft, 3.25ft and 0ft, are at the higher end for observed reliable longitudinal macrodispersivities and beyond any reliable transverse dispersivity values$^{4,5}$.

Find below a comparison of the model with the two sets of hydrogeological parameters, we redefine the respective data input class and create a new model instance.

In [27]:

### define new set of hydrological parameters with alternative dispersivities
hydro_disp = mbt.HydrologicalParameters(
    h_conductivity = 0.011 / 100 * 3600 * 24, #m/d
    h_gradient = 0.003, # m/m
    porosity = 0.3,
    alpha_x = 13.3 /ft, #initially used dispersivity (before calibration)
    alpha_y = 1.3/ft, #initially used dispersivity (before calibration)
    alpha_z = 0
)

### define new set of hydrological parameters with alternative dispersivities hydro_disp = mbt.HydrologicalParameters( h_conductivity = 0.011 / 100 * 3600 * 24, #m/d h_gradient = 0.003, # m/m porosity = 0.3, alpha_x = 13.3 /ft, #initially used dispersivity (before calibration) alpha_y = 1.3/ft, #initially used dispersivity (before calibration) alpha_z = 0 )

Rerun models with alternative dispersivities

In [28]:

mbt_model.hydrological_parameters = hydro_disp
mbt_keesler_decay_disp = mbt_model.run()

mbt_model.instant_reaction(ea_concentrations)
mbt_keesler_instant_disp = mbt_model.run()

mbt_model = mbt.Mibitrans(hydro_disp,att,source,model)
mbt_keesler_nodecay_disp = mbt_model.run()

mbt_model.hydrological_parameters = hydro_disp mbt_keesler_decay_disp = mbt_model.run() mbt_model.instant_reaction(ea_concentrations) mbt_keesler_instant_disp = mbt_model.run() mbt_model = mbt.Mibitrans(hydro_disp,att,source,model) mbt_keesler_nodecay_disp = mbt_model.run()

In [29]:

mbt_keesler_decay.centerline(time=time_point_2, label="Linear decay: Mbt", color="C0", linestyle="-",lw = 2)
mbt_keesler_instant.centerline(time=time_point_2, label="Inst reaction: Mbt", color="C2", linestyle="-",lw = 2)
mbt_keesler_nodecay.centerline(time=time_point_2, label="No decay: Mbt", color="C1",lw = 2)

mbt_keesler_decay_disp.centerline(time=time_point_2, label=r"Mbt$^*$", color="C0", linestyle=":",lw = 2)
mbt_keesler_instant_disp.centerline(time=time_point_2, label=r"Mbt$^*$", color="C2", linestyle=":",lw = 2)
mbt_keesler_nodecay_disp.centerline(time=time_point_2, label=r"Mbt$^*$", color="C1", linestyle=":",lw = 2)

bio_keesler_decay.centerline(time=time_point_2, label="Bioscreen", color="C0", linestyle="--",lw = 2)
bio_keesler_instant.centerline(time=time_point_2, label="Bioscreen", color="C2", linestyle="--",lw = 2)
plt.scatter(field_data_x, field_data_c, linestyle="--", label="field data", color="black",s = 50,zorder = 5)

plt.title(r'Effect of dispersivity$^*$ ($\alpha_L=13.3$ft, $\alpha_T=1.3$ft), $t={}$y'.format(time_point_2/365))
plt.legend(ncol = 3)
plt.xlim(-1,90)

mbt_keesler_decay.centerline(time=time_point_2, label="Linear decay: Mbt", color="C0", linestyle="-",lw = 2) mbt_keesler_instant.centerline(time=time_point_2, label="Inst reaction: Mbt", color="C2", linestyle="-",lw = 2) mbt_keesler_nodecay.centerline(time=time_point_2, label="No decay: Mbt", color="C1",lw = 2) mbt_keesler_decay_disp.centerline(time=time_point_2, label=r"Mbt$^*$", color="C0", linestyle=":",lw = 2) mbt_keesler_instant_disp.centerline(time=time_point_2, label=r"Mbt$^*$", color="C2", linestyle=":",lw = 2) mbt_keesler_nodecay_disp.centerline(time=time_point_2, label=r"Mbt$^*$", color="C1", linestyle=":",lw = 2) bio_keesler_decay.centerline(time=time_point_2, label="Bioscreen", color="C0", linestyle="--",lw = 2) bio_keesler_instant.centerline(time=time_point_2, label="Bioscreen", color="C2", linestyle="--",lw = 2) plt.scatter(field_data_x, field_data_c, linestyle="--", label="field data", color="black",s = 50,zorder = 5) plt.title(r'Effect of dispersivity$^*$ ($\alpha_L=13.3$ft, $\alpha_T=1.3$ft), $t={}$y'.format(time_point_2/365)) plt.legend(ncol = 3) plt.xlim(-1,90)

Out[29]:

(-1.0, 90.0)

Interpration

We see significant differences of the effect of dispersivities for the different degradation models. The results of the instant reaction model change significantly with increasing dispersivity. The match to plume length is improved. However, the instant reaction model still does not fit well to concentration observation closer to the source.

Changing the dispersivity in the linear decay model shows very little change in the modelled results. While differences between the Mibitrans and the Bioscreen model increase, the overall effect of increased dispersivity is small.

They are even smaller then differences based on the model choice, i.e. Mibitrans model representing the exact solution to the transport problem and the Bioscreen model that contains equation truncation and simplifications. The first matches the concentration at about $10$m better, while the latter is a better fit to the concentration around $20$m. Given measurement uncertainty, both models provide a reasonable fit under the use of a linear decay model with the fitted degradation rate.