Benchmarking to BIOSCREEN¶

This notebook provides a comparison of concentrations calculated with the mibitrans package and concentrations calculated with the Excel-based distributions BIOSCREEN$^{1,2}$ and BIOSCREEN-AT$^{3}$. The notebooks benchmarks mibitrans showing that equation implementations are correct.

Authors: Jorrit Bakker, Alraune Zech

Background¶

BIOSCREEN-AT$^{3}$ implements the exact solution of the transport equations for the general setting of the package, following the solution of Wexler, 1992$^4$. Calculations require numerical integration. The model class Mibitrans of the mibitrans package implements the same exact equation in python for contaminant transport modelling.

BIOSCREEN$^{1,2}$ implements an approximate, but fully closed-form, equation as solution of the contaminant transport, based on the empirical model of Domenico, 1987$^5$. West et al., 2007$^6$ discussed in detail the approximations and their effects. The same equation is implemented in the mibitrans package as Bioscreen model class. Mostly for comparison with BIOSCREEN.

BIOSCREEN example dataset and BIOSCREEN results¶

We use for comparison the set of transport and model parameters, those provided as BIOSCREEN example dataset in BIOSCREEN. We ran both, BIOSCREEN and BIOSCREEN-AT in Excel to retreive concentrations for comparison. Values are listed below in arrays.

References¶

$^1$ Newell, C. J., R. K. McLeod, and J. R. Gonzales (1996), BIOSCREEN Natural Attenuation Decision Support System User’s Manual Version 1.3, Tech. Rep. EPA/600/R-96/087, US EPA

$^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$Karanovic, M., C. J. Neville, and C. B. Andrews,(2007), BIOSCREEN-AT: BIOSCREEN with an exact analytical solution, Groundwater, 45 (2), 242–245.

$^4$ Wexler, E. J., (1992), Analytical solutions for one-, two-, and three-dimensional solute transport in ground-water systems with uniform flow, Tech. Rep. 03-B7, U.S. G.P.O.

$^5$ Domenico, P. A., 1987, An analytical model for multidimensional transport of a decaying contaminant species, Journal of Hydrology, 91 (1), 49–58.

$^6$ West, M. R., B. H. Kueper, and M. J. Ungs, (2007), On the use and error of approximation in the Domenico (1987) solution, Groundwater, 45 (2), 126–135.

In [1]:

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

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

In [2]:

ft = 3.281 # factor to convert ft to m

ft = 3.281 # factor to convert ft to m

Specify Data and Settings for Benchmarking¶

Define BIOSCREEN example dataset¶

In [3]:

bio_hydro = mbt.HydrologicalParameters(
    velocity=1609.081 / ft / 365, #[m/d]
    porosity=0.25, #[-]
    alpha_x=28.50253 / ft, #[m]
    alpha_y=2.85025 / 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.0000000001, 125.0000000001]) / 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=1450 / ft, #[m]
    model_width=300 / 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_hydro = mbt.HydrologicalParameters( velocity=1609.081 / ft / 365, #[m/d] porosity=0.25, #[-] alpha_x=28.50253 / ft, #[m] alpha_y=2.85025 / 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.0000000001, 125.0000000001]) / 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=1450 / ft, #[m] model_width=300 / 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] )

Load simulation results from BIOSCREEN¶

We ran BIOSCREEN with the specified parameters and stored calculated values in a json-file (mibitrans.data.example_data.json). It is accessed through the BioscreenData class. We define the observation points in $x$, $y$ and $t$ as separate arrays.

In [4]:

bioscreen_data = mbt.BioscreenData()

bioscreen_data = mbt.BioscreenData()

In [5]:

bio_x = np.array([0,145,290,435,580,725,870,1015,1160,1305,1450])/ft
bio_y = np.array([-150,-75,0,75,150])/ft
bio_t = np.array([0.5*365, 365, 2*365, 5*365])

bio_x = np.array([0,145,290,435,580,725,870,1015,1160,1305,1450])/ft bio_y = np.array([-150,-75,0,75,150])/ft bio_t = np.array([0.5*365, 365, 2*365, 5*365])

BIOSCREEN-AT example dataset¶

Input parameters are the default parameters when opening up BIOSCREEN-AT. All imperial units are converted to consistent metric units.

In [6]:

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=1500 / ft, #[m]
    model_width=300 / ft, #[m]
    model_time=5 * 365, #[d]
    dx=5 / 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]
)

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=1500 / ft, #[m] model_width=300 / ft, #[m] model_time=5 * 365, #[d] dx=5 / 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] )

Load Simulation results from BIOSCREEN-AT¶

We ran BIOSCREEN with the specified parameters and stored calculated values in a json-file (mibitrans.data.example_data.json). It is accessed through the BioscreenATData class. We define the observation points in $x$, $y$ and $t$ as separate arrays.

In [7]:

bioscreenAT_data = mbt.BioscreenATData()

bioscreenAT_data = mbt.BioscreenATData()

In [8]:

bat_x = np.array([0, 150, 300, 450, 600, 750, 900, 1050, 1200, 1350, 1500])/ft
bat_y = np.array([-150,-75,0,75,150])/ft
bat_t = np.array([365, 2*365, 3*365, 4*365, 5*365])

bat_x = np.array([0, 150, 300, 450, 600, 750, 900, 1050, 1200, 1350, 1500])/ft bat_y = np.array([-150,-75,0,75,150])/ft bat_t = np.array([365, 2*365, 3*365, 4*365, 5*365])

Comparison with BIOSCREEN¶

Making use of Example data of BIOSCREEN.

BIOSCREEN vs Bioscreen-model class of mibitrans¶

Comparison for all three options of biodegradation:

• no decay
• linear decay
• instantaneous reaction

In [9]:

#bio_bio; Bioscreen model with BIOSCREEN input
bio_bio_object = mbt.Bioscreen(bio_hydro, bio_att, bio_source, bio_model)

#nd; no decay
bio_bio_nd_results = bio_bio_object.run()

#dc; decay
bio_bio_object.attenuation_parameters.half_life = 0.1 * 365 #[d]
bio_bio_dc_results = bio_bio_object.run()

#inst; instant reaction
bio_bio_object.instant_reaction(bio_ea)
bio_bio_inst_results = bio_bio_object.run()

#bio_bio; Bioscreen model with BIOSCREEN input bio_bio_object = mbt.Bioscreen(bio_hydro, bio_att, bio_source, bio_model) #nd; no decay bio_bio_nd_results = bio_bio_object.run() #dc; decay bio_bio_object.attenuation_parameters.half_life = 0.1 * 365 #[d] bio_bio_dc_results = bio_bio_object.run() #inst; instant reaction bio_bio_object.instant_reaction(bio_ea) bio_bio_inst_results = bio_bio_object.run()

Note that BIOSCREEN only provides limited resolution which causes visual differences in the results, particularly for instantaneous reaction model. Actually calculated values (indicated with markers) are identical.

In [10]:

cmap = plt.get_cmap('Paired')
colors = cmap.colors[::2]
colors_bio = cmap.colors[1::2]

cmap = plt.get_cmap('Paired') colors = cmap.colors[::2] colors_bio = cmap.colors[1::2]

In [11]:

for i,ti in enumerate(bio_t):
    mbt.centerline(bio_bio_nd_results, time=ti,
                   color=colors[i],lw = 4,
                   label="Bioscreen")
for i,ti in enumerate(bio_t):
    plt.scatter(bio_x, bioscreen_data.nodecay[i, 2, :],
             color=colors_bio[i], marker='+', s=60, zorder=3,
             label=f"BIOSCREEN, t={ti/365}y")

plt.title("BIOSCREEN vs. Bioscreen class without decay")
plt.legend(ncols=2, loc='upper right', fontsize=9)
plt.ylim([0,11])
plt.show()

for i,ti in enumerate(bio_t): mbt.centerline(bio_bio_nd_results, time=ti, color=colors[i],lw = 4, label="Bioscreen") for i,ti in enumerate(bio_t): plt.scatter(bio_x, bioscreen_data.nodecay[i, 2, :], color=colors_bio[i], marker='+', s=60, zorder=3, label=f"BIOSCREEN, t={ti/365}y") plt.title("BIOSCREEN vs. Bioscreen class without decay") plt.legend(ncols=2, loc='upper right', fontsize=9) plt.ylim([0,11]) plt.show()

In [12]:

for i,ti in enumerate(bio_t):
    mbt.centerline(bio_bio_dc_results, time=ti,lw = 4,
                   color=colors[i], label="Bioscreen")

for i,ti in enumerate(bio_t):
    plt.scatter(bio_x, bioscreen_data.lineardecay[i, 2, :],
             color=colors_bio[i], marker = '+', s=60, zorder=3,
             label=f"BIOSCREEN, t={ti/365}y")

plt.title("BIOSCREEN vs Bioscreen class with linear decay")
plt.legend(ncols = 2,loc = 'upper right')
plt.show()

for i,ti in enumerate(bio_t): mbt.centerline(bio_bio_dc_results, time=ti,lw = 4, color=colors[i], label="Bioscreen") for i,ti in enumerate(bio_t): plt.scatter(bio_x, bioscreen_data.lineardecay[i, 2, :], color=colors_bio[i], marker = '+', s=60, zorder=3, label=f"BIOSCREEN, t={ti/365}y") plt.title("BIOSCREEN vs Bioscreen class with linear decay") plt.legend(ncols = 2,loc = 'upper right') plt.show()

In [13]:

for i,ti in enumerate(bio_t):
    mbt.centerline(bio_bio_inst_results, time=ti,
                   color=colors[i],lw = 4,
                   label="Bioscreen")

for i,ti in enumerate(bio_t):
    plt.scatter(bio_x, bioscreen_data.instant[i, 2, :],
             color=colors_bio[i], marker = '+', s=60, zorder=3,
             label=f"BIOSCREEN, t={ti/365}y")

plt.title("BIOSCREEN vs Bioscreen class with instant reaction")
plt.legend(ncols = 2,loc = 'upper right',fontsize=9)
plt.show()

for i,ti in enumerate(bio_t): mbt.centerline(bio_bio_inst_results, time=ti, color=colors[i],lw = 4, label="Bioscreen") for i,ti in enumerate(bio_t): plt.scatter(bio_x, bioscreen_data.instant[i, 2, :], color=colors_bio[i], marker = '+', s=60, zorder=3, label=f"BIOSCREEN, t={ti/365}y") plt.title("BIOSCREEN vs Bioscreen class with instant reaction") plt.legend(ncols = 2,loc = 'upper right',fontsize=9) plt.show()

BIOSCREEN vs Anatrans-model class of mibitrans¶

Comparison for all three options of biodegradation:

• no decay
• linear decay
• instantaneous reaction

In [14]:

#bio_ana; Anatrans model with BIOSCREEN input
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()

#bio_ana; Anatrans model with BIOSCREEN input 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()

Note that BIOSCREEN only provides limited resolution (indicated with markers). Differences in results, particularly for instantaneous reactin model are caused by that.

In [15]:

time_bio = bio_t
for i,ti in enumerate(bio_t):
    mbt.centerline(bio_ana_nd_results, time=ti,
                   color=colors[i],lw = 4,
                   label="Anatrans")

for i,ti in enumerate(bio_t):
    plt.scatter(bio_x, bioscreen_data.nodecay[i, 2, :],
             color=colors_bio[i], marker = '+', s=60, zorder=3,
             label=f"BIOSCREEN, t={ti/365}y")

plt.title("BIOSCREEN vs Anatrans models without decay")
plt.legend(ncols=2,loc = 'upper right',fontsize=9)
plt.show()

time_bio = bio_t for i,ti in enumerate(bio_t): mbt.centerline(bio_ana_nd_results, time=ti, color=colors[i],lw = 4, label="Anatrans") for i,ti in enumerate(bio_t): plt.scatter(bio_x, bioscreen_data.nodecay[i, 2, :], color=colors_bio[i], marker = '+', s=60, zorder=3, label=f"BIOSCREEN, t={ti/365}y") plt.title("BIOSCREEN vs Anatrans models without decay") plt.legend(ncols=2,loc = 'upper right',fontsize=9) plt.show()

In [16]:

for i,ti in enumerate(bio_t):
    mbt.centerline(bio_ana_dc_results, time=ti,
                   color=colors[i],lw = 4,
                   label="Anatrans")

for i,ti in enumerate(bio_t):
    plt.scatter(bio_x, bioscreen_data.lineardecay[i, 2, :],
             color=colors_bio[i], marker = '+', s=60, zorder=3,
             label=f"BIOSCREEN, t={ti/365}y")

plt.title("BIOSCREEN vs Anatrans models with linear decay")
plt.legend(ncols =2,loc = 'upper right')
plt.show()

for i,ti in enumerate(bio_t): mbt.centerline(bio_ana_dc_results, time=ti, color=colors[i],lw = 4, label="Anatrans") for i,ti in enumerate(bio_t): plt.scatter(bio_x, bioscreen_data.lineardecay[i, 2, :], color=colors_bio[i], marker = '+', s=60, zorder=3, label=f"BIOSCREEN, t={ti/365}y") plt.title("BIOSCREEN vs Anatrans models with linear decay") plt.legend(ncols =2,loc = 'upper right') plt.show()

In [17]:

for i,ti in enumerate(bio_t):
    mbt.centerline(bio_ana_inst_results, time=ti,
                   color=colors[i],lw = 4,
                   label="Anatrans")

for i,ti in enumerate(bio_t):
    plt.scatter(bio_x, bioscreen_data.instant[i, 2, :],
             color=colors_bio[i], marker = '+', s=60, zorder=3,
             label=f"BIOSCREEN, t={ti/365}y")

plt.title("BIOSCREEN vs Anatrans models with instant reaction")
plt.legend(ncols =2,loc = 'upper right',fontsize=9)
plt.show()

for i,ti in enumerate(bio_t): mbt.centerline(bio_ana_inst_results, time=ti, color=colors[i],lw = 4, label="Anatrans") for i,ti in enumerate(bio_t): plt.scatter(bio_x, bioscreen_data.instant[i, 2, :], color=colors_bio[i], marker = '+', s=60, zorder=3, label=f"BIOSCREEN, t={ti/365}y") plt.title("BIOSCREEN vs Anatrans models with instant reaction") plt.legend(ncols =2,loc = 'upper right',fontsize=9) plt.show()

Comparison with BIOSCREEN-AT¶

Making use of Example data of BIOSCREEN-AT.

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.

BIOSCREEN-AT vs Mibitrans-model class of mibitrans¶

Comparison for two options of biodegradation:

• no decay
• linear decay

Note that BIOSCREEN-AT does not have an implementation of the instantaneous reaction model for degradation (despite BIOSCREEN-AT taking in instant reaction parameters). Consequently, there is no comparison.

In [18]:

# 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()

# 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()

BIOSCREEN-AT (just as BIOSCREEN) only provides limited resolution which causes visual differences in the results. Actually calculated values (indicated with markers) of BIOSCREEN-AT and the Mibitrans-model class of mibitrans are identical.

In [19]:

for i,ti in enumerate(bat_t):
    mbt.centerline(bat_mbt_nd_results, time=ti,
                   color=colors[i],lw = 4,
                   label="Mibitrans")
for i,ti in enumerate(bat_t):
    plt.scatter(bat_x, bioscreenAT_data.nodecay[i, 2, :],
             color=colors_bio[i], marker = '+', s=60, zorder=3,
             label=f"BIOSCREEN-AT, t={ti/365}y")

plt.title("BIOSCREEN-AT vs Mibitrans models with no decay")
plt.legend(ncols =2,loc = 'upper right')
plt.show()

for i,ti in enumerate(bat_t): mbt.centerline(bat_mbt_nd_results, time=ti, color=colors[i],lw = 4, label="Mibitrans") for i,ti in enumerate(bat_t): plt.scatter(bat_x, bioscreenAT_data.nodecay[i, 2, :], color=colors_bio[i], marker = '+', s=60, zorder=3, label=f"BIOSCREEN-AT, t={ti/365}y") plt.title("BIOSCREEN-AT vs Mibitrans models with no decay") plt.legend(ncols =2,loc = 'upper right') plt.show()

In [20]:

bat_mbt_dc_results.plume_3d()

bat_mbt_dc_results.plume_3d()

Out[20]:

<Axes3D: title={'center': 'Plume plot of Mibitrans linear model, at t = 1825 days.'}, xlabel='Distance from source (m)', ylabel='Distance from plume center (m)', zlabel='Concentration [g/$m^{3}$]'>

In [21]:

for i,ti in enumerate(bat_t):
    mbt.centerline(bat_mbt_dc_results, time=ti,
                   color=colors[i],lw = 4,
                   label="Mibitrans")
for i,ti in enumerate(bat_t):
    plt.scatter(bat_x, bioscreenAT_data.lineardecay[i, 2, :],
             color=colors_bio[i], marker = '+', s=60, zorder=3,
             label=f"BIOSCREEN-AT, t={ti/365}y")

plt.title("BIOSCREEN-AT vs Mibitrans models with linear decay")
plt.legend(ncols =2,loc = 'upper right')
plt.show()

for i,ti in enumerate(bat_t): mbt.centerline(bat_mbt_dc_results, time=ti, color=colors[i],lw = 4, label="Mibitrans") for i,ti in enumerate(bat_t): plt.scatter(bat_x, bioscreenAT_data.lineardecay[i, 2, :], color=colors_bio[i], marker = '+', s=60, zorder=3, label=f"BIOSCREEN-AT, t={ti/365}y") plt.title("BIOSCREEN-AT vs Mibitrans models with linear decay") plt.legend(ncols =2,loc = 'upper right') plt.show()

Quantify differences¶

In [22]:

###  Functions for quantitative comparison
from mibitrans.analysis import differences as diff

### Functions for quantitative comparison from mibitrans.analysis import differences as diff

Redefine model dimension for mibitrans models to match data from BIOSCREEN and BIOSCREEN-AT¶

In [23]:

bio_model_comp = mbt.ModelParameters(
    model_length=1450 / ft, #[m]
    model_width=300 / ft, #[m]
    model_time=5 * 365, #[d]
    dx=145 / ft, #[m]
    dy=75 / ft, #[m]
    dt=365/2, #[d]
)

bat_model_comp = mbt.ModelParameters(
    model_length=1500 / ft, #[m]
    model_width=300 / ft, #[m]
    model_time=5 * 365, #[d]
    dx=150 / ft, #[m]
    dy=75 / ft, #[m]
    dt=365, #[d]
)

bio_model_comp = mbt.ModelParameters( model_length=1450 / ft, #[m] model_width=300 / ft, #[m] model_time=5 * 365, #[d] dx=145 / ft, #[m] dy=75 / ft, #[m] dt=365/2, #[d] ) bat_model_comp = mbt.ModelParameters( model_length=1500 / ft, #[m] model_width=300 / ft, #[m] model_time=5 * 365, #[d] dx=150 / ft, #[m] dy=75 / ft, #[m] dt=365, #[d] )

In [24]:

bio_object_comp = mbt.Bioscreen(bio_hydro, bio_att, bio_source, bio_model_comp)
bio_results_comp_nodecay = bio_object_comp.run()

bio_object_comp.attenuation_parameters.half_life = 0.1 * 365
bio_results_comp_linear= bio_object_comp.run()

bio_object_comp.instant_reaction(bio_ea)
bio_results_comp_instant = bio_object_comp.run()

mbt_object_comp = mbt.Mibitrans(bat_hydro, bat_att, bat_source, bat_model_comp)
mbt_results_comp_nodecay = mbt_object_comp.run()

mbt_object_comp.attenuation_parameters.half_life = 10*365
mbt_results_comp_linear = mbt_object_comp.run()

bio_object_comp = mbt.Bioscreen(bio_hydro, bio_att, bio_source, bio_model_comp) bio_results_comp_nodecay = bio_object_comp.run() bio_object_comp.attenuation_parameters.half_life = 0.1 * 365 bio_results_comp_linear= bio_object_comp.run() bio_object_comp.instant_reaction(bio_ea) bio_results_comp_instant = bio_object_comp.run() mbt_object_comp = mbt.Mibitrans(bat_hydro, bat_att, bat_source, bat_model_comp) mbt_results_comp_nodecay = mbt_object_comp.run() mbt_object_comp.attenuation_parameters.half_life = 10*365 mbt_results_comp_linear = mbt_object_comp.run()

In [25]:

bio_time_index = [0,1,3,-1]

bio_time_index = [0,1,3,-1]

In [26]:

#Quick check if model dimensions are similar to those used to obtain BIOSCREEN results
print(np.sum(abs(bio_x - bio_results_comp_instant.x)))
print(np.sum(abs(bio_y - bio_results_comp_instant.y)))
print(np.sum(abs(bio_t - bio_results_comp_instant.t[bio_time_index])))
print(np.sum(abs(bat_x - mbt_results_comp_linear.x)))
print(np.sum(abs(bat_y - mbt_results_comp_linear.y)))
print(np.sum(abs(bat_t - mbt_results_comp_linear.t)))

#Quick check if model dimensions are similar to those used to obtain BIOSCREEN results print(np.sum(abs(bio_x - bio_results_comp_instant.x))) print(np.sum(abs(bio_y - bio_results_comp_instant.y))) print(np.sum(abs(bio_t - bio_results_comp_instant.t[bio_time_index]))) print(np.sum(abs(bat_x - mbt_results_comp_linear.x))) print(np.sum(abs(bat_y - mbt_results_comp_linear.y))) print(np.sum(abs(bat_t - mbt_results_comp_linear.t)))

8.526512829121202e-14
3.552713678800501e-15
0.0
1.4210854715202004e-13
3.552713678800501e-15
0

In [27]:

max_source_conc = bio_source.source_zone_concentration[0]

mad_bio_bio_nd = diff.mean_absolute_difference(
    bioscreen_data.nodecay / max_source_conc, bio_results_comp_nodecay.relative_cxyt[bio_time_index,:,:]
)

mad_bio_bio_lindec = diff.mean_absolute_difference(
    bioscreen_data.lineardecay / max_source_conc, bio_results_comp_linear.relative_cxyt[bio_time_index,:,:]
)

mad_bio_bio_instant = diff.mean_absolute_difference(
    bioscreen_data.instant / max_source_conc, bio_results_comp_instant.relative_cxyt[bio_time_index,:,:]
)

max_source_conc = bio_source.source_zone_concentration[0] mad_bio_bio_nd = diff.mean_absolute_difference( bioscreen_data.nodecay / max_source_conc, bio_results_comp_nodecay.relative_cxyt[bio_time_index,:,:] ) mad_bio_bio_lindec = diff.mean_absolute_difference( bioscreen_data.lineardecay / max_source_conc, bio_results_comp_linear.relative_cxyt[bio_time_index,:,:] ) mad_bio_bio_instant = diff.mean_absolute_difference( bioscreen_data.instant / max_source_conc, bio_results_comp_instant.relative_cxyt[bio_time_index,:,:] )

In [28]:

#print(abs(BIOSCREEN_data_instant - bio_results_comp_instant.cxyt[bio_time_index,:,:]))

#print(abs(BIOSCREEN_data_instant - bio_results_comp_instant.cxyt[bio_time_index,:,:]))

In [29]:

print("Mean absolute difference in relative concentration between BIOSCREEN "
      "and Bioscreen model class with no decay:",mad_bio_bio_nd)
print("Mean absolute difference in relative concentration between BIOSCREEN "
      "and Bioscreen model class with linear decay:",mad_bio_bio_lindec)
print("Mean absolute difference in relative concentration between BIOSCREEN "
      "and Bioscreen model class with instant reaction:",mad_bio_bio_instant)

print("Mean absolute difference in relative concentration between BIOSCREEN " "and Bioscreen model class with no decay:",mad_bio_bio_nd) print("Mean absolute difference in relative concentration between BIOSCREEN " "and Bioscreen model class with linear decay:",mad_bio_bio_lindec) print("Mean absolute difference in relative concentration between BIOSCREEN " "and Bioscreen model class with instant reaction:",mad_bio_bio_instant)

Mean absolute difference in relative concentration between BIOSCREEN and Bioscreen model class with no decay: 2.4880399540013646e-06
Mean absolute difference in relative concentration between BIOSCREEN and Bioscreen model class with linear decay: 8.033000606564948e-06
Mean absolute difference in relative concentration between BIOSCREEN and Bioscreen model class with instant reaction: 6.444204434847127e-05

In [30]:

max_source_conc = bat_source.source_zone_concentration[0]

mad_bat_mibi_nd = diff.mean_absolute_difference(
    bioscreenAT_data.nodecay / max_source_conc,
    mbt_results_comp_nodecay.relative_cxyt
)

mad_bat_mibi_lindec = diff.mean_absolute_difference(
    bioscreenAT_data.lineardecay / max_source_conc,
    mbt_results_comp_linear.relative_cxyt
)

max_source_conc = bat_source.source_zone_concentration[0] mad_bat_mibi_nd = diff.mean_absolute_difference( bioscreenAT_data.nodecay / max_source_conc, mbt_results_comp_nodecay.relative_cxyt ) mad_bat_mibi_lindec = diff.mean_absolute_difference( bioscreenAT_data.lineardecay / max_source_conc, mbt_results_comp_linear.relative_cxyt )

In [31]:

print("Mean absolute difference in relative concentration between BIOSCREEN-AT "
      "and Mibitrans model class with no decay:",mad_bat_mibi_nd)
print("Mean absolute difference in relative concentration between BIOSCREEN-AT "
      "and Mibitrans model class with linear decay:",mad_bat_mibi_lindec)

print("Mean absolute difference in relative concentration between BIOSCREEN-AT " "and Mibitrans model class with no decay:",mad_bat_mibi_nd) print("Mean absolute difference in relative concentration between BIOSCREEN-AT " "and Mibitrans model class with linear decay:",mad_bat_mibi_lindec)

Mean absolute difference in relative concentration between BIOSCREEN-AT and Mibitrans model class with no decay: 1.0600970801086523e-05
Mean absolute difference in relative concentration between BIOSCREEN-AT and Mibitrans model class with linear decay: 1.0865108469319337e-05