Comparison with BIOSCREEN
import matplotlib.pyplot as plt
import numpy as np
import mibitrans as mbt
ft = 3.281 # factor to convert ft to m
This notebook compares the examples present in BIOSCREEN and BIOSCREEN-AT with the similar models (Bioscreen/Anatrans and Mibitrans respectively) in the mibitrans package.
BIOSCREEN example dataset¶
bio_hydro = mbt.HydrologicalParameters(
velocity=1609.1 / ft / 365, #[m/d]
porosity=0.25, #[-]
alpha_x=28.5 / ft, #[m]
alpha_y=2.9 / ft, #[m]
alpha_z=0 / ft, #[m]
)
bio_att = mbt.AttenuationParameters(
bulk_density=1700, #[kg/m3]
partition_coefficient=0.038, #[m3/kg]
fraction_organic_carbon=0.0008, #[-]
half_life=0, #[d]
)
bio_source = mbt.SourceParameters(
source_zone_boundary=np.array([50, 75, 125]) / ft, #[m]
source_zone_concentration=np.array([9, 2.8, 0.07]), #[g/m3]
depth=10 / ft, #[m]
total_mass="inf", #[g]
)
bio_model = mbt.ModelParameters(
model_length=1500 / ft, #[m]
model_width=400 / ft, #[m]
model_time=5 * 365, #[d]
dx=5 / ft, #[m]
dy=5 / ft, #[m]
dt=365/2, #[d]
)
bio_ea = mbt.ElectronAcceptors(
delta_oxygen=5.78, #[g/m3]
delta_nitrate=17, #[g/m3]
ferrous_iron=11.3, #[g/m3]
delta_sulfate=100, #[g/m3]
methane=0.414, #[g/m3]
)
#bio_bio; BIOSCREEN input with Bioscreen model
bio_bio_object = mbt.Bioscreen(bio_hydro, bio_att, bio_source, bio_model)
#nd; no decay
bio_bio_nd_results = bio_bio_object.run()
bio_bio_object.attenuation_parameters.half_life = 0.1 * 365 #[d]
#dc; decay
bio_bio_dc_results = bio_bio_object.run()
bio_bio_object.instant_reaction(bio_ea)
#inst; instant reaction
bio_bio_inst_results = bio_bio_object.run()
#bio_ana; BIOSCREEN input with Anatrans model
bio_ana_object = mbt.Anatrans(bio_hydro, bio_att, bio_source, bio_model)
bio_ana_nd_results = bio_ana_object.run()
bio_ana_object.attenuation_parameters.half_life = 0.1 * 365 #[d]
bio_ana_dc_results = bio_ana_object.run()
bio_ana_object.instant_reaction(bio_ea)
bio_ana_inst_results = bio_ana_object.run()
time_bio = [365*0.5, 365, 2*365, 5*365] #[d]
BIOSCREEN_data_nodecay = np.array([
[9.000, 8.435, 7.267, 5.900, 4.091, 2.151, 0.780, 0.183, 0.027, 0.002, 0.000], # 0.5y
[9.000, 8.467, 7.465, 6.680, 6.070, 5.548, 5.006, 4.327, 3.441, 2.410, 1.432], # 1y
[9.000, 8.467, 7.466, 6.684, 6.089, 5.624, 5.250, 4.940, 4.679, 4.454, 4.256], # 2y
[9.000, 8.467, 7.466, 6.684, 6.089, 5.624, 5.250, 4.940, 4.679, 4.455, 4.260] # 5y
])
BIOSCREEN_data_lineardecay = np.array([
[9.000, 4.348, 1.965, 0.889, 0.386, 0.146, 0.043, 0.009, 0.001, 0.000, 0.000],
[9.000, 4.348, 1.969, 0.905, 0.424, 0.201, 0.096, 0.046, 0.022, 0.010, 0.004],
[9.000, 4.348, 1.969, 0.905, 0.424, 0.201, 0.096, 0.047, 0.023, 0.011, 0.005],
[9.000, 4.348, 1.969, 0.905, 0.424, 0.201, 0.096, 0.047, 0.023, 0.011, 0.005]
])
BIOSCREEN_data_instant = np.array([
[9.000, 8.332, 6.473, 2.364, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000],
[9.000, 8.466, 7.405, 6.331, 5.165, 3.757, 1.724, 0.000, 0.000, 0.000, 0.000],
[9.000, 8.466, 7.407, 6.350, 5.268, 4.192, 3.152, 2.167, 1.244, 0.378, 0.000],
[9.000, 8.466, 7.407, 6.350, 5.268, 4.192, 3.152, 2.168, 1.245, 0.385, 0.000]
])
bio_x = np.array([0, 145, 290, 435, 580, 725, 870, 1015, 1160, 1305, 1450]) / ft #[m]
colors = ["tomato", "sandybrown", "khaki", "greenyellow", "lightgreen"]
colors_bio = ["red", "orange", "gold", "yellowgreen", "green"]
for i in range(len(time_bio)):
mbt.centerline(bio_bio_nd_results, time=time_bio[i], color=colors[i], label=f"Bioscreen, t={time_bio[i]}d")
plt.plot(bio_x, BIOSCREEN_data_nodecay[i, :], color=colors_bio[i],
linestyle=":", label=f"BIOSCREEN, t={time_bio[i]}d")
plt.title("Comparison of BIOSCREEN and Bioscreen models without decay")
plt.legend()
plt.show()
for i in range(len(time_bio)):
mbt.centerline(bio_bio_dc_results, time=time_bio[i], color=colors[i], label=f"Bioscreen, t={time_bio[i]}d")
plt.plot(bio_x, BIOSCREEN_data_lineardecay[i, :], color=colors_bio[i],
linestyle=":", label=f"BIOSCREEN, t={time_bio[i]}d")
plt.title("Comparison of BIOSCREEN and Bioscreen models with linear decay")
plt.legend()
plt.show()
for i in range(len(time_bio)):
mbt.centerline(bio_bio_inst_results, time=time_bio[i], color=colors[i], label=f"Bioscreen, t={time_bio[i]}d")
plt.plot(bio_x, BIOSCREEN_data_instant[i, :], color=colors_bio[i],
linestyle=":", label=f"BIOSCREEN, t={time_bio[i]}d")
plt.title("Comparison of BIOSCREEN and Bioscreen models with instant reaction")
plt.legend()
plt.show()
colors = ["tomato", "sandybrown", "khaki", "greenyellow", "lightgreen"]
colors_bio = ["red", "orange", "gold", "yellowgreen", "green"]
for i in range(len(time_bio)):
mbt.centerline(bio_ana_nd_results, time=time_bio[i], color=colors[i], label=f"Anatrans, t={time_bio[i]}d")
plt.plot(bio_x, BIOSCREEN_data_nodecay[i, :], color=colors_bio[i],
linestyle=":", label=f"BIOSCREEN, t={time_bio[i]}d")
plt.title("Comparison of BIOSCREEN and Anatrans models without decay")
plt.legend()
plt.show()
for i in range(len(time_bio)):
mbt.centerline(bio_ana_dc_results, time=time_bio[i], color=colors[i], label=f"Anatrans, t={time_bio[i]}")
plt.plot(bio_x, BIOSCREEN_data_lineardecay[i, :], color=colors_bio[i],
linestyle=":", label=f"BIOSCREEN, t={time_bio[i]}d")
plt.title("Comparison of BIOSCREEN and Anatrans models with linear decay")
plt.legend()
plt.show()
for i in range(len(time_bio)):
mbt.centerline(bio_ana_inst_results, time=time_bio[i], color=colors[i], label=f"Anatrans, t={time_bio[i]}")
plt.plot(bio_x, BIOSCREEN_data_instant[i, :], color=colors_bio[i],
linestyle=":", label=f"BIOSCREEN, t={time_bio[i]}d")
plt.title("Comparison of BIOSCREEN and Anatrans models with instant reaction")
plt.legend()
plt.show()
Input parameters¶
Input parameters are the default parameters when opening up BIOSCREEN-AT. All imperial units are converted to consistent metric units.
bat_hydro = mbt.HydrologicalParameters(
velocity=335.2 / ft / 365, #[m/d]
porosity=0.25, #[-]
alpha_x=28.887 / ft, #[m]
alpha_y=2.889 / ft, #[m]
alpha_z=0.289 / ft, #[m]
)
bat_att = mbt.AttenuationParameters(
retardation=1.20672, #[-]
half_life=0, #[d]
)
bat_source = mbt.SourceParameters(
source_zone_boundary=np.array([50 / ft]), #[m]
source_zone_concentration=np.array([9]), #[g/m3]
depth=10 / ft, #[m]
total_mass="inf", #[g]
)
bat_model = mbt.ModelParameters(
model_length=2500 / ft, #[m]
model_width=500 / ft, #[m]
model_time=5 * 365, #[d]
dx=20 / ft, #[m]
dy=5 / ft, #[m]
dt=365, #[d]
)
bat_ea = mbt.ElectronAcceptors(
delta_oxygen=5.78, #[g/m3]
delta_nitrate=17, #[g/m3]
ferrous_iron=11.3, #[g/m3]
delta_sulfate=100, #[g/m3]
methane=0.414, #[g/m3]
)
Running model¶
As the exact solution uses an integral, run time is longer than that of the Bioscreen solution, depending on model resolution. Using verbose when running the model is recommended, to keep track of the calculation process.
# bat_mbt; BIOSCREEN-AT data with Mibitrans model
bat_mbt_object = mbt.Mibitrans(
hydrological_parameters=bat_hydro,
attenuation_parameters=bat_att,
source_parameters=bat_source,
model_parameters=bat_model,
)
bat_mbt_nd_results = bat_mbt_object.run()
bat_mbt_object.attenuation_parameters.half_life = 365 * 10 #[d]
bat_mbt_dc_results = bat_mbt_object.run()
bat_mbt_object.instant_reaction(electron_acceptors=bat_ea)
bat_mbt_inst_results = bat_mbt_object.run()
Visualize¶
bat_mbt_nd_results.centerline(color="red", label="No decay")
bat_mbt_dc_results.centerline(color="blue", label="Linear decay")
bat_mbt_inst_results.centerline(color="green", label="Instant reaction")
plt.legend()
plt.show()
Running Bioscreen solution for comparison
Comparing to BIOSCREEN-AT¶
bat_x = np.array([0, 250, 500, 750, 1000, 1250, 1500, 1750, 2000, 2250, 2500]) / ft #[m]
BIOSCREENAT_data_nodecay = np.array([
[9.000, 3.635, 0.236, 0.001, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000],
[9.000, 4.649, 2.144, 0.452, 0.022, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000],
[9.000, 4.696, 2.707, 1.505, 0.493, 0.064, 0.003, 0.000, 0.000, 0.000, 0.000],
[9.000, 4.699, 2.767, 1.872, 1.161, 0.483, 0.104, 0.010, 0.000, 0.000, 0.000],
[9.000, 4.699, 2.772, 1.932, 1.423, 0.947, 0.459, 0.134, 0.021, 0.002, 0.000],
])
BIOSCREENAT_data_lineardecay = np.array([
[9.000, 3.479, 0.222, 0.001, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000],
[9.000, 4.407, 1.944, 0.401, 0.019, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000],
[9.000, 4.448, 2.423, 1.289, 0.412, 0.053, 0.002, 0.000, 0.000, 0.000, 0.000],
[9.000, 4.450, 2.471, 1.580, 0.938, 0.380, 0.081, 0.008, 0.000, 0.000, 0.000],
[9.000, 4.450, 2.475, 1.624, 1.131, 0.721, 0.340, 0.097, 0.015, 0.001, 0.000],
])
colors = ["tomato", "sandybrown", "khaki", "greenyellow", "lightgreen"]
colors_bio = ["red", "orange", "gold", "yellowgreen", "green"]
for i in range(len(bat_mbt_nd_results.t)):
mbt.centerline(bat_mbt_nd_results, time=bat_mbt_object.t[i],
color=colors[i], label=f"Mibitrans, t={bat_mbt_object.t[i]}")
plt.plot(bat_x, BIOSCREENAT_data_nodecay[i, :], color=colors_bio[i],
linestyle=":", label=f"bio_at, t={bat_mbt_object.t[i]}")
plt.legend()
plt.show()
for i in range(len(bat_mbt_dc_results.t)):
mbt.centerline(bat_mbt_dc_results, time=bat_mbt_object.t[i],
color=colors[i], label=f"Mibitrans, t={bat_mbt_object.t[i]}")
plt.plot(bat_x, BIOSCREENAT_data_lineardecay[i, :], color=colors_bio[i],
linestyle=":", label=f"bio_at, t={bat_mbt_object.t[i]}")
plt.legend()
plt.show()
mibitrans output corresponds very well with the BIOSCREEN-AT output, and has higher resolution. Note that despite BIOSCREEN-AT taking in instant reaction parameters, it does not have the option to calculate biodegradation this way, therefore, instant reaction model can not be compared.