In [1]:

import matplotlib.pyplot as plt
import numpy as np
import mibitrans as mbt
from mibitrans.analysis import differences as diff

import matplotlib.pyplot as plt import numpy as np import mibitrans as mbt from mibitrans.analysis import differences as diff

Exact solution versus untruncated solution¶

The mibitrans solution uses the exact solution of the 3D transport ADE, as described in Wexler et al. (1992). The anatrans solution approximates this solution by the method described in Domenico (1987), but without truncating the x part of the equation. Both are adapted for multiple source zones and source depletion, see equations below.

Anatrans solution¶

\begin{align}\tag{1} C(x, y, t) &= \sum_{i=1}^{n}\left\{ \frac{C^*_{0,i}}{8} \exp \left(-\gamma_s t\right) \right. \\ &\quad \cdot \left( \exp \left[ \frac{x\left(1-P\right)}{2\alpha_x}\right] \cdot \operatorname{erfc} \left[ \frac{x - Pvt}{2\sqrt{\alpha_x vt }} \right] \right.\\ &\quad \: + \left. \exp \left[ \frac{x\left(1+P\right)}{2\alpha_x}\right] \cdot \operatorname{erfc} \left[ \frac{x + Pvt}{2\sqrt{\alpha_x vt }} \right] \right) \\ &\quad \cdot \left( \operatorname{erf} \left[ \frac{y + Y_i/2}{2\sqrt{\alpha_y x}} \right] - \operatorname{erf} \left[ \frac{y - Y_i/2}{2\sqrt{\alpha_y x)}} \right] \right) \\ &\quad \cdot \left. \left( \operatorname{erf} \left[ \frac{Z}{2\sqrt{\alpha_z x)}} \right] - \operatorname{erf} \left[ \frac{-Z}{2\sqrt{\alpha_z x}} \right] \right) \right\} \\ & \text{with} \quad P = \sqrt{1+4\left(\mu - \gamma_s \right) \alpha_x/v} \end{align}

Mibitrans solution¶

$$ \begin{equation}\tag{2} \begin{aligned} C(x,y,t) &= \sum_{i=1}^{n}\left(C^*_{0,i}\frac{x}{8\sqrt{\pi \alpha_{x}\frac{v}{R}}}\exp(-\gamma_s t) \right. \\ &\quad \cdot \int_{0}^{t}\left[\frac{1}{\tau^{\frac{3}{2}}} \exp\left((\gamma_s - \mu)\tau - \frac{(x-\frac{v\tau}{R})^2}{4\alpha_{x}\frac{v\tau}{R}}\right) \right. \\ &\quad \cdot \left\{\operatorname{erfc}\left(\frac{y-Y_i / 2}{2 \sqrt{\alpha_{y}\frac{v\tau}{R}}}\right)-\operatorname{erfc}\left(\frac{y+Y_i / 2}{2 \sqrt{\alpha_{y}\frac{v\tau}{R}}}\right) \right\} \\ &\quad \left. \left. \cdot \left\{\operatorname{erfc}\left(\frac{-Z}{2 \sqrt{\alpha_{z}\frac{v\tau}{R}}}\right)-\operatorname{erfc}\left(\frac{Z}{2 \sqrt{\alpha_{z}\frac{v\tau}{R}}}\right) \right\}\right] d\tau \right) \end{aligned} \end{equation} $$

The approximation (see notebook example_validity_substitution, for more information) for the transverse dispersivity terms in the Anatrans solution introduces an error, of which the size depends on parameter choices. For the Domenico (1987) (Bioscreen in this package) solution specifically, the error size is dependent on the Peclet number defined as $Pe = \frac{vx_R}{D_x}$, where $x_R$ is the characteristic travel distance (Guyonnet and Neville, 2004). $D_x$ is the longitudinal dispersion coefficient, defined as $D_{eff} + \alpha_x v$. If diffusion is considered to be negligible, the expression becomes, $Pe = \frac{x_R}{\alpha_x}$, and thus is mainly dependent on the dispersivity. For the Anatrans solution however, the longitudinal component of the equation is nearly identical to that of the Mibitrans solution. Thus, if $\alpha_y$ and $\alpha_z = 0$, the Mibitrans and Anatrans solutions resolve to the same concentration distribution, which is shown below.

In [2]:

hydro = mbt.HydrologicalParameters(
    velocity=0.277, # Flow velocity [m/d]
    porosity=0.25,  # Effective soil porosity [-]
    alpha_x=10,     # Longitudinal dispersivity, in [m]
    alpha_y=0,      # Transverse horizontal dispersivity, in [m]
    alpha_z=0,    # Transverse vertical dispersivity, in [m]
    diffusion=0, # Molecular diffusion, in [m2/day]
)

att = mbt.AttenuationParameters(
    retardation=1,
    # Contaminant half life, in [days]
    half_life=0
)

source = mbt.SourceParameters(
    source_zone_boundary=np.array([10]),
    source_zone_concentration=np.array([11]),
    depth=2.5,
    total_mass="inf"
)

model = mbt.ModelParameters(
    # Model extent in the longitudinal (x) direction in [m].
    model_length = 800,
    # Model extent in the transverse horizontal (y) direction in [m].
    model_width = 150,
    # Model duration in [days].
    model_time = 10*365,
    # Model grid discretization step size in the longitudinal (x) direction, in [m].
    dx = 2,
    # Model grid discretization step size in the transverse horizontal (y) direction, in [m].
    dy = 1,
    # Model time discretization step size, in [days]
    dt = 365 / 5
)

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

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

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

hydro = mbt.HydrologicalParameters( velocity=0.277, # Flow velocity [m/d] porosity=0.25, # Effective soil porosity [-] alpha_x=10, # Longitudinal dispersivity, in [m] alpha_y=0, # Transverse horizontal dispersivity, in [m] alpha_z=0, # Transverse vertical dispersivity, in [m] diffusion=0, # Molecular diffusion, in [m2/day] ) att = mbt.AttenuationParameters( retardation=1, # Contaminant half life, in [days] half_life=0 ) source = mbt.SourceParameters( source_zone_boundary=np.array([10]), source_zone_concentration=np.array([11]), depth=2.5, total_mass="inf" ) model = mbt.ModelParameters( # Model extent in the longitudinal (x) direction in [m]. model_length = 800, # Model extent in the transverse horizontal (y) direction in [m]. model_width = 150, # Model duration in [days]. model_time = 10*365, # Model grid discretization step size in the longitudinal (x) direction, in [m]. dx = 2, # Model grid discretization step size in the transverse horizontal (y) direction, in [m]. dy = 1, # Model time discretization step size, in [days] dt = 365 / 5 ) mbt_object = mbt.Mibitrans(hydro, att, source, model) mbt_results_longitudinal = mbt_object.run() ana_object = mbt.Anatrans(hydro, att, source, model) ana_results_longitudinal = ana_object.run() bio_object = mbt.Bioscreen(hydro, att, source, model) bio_results_longitudinal = bio_object.run()

In [3]:

times=np.array([365, 2*365, 4*365, 6*365])
colors_ana = ["darkblue", "blue", "royalblue", "cornflowerblue", "skyblue"]
colors_mbt = ["darkgreen", "forestgreen", "limegreen", "greenyellow", "palegreen"]
colors_bio = ["maroon", "red", "orangered", "tomato", "coral"]

for i in range(len(times)):
    mbt_results_longitudinal.centerline(
        time=times[i],
        color=colors_mbt[i],
        lw=3,
        alpha=0.7,
        label=f"mbt, $t={times[i]}d$"
    )
    ana_results_longitudinal.centerline(
        time=times[i],
        color=colors_ana[i],
        lw=2,
        linestyle="--",
        label=f"ana, $t={times[i]}$d"
    )
    bio_results_longitudinal.centerline(
        time=times[i],
        color=colors_bio[i],
        lw=2,
        linestyle=":",
        label=f"bio, $t={times[i]}$d"
    )
    plt.legend(bbox_to_anchor=(1,1))
    plt.title("Temporal evolution of plume centerline for various mibitrans models")

plt.show()

times=np.array([365, 2*365, 4*365, 6*365]) colors_ana = ["darkblue", "blue", "royalblue", "cornflowerblue", "skyblue"] colors_mbt = ["darkgreen", "forestgreen", "limegreen", "greenyellow", "palegreen"] colors_bio = ["maroon", "red", "orangered", "tomato", "coral"] for i in range(len(times)): mbt_results_longitudinal.centerline( time=times[i], color=colors_mbt[i], lw=3, alpha=0.7, label=f"mbt, $t={times[i]}d$" ) ana_results_longitudinal.centerline( time=times[i], color=colors_ana[i], lw=2, linestyle="--", label=f"ana, $t={times[i]}$d" ) bio_results_longitudinal.centerline( time=times[i], color=colors_bio[i], lw=2, linestyle=":", label=f"bio, $t={times[i]}$d" ) plt.legend(bbox_to_anchor=(1,1)) plt.title("Temporal evolution of plume centerline for various mibitrans models") plt.show()

As shown, the Bioscreen solution has a slight error compared to the Anatrans and Mibitrans solutions, caused by the truncation in the Domenico (1987) solution. Thus, difference between Anatrans and Mibitrans is caused by the transverse terms.

In [4]:

hydro = mbt.HydrologicalParameters(
    velocity=0.277, # Flow velocity [m/d]
    porosity=0.25,  # Effective soil porosity [-]
    alpha_x=10,     # Longitudinal dispersivity, in [m]
    alpha_y=0.5,      # Transverse horizontal dispersivity, in [m]
    alpha_z=0.001,    # Transverse vertical dispersivity, in [m]
    diffusion=0, # Molecular diffusion, in [m2/day]
)

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

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

hydro = mbt.HydrologicalParameters( velocity=0.277, # Flow velocity [m/d] porosity=0.25, # Effective soil porosity [-] alpha_x=10, # Longitudinal dispersivity, in [m] alpha_y=0.5, # Transverse horizontal dispersivity, in [m] alpha_z=0.001, # Transverse vertical dispersivity, in [m] diffusion=0, # Molecular diffusion, in [m2/day] ) mbt_object = mbt.Mibitrans(hydro, att, source, model) mbt_results = mbt_object.run() ana_object = mbt.Anatrans(hydro, att, source, model) ana_results = ana_object.run()

In [5]:

times=np.array([365, 2*365, 4*365, 6*365])

for i in range(len(times)):
    mbt_results.centerline(
        time=times[i],
        color=colors_mbt[i],
        lw=3,
        label=f"mbt, $t={times[i]/365}y$"
    )
    ana_results.centerline(
        time=times[i],
        color=colors_ana[i],
        lw=2,
        alpha=0.8,
        linestyle="--",
        label=f"ana, $t={times[i]/365}$y"
    )
plt.title("Differences between Mibitrans and Anatrans models in plume centerline")
plt.legend()
plt.show()

times=np.array([365, 2*365, 4*365, 6*365]) for i in range(len(times)): mbt_results.centerline( time=times[i], color=colors_mbt[i], lw=3, label=f"mbt, $t={times[i]/365}y$" ) ana_results.centerline( time=times[i], color=colors_ana[i], lw=2, alpha=0.8, linestyle="--", label=f"ana, $t={times[i]/365}$y" ) plt.title("Differences between Mibitrans and Anatrans models in plume centerline") plt.legend() plt.show()

In [6]:

times=np.array([1*365, 2*365, 3*365, 4*365, 5*365])

for i in range(len(times)):
    mbt_results.transverse(
        x_position = 300,
        time=times[i],
        color=colors_mbt[i],
        lw=3,
        label=f"mbt, $t={times[i]/365}y$"
    )
    ana_results.transverse(
        x_position = 300,
        time=times[i],
        color=colors_ana[i],
        lw=2,
        linestyle="--",
        label=f"ana, $t={times[i]/365}y$"
    )

plt.ylim((0, 5))
plt.xlim((-50,50))
plt.title("Differences between Mibitrans and Anatrans models transverse to plume")
plt.legend(bbox_to_anchor=(1,1))
plt.show()

times=np.array([1*365, 2*365, 3*365, 4*365, 5*365]) for i in range(len(times)): mbt_results.transverse( x_position = 300, time=times[i], color=colors_mbt[i], lw=3, label=f"mbt, $t={times[i]/365}y$" ) ana_results.transverse( x_position = 300, time=times[i], color=colors_ana[i], lw=2, linestyle="--", label=f"ana, $t={times[i]/365}y$" ) plt.ylim((0, 5)) plt.xlim((-50,50)) plt.title("Differences between Mibitrans and Anatrans models transverse to plume") plt.legend(bbox_to_anchor=(1,1)) plt.show()

In [7]:

x_positions = [75, 150, 300, 600]

for i in range(len(x_positions)):
    mbt_results.breakthrough(
        x_position=x_positions[i],
        color=colors_mbt[i],
        lw=3,
        label=f"mbt, $x={x_positions[i]}m$"
    )
    ana_results.breakthrough(
        x_position=x_positions[i],
        color=colors_ana[i],
        lw=2,
        linestyle="--",
        label=f"ana, $x={x_positions[i]}m$"
    )

plt.title("Breakthrough curve of Mibitrans and Anatrans model, at various x locations.")
plt.legend(bbox_to_anchor=(1,1))
plt.show()

x_positions = [75, 150, 300, 600] for i in range(len(x_positions)): mbt_results.breakthrough( x_position=x_positions[i], color=colors_mbt[i], lw=3, label=f"mbt, $x={x_positions[i]}m$" ) ana_results.breakthrough( x_position=x_positions[i], color=colors_ana[i], lw=2, linestyle="--", label=f"ana, $x={x_positions[i]}m$" ) plt.title("Breakthrough curve of Mibitrans and Anatrans model, at various x locations.") plt.legend(bbox_to_anchor=(1,1)) plt.show()

Under these transport parameters, the Anatrans model consistently underestimates the contaminant concentrations at the plume centerline. The transverse distribution shows that Anatrans as increased concentrations compared to Mibitrans at the fringes of the plume. Thus, the approximation used in the Anatrans solution overestimates the effect of transverse dispersion on the concentration distribution.

Comparing model differences for various flow velocities and dispersivities¶

Assuming that diffusion is negligible, the flow velocity is of little influence on the differences between the Anatrans and Mibitrans models.

In [8]:

hydro = mbt.HydrologicalParameters(
    velocity=0.1, # Flow velocity [m/d]
    porosity=0.25,  # Effective soil porosity [-]
    alpha_x=5,     # Longitudinal dispersivity, in [m]
    alpha_y=0.5,      # Transverse horizontal dispersivity, in [m]
    alpha_z=0.001,    # Transverse vertical dispersivity, in [m]
    diffusion=0, # Molecular diffusion, in [m2/day]
)

mbt_object = mbt.Mibitrans(hydro, att, source, model)
mbt_object.hydrological_parameters.velocity = 0.1
mbt_results_v01 = mbt_object.run()

ana_object = mbt.Anatrans(hydro, att, source, model)
ana_object.hydrological_parameters.velocity = 0.1
ana_results_v01 = ana_object.run()

mbt_object.hydrological_parameters.velocity = 0.2
mbt_results_v02 = mbt_object.run()

ana_object.hydrological_parameters.velocity = 0.2
ana_results_v02 = ana_object.run()

hydro = mbt.HydrologicalParameters( velocity=0.1, # Flow velocity [m/d] porosity=0.25, # Effective soil porosity [-] alpha_x=5, # Longitudinal dispersivity, in [m] alpha_y=0.5, # Transverse horizontal dispersivity, in [m] alpha_z=0.001, # Transverse vertical dispersivity, in [m] diffusion=0, # Molecular diffusion, in [m2/day] ) mbt_object = mbt.Mibitrans(hydro, att, source, model) mbt_object.hydrological_parameters.velocity = 0.1 mbt_results_v01 = mbt_object.run() ana_object = mbt.Anatrans(hydro, att, source, model) ana_object.hydrological_parameters.velocity = 0.1 ana_results_v01 = ana_object.run() mbt_object.hydrological_parameters.velocity = 0.2 mbt_results_v02 = mbt_object.run() ana_object.hydrological_parameters.velocity = 0.2 ana_results_v02 = ana_object.run()

In [9]:

mbt_results_v01.centerline(time=3*365, color=colors_mbt[0], lw=3, alpha=0.7, label="mbt v=0.1, t=3y")
mbt_results_v02.centerline(time=3*365, color=colors_mbt[1], lw=3, alpha=0.7, label="mbt v=0.2, t=3y")
mbt_results_v01.centerline(time=6*365, color=colors_mbt[2], lw=2, linestyle="--", label="mbt v=0.1, t=6y")
mbt_results_v02.centerline(time=6*365, color=colors_mbt[3], lw=3, linestyle="--", label="mbt v=0.2, t=6y")
ana_results_v01.centerline(time=3*365, color=colors_ana[0], lw=3, alpha=0.7, label="ana v=0.1, t=3y")
ana_results_v02.centerline(time=3*365, color=colors_ana[1], lw=3, alpha=0.7, label="ana v=0.2, t=3y")
ana_results_v01.centerline(time=6*365, color=colors_ana[3], lw=2, linestyle="--", label="ana v=0.1, t=6y")
ana_results_v02.centerline(time=6*365, color=colors_ana[4], lw=3, linestyle="--", label="ana v=0.2, t=6y")

plt.xlim(-10, 600)
plt.title("Differences between Mibitrans and Anatrans for various flow velocities")
plt.legend()
plt.show()

mbt_results_v01.centerline(time=3*365, color=colors_mbt[0], lw=3, alpha=0.7, label="mbt v=0.1, t=3y") mbt_results_v02.centerline(time=3*365, color=colors_mbt[1], lw=3, alpha=0.7, label="mbt v=0.2, t=3y") mbt_results_v01.centerline(time=6*365, color=colors_mbt[2], lw=2, linestyle="--", label="mbt v=0.1, t=6y") mbt_results_v02.centerline(time=6*365, color=colors_mbt[3], lw=3, linestyle="--", label="mbt v=0.2, t=6y") ana_results_v01.centerline(time=3*365, color=colors_ana[0], lw=3, alpha=0.7, label="ana v=0.1, t=3y") ana_results_v02.centerline(time=3*365, color=colors_ana[1], lw=3, alpha=0.7, label="ana v=0.2, t=3y") ana_results_v01.centerline(time=6*365, color=colors_ana[3], lw=2, linestyle="--", label="ana v=0.1, t=6y") ana_results_v02.centerline(time=6*365, color=colors_ana[4], lw=3, linestyle="--", label="ana v=0.2, t=6y") plt.xlim(-10, 600) plt.title("Differences between Mibitrans and Anatrans for various flow velocities") plt.legend() plt.show()

The curves for $v=0.1m/s$ at $t=6y$ and $v=0.2m/s$ at $t=3y$ of Mibitrans are expected to overlap. Diffusion $=0$ and thus with twice as low velocity, it takes twice as long for the advective front to get into the same position. The distributions for Anatrans also show this overlap, implying that regardless of flow velocity, the difference between Anatrans and Mibitrans remain the same.

However, changing the transverse dispersivity does change the differences between the models, as is shown below.

In [10]:

alpha_x_list = [10,5,2,1]
# Varying transverse horizontal dispersivities
alpha_y_list = [0.5, 0.1, 0.05, 0.01]
# And corresponding transverse vertical dispersivities
alpha_z_list = [0.01, 0.005, 0.002, 0.0005]
output_mbt = []
output_ana = []

mbt_object_alpha = mbt.Mibitrans(hydro, att, source, model)
ana_object_alpha = mbt.Anatrans(hydro, att, source, model)

for i in range(len(alpha_y_list)):
    hydro = mbt.HydrologicalParameters(
        velocity=0.277,             # [m/d]
        porosity=0.25,              # [-]
        alpha_x=5,#alpha_x_list[i],                 # [m]
        alpha_y=alpha_y_list[i],    # [m]
        alpha_z=alpha_z_list[i],    # [m]
        diffusion=0,                # [m2/day]
    )

    mbt_object_alpha.hydrological_parameters = hydro
    mbt_results_alpha = mbt_object_alpha.run()

    ana_object_alpha.hydrological_parameters = hydro
    ana_results_alpha = ana_object_alpha.run()

    output_mbt.append(mbt_results_alpha)
    output_ana.append(ana_results_alpha)

alpha_x_list = [10,5,2,1] # Varying transverse horizontal dispersivities alpha_y_list = [0.5, 0.1, 0.05, 0.01] # And corresponding transverse vertical dispersivities alpha_z_list = [0.01, 0.005, 0.002, 0.0005] output_mbt = [] output_ana = [] mbt_object_alpha = mbt.Mibitrans(hydro, att, source, model) ana_object_alpha = mbt.Anatrans(hydro, att, source, model) for i in range(len(alpha_y_list)): hydro = mbt.HydrologicalParameters( velocity=0.277, # [m/d] porosity=0.25, # [-] alpha_x=5,#alpha_x_list[i], # [m] alpha_y=alpha_y_list[i], # [m] alpha_z=alpha_z_list[i], # [m] diffusion=0, # [m2/day] ) mbt_object_alpha.hydrological_parameters = hydro mbt_results_alpha = mbt_object_alpha.run() ana_object_alpha.hydrological_parameters = hydro ana_results_alpha = ana_object_alpha.run() output_mbt.append(mbt_results_alpha) output_ana.append(ana_results_alpha)

In [11]:

for i in range(len(alpha_y_list)):
    output_mbt[i].centerline(
        time=5*365,
        color=colors_mbt[i],
        lw=2.5,
        label=fr"mbt, $\alpha_y={alpha_y_list[i]}m, \alpha_z={alpha_z_list[i]}m$"
    )
    output_ana[i].centerline(
        time=5*365,
        color=colors_ana[i],
        linestyle="--",
        lw=2,
        label=fr"ana, $\alpha_y={alpha_y_list[i]}m, \alpha_z={alpha_z_list[i]}m$"
    )
plt.title("Plume centerline of Mibitrans and Anatrans models for various transverse dispersivities")
plt.legend(bbox_to_anchor=(1,1))
plt.xlim((-10, 700))

for i in range(len(alpha_y_list)): output_mbt[i].centerline( time=5*365, color=colors_mbt[i], lw=2.5, label=fr"mbt, $\alpha_y={alpha_y_list[i]}m, \alpha_z={alpha_z_list[i]}m$" ) output_ana[i].centerline( time=5*365, color=colors_ana[i], linestyle="--", lw=2, label=fr"ana, $\alpha_y={alpha_y_list[i]}m, \alpha_z={alpha_z_list[i]}m$" ) plt.title("Plume centerline of Mibitrans and Anatrans models for various transverse dispersivities") plt.legend(bbox_to_anchor=(1,1)) plt.xlim((-10, 700))

Out[11]:

(-10.0, 700.0)

In [12]:

x_pos = 250

for i in range(len(alpha_y_list)):
    output_mbt[i].transverse(
        x_position=x_pos,
        color=colors_mbt[i],
        lw=2.5,
        label=fr"mbt, $\alpha_y={alpha_y_list[i]}m, \alpha_z={alpha_z_list[i]}m$"
    )
    output_ana[i].transverse(
        x_position=x_pos,
        color=colors_ana[i],
        linestyle=":",
        lw=2,
        label=fr"ana, $\alpha_y={alpha_y_list[i]}m, \alpha_z={alpha_z_list[i]}m$"
    )

plt.legend(bbox_to_anchor=(1,1))
plt.title("Transverse plot of Mibitrans and Anatrans models for various transverse dispersivities")
plt.xlim(-60,60)
plt.ylim(-0.2,11.5)
plt.show()

x_pos = 250 for i in range(len(alpha_y_list)): output_mbt[i].transverse( x_position=x_pos, color=colors_mbt[i], lw=2.5, label=fr"mbt, $\alpha_y={alpha_y_list[i]}m, \alpha_z={alpha_z_list[i]}m$" ) output_ana[i].transverse( x_position=x_pos, color=colors_ana[i], linestyle=":", lw=2, label=fr"ana, $\alpha_y={alpha_y_list[i]}m, \alpha_z={alpha_z_list[i]}m$" ) plt.legend(bbox_to_anchor=(1,1)) plt.title("Transverse plot of Mibitrans and Anatrans models for various transverse dispersivities") plt.xlim(-60,60) plt.ylim(-0.2,11.5) plt.show()

In [13]:

x_pos = 250

for i in range(len(alpha_y_list)):
    output_mbt[i].breakthrough(
        x_position=x_pos,
        color=colors_mbt[i],
        lw=2.5,
        label=fr"mbt, $\alpha_y={alpha_y_list[i]}m, \alpha_z={alpha_z_list[i]}m$"
    )
    output_ana[i].breakthrough(
        x_position=x_pos,
        color=colors_ana[i],
        linestyle="--",
        lw=2,
        label=fr"ana, $\alpha_y={alpha_y_list[i]}m, \alpha_z={alpha_z_list[i]}m$"
    )

plt.legend(bbox_to_anchor=(1.01,0.6))
plt.xlim(-10,2000)
plt.ylim(-0.2,12)
plt.title("BTC of Mibitrans and Anatrans models for various transverse dispersivities")
plt.show()

x_pos = 250 for i in range(len(alpha_y_list)): output_mbt[i].breakthrough( x_position=x_pos, color=colors_mbt[i], lw=2.5, label=fr"mbt, $\alpha_y={alpha_y_list[i]}m, \alpha_z={alpha_z_list[i]}m$" ) output_ana[i].breakthrough( x_position=x_pos, color=colors_ana[i], linestyle="--", lw=2, label=fr"ana, $\alpha_y={alpha_y_list[i]}m, \alpha_z={alpha_z_list[i]}m$" ) plt.legend(bbox_to_anchor=(1.01,0.6)) plt.xlim(-10,2000) plt.ylim(-0.2,12) plt.title("BTC of Mibitrans and Anatrans models for various transverse dispersivities") plt.show()

Differences are better visualized as the relative and/or absolute difference between the two models.

In [14]:

colors = ["maroon", "red", "orangered", "coral"]
time_step = len(mbt_object.t)//2 - 1
y_position = len(mbt_object.y)//2
print("plotting at t =", mbt_object.t[time_step] / 365, "years and a y-position =", mbt_object.y[y_position], "m")

for i in range(len(alpha_y_list)):
    diff_ana_mbt = (output_ana[i].cxyt - output_mbt[i].cxyt) / output_mbt[i].cxyt
    plt.plot(
        mbt_object.x,
        diff_ana_mbt[time_step,y_position,:],
        color=colors[i],
        lw=2,
        label=fr"$\alpha_y={alpha_y_list[i]}m, \alpha_z={alpha_z_list[i]}m$"
    )

plt.xlabel("Distance from source [m]")
plt.ylabel("Relative difference")
plt.title("Relative difference between Mibitrans and Anatrans models")
plt.legend()
plt.show()

for i in range(len(alpha_y_list)):
    diff_ana_mbt = output_ana[i].cxyt - output_mbt[i].cxyt
    plt.plot(
        mbt_object.x,
        diff_ana_mbt[time_step,y_position,:],
        color=colors[i],
        lw=2,
        label=fr"$\alpha_y={alpha_y_list[i]}m, \alpha_z={alpha_z_list[i]}m$"
    )

plt.xlabel("Distance from source [m]")
plt.ylabel("Absolute difference")
plt.title("Absolute difference between Mibitrans and Anatrans models")
plt.legend()
plt.show()

colors = ["maroon", "red", "orangered", "coral"] time_step = len(mbt_object.t)//2 - 1 y_position = len(mbt_object.y)//2 print("plotting at t =", mbt_object.t[time_step] / 365, "years and a y-position =", mbt_object.y[y_position], "m") for i in range(len(alpha_y_list)): diff_ana_mbt = (output_ana[i].cxyt - output_mbt[i].cxyt) / output_mbt[i].cxyt plt.plot( mbt_object.x, diff_ana_mbt[time_step,y_position,:], color=colors[i], lw=2, label=fr"$\alpha_y={alpha_y_list[i]}m, \alpha_z={alpha_z_list[i]}m$" ) plt.xlabel("Distance from source [m]") plt.ylabel("Relative difference") plt.title("Relative difference between Mibitrans and Anatrans models") plt.legend() plt.show() for i in range(len(alpha_y_list)): diff_ana_mbt = output_ana[i].cxyt - output_mbt[i].cxyt plt.plot( mbt_object.x, diff_ana_mbt[time_step,y_position,:], color=colors[i], lw=2, label=fr"$\alpha_y={alpha_y_list[i]}m, \alpha_z={alpha_z_list[i]}m$" ) plt.xlabel("Distance from source [m]") plt.ylabel("Absolute difference") plt.title("Absolute difference between Mibitrans and Anatrans models") plt.legend() plt.show()

plotting at t = 5.0 years and a y-position = 0.0 m

/tmp/ipykernel_2476/1395217995.py:7: RuntimeWarning: divide by zero encountered in divide
  diff_ana_mbt = (output_ana[i].cxyt - output_mbt[i].cxyt) / output_mbt[i].cxyt
/tmp/ipykernel_2476/1395217995.py:7: RuntimeWarning: invalid value encountered in divide
  diff_ana_mbt = (output_ana[i].cxyt - output_mbt[i].cxyt) / output_mbt[i].cxyt

In [15]:

time_step = len(mbt_object.t)//2 - 1

fig, ax = plt.subplots(2,2, sharex=True, sharey=True)
plot_list = []
i=0

for j in range(len(ax[:,0])):
    for k in range(len(ax[0,:])):
        diff_ana_mbt = output_ana[i].relative_cxyt - output_mbt[i].relative_cxyt
        plot_list.append(ax[j,k].pcolormesh(
            mbt_object.x,
            mbt_object.y,
            diff_ana_mbt[time_step,:,:],
            vmin=-0.15,
            vmax=0.15,
            cmap="seismic",
            shading="gouraud"
        ))
        ax[j,k].set_title(fr"$\alpha_y={alpha_y_list[i]}m, \alpha_z={alpha_z_list[i]}m$", fontsize=9)
        # To show differences more clearly, x-limit is 1/4 of model extent
        #ax[j,k].set_xlim(-1, mbt_object.x[-1]//4)
        # # Also limit y-extent
        #ax[j,k].set_ylim(mbt_object.y[len(mbt_object.y)//4], mbt_object.y[-len(mbt_object.y)//4])
        i+=1

fig.colorbar(plot_list[0], ax=ax, orientation="vertical", label="Absolute difference in relative concentration")
fig.supxlabel("Distance from source [m]", fontsize=10)
fig.supylabel("y-position [m]", fontsize=10)
fig.suptitle("Difference between Mibitrans and Anatrans model for various dispersivities", fontsize=11)
plt.savefig("mbt_ana_comp_cmesh.png", dpi=300)

time_step = len(mbt_object.t)//2 - 1 fig, ax = plt.subplots(2,2, sharex=True, sharey=True) plot_list = [] i=0 for j in range(len(ax[:,0])): for k in range(len(ax[0,:])): diff_ana_mbt = output_ana[i].relative_cxyt - output_mbt[i].relative_cxyt plot_list.append(ax[j,k].pcolormesh( mbt_object.x, mbt_object.y, diff_ana_mbt[time_step,:,:], vmin=-0.15, vmax=0.15, cmap="seismic", shading="gouraud" )) ax[j,k].set_title(fr"$\alpha_y={alpha_y_list[i]}m, \alpha_z={alpha_z_list[i]}m$", fontsize=9) # To show differences more clearly, x-limit is 1/4 of model extent #ax[j,k].set_xlim(-1, mbt_object.x[-1]//4) # # Also limit y-extent #ax[j,k].set_ylim(mbt_object.y[len(mbt_object.y)//4], mbt_object.y[-len(mbt_object.y)//4]) i+=1 fig.colorbar(plot_list[0], ax=ax, orientation="vertical", label="Absolute difference in relative concentration") fig.supxlabel("Distance from source [m]", fontsize=10) fig.supylabel("y-position [m]", fontsize=10) fig.suptitle("Difference between Mibitrans and Anatrans model for various dispersivities", fontsize=11) plt.savefig("mbt_ana_comp_cmesh.png", dpi=300)

Develop note: Below trying out if plotting this logarithmically would show clearer results. While they definitely look nice, I am not sure if they are useful yet.

In [16]:

time_step = len(mbt_object.t)//2 - 1
fig, ax = plt.subplots(2,2, sharex=True, sharey=True)
plot_list = []
i=0

for j in range(len(ax[:,0])):
    for k in range(len(ax[0,:])):
        diff_ana_mbt = output_ana[i].relative_cxyt - output_mbt[i].relative_cxyt
        plot_list.append(ax[j,k].pcolormesh(
            mbt_object.x,
            mbt_object.y,
            np.log(abs(diff_ana_mbt[time_step,:,:])),
            vmin=-20,
            vmax=0,
            cmap="inferno",
        ))
        ax[j,k].set_title(fr"$\alpha_y={alpha_y_list[i]}m, \alpha_z={alpha_z_list[i]}m$", fontsize=9)
        # To show differences more clearly, x-limit is 1/4 of model extent
        ax[j,k].set_xlim(-1, mbt_object.x[-1]//4)
        # # Also limit y-extent
        ax[j,k].set_ylim(mbt_object.y[len(mbt_object.y)//5], mbt_object.y[-len(mbt_object.y)//5])
        i+=1

fig.colorbar(plot_list[0], ax=ax, orientation="vertical", label="Absolute difference in relative concentration")
fig.supxlabel("Distance from source [m]", fontsize=10)
fig.supylabel("y-position [m]", fontsize=10)
fig.suptitle("Difference between Mibitrans and Anatrans model for various dispersivities", fontsize=11)
plt.show()

time_step = len(mbt_object.t)//2 - 1 fig, ax = plt.subplots(2,2, sharex=True, sharey=True) plot_list = [] i=0 for j in range(len(ax[:,0])): for k in range(len(ax[0,:])): diff_ana_mbt = output_ana[i].relative_cxyt - output_mbt[i].relative_cxyt plot_list.append(ax[j,k].pcolormesh( mbt_object.x, mbt_object.y, np.log(abs(diff_ana_mbt[time_step,:,:])), vmin=-20, vmax=0, cmap="inferno", )) ax[j,k].set_title(fr"$\alpha_y={alpha_y_list[i]}m, \alpha_z={alpha_z_list[i]}m$", fontsize=9) # To show differences more clearly, x-limit is 1/4 of model extent ax[j,k].set_xlim(-1, mbt_object.x[-1]//4) # # Also limit y-extent ax[j,k].set_ylim(mbt_object.y[len(mbt_object.y)//5], mbt_object.y[-len(mbt_object.y)//5]) i+=1 fig.colorbar(plot_list[0], ax=ax, orientation="vertical", label="Absolute difference in relative concentration") fig.supxlabel("Distance from source [m]", fontsize=10) fig.supylabel("y-position [m]", fontsize=10) fig.suptitle("Difference between Mibitrans and Anatrans model for various dispersivities", fontsize=11) plt.show()

/tmp/ipykernel_2476/568875887.py:12: RuntimeWarning: divide by zero encountered in log
  np.log(abs(diff_ana_mbt[time_step,:,:])),

In [17]:

for i in range(len(output_mbt)):
    plt.plot(output_mbt[i].t, diff.mean_absolute_difference(output_mbt[i].relative_cxyt, output_ana[i].relative_cxyt,
                                                          mean_axis=(1,2), concentration_cutoff=1e-5))
for i in range(len(output_mbt)):
    print(diff.mean_absolute_difference(output_mbt[i].cxyt, output_ana[i].cxyt, concentration_cutoff=1e-5))

for i in range(len(output_mbt)): plt.plot(output_mbt[i].t, diff.mean_absolute_difference(output_mbt[i].relative_cxyt, output_ana[i].relative_cxyt, mean_axis=(1,2), concentration_cutoff=1e-5)) for i in range(len(output_mbt)): print(diff.mean_absolute_difference(output_mbt[i].cxyt, output_ana[i].cxyt, concentration_cutoff=1e-5))

0.021556768635891243
0.02085333697799966
0.019325712944400733
0.013618783450521714

In [18]:

fig, ax = plt.subplots(2,2, sharex=True, sharey=True)
plot_list = []
i=0

for j in range(len(ax[:,0])):
    for k in range(len(ax[0,:])):
        #diff_ana_mbt = output_ana[i].relative_cxyt - output_mbt[i].relative_cxyt
        plot_list.append(ax[j,k].pcolormesh(
            mbt_object.x,
            mbt_object.t,
            diff.mean_absolute_difference(output_mbt[i].relative_cxyt, output_ana[i].relative_cxyt,
                                        mean_axis=1, concentration_cutoff=0),
            vmin=0,
            vmax=0.01,
            cmap="plasma",
            shading="gouraud"
        ))
        ax[j,k].set_title(fr"$\alpha_y={alpha_y_list[i]}m, \alpha_z={alpha_z_list[i]}m$", fontsize=9)
        # To show differences more clearly, x-limit is 1/4 of model extent
        #ax[j,k].set_xlim(-1, mbt_object.x[-1]//4)
        # # Also limit y-extent
        #ax[j,k].set_ylim(mbt_object.y[len(mbt_object.y)//5], mbt_object.y[-len(mbt_object.y)//5])
        i+=1

fig.colorbar(plot_list[0], ax=ax, orientation="vertical", label="Mean absolute difference in relative concentration")
fig.supxlabel("Distance from source [m]", fontsize=10)
fig.supylabel("time [d]", fontsize=10)
fig.suptitle("Difference between Mibitrans and Anatrans model for various dispersivities", fontsize=11)
plt.show()

fig, ax = plt.subplots(2,2, sharex=True, sharey=True) plot_list = [] i=0 for j in range(len(ax[:,0])): for k in range(len(ax[0,:])): #diff_ana_mbt = output_ana[i].relative_cxyt - output_mbt[i].relative_cxyt plot_list.append(ax[j,k].pcolormesh( mbt_object.x, mbt_object.t, diff.mean_absolute_difference(output_mbt[i].relative_cxyt, output_ana[i].relative_cxyt, mean_axis=1, concentration_cutoff=0), vmin=0, vmax=0.01, cmap="plasma", shading="gouraud" )) ax[j,k].set_title(fr"$\alpha_y={alpha_y_list[i]}m, \alpha_z={alpha_z_list[i]}m$", fontsize=9) # To show differences more clearly, x-limit is 1/4 of model extent #ax[j,k].set_xlim(-1, mbt_object.x[-1]//4) # # Also limit y-extent #ax[j,k].set_ylim(mbt_object.y[len(mbt_object.y)//5], mbt_object.y[-len(mbt_object.y)//5]) i+=1 fig.colorbar(plot_list[0], ax=ax, orientation="vertical", label="Mean absolute difference in relative concentration") fig.supxlabel("Distance from source [m]", fontsize=10) fig.supylabel("time [d]", fontsize=10) fig.suptitle("Difference between Mibitrans and Anatrans model for various dispersivities", fontsize=11) plt.show()

In [19]:

x_disp = np.array([10,5,2,1,0.5])
y_disp = np.array([0.5,0.2,0.1,0.05,0.02])
z_disp = np.array([0.01, 0.005, 0.002, 0.001, 0])

att_loop = mbt.AttenuationParameters(
    retardation=1,
    # Contaminant half life, in [days]
    half_life=5*365
)

source_loop = mbt.SourceParameters(
    source_zone_boundary=np.array([10]),
    source_zone_concentration=np.array([11]),
    depth=2.5,
    total_mass=1000000
)

model_loop = mbt.ModelParameters(
    # Model extent in the longitudinal (x) direction in [m].
    model_length = 1100,
    # Model extent in the transverse horizontal (y) direction in [m].
    model_width = 150,
    # Model duration in [days].
    model_time = 10*365,
    # Model grid discretization step size in the longitudinal (x) direction, in [m].
    dx = 2,
    # Model grid discretization step size in the transverse horizontal (y) direction, in [m].
    dy = 1,
    # Model time discretization step size, in [days]
    dt = 365
)

list_mbt = [[0]*len(x_disp) for _ in range(len(y_disp))]
list_ana = [[0]*len(x_disp) for _ in range(len(y_disp))]

mbt_object_disp = mbt.Mibitrans(hydro, att_loop, source_loop, model_loop)
ana_object_disp = mbt.Anatrans(hydro, att_loop, source_loop, model_loop)

# Loops for a total of 25 runs, will take about 1min.
for i in range(len(x_disp)):
    for j in range(len(y_disp)):
        hydro_disp = mbt.HydrologicalParameters(
            velocity=0.2,             # [m/d]
            porosity=0.25,              # [-]
            alpha_x=x_disp[i],#alpha_x_list[i],                 # [m]
            alpha_y=y_disp[j],    # [m]
            alpha_z=z_disp[j],    # [m]
            diffusion=0,                # [m2/day]
        )
        print(i,j) # Print to track progress
        mbt_object_disp.hydrological_parameters = hydro_disp
        list_mbt[j][i] = mbt_object_disp.run()
        ana_object_disp.hydrological_parameters = hydro_disp
        list_ana[j][i] = ana_object_disp.run()

x_disp = np.array([10,5,2,1,0.5]) y_disp = np.array([0.5,0.2,0.1,0.05,0.02]) z_disp = np.array([0.01, 0.005, 0.002, 0.001, 0]) att_loop = mbt.AttenuationParameters( retardation=1, # Contaminant half life, in [days] half_life=5*365 ) source_loop = mbt.SourceParameters( source_zone_boundary=np.array([10]), source_zone_concentration=np.array([11]), depth=2.5, total_mass=1000000 ) model_loop = mbt.ModelParameters( # Model extent in the longitudinal (x) direction in [m]. model_length = 1100, # Model extent in the transverse horizontal (y) direction in [m]. model_width = 150, # Model duration in [days]. model_time = 10*365, # Model grid discretization step size in the longitudinal (x) direction, in [m]. dx = 2, # Model grid discretization step size in the transverse horizontal (y) direction, in [m]. dy = 1, # Model time discretization step size, in [days] dt = 365 ) list_mbt = [[0]*len(x_disp) for _ in range(len(y_disp))] list_ana = [[0]*len(x_disp) for _ in range(len(y_disp))] mbt_object_disp = mbt.Mibitrans(hydro, att_loop, source_loop, model_loop) ana_object_disp = mbt.Anatrans(hydro, att_loop, source_loop, model_loop) # Loops for a total of 25 runs, will take about 1min. for i in range(len(x_disp)): for j in range(len(y_disp)): hydro_disp = mbt.HydrologicalParameters( velocity=0.2, # [m/d] porosity=0.25, # [-] alpha_x=x_disp[i],#alpha_x_list[i], # [m] alpha_y=y_disp[j], # [m] alpha_z=z_disp[j], # [m] diffusion=0, # [m2/day] ) print(i,j) # Print to track progress mbt_object_disp.hydrological_parameters = hydro_disp list_mbt[j][i] = mbt_object_disp.run() ana_object_disp.hydrological_parameters = hydro_disp list_ana[j][i] = ana_object_disp.run()

0 0

0 1

0 2

0 3

0 4

1 0

1 1

1 2

1 3

1 4

2 0

2 1

2 2

2 3

2 4

3 0

3 1

3 2

3 3

3 4

4 0

4 1

4 2

4 3

4 4

Plot to see if plume extends reasonably in model boundary

Sensitivity plots more explored and explained in the model_sensitivities.ipynb notebook.

In [20]:

colors=["red", "orange", "green", "blue", "darkorchid"]
linestyles=["-", "--", ":", "-.", "--"]
for i in range(len(x_disp)):
    for j in range(len(y_disp)):
        list_mbt[j][i].centerline(color=colors[i], linestyle=linestyles[j])
plt.title(r"Collection of plume centerlines with varying combination of $\alpha_x$, $\alpha_y$ and $\alpha_z$")
plt.show()
for i in range(len(x_disp)):
    for j in range(len(y_disp)):
        list_mbt[j][i].transverse(x_position=600, color=colors[i], linestyle=linestyles[j])
plt.ylim((0,4))
plt.title(r"Collection of transverse plots with varying combination of $\alpha_x$, $\alpha_y$ and $\alpha_z$")
plt.show()

colors=["red", "orange", "green", "blue", "darkorchid"] linestyles=["-", "--", ":", "-.", "--"] for i in range(len(x_disp)): for j in range(len(y_disp)): list_mbt[j][i].centerline(color=colors[i], linestyle=linestyles[j]) plt.title(r"Collection of plume centerlines with varying combination of $\alpha_x$, $\alpha_y$ and $\alpha_z$") plt.show() for i in range(len(x_disp)): for j in range(len(y_disp)): list_mbt[j][i].transverse(x_position=600, color=colors[i], linestyle=linestyles[j]) plt.ylim((0,4)) plt.title(r"Collection of transverse plots with varying combination of $\alpha_x$, $\alpha_y$ and $\alpha_z$") plt.show()

In [21]:

diff_array = np.zeros((len(y_disp), len(x_disp)))
for i in range(len(x_disp)):
    for j in range(len(y_disp)):
        diff_array[j,i] = diff.mean_absolute_difference(list_mbt[j][i].relative_cxyt, list_ana[j][i].relative_cxyt,
                                                      concentration_cutoff=1e-5)*100
print(diff_array)

y_ticks = []
for i in range(len(y_disp)):
    y_ticks.append(f"({y_disp[i]}, {z_disp[i]})")

fig, ax = plt.subplots()
im = ax.imshow(diff_array, cmap="viridis")
ax.set_xticks(range(len(x_disp)), labels=x_disp)
ax.set_yticks(range(len(y_ticks)), labels=y_ticks)

for i in range(len(x_disp)):
    for j in range(len(y_disp)):
        text = ax.text(i, j, np.round(diff_array[j, i], decimals=3),
                       ha="center", va="center", color="w")

ax.set_xlabel(r"Longitudinal dispersivity ($\alpha_x$) [m]")
ax.set_ylabel(r"Transverse dispersivity ($\alpha_y, \alpha_z$) [m]")
fig.colorbar(im, label="mean absolute difference in relative concentration [%]")
plt.savefig("mbt_ana_comp.png", dpi =300)

diff_array = np.zeros((len(y_disp), len(x_disp))) for i in range(len(x_disp)): for j in range(len(y_disp)): diff_array[j,i] = diff.mean_absolute_difference(list_mbt[j][i].relative_cxyt, list_ana[j][i].relative_cxyt, concentration_cutoff=1e-5)*100 print(diff_array) y_ticks = [] for i in range(len(y_disp)): y_ticks.append(f"({y_disp[i]}, {z_disp[i]})") fig, ax = plt.subplots() im = ax.imshow(diff_array, cmap="viridis") ax.set_xticks(range(len(x_disp)), labels=x_disp) ax.set_yticks(range(len(y_ticks)), labels=y_ticks) for i in range(len(x_disp)): for j in range(len(y_disp)): text = ax.text(i, j, np.round(diff_array[j, i], decimals=3), ha="center", va="center", color="w") ax.set_xlabel(r"Longitudinal dispersivity ($\alpha_x$) [m]") ax.set_ylabel(r"Transverse dispersivity ($\alpha_y, \alpha_z$) [m]") fig.colorbar(im, label="mean absolute difference in relative concentration [%]") plt.savefig("mbt_ana_comp.png", dpi =300)

[[0.27468921 0.17536177 0.09131589 0.05441245 0.03277492]
 [0.26769541 0.17151871 0.09029031 0.0546924  0.03385815]
 [0.25320083 0.16252541 0.08674549 0.05371712 0.03448358]
 [0.22879498 0.14748755 0.08048884 0.05144208 0.03464666]
 [0.18915283 0.12339181 0.0701624  0.04746385 0.03457523]]

In [22]:

diff_array = np.zeros((len(y_disp), len(x_disp)))
for i in range(len(x_disp)):
    for j in range(len(y_disp)):
        diff_array[j,i] = diff.mean_relative_difference(list_mbt[j][i].relative_cxyt, list_ana[j][i].relative_cxyt,
                                                      concentration_cutoff=5e-3)
print(diff_array)

y_ticks = []
for i in range(len(y_disp)):
    y_ticks.append(f"({y_disp[i]}, {z_disp[i]})")

fig, ax = plt.subplots(figsize=(6,6))
im = ax.imshow(diff_array, cmap="viridis")
ax.set_xticks(range(len(x_disp)), labels=x_disp)
ax.set_yticks(range(len(y_ticks)), labels=y_ticks)

for i in range(len(x_disp)):
    for j in range(len(y_disp)):
        text = ax.text(i, j, np.round(diff_array[j, i], decimals=3),
                       ha="center", va="center", color="w")

ax.set_xlabel(r"Longitudinal dispersivity ($\alpha_x$) [m]")
ax.set_ylabel(r"Transverse dispersivity ($\alpha_y, \alpha_z$) [m]")
fig.colorbar(im, label="mean relative difference in relative concentration")

diff_array = np.zeros((len(y_disp), len(x_disp))) for i in range(len(x_disp)): for j in range(len(y_disp)): diff_array[j,i] = diff.mean_relative_difference(list_mbt[j][i].relative_cxyt, list_ana[j][i].relative_cxyt, concentration_cutoff=5e-3) print(diff_array) y_ticks = [] for i in range(len(y_disp)): y_ticks.append(f"({y_disp[i]}, {z_disp[i]})") fig, ax = plt.subplots(figsize=(6,6)) im = ax.imshow(diff_array, cmap="viridis") ax.set_xticks(range(len(x_disp)), labels=x_disp) ax.set_yticks(range(len(y_ticks)), labels=y_ticks) for i in range(len(x_disp)): for j in range(len(y_disp)): text = ax.text(i, j, np.round(diff_array[j, i], decimals=3), ha="center", va="center", color="w") ax.set_xlabel(r"Longitudinal dispersivity ($\alpha_x$) [m]") ax.set_ylabel(r"Transverse dispersivity ($\alpha_y, \alpha_z$) [m]") fig.colorbar(im, label="mean relative difference in relative concentration")

[[0.10701172 0.0590436  0.02596625 0.01355647 0.00710582]
 [0.10308146 0.05644955 0.02469533 0.01288159 0.00673698]
 [0.09583099 0.05184132 0.02241242 0.01176208 0.00606193]
 [0.08279044 0.04464152 0.01930426 0.01012225 0.00535351]
 [0.0645975  0.03443319 0.01465069 0.00782025 0.00415009]]

Out[22]:

<matplotlib.colorbar.Colorbar at 0x7fd9d87de710>

In [23]:

x_disp_3 = np.array([10,3,1])
y_disp_3 = np.array([0.2,0.02,0.002])
z_disp_3 = np.array([0.05, 0.005, 0.0005, 0])

att_loop = mbt.AttenuationParameters(
    retardation=1,
    # Contaminant half life, in [days]
    half_life=5*365
)

source_loop = mbt.SourceParameters(
    source_zone_boundary=np.array([10]),
    source_zone_concentration=np.array([11]),
    depth=2.5,
    total_mass=1000000
)

model_loop = mbt.ModelParameters(
    # Model extent in the longitudinal (x) direction in [m].
    model_length = 1100,
    # Model extent in the transverse horizontal (y) direction in [m].
    model_width = 120,
    # Model duration in [days].
    model_time = 10*365,
    # Model grid discretization step size in the longitudinal (x) direction, in [m].
    dx = 2,
    # Model grid discretization step size in the transverse horizontal (y) direction, in [m].
    dy = 1,
    # Model time discretization step size, in [days]
    dt = 365
)

list_mbt_3 = [[[0]*len(x_disp_3) for _ in range(len(y_disp_3))] for _ in range(len(z_disp_3))]
list_ana_3 = [[[0]*len(x_disp_3) for _ in range(len(y_disp_3))] for _ in range(len(z_disp_3))]

mbt_object_disp = mbt.Mibitrans(hydro, att_loop, source_loop, model_loop)
ana_object_disp = mbt.Anatrans(hydro, att_loop, source_loop, model_loop)

# Loops for a total of 25 runs, will take about 1min.
for i in range(len(x_disp_3)):
    for j in range(len(y_disp_3)):
        for k in range(len(z_disp_3)):
            hydro_disp = mbt.HydrologicalParameters(
                velocity=0.2,             # [m/d]
                porosity=0.25,              # [-]
                alpha_x=x_disp_3[i],#alpha_x_list[i],                 # [m]
                alpha_y=y_disp_3[j],    # [m]
                alpha_z=z_disp_3[k],    # [m]
                diffusion=0,                # [m2/day]
            )
            print(i,j,k) # Print to track progress
            mbt_object_disp.hydrological_parameters = hydro_disp
            list_mbt_3[k][j][i] = mbt_object_disp.run()
            ana_object_disp.hydrological_parameters = hydro_disp
            list_ana_3[k][j][i] = ana_object_disp.run()

x_disp_3 = np.array([10,3,1]) y_disp_3 = np.array([0.2,0.02,0.002]) z_disp_3 = np.array([0.05, 0.005, 0.0005, 0]) att_loop = mbt.AttenuationParameters( retardation=1, # Contaminant half life, in [days] half_life=5*365 ) source_loop = mbt.SourceParameters( source_zone_boundary=np.array([10]), source_zone_concentration=np.array([11]), depth=2.5, total_mass=1000000 ) model_loop = mbt.ModelParameters( # Model extent in the longitudinal (x) direction in [m]. model_length = 1100, # Model extent in the transverse horizontal (y) direction in [m]. model_width = 120, # Model duration in [days]. model_time = 10*365, # Model grid discretization step size in the longitudinal (x) direction, in [m]. dx = 2, # Model grid discretization step size in the transverse horizontal (y) direction, in [m]. dy = 1, # Model time discretization step size, in [days] dt = 365 ) list_mbt_3 = [[[0]*len(x_disp_3) for _ in range(len(y_disp_3))] for _ in range(len(z_disp_3))] list_ana_3 = [[[0]*len(x_disp_3) for _ in range(len(y_disp_3))] for _ in range(len(z_disp_3))] mbt_object_disp = mbt.Mibitrans(hydro, att_loop, source_loop, model_loop) ana_object_disp = mbt.Anatrans(hydro, att_loop, source_loop, model_loop) # Loops for a total of 25 runs, will take about 1min. for i in range(len(x_disp_3)): for j in range(len(y_disp_3)): for k in range(len(z_disp_3)): hydro_disp = mbt.HydrologicalParameters( velocity=0.2, # [m/d] porosity=0.25, # [-] alpha_x=x_disp_3[i],#alpha_x_list[i], # [m] alpha_y=y_disp_3[j], # [m] alpha_z=z_disp_3[k], # [m] diffusion=0, # [m2/day] ) print(i,j,k) # Print to track progress mbt_object_disp.hydrological_parameters = hydro_disp list_mbt_3[k][j][i] = mbt_object_disp.run() ana_object_disp.hydrological_parameters = hydro_disp list_ana_3[k][j][i] = ana_object_disp.run()

0 0 0

0 0 1

0 0 2

0 0 3

0 1 0

0 1 1

0 1 2

0 1 3

0 2 0

0 2 1

0 2 2

0 2 3

1 0 0

1 0 1

1 0 2

1 0 3

1 1 0

1 1 1

1 1 2

1 1 3

1 2 0

1 2 1

1 2 2

1 2 3

2 0 0

2 0 1

2 0 2

2 0 3

2 1 0

2 1 1

2 1 2

2 1 3

2 2 0

2 2 1

2 2 2

2 2 3

In [24]:

list_mbt_3[0][0][0].centerline()
list_mbt_3[0][-1][0].centerline()
list_ana_3[0][0][0].centerline()
list_ana_3[0][-1][0].centerline()
plt.show()
xpos = 100
time= 2*365
list_mbt_3[0][0][1].transverse(x_position=xpos, time=time, label=r"mbt, $\alpha_y=0.2$")
list_mbt_3[0][-1][1].transverse(x_position=xpos, time=time, label=r"mbt, $\alpha_y=0.002$")
list_ana_3[0][0][1].transverse(x_position=xpos, time=time, label=r"ana, $\alpha_y=0.2$")
list_ana_3[0][-1][1].transverse(x_position=xpos, time=time, label=r"ana, $\alpha_y=0.002$")
plt.legend()
plt.ylim(0,6)

list_mbt_3[0][0][0].centerline() list_mbt_3[0][-1][0].centerline() list_ana_3[0][0][0].centerline() list_ana_3[0][-1][0].centerline() plt.show() xpos = 100 time= 2*365 list_mbt_3[0][0][1].transverse(x_position=xpos, time=time, label=r"mbt, $\alpha_y=0.2$") list_mbt_3[0][-1][1].transverse(x_position=xpos, time=time, label=r"mbt, $\alpha_y=0.002$") list_ana_3[0][0][1].transverse(x_position=xpos, time=time, label=r"ana, $\alpha_y=0.2$") list_ana_3[0][-1][1].transverse(x_position=xpos, time=time, label=r"ana, $\alpha_y=0.002$") plt.legend() plt.ylim(0,6)

Out[24]:

(0.0, 6.0)

In [25]:

##### Plotting preferences #####
cutoff = 0
as_percentage = True
value_decimals = 4
figuresize = (7,5)
colormap = "viridis"
second_y_axis_location = -0.24
x_label = r"Longitudinal dispersivity ($\alpha_x$) [m]"
first_y_label = r"Transverse vertical dispersivity $\alpha_z$ [m]"
second_y_label = r"Transverse horizontal dispersivity $\alpha_y$ [m]"
colorbar_label = "mean absolute difference in relative concentration [%]"
#################################

diff_array_3 = np.zeros((len(y_disp_3) * len(z_disp_3), len(x_disp_3)))
for i in range(len(x_disp_3)):
    for j in range(len(y_disp_3)):
        for k in range(len(z_disp_3)):
            diff_array_3[len(z_disp_3)*j+k,i] = diff.mean_absolute_difference(
                list_mbt_3[k][j][i].relative_cxyt, list_ana_3[k][j][i].relative_cxyt,
                concentration_cutoff=cutoff
            )

if as_percentage:
    diff_array_3 = diff_array_3 * 100

z_ticks_3 = []
for i in range(len(y_disp_3)):
    for j in range(len(z_disp_3)):
        z_ticks_3.append(f"{z_disp_3[j]}")

fig, ax = plt.subplots(figsize=figuresize)
imshow_aspect = str(1/len(z_disp_3))
im = ax.imshow(
    diff_array_3,
    cmap=colormap,
    aspect=imshow_aspect
)
v2 = ax.secondary_yaxis(second_y_axis_location)
middle_z_disp = len(z_disp_3)/2-0.5

second_y_tick_location = np.linspace(middle_z_disp, (len(y_disp_3)-1)*len(z_disp_3) + middle_z_disp, len(y_disp_3))

ax.set_xticks(range(len(x_disp_3)), labels=x_disp_3)
ax.set_yticks(range(len(z_ticks_3)), labels=z_ticks_3)
v2.set_yticks(second_y_tick_location, labels=y_disp_3)

ax.hlines(second_y_tick_location[:-1]+len(z_disp_3)/2, -0.5, len(x_disp_3)-0.5, color="white", lw=2)
ax.vlines(
    np.arange(0.5,len(x_disp_3)-1,1), -0.5, len(y_disp_3)*len(z_disp_3)-0.5,
    color="white",
    lw=2
)

for i in range(len(x_disp_3)):
    for j in range(len(y_disp_3)*len(z_disp_3)):
        text3 = ax.text(
            i, j, np.round(diff_array_3[j, i], decimals=value_decimals),
            ha="center",
            va="center",
            color="w"
        )

ax.set_xlabel(x_label)
ax.set_ylabel(first_y_label)
v2.set_ylabel(second_y_label)
fig.colorbar(im, label=colorbar_label)
plt.savefig("mbt_ana_comp_dual_axis.png", dpi=500, bbox_inches="tight")
plt.show()

##### Plotting preferences ##### cutoff = 0 as_percentage = True value_decimals = 4 figuresize = (7,5) colormap = "viridis" second_y_axis_location = -0.24 x_label = r"Longitudinal dispersivity ($\alpha_x$) [m]" first_y_label = r"Transverse vertical dispersivity $\alpha_z$ [m]" second_y_label = r"Transverse horizontal dispersivity $\alpha_y$ [m]" colorbar_label = "mean absolute difference in relative concentration [%]" ################################# diff_array_3 = np.zeros((len(y_disp_3) * len(z_disp_3), len(x_disp_3))) for i in range(len(x_disp_3)): for j in range(len(y_disp_3)): for k in range(len(z_disp_3)): diff_array_3[len(z_disp_3)*j+k,i] = diff.mean_absolute_difference( list_mbt_3[k][j][i].relative_cxyt, list_ana_3[k][j][i].relative_cxyt, concentration_cutoff=cutoff ) if as_percentage: diff_array_3 = diff_array_3 * 100 z_ticks_3 = [] for i in range(len(y_disp_3)): for j in range(len(z_disp_3)): z_ticks_3.append(f"{z_disp_3[j]}") fig, ax = plt.subplots(figsize=figuresize) imshow_aspect = str(1/len(z_disp_3)) im = ax.imshow( diff_array_3, cmap=colormap, aspect=imshow_aspect ) v2 = ax.secondary_yaxis(second_y_axis_location) middle_z_disp = len(z_disp_3)/2-0.5 second_y_tick_location = np.linspace(middle_z_disp, (len(y_disp_3)-1)*len(z_disp_3) + middle_z_disp, len(y_disp_3)) ax.set_xticks(range(len(x_disp_3)), labels=x_disp_3) ax.set_yticks(range(len(z_ticks_3)), labels=z_ticks_3) v2.set_yticks(second_y_tick_location, labels=y_disp_3) ax.hlines(second_y_tick_location[:-1]+len(z_disp_3)/2, -0.5, len(x_disp_3)-0.5, color="white", lw=2) ax.vlines( np.arange(0.5,len(x_disp_3)-1,1), -0.5, len(y_disp_3)*len(z_disp_3)-0.5, color="white", lw=2 ) for i in range(len(x_disp_3)): for j in range(len(y_disp_3)*len(z_disp_3)): text3 = ax.text( i, j, np.round(diff_array_3[j, i], decimals=value_decimals), ha="center", va="center", color="w" ) ax.set_xlabel(x_label) ax.set_ylabel(first_y_label) v2.set_ylabel(second_y_label) fig.colorbar(im, label=colorbar_label) plt.savefig("mbt_ana_comp_dual_axis.png", dpi=500, bbox_inches="tight") plt.show()

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.

Guyonnet, D., & Neville, C., Dimensionless analysis of two analytical solutions for 3-D solute transport in groundwater. Journal of Contaminant Hydrology, 75(1), 141–153, https://doi.org/10.1016/j.jconhyd.2004.06.004, 2004.