Basic use of mibitrans¶

This example notebook presents the features implemented in the mibitrans package, and explains how to use them. The parameters used here are example field data in reasonable orders of magnitude, but do not represent any field site in paticular.

Authors: Jorrit Bakker, Alraune Zech

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

Input by dataclasses¶

mibitrans use dataclasses to handle data input. Input is grouped into four categories, each having a separate input dataclass:

• hydrological parameters
• attenuation parameters
• source parameters
• model parameters

Those are required for a regular transport model. An instantaneous reaction model additionally needs electron acceptor concentrations, this is discussed in a later section.

Units of parameters have default values, consistent between data input categories. Input values need to be expressed accordingly to avoid mishaps. It is possible to provide input in other units, but this requires the user to ensure that they are consistent throughout the entire modelling process.

Hydrological parameters¶

The data class HydrologicalParameters contains parameters that are inherent to the aquifer properties and relevant for the conservative transport processes of advection and dispersion: flow velocity, porosity, dispersivity and diffusion.

In [2]:

hydro = mbt.HydrologicalParameters(
    velocity=0.1,   # Groundwater flow velocity, in [m/day]
    porosity=0.3,   # Effective soil porosity [-]
    alpha_x=3,      # Longitudinal dispersivity, in [m]
    alpha_y=0.02,   # Transverse horizontal dispersivity, in [m]
    alpha_z=0.005,  # Transverse vertical dispersivity, in [m]
    diffusion=1e-5  # Molecular diffusion, in [m2/day]
)

hydro = mbt.HydrologicalParameters( velocity=0.1, # Groundwater flow velocity, in [m/day] porosity=0.3, # Effective soil porosity [-] alpha_x=3, # Longitudinal dispersivity, in [m] alpha_y=0.02, # Transverse horizontal dispersivity, in [m] alpha_z=0.005, # Transverse vertical dispersivity, in [m] diffusion=1e-5 # Molecular diffusion, in [m2/day] )

Alternatively, hydraulic conductivity and hydraulic gradient can be provided. Flow velocity is then calculated from those.

In [3]:

hydro = mbt.HydrologicalParameters(
    h_gradient=0.005,   # Hydraulic gradient [-]
    h_conductivity=6,   # Hydraulic conductivity [m/day]
    porosity=0.3,       # Effective soil porosity [-]
    alpha_x=3,          # Longitudinal dispersivity, in [m]
    alpha_y=0.02,       # Transverse horizontal dispersivity, in [m]
    alpha_z=0.005,      # Transverse vertical dispersivity, in [m]
)

hydro = mbt.HydrologicalParameters( h_gradient=0.005, # Hydraulic gradient [-] h_conductivity=6, # Hydraulic conductivity [m/day] porosity=0.3, # Effective soil porosity [-] alpha_x=3, # Longitudinal dispersivity, in [m] alpha_y=0.02, # Transverse horizontal dispersivity, in [m] alpha_z=0.005, # Transverse vertical dispersivity, in [m] )

Check the value of an input parameter by calling the dataclass with the corresponding attribute name.

In [4]:

print("The calculated groundwater flow velocity is: ", hydro.velocity, r"m/d")

print(f"The dispersivity values are: {hydro.alpha_x}m, {hydro.alpha_y}m and {hydro.alpha_z}m, "
      "for the x, y and z directions respectively.")

print("The calculated groundwater flow velocity is: ", hydro.velocity, r"m/d") print(f"The dispersivity values are: {hydro.alpha_x}m, {hydro.alpha_y}m and {hydro.alpha_z}m, " "for the x, y and z directions respectively.")

The calculated groundwater flow velocity is:  0.1 m/d
The dispersivity values are: 3m, 0.02m and 0.005m, for the x, y and z directions respectively.

If transverse vertical dispersivity, $\alpha_z$ and/or diffusion are not given, they default to $0$.

Note that only the Mibitrans model includes diffusion in its calculations. Anatrans and Bioscreen do not, due to assumptions taken in their derivation.

Attenuation parameters¶

The data class AttenuationParameters handles parameters related to attenuation by adsorption and degradation through linear decay.

Attenuation: Retardation¶

Retardation can be simply given as a value $R \ge 1$.

In [5]:

att = mbt.AttenuationParameters(
    retardation=1.3 # Retardation factor for transported contaminants [-]
)

att = mbt.AttenuationParameters( retardation=1.3 # Retardation factor for transported contaminants [-] )

Alternatively, the retardation coefficient $R$ can be calculated from soil bulk density, partition coefficient and the soil organic carbon fraction - when these parameters are provided. Note that the calculation of the retardation factor also requires porosity, which has to be provided in the HydrologicalParameters class. The calculation of retardation is automatically performed through the function calculate_retardation() when running the model. But this function can also be manually called.

In [6]:

att = mbt.AttenuationParameters(
    bulk_density=1.7,               # Soil bulk density in [g/m^3]
    partition_coefficient=53,       # Partition coefficient  in [m^3/g]
    fraction_organic_carbon=1e-3    # Fraction of organic material in soil [-]
)

# Calculate retardation factor beforehand (provide porosity):
att.calculate_retardation(porosity=hydro.porosity)

# And then check the value:
print("The calculated retardation value is: ", att.retardation)

att = mbt.AttenuationParameters( bulk_density=1.7, # Soil bulk density in [g/m^3] partition_coefficient=53, # Partition coefficient in [m^3/g] fraction_organic_carbon=1e-3 # Fraction of organic material in soil [-] ) # Calculate retardation factor beforehand (provide porosity): att.calculate_retardation(porosity=hydro.porosity) # And then check the value: print("The calculated retardation value is: ", att.retardation)

The calculated retardation value is:  1.3003333333333333

Attenuation: degradation parameters¶

The basic concept of modelling degradation (e.g. as result of biological activity) is by a linear decay model, independent of any geochemical and biological settings. Input for the linear decay model is either the contaminant decay rate $\mu$ or half life $t_{1/2}$.

Both options to provide these parameters are implemented. The other corresponding parameter is automatically calculated.

In [7]:

att = mbt.AttenuationParameters(
    half_life=365/2,         # Contaminant half life, in [days]
)
print("The internally decay rate (1/days) is: ", att.decay_rate)

att = mbt.AttenuationParameters(
    decay_rate=3.798e-3     # Contaminant first order decay rate in [1/days]
)
print("The internally calculated contaminant half life (days) is: ", att.half_life)

att = mbt.AttenuationParameters( half_life=365/2, # Contaminant half life, in [days] ) print("The internally decay rate (1/days) is: ", att.decay_rate) att = mbt.AttenuationParameters( decay_rate=3.798e-3 # Contaminant first order decay rate in [1/days] ) print("The internally calculated contaminant half life (days) is: ", att.half_life)

The internally decay rate (1/days) is:  0.0037980667427942208
The internally calculated contaminant half life (days) is:  182.50320709845846

A complete set of attenuation parameters (for the examples below) would look like:

In [8]:

att = mbt.AttenuationParameters(
    bulk_density=1.7,               # Soil bulk density in [g/m^3]
    partition_coefficient=53,       # Partition coefficient in [m^3/g]
    fraction_organic_carbon=1e-3,   # Fraction of organic material in soil [-]
    half_life=365/2,                 # Contaminant half life, in [days]
)

att = mbt.AttenuationParameters( bulk_density=1.7, # Soil bulk density in [g/m^3] partition_coefficient=53, # Partition coefficient in [m^3/g] fraction_organic_carbon=1e-3, # Fraction of organic material in soil [-] half_life=365/2, # Contaminant half life, in [days] )

Source parameters¶

The data class SourceParameters handles all parameters related to the contaminant source. Conceptually, the source is treated as a seperate phase, which dissolves into the groundwater over time. mibitrans also makes use of source superposition that reflects the typical field observation that the source expands over a certain width with decreasing concentrations towards the margins. Input parameters are thus:

• source depth depth: extension of contaminant source in depth ($z$-direction); is assumed constant over the entire width of the source
• total mass total_mass: total amound of contaminant source in NAPL phase; can be a specific value or infinite.
• width of the different horizontal source zones source_zone_boundary: defined as array of boundary locations relative to source center
• the contaminant concentrations within each source zone source_zone_concentration: array of similar shape as source_zone_boundary as it specifies the contaminant concentration released to the groundwater within each of the source zones.

First an example for a simple source zone with no separate zones having a width of $20$m, depth of $3m$, a continuous input (infinite source mass) and a release concentration of $5$ g/m^3:

In [9]:

source = mbt.SourceParameters(
    source_zone_boundary=np.array([20]),     # Source zone extent, in [m]
    source_zone_concentration=np.array([5]), # Source zone concentrations, in [g/m^3]
    depth=2,                                 # Source depth extent, in [m]
    total_mass="infinite"                    # Source mass, in [g] or as np.inf"
)

source = mbt.SourceParameters( source_zone_boundary=np.array([20]), # Source zone extent, in [m] source_zone_concentration=np.array([5]), # Source zone concentrations, in [g/m^3] depth=2, # Source depth extent, in [m] total_mass="infinite" # Source mass, in [g] or as np.inf" )

Source zonation¶

The source is assumed to be symmetrical in $y$-direction with the center located at $y=0$. The width of the source is divided into $n$ zones, where $n$ is a free choice of the user. The array source_zone_boundary defines the boundaries between zones. The release concentration within each zone is defined in the array source_zone_concentration. Release concentrations decrease from the center to the fringes. The source depth is assumed to be constant over the entire width.

Second example: Consider a source of total width of $40$m. It has a very high concentration in the source center at $y \in [-2m,2m]$, a lower concentration up to the distance of $\pm 11$ m (specifically for $y \in [-11m,-2m]$ and $y \in [2m,11m]$) and the lowest concentration from $\pm 11m$ up to $\pm 20m$. This would be entered as array of $[2,11,20]$ for source_zone_boundary and an array of the same shape for the corresponding source zone concentrations.

The shape of the source can be easily visualized with source.visualize()

In [10]:

source = mbt.SourceParameters(
    source_zone_boundary=np.array([5, 10, 20]),       #[m]
    source_zone_concentration=np.array([10, 4, 0.5]), #[g/m3]
    depth=2,                                        #[m]
    total_mass=np.inf                               #[g]
)

source = mbt.SourceParameters( source_zone_boundary=np.array([5, 10, 20]), #[m] source_zone_concentration=np.array([10, 4, 0.5]), #[g/m3] depth=2, #[m] total_mass=np.inf #[g] )

Visualize what the source zone looks like to check your input:

In [11]:

source.visualize_source_zone()
plt.show()

source.visualize_source_zone() plt.show()

By giving the total solid-phase contaminant source mass (total_mass) a finite amount, the source will not stay constant in strength but deplete over time. Consequently during simulation, the source mass and the source zone concentrations will diminish over time. This aspect is implemented as source depletion at the rate $\gamma_s$ in $1/day$. The rate $\gamma_s$ depends on the flow velocity, size of the source zone and source zone concentrations. When the source mass is considered infinite, contaminant source concentrations will not diminish over time (i.e. $\gamma_s = 0$). For more information about source depletion, see the theory section of the documentation.

In [12]:

source = mbt.SourceParameters(
    source_zone_boundary=np.array([2, 5, 10, 20]),      #[m]
    source_zone_concentration=np.array([15, 10, 4, 0.5]),  #[g/m3]
    depth=2,                                            #[m]
    total_mass=5e6                                      #[g]
)

source = mbt.SourceParameters( source_zone_boundary=np.array([2, 5, 10, 20]), #[m] source_zone_concentration=np.array([15, 10, 4, 0.5]), #[g/m3] depth=2, #[m] total_mass=5e6 #[g] )

The source depletion can be visualized using source.visualize_source_depletion(). And requires additional input of the hydrological parameters to calculate. For the instant reaction model, electron acceptor concentrations and (optionally) utilization factors are also needed. The input method for the latter two parameters is the same as described in the 'Instant reaction model' section below.

In [13]:

source.visualize_source_depletion(hydro)
plt.show()

source.visualize_source_depletion(hydro) plt.show()

Model parameters¶

The data class ModelParameters handles all parameters related to the spatial and temporal settings of the domain in which the contaminant plume should be calculated, including model dimensions and discretization parameters:

• model_length is the domain extent in the (x) direction parallel to the groundwater flow direction.
• model_width is the domain extent of the model perpendicular (y) to the groundwater flow direction.
• model_time is the model duration, i.e. the final time of plume calculation.
• dx and dy handle step size of the grid at which values are calculated.
• dt handles time step sizes, i.e. time steps at which values are calculated.

The user should ensure that the source zone fits inside of the given model width. If step sizes are undefined, default values are chosen as ratios to model dimensions with model_length$/100$, model_width$/50$ and model_time$/10$.

In [14]:

model = mbt.ModelParameters(
    model_length=220,       # Model extent in the x direction in [m].
    model_width=60,         # Model extent in the y direction in [m].
    model_time=5 * 365,     # Model duration in [days].
    dx=1,                   # Model discretization step size in x direction, in [m].
    dy=0.5,                 # Model discretization step size in y direction, in [m].
    dt=365/5                # Model time discretization step size, in [days]
)

model = mbt.ModelParameters( model_length=220, # Model extent in the x direction in [m]. model_width=60, # Model extent in the y direction in [m]. model_time=5 * 365, # Model duration in [days]. dx=1, # Model discretization step size in x direction, in [m]. dy=0.5, # Model discretization step size in y direction, in [m]. dt=365/5 # Model time discretization step size, in [days] )

Input checking¶

All data classes for parameter input are checked. It evaluates if the required input parameters are present, if they are of the correct data type and if they are in the expected domain. We provide here some examples of data input that fails the checks.

Uncomment the cells below to see results of the input check - Cells are intended to raise an error!

In [15]:

# bad_hydro = mbt.HydrologicalParameters(velocity = 1)

# bad_hydro = mbt.HydrologicalParameters(velocity = 1)

Not all required parameters are specified; HydrologicalParameters needs the porosity, longitudinal dispersivity and transverse horizontal dispersivity as well.

In [16]:

# bad_source = mbt.SourceParameters(
#     source_zone_boundary = np.array([1,2,3]),
#     source_zone_concentration = "this is a string, not an array",
#     depth = 10,
#     total_mass = "inf",
# )

# bad_source = mbt.SourceParameters( # source_zone_boundary = np.array([1,2,3]), # source_zone_concentration = "this is a string, not an array", # depth = 10, # total_mass = "inf", # )

The given input datatype is not what is expected; source_zone_concentration expects a numpy array of the same length as the array given in source_zone_boundary.

In [17]:

# bad_att = mbt.AttenuationParameters(
#     retardation = 0.1
# )

# bad_att = mbt.AttenuationParameters( # retardation = 0.1 # )

The parameter has a value outside its valid domain; retardation should have a value $R\geq1$.

Transport simulation¶

Choice of transport solution¶

The package mibitrans offers three analytical solutions to the contaminant transport equation ADE (with advection, dispersion, diffusion, retardation and decay) under uniform groundwater flow in 3D and constant source supply. All three solutions implemented different equations:

• Mibitrans solution is the full exact solution to the ADE under the specified initial and boundary conditions (Wexler, 1992; Karanovic 2007). The Mibitrans solution requires numerical integration which causes a higher computational effort and duration of calculation.
• Anatrans solution is derived based on a small simplification for transport in transverse directions. It provides a fully analytical (i.e. does not require numerical integration), but approximate solution to the ADE. It is interesting for screening given speed up of calculation, particularly when concentrations in transverse direction are of minor interest.
• Bioscreen solution makes use of the approximate analytical solution introduced in BIOSCREEN (Newell et al., 1996,1997) building up on the solution of Domenico (1987) in various adaptions. It is a further simplification of the equation used in Anatrans, introducing a mathematical truncation affecting correctness of calculated contaminant concentration in flow direction, particularly for Peclet numbers lower that 100.

Each solution is implemented as separate model class, located in mibitrans.transport.models. A model is first initialized through making the model object and passing on the model parameters defined in the input data classes$^1$. During model initialization, the model grid is generated and some prior parameter are calculated. After initialization, no transport results have been calculated yet.

Contaminant concentrations as result of transport are calculated when calling the model instance with its class method run. This generates a results class object, which contains the concentration distribution as a 3D numpy array (indexed as $C[t,y,x]$) and the parameters used for this specific run.

$^1$Note that the input argument names for the dataclass follow snake_case instead of PascalCase (e.g. hydrological_paramaters takes in the data class HydrologicalParameters) following Python conventions.

Mibitrans model¶

As the namesake model of this package, it is the recommended model to use. It is the fully exact solution to the contaminant transport equation under uniform flow and a constant contaminant source while the other two models introduce a margin of error as consequence of taken assumptions. The exact nature of these differences as well as the underlying equation can be found in the package documentation.

However, since this model requires evaluation of an integral, computation time might be longer, depending on model discretization.

In [18]:

# Initializing the model object
mbt_object = mbt.Mibitrans(
    hydrological_parameters=hydro,
    attenuation_parameters=att,
    source_parameters=source,
    model_parameters=model,
    verbose=False #Set to True to follow model calculation steps
)

# Initializing the model object mbt_object = mbt.Mibitrans( hydrological_parameters=hydro, attenuation_parameters=att, source_parameters=source, model_parameters=model, verbose=False #Set to True to follow model calculation steps )

In [19]:

# Can use model object to request input and various attributes
print("Retarded flow velocity: ", mbt_object.rv, r"r$m/day$")
print(len(mbt_object.x), "x-steps,", len(mbt_object.y), "y-steps,", len(mbt_object.t), "t-steps")
print("Input source mass: ", mbt_object.source_parameters.total_mass, "g")

# Can use model object to request input and various attributes print("Retarded flow velocity: ", mbt_object.rv, r"r$m/day$") print(len(mbt_object.x), "x-steps,", len(mbt_object.y), "y-steps,", len(mbt_object.t), "t-steps") print("Input source mass: ", mbt_object.source_parameters.total_mass, "g")

Retarded flow velocity:  0.07690335811330429 r$m/day$
221 x-steps, 121 y-steps, 25 t-steps
Input source mass:  5000000.0 g

In [20]:

# Run the model once initialized and obtain the results object
mbt_results = mbt_object.run()

# Run the model once initialized and obtain the results object mbt_results = mbt_object.run()

Model results¶

Contaminant concentrations and the input parameters used for this specific run are returned as result object (here mbt_results) after running the model. This ensures that it is always clear which conditions correspond with which result, even when re-running the model.

The model results object does not only contain all data, but also offers features for visualization and analysis, such as plotting and generating a mass balance. This is further elaborated upon in the Visualization and Mass Balance sections below.

In [21]:

model_cxyt = mbt_results.cxyt # Concentration distribution in cxyt class-attribute

model_cxyt = mbt_results.cxyt # Concentration distribution in cxyt class-attribute

cxyt is indexed as [time, y-position, x-position].

Get the concentration at the last time step, in the center of the plume for all $x$ via:

In [22]:

plume_center = model_cxyt[-1, len(mbt_results.y)//2, :]
print("The first 3 concentrations at the plume center are: ", plume_center[:3], "g/m3")

plume_center = model_cxyt[-1, len(mbt_results.y)//2, :] print("The first 3 concentrations at the plume center are: ", plume_center[:3], "g/m3")

The first 3 concentrations at the plume center are:  [14.94425883 14.3056162  13.69409429] g/m3

Access information of parameters of model run from results object via:

In [23]:

print("The flow velocity of the model was:", mbt_results.hydrological_parameters.velocity, "m/d")

print("The flow velocity of the model was:", mbt_results.hydrological_parameters.velocity, "m/d")

The flow velocity of the model was: 0.1 m/d

Use build-in visualization methods of the results object:

(more details in the visualization section)

In [24]:

mbt_results.centerline()
plt.show()

mbt_results.centerline() plt.show()

Anatrans model¶

The Anatrans model implements a fully analytical equation speeding calculation up as no numerical integration is needed as in the Mibitrans model. Anatrans is fully exact for transport in horizontal ($x$) direction of the groundwater flow and takes assumption for transport in vertical directions ($y$ and $z$). The underlying equation can be found in the package documentation. Differences in concentrations in vertical directions introduced through the approximation are discussed in the notebook example_mibitrans_anatrans.ipynb.

In [25]:

ana_object = mbt.Anatrans(
    hydrological_parameters=hydro,
    attenuation_parameters=att,
    source_parameters=source,
    model_parameters=model
)

ana_results = ana_object.run()
ana_results.centerline()
plt.show()

ana_object = mbt.Anatrans( hydrological_parameters=hydro, attenuation_parameters=att, source_parameters=source, model_parameters=model ) ana_results = ana_object.run() ana_results.centerline() plt.show()

Bioscreen model¶

This model is an exact implementation of the transport equations implemented in the BIOSCREEN screening model of Newell et al. (1997), which is based on the Domenico (1987) analytical model. Using a truncated version of the equation used in the Anatrans model. This model is implemented as a method of comparison with the original BIOSCREEN software. And is included for legacy reasons, since it is the first model implemented in the mibitrans package, serving as a basis for the other models. However, caution should be taken when using this model, since a varying error is introduced by using the truncated analytical solution. The error is most prominent for shorter times and distances from the source, and depends on the ratio of flow velocity and longitudinal dispersivity. For modelling, the Anatrans (untruncated approximate solution) and Mibitrans (exact analytical solution) models are recommended instead.

In [26]:

bio_object = mbt.Bioscreen(
    hydrological_parameters=hydro,
    attenuation_parameters=att,
    source_parameters=source,
    model_parameters=model
)
bio_results = bio_object.run()
bio_results.centerline()
plt.show()

bio_object = mbt.Bioscreen( hydrological_parameters=hydro, attenuation_parameters=att, source_parameters=source, model_parameters=model ) bio_results = bio_object.run() bio_results.centerline() plt.show()

Instant reaction model¶

The mibitrans package offers another option for modeling degradation of contaminants during transport (beside the linear decay), which is the instant reaction model. If geochemical field site data (e.g. from observations) is available, contaminant transport and biodegradation can be modeled by usage of the instantaneous reaction model coupled with transport. This model superimposes the effect of biodegradation reactions that are limited by electron acceptor availability to the conservative analytical transport model. Details on the theoretical background of the method are provided in the documentation, as are the details on this type of reaction model.

The data class ElectronAcceptors handles the concentrations of electron acceptors (EA) available for degradation (oxygen, nitrate, sulfate) and metabolic by-products of the degradation reactions (ferrous iron and methane). For the EA concentrations (delta_XX), the available amount is assumed to be the difference between background oxygen and current concentration in groundwater. For metabolic by-products the current concentrations is used.

The data class UtilizationFactors takes the utilization factor (UF) for each respective electron acceptor. These factors relate EA concentrations to the amount of degraded contaminants based on stoichiometric relations in chemical reactions. E.g. the UF of oxygen determines the mass of oxygen consumed per mass of biodegraded contaminant. For ferrous iron it identifies the mass of ferrous iron generated per mass of biodegraded contaminant.

These data classes are input for the instant_reaction class method of the transport class (available for Mibitrans, Anatrans and Bioscreen). Electron acceptor concentrations are required as input. For the utilization factor, default values for a BTEX mixture are used, unless different input is provided.

In [27]:

# For streamlined input, use the ElectronAcceptor dataclass;
ea = mbt.ElectronAcceptors(
    delta_oxygen=9,     # Available oxygen in groundwater, in [g/m^3]
    delta_nitrate=6,    # Available nitrate in groundwater, in [g/m^3]
    delta_sulfate=5,    # Available sulfate in groundwater, in [g/m^3]
    ferrous_iron=8,    # Ferrous iron concentration in groundwater, in [g/m^3]
    methane=4          # Methane concentration in groundwater, in [g/m^3]
)
mbt_object.instant_reaction(electron_acceptors=ea)

# For streamlined input, use the ElectronAcceptor dataclass; ea = mbt.ElectronAcceptors( delta_oxygen=9, # Available oxygen in groundwater, in [g/m^3] delta_nitrate=6, # Available nitrate in groundwater, in [g/m^3] delta_sulfate=5, # Available sulfate in groundwater, in [g/m^3] ferrous_iron=8, # Ferrous iron concentration in groundwater, in [g/m^3] methane=4 # Methane concentration in groundwater, in [g/m^3] ) mbt_object.instant_reaction(electron_acceptors=ea)

In [28]:

# Can adapt utilization factors if needed
uf = mbt.UtilizationFactor(
    util_oxygen=2,        # utilization factor of oxygen, in [g/g]
    util_nitrate=1,       # utilization factor of nitrate, in [g/g]
    util_sulfate=3,       # utilization factor of sulfate, in [g/g]
    util_ferrous_iron=4,  # utilization factor of ferrous iron, in [g/g]
    util_methane=5,       # utilization factor of methane, in [g/g]
)
mbt_object.instant_reaction(electron_acceptors=ea, utilization_factor=uf)

# Can adapt utilization factors if needed uf = mbt.UtilizationFactor( util_oxygen=2, # utilization factor of oxygen, in [g/g] util_nitrate=1, # utilization factor of nitrate, in [g/g] util_sulfate=3, # utilization factor of sulfate, in [g/g] util_ferrous_iron=4, # utilization factor of ferrous iron, in [g/g] util_methane=5, # utilization factor of methane, in [g/g] ) mbt_object.instant_reaction(electron_acceptors=ea, utilization_factor=uf)

Alternatively, electron acceptors and/or utilization factors can be entered as a dictionary:

In [29]:

# Alternatively, electron acceptors and/or utilization factors can be entered as a dictionary
mbt_object.instant_reaction(
    electron_acceptors={
        "delta_oxygen":9,        # Available oxygen concentration in groundwater, in [g/m^3]
        "delta_nitrate":6,       # Available nitrate concentration in groundwater, in [g/m^3]
        "delta_sulfate":5,       # Available sulfate concentration in groundwater, in [g/m^3]
        "ferrous_iron":8,        # Ferrous iron concentration in groundwater, in [g/m^3]
        "methane":4,             # Methane concentration in groundwater, in [g/m^3]
    }
)

# Alternatively, electron acceptors and/or utilization factors can be entered as a dictionary mbt_object.instant_reaction( electron_acceptors={ "delta_oxygen":9, # Available oxygen concentration in groundwater, in [g/m^3] "delta_nitrate":6, # Available nitrate concentration in groundwater, in [g/m^3] "delta_sulfate":5, # Available sulfate concentration in groundwater, in [g/m^3] "ferrous_iron":8, # Ferrous iron concentration in groundwater, in [g/m^3] "methane":4, # Methane concentration in groundwater, in [g/m^3] } )

Or as a list:

In [30]:

mbt_object.instant_reaction(
    electron_acceptors=[9, 6, 8, 5, 4],
    utilization_factor=[3.14, 4.9, 21.8, 4.7, 0.78]
)

mbt_object.instant_reaction( electron_acceptors=[9, 6, 8, 5, 4], utilization_factor=[3.14, 4.9, 21.8, 4.7, 0.78] )

Note that using instant_reaction also resets the utilization factor to default:

In [31]:

print("electron acceptor concentrations: ", mbt_object.electron_acceptors)
print("electron acceptor utilization factors: ", mbt_object.utilization_factor)

print("electron acceptor concentrations: ", mbt_object.electron_acceptors) print("electron acceptor utilization factors: ", mbt_object.utilization_factor)

electron acceptor concentrations:  {'delta_oxygen': 9, 'delta_nitrate': 6, 'ferrous_iron': 8, 'delta_sulfate': 5, 'methane': 4}
electron acceptor utilization factors:  {'util_oxygen': 3.14, 'util_nitrate': 4.9, 'util_ferrous_iron': 21.8, 'util_sulfate': 4.7, 'util_methane': 0.78}

Running the model now would provide results using the instantaneous reaction model. You can switch between linear degradation and instant reaction using the model object.

Switch back from between instant reaction to linear decay model:

In [32]:

# Switch between instant reaction and linear model
mbt_object.mode = "linear"
print("The current model mode is: ", mbt_object.mode)

# Switch between instant reaction and linear model mbt_object.mode = "linear" print("The current model mode is: ", mbt_object.mode)

The current model mode is:  linear

Switch back to instant reaction as follows:

In [33]:

mbt_object.mode = "instant_reaction"
print("The current model mode is: ", mbt_object.mode)

mbt_object.mode = "instant_reaction" print("The current model mode is: ", mbt_object.mode)

The current model mode is:  instant_reaction

Sample¶

Use the sample method to get the concentration at a specific location and time.

In [34]:

# Instead of distribution, sample a specific location at any point along the plume and any point in time
concentration = mbt_object.sample(
    x_position=20,
    y_position=0,
    time=365)
print("The concentration at 20m downstream from the source after 1 year is:", concentration, r"g/m^3")

# Instead of distribution, sample a specific location at any point along the plume and any point in time concentration = mbt_object.sample( x_position=20, y_position=0, time=365) print("The concentration at 20m downstream from the source after 1 year is:", concentration, r"g/m^3")

The concentration at 20m downstream from the source after 1 year is: 10.251893081117473 g/m^3

Visualization¶

The mibitrans package has various ways to visualize and output the model results, including 1D plotting and multidimensional plots. Other options to retrieve information from modelling results are output of sample data and calculations of the mass balance. Find examples for each feature below.

First, make some different model objects to use in the plotting, takes a couple of seconds to run.

In [35]:

#Mibitrans model
mbt_object = mbt.Mibitrans(hydro, att, source, model)
mbt_lineardecay = mbt_object.run()

mbt_object.attenuation_parameters.decay_rate = 0
mbt_nodecay = mbt_object.run()

mbt_object.instant_reaction(electron_acceptors=ea)
mbt_instant = mbt_object.run()

#Anatrans model
ana_object = mbt.Anatrans(hydro, att, source, model)
ana_lineardecay = ana_object.run()

ana_object.attenuation_parameters.decay_rate = 0
ana_nodecay = ana_object.run()

ana_object.instant_reaction(electron_acceptors=ea)
ana_instant =  ana_object.run()

#Bioscreen model
bio_object = mbt.Bioscreen(hydro, att, source, model)
bio_lineardecay = bio_object.run()

bio_object.attenuation_parameters.decay_rate = 0
bio_nodecay = bio_object.run()

bio_object.instant_reaction(electron_acceptors=ea)
bio_instant = bio_object.run()

#Mibitrans model mbt_object = mbt.Mibitrans(hydro, att, source, model) mbt_lineardecay = mbt_object.run() mbt_object.attenuation_parameters.decay_rate = 0 mbt_nodecay = mbt_object.run() mbt_object.instant_reaction(electron_acceptors=ea) mbt_instant = mbt_object.run() #Anatrans model ana_object = mbt.Anatrans(hydro, att, source, model) ana_lineardecay = ana_object.run() ana_object.attenuation_parameters.decay_rate = 0 ana_nodecay = ana_object.run() ana_object.instant_reaction(electron_acceptors=ea) ana_instant = ana_object.run() #Bioscreen model bio_object = mbt.Bioscreen(hydro, att, source, model) bio_lineardecay = bio_object.run() bio_object.attenuation_parameters.decay_rate = 0 bio_nodecay = bio_object.run() bio_object.instant_reaction(electron_acceptors=ea) bio_instant = bio_object.run()

Centerline plot¶

It shows the concentration distribution $C(x, y_s, t_s)$ in the longitudinal ($x$) direction, at a specific vertical plume location $y_s$ and specific time $t_s$. By default, the fixed vertical location is the center of the plume ($y_s=0$) and the default time is the last time step model_time as defined in the data class ModelParameters. The $y_s$-position along the plume and point of time $t_s$ can be adapted with function arguments.

Plotting can be performed either as class method, i.e. model_object.centerline(), or by calling the centerline function with the model object as argument, i.e. centerline(model_object).

In [36]:

# Pass result object to plotting function, using default time
mbt.centerline(mbt_nodecay)
plt.show()

# Pass result object to plotting function, using default time mbt.centerline(mbt_nodecay) plt.show()

Plume location¶

Define arguments time and y_position to plot somewhere of the plume center and at a specific time. If the time is not identical to a specific time step defined during time stepping (determined by dt in the data class ModelParameters), the closest point in time is identified and the concentration profile at that time is provided.

In [37]:

mbt.centerline(mbt_nodecay, time=5*365, y_position=7)
plt.show()

mbt.centerline(mbt_nodecay, time=5*365, y_position=7) plt.show()

Plot modification¶

Plot settings can be changed to adapt appearance and combine plots.

The centerline function (transverse and breakthrough as well) accepts a list of multiple models and can plot them at the same time. It further accepts a list of names to use for the plot legend.

In [38]:

mbt.centerline([mbt_nodecay, mbt_lineardecay, mbt_instant],
               time=5*365,
               legend_names=["no degradation", "linear decay", "instant reaction"])
plt.show()

mbt.centerline([mbt_nodecay, mbt_lineardecay, mbt_instant], time=5*365, legend_names=["no degradation", "linear decay", "instant reaction"]) plt.show()

Change (axis) title, figure size, plot extent, etc.¶

You can modify the plot with the respective commands for modifying lineplots (plt.plot()).

In [39]:

mbt.centerline([mbt_nodecay, ana_nodecay, bio_nodecay],time=3*365)
plt.title("Better title than the one generated automatically")
plt.xlabel("I have changed!")
plt.ylabel("And so have I")
plt.legend(["Mibitrans", "Anatrans", "Bioscreen"])
plt.xlim((0, 150))
plt.show()

mbt.centerline([mbt_nodecay, ana_nodecay, bio_nodecay],time=3*365) plt.title("Better title than the one generated automatically") plt.xlabel("I have changed!") plt.ylabel("And so have I") plt.legend(["Mibitrans", "Anatrans", "Bioscreen"]) plt.xlim((0, 150)) plt.show()

Calling each centerline function separately per model gives more control over plot layout. Keyword arguments for plt.plot can be passed on through the function.

In [40]:

time = 5 #[years]
y0, y1, y2, y3 = 0, 6, 9, 20
mbt.centerline(mbt_nodecay, time = time*365, y_position = y0, linestyle="-", color="blue", label = f"y={y0}m")
mbt.centerline(mbt_nodecay, time = time*365, y_position = y1, linestyle="--", color="orange", label = f"y={y1}m")
mbt.centerline(mbt_nodecay, time = time*365, y_position = y2, linestyle="-.", color="green", label = f"y={y2}m")
mbt.centerline(mbt_nodecay, time = time*365, y_position = y3, linestyle=":", color="red", label = f"y={y3}m")
plt.legend()
plt.title("Mibitrans with no degradation at various plume locations $y$",fontsize =10)
plt.xlabel("Plume travel distance $x$ in m")
plt.ylabel(f"Concentration distribution $C(x,y,t={time}y) [g/m^3]$")
plt.show()

time = 5 #[years] y0, y1, y2, y3 = 0, 6, 9, 20 mbt.centerline(mbt_nodecay, time = time*365, y_position = y0, linestyle="-", color="blue", label = f"y={y0}m") mbt.centerline(mbt_nodecay, time = time*365, y_position = y1, linestyle="--", color="orange", label = f"y={y1}m") mbt.centerline(mbt_nodecay, time = time*365, y_position = y2, linestyle="-.", color="green", label = f"y={y2}m") mbt.centerline(mbt_nodecay, time = time*365, y_position = y3, linestyle=":", color="red", label = f"y={y3}m") plt.legend() plt.title("Mibitrans with no degradation at various plume locations $y$",fontsize =10) plt.xlabel("Plume travel distance $x$ in m") plt.ylabel(f"Concentration distribution $C(x,y,t={time}y) [g/m^3]$") plt.show()

This also works with the use of the inherent plotting class method of the result objects. Loop to make combining plots even easier.

In [41]:

linestyles = ["-", "--", "-.", ":", "-"]
cmap = plt.get_cmap('Dark2')
colors = cmap.colors   # list of RGB tuples

for i,time in enumerate(range(1,6)):
    ana_instant.centerline(time = time*365, color=colors[i], linestyle=linestyles[i], label=f"t={time}y")
plt.title("Anatrans instant reaction model, at various times")
plt.xlim((0, 150))
plt.legend()
plt.show()

linestyles = ["-", "--", "-.", ":", "-"] cmap = plt.get_cmap('Dark2') colors = cmap.colors # list of RGB tuples for i,time in enumerate(range(1,6)): ana_instant.centerline(time = time*365, color=colors[i], linestyle=linestyles[i], label=f"t={time}y") plt.title("Anatrans instant reaction model, at various times") plt.xlim((0, 150)) plt.legend() plt.show()

Transverse distribution¶

The plotting routine transverse visualizes the concentration distribution $C(x_s, y, t_s)$ at a plume cross section, i.e. in horizontal transverse direction at a specified distance $x_s$ and time $t_s$. By default the time is the last time step model_time as defined in the data class ModelParameters.

Plotting can be performed either as class method, i.e. model_object.transverse(), or by calling the function with the model object as argument, i.e. transverse(model_object).

In [42]:

mbt.transverse(mbt_lineardecay, x_position=50)
plt.show()

mbt.transverse(mbt_lineardecay, x_position=50) plt.show()

Plot modification¶

Modifications to transverse plots can be done similarly to the centerline plots.

In [43]:

mbt.transverse(ana_nodecay, x_position=20, time=365, linestyle="--", color="green", label="no degradation")
mbt.transverse(ana_instant, x_position=20, time=365, linestyle="-.", color="red", label="instant reaction")
plt.title("No degradation and instant reaction Domenico model at x=80m, t=6years")
plt.legend()
plt.show()

mbt.transverse(ana_nodecay, x_position=20, time=365, linestyle="--", color="green", label="no degradation") mbt.transverse(ana_instant, x_position=20, time=365, linestyle="-.", color="red", label="instant reaction") plt.title("No degradation and instant reaction Domenico model at x=80m, t=6years") plt.legend() plt.show()

In [44]:

mbt_instant.transverse(5, label="x = 5m", color="darkolivegreen")
mbt_instant.transverse(20, label="x = 20m", color="yellowgreen")
mbt_instant.transverse(50, label="x = 100m", color="greenyellow")
mbt_instant.transverse(100, label="x = 300m", color="palegreen")
plt.xlim([-20,20])
plt.legend()
plt.title("Plume cross section for Mibitrans inst reaction model, $t=20$y", fontsize = 10)
plt.show()

mbt_instant.transverse(5, label="x = 5m", color="darkolivegreen") mbt_instant.transverse(20, label="x = 20m", color="yellowgreen") mbt_instant.transverse(50, label="x = 100m", color="greenyellow") mbt_instant.transverse(100, label="x = 300m", color="palegreen") plt.xlim([-20,20]) plt.legend() plt.title("Plume cross section for Mibitrans inst reaction model, $t=20$y", fontsize = 10) plt.show()

Breakthrough curve¶

Breakthrough curve shows the concentration over time for a specified point in the contaminant plume. Specifically, breakthrough visualizes the concentration distribution $C(x_s, y_s, t)$ at the plume location ($x_s$,$y_s$) over time. By default, the fixed vertical location is the center of the plume ($y_s=0$). Total displayed time corresponds to the last time step model_time as defined in the data class ModelParameters and intermediate points reflect the defined time stepping dt.

Plotting can be performed either as class method, i.e. model_object.breakthrough(), or by calling the function with the model object as argument, i.e. breakthrough(model_object). In the latter form, it can accept a list of models to be displayed in a single plot. Calling each function separately per model gives more control over plot layout. Keyword arguments for plt.plot can be passed on through the function.

In [45]:

# Concentration distribution can also be plotted in the transverse direction
mbt.breakthrough([bio_nodecay, bio_lineardecay, bio_instant], x_position=80,
             legend_names=["no degradation", "linear decay", "instant reaction"])
plt.title("Compare degredation options for Bioscreen model, BTC at x=80m",fontsize = 10)
plt.show()

# Concentration distribution can also be plotted in the transverse direction mbt.breakthrough([bio_nodecay, bio_lineardecay, bio_instant], x_position=80, legend_names=["no degradation", "linear decay", "instant reaction"]) plt.title("Compare degredation options for Bioscreen model, BTC at x=80m",fontsize = 10) plt.show()

In [46]:

mbt_nodecay.breakthrough(x_position=50,marker = 'o',mec = 'k', color="red")
plt.show()

mbt_nodecay.breakthrough(x_position=50,marker = 'o',mec = 'k', color="red") plt.show()

Two-dimensional plume¶

The function plume_2d visualizes the concentration distribution $C(x,y,t_s)$ at a desired point in time $t_s$ as color-coded 2D image. By default, concentration distribution of the last time step model_time is shown, as defined in the data class ModelParameters. Total displayed dimensions corresponds to the domain parameters defined in the data class ModelParameters.

Plotting can be performed either as class method, i.e. model_object.plume_2d(), or by calling the function with the model object as argument, i.e. plume_2d(model_object). Keyword arguments for plt.pcolormesh (on which plume_2d is based) can be passed on through the function.

In [47]:

mbt.plume_2d(mbt_nodecay)
plt.show()

mbt.plume_2d(mbt_nodecay) plt.show()

Plot modification: You can pass keyword arguments of plt.colormesh() for plot modification and use regular commands of matplotlib.pyplot for plot modification.

In [48]:

# Function passes plt.colormesh keyword arguments
mbt.plume_2d(mbt_lineardecay, cmap="plasma")
plt.ylim([-15,15])
plt.show()

# Function passes plt.colormesh keyword arguments mbt.plume_2d(mbt_lineardecay, cmap="plasma") plt.ylim([-15,15]) plt.show()

Alternatively, it can be accessed through the class method:

In [49]:

# Once again also can be accessed through the class method
ana_instant.plume_2d(time=5*365, cmap="Spectral", shading="gouraud")
plt.title("Anatrans instant reaction model, at t=5y")
plt.show()

# Once again also can be accessed through the class method ana_instant.plume_2d(time=5*365, cmap="Spectral", shading="gouraud") plt.title("Anatrans instant reaction model, at t=5y") plt.show()

3-dimensional plume¶

The concentration distribution can be visualized as plume in 3D view using plume_3d. Specifically, plume_3d shows the distribution $C(x,y,t_s)$ at desired point in time $t_s$ in a 3D plot, with the vertical axis representing the concentration. By default, $t_s$ is the last time step. Total displayed dimensions corresponds to the domain parameters defined in the data class ModelParameters.

Plotting can be performed either as class method, i.e. model_object.plume_3d(), or by calling the function with the model object as argument, i.e. plume_3d(model_object). Keyword arguments for plt.plot_surface (on which plume_3d is based) can be passed on through the function.

In [50]:

# Plot the x and y concentration distribution for no degradation decay model, uses plot_surface
ana_nodecay.plume_3d(time=5*365)
plt.title("Contaminant plume with Anatrans no degradation  model, t = 5 years")
plt.show()

# Plot the x and y concentration distribution for no degradation decay model, uses plot_surface ana_nodecay.plume_3d(time=5*365) plt.title("Contaminant plume with Anatrans no degradation model, t = 5 years") plt.show()

Alternatively, by calling the function with the model object and passing plt.plot_surface() keyword arguments:

In [51]:

# Function passes plot_surface keyword arguments
mbt.plume_3d(ana_lineardecay, cmap="summer")
plt.title("Contaminant plume with linear degradation Anatrans model, t = 5 years")
plt.show()

# Function passes plot_surface keyword arguments mbt.plume_3d(ana_lineardecay, cmap="summer") plt.title("Contaminant plume with linear degradation Anatrans model, t = 5 years") plt.show()

Adapting view point of plot¶

Functionality of plume_3d allows modification to the view point of the plot similarly to plt.plot_surface. Call the function and store the returned ax object. Adaptions are performed on ax as outlined below. See the documentation of plt.plot_surface for more information.

In [52]:

ax = ana_instant.plume_3d(cmap="cividis")
plt.title("Contaminant plume with instant reaction Anatrans model, t = 5 years")
ax.view_init(elev=15, azim=340)
plt.show()

ax = ana_instant.plume_3d(cmap="cividis") plt.title("Contaminant plume with instant reaction Anatrans model, t = 5 years") ax.view_init(elev=15, azim=340) plt.show()

Animate plots¶

All plots have the option to be animated, producing an animation that can be saved as file. Multiple models can be combined in a single animation, but all models should have the same spatial and temporal discretization.

Animations unfortunately do not work well in Jupyter Notebooks, in combination with other plots. For examples how to make animations in mibitrans, see animate_plume_plots.py

Mass balance¶

The mibitrans package provides the feature to calculate mass balances, specifically mass within the source and mass within the plume. It operates as method mass_balance on the results object of the model of choice (in our example we use the results of the Mibitrans model with the instantaneous reaction model, mbt_instant.) Each of the mass balance elements is a class property.

Examples of mass balance elements are:

• plume mass for given point in time mb.plume_mass
• source mass for given time mb.source_mass
• contaminant mass that has left the source zone mb.delta_source
• degraded mass over time mb.degraded_mass
• change in electron acceptor mass due to biodegradation in instantaneous reaction model mb.electron_acceptor_change

Let's calculate the mass balance for a specified point in time:

In [53]:

# Get mass balance for a specified point in time
mb = mbt_lineardecay.mass_balance(time=5*365)
print(mb.plume_mass)

# Get mass balance for a specified point in time mb = mbt_lineardecay.mass_balance(time=5*365) print(mb.plume_mass)

2272.221652649384

Call the mass balance object to change the time of the mass balance. If time is set to all (or call is left empty), mass balance will instead show all time steps.

In [54]:

mb(time="all")
print(mb.plume_mass)

mb(time="all") print(mb.plume_mass)

[ 691.4787647  1090.64208574 1383.13013161 1602.33546849 1767.65749881
 1892.59633504 1987.05717549 2058.4508047  2112.36879416 2153.04629171
 2183.69601949 2206.75585543 2224.07505259 2237.05539427 2246.75858639
 2253.98806935 2259.35133886 2263.30738137 2266.20273712 2268.29887981
 2269.7929741  2270.8335929  2271.53260962 2271.97419823 2272.22165265]

These mass balance elements can be used for various plots. Before doing so, let's also calculate the mass balance for the no-decay degredation option for comparison in the plot.

In [55]:

mb_nd = mbt_nodecay.mass_balance()

mb_nd = mbt_nodecay.mass_balance()

In [56]:

plt.plot(mb.t, mb_nd.plume_mass, color="C1",lw = 4, label=r"plume mass, no decay")
plt.plot(mb.t, mb.plume_mass, color="C0", lw = 3, label=r"plume mass, decay")
plt.plot(mb.t, mb.degraded_mass, color="C0",lw = 3,ls = ':', label=r"degraded mass, decay")
plt.plot(mb.t, mb.plume_mass + mb.degraded_mass, color="C0",lw = 2, linestyle="--", label=r"plume+degraded mass, decay")
plt.xlabel("Time [days]")
plt.ylabel("Contaminant mass [g]")
plt.title("Mass balance over time")
plt.legend()
plt.show()

plt.plot(mb.t, mb_nd.plume_mass, color="C1",lw = 4, label=r"plume mass, no decay") plt.plot(mb.t, mb.plume_mass, color="C0", lw = 3, label=r"plume mass, decay") plt.plot(mb.t, mb.degraded_mass, color="C0",lw = 3,ls = ':', label=r"degraded mass, decay") plt.plot(mb.t, mb.plume_mass + mb.degraded_mass, color="C0",lw = 2, linestyle="--", label=r"plume+degraded mass, decay") plt.xlabel("Time [days]") plt.ylabel("Contaminant mass [g]") plt.title("Mass balance over time") plt.legend() plt.show()

Here, the entire plume fits neatly in the model domain. However, the mass balance does not give representative results if the contaminant plume does not fit in the model domain (i.e. at the boundaries of the domain, concentrations are more than 1% of the maximum initial source concentration), since the mass balance is only calculated inside the model domain. In such case, creating the mass balance object will display a warning.

Sources¶

Domenico, P. A., An analytical model for multidimensional transport of a decaying contaminant species, Journal of Hydrology, 91 (1), 49–58, doi:10.1016/0022-1694(87)90127-2, 1987.

Karanovic, M., C. J. Neville, and C. B. Andrews, BIOSCREEN-AT: BIOSCREEN with an exact analytical solution, Groundwater, 45 (2), 242–245, doi:10.1111/j.1745-6584.2006.00296.x, 2007.

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

Newell, C. J., R. K. McLeod, and J. R. Gonzales, BIOSCREEN natural attenuation decision support system version 1.4 revisions, Tech. rep., U.S. EPA, 1997.

Wexler, E. J., 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. ; doi:10.3133/twri03B7, 1992.