Bemidji oil pipeline burst
import matplotlib.pyplot as plt
import numpy as np
import mibitrans as mbt
Field site example: Bemidji oil pipeline burst¶
This example explores a well studied contaminated field site near Bemidji (Minnesota, USA), and aims to find the plume length under steady state conditions.
Site description¶
Near the town of Bemidji, a crude-oil pipeline burst in 1979, releasing approximately $1700 m^3$ of crude oil. After initial remediation efforts, about $400 m^3$ remained. The subsurface consists of glacial outwash deposits, comprised of poorly sorted sandy gravel, gravely sand and sand interbedded with thin silt layers. This aquifer is underlain by a low permeable till layer [1]. A bromide tracer test performed in 1997 between two wells ($1.6m$ apart) calculated a flow velocity of $0.06m/day$ and an estimated longitudinal dispersivity of $0.15m$ [2]. The flow velocity is close to the range ($0.004$ to $0.056 m/day$), based on observed average porosity ($0.38$), hydraulic gradients and hydraulic conductivities. As the tracer test was performed over a short distance, any macro-dispersivity is not included in the estimate of $\alpha_l$. Therefore, the actual longitudinal dispersivity is expected to be higher. No information is given on transverse horizontal or transverse vertical dispersivities. Therefore, it is assumed that $\alpha_{th} = \alpha_l / 10$ and $\alpha_{tv} = \alpha_l / 100$. Although these estimates are crude, they correspond with reasonable orders of magnitude for the respective parameters [3].
hydro = mbt.HydrologicalParameters(
velocity=0.06, #[m/d]
porosity=0.38, #[-]
alpha_x=0.15, #[m]
alpha_y=0.015, #[m]
alpha_z=0.0015 #[m]
)
For biodegradation, the instant reaction approach is used. Therefore, the decay rate is not relevant in this scenario. As the example looks at steady state, contaminant retardation is not relevant either, since this only influences the time until steady state, not the steady state itself.
Here, only one electron acceptor is considered; oxygen, with a concentration of $8 g/m^3$ and the default utilization factor of $3.14 g/g$ [4].
att = mbt.AttenuationParameters(
decay_rate=0, #[1/day]
retardation=1 #[-]
)
ea = mbt.ElectronAcceptors(
delta_oxygen=8, #[g/m3]
delta_nitrate=0,
ferrous_iron=0,
delta_sulfate=0,
methane=0
)
No source depletion is considered, since that is incompatible with steady state conditions. Initial source concentrations are reported to be in the range of $4-7 g/m^3$ [5]. Here, the (approximately) same concentration ($6 g/m^3$) and source dimensions ($1m$ radius) are used as a recent study modelling the same site [4], which in turn are based on field site literature [1].
source = mbt.SourceParameters(
source_zone_boundary=np.array([1]), #[m]
source_zone_concentration=np.array([6]), #[g/m3]
total_mass=np.inf, #[g]
depth=1, #[m]
)
model = mbt.ModelParameters(
model_length=200, #[m]
model_width=10, #[m]
model_time=10*365, #[days]
dx=1, #[m]
dy=0.1, #[m]
dt=365/5, #[days]
)
Modelling¶
Use Mibitrans solution to calculate the concentration distribution.
mbt_obj = mbt.Mibitrans(hydro, att, source, model)
mbt_obj.instant_reaction(ea)
results_mbt = mbt_obj.run()
colors=["lightgreen", "greenyellow", "khaki", "goldenrod", "orangered"]
plot_times=[10*365, 7*365, 5*365, 3*365, 1*365]
lw = [2, 2.25, 2.5, 2.75, 3]
linestyle=["-", "--", "--", ":", ":"]
for i in range(len(plot_times)):
results_mbt.centerline(time=plot_times[i], color=colors[i], linestyle=linestyle[i],
lw=lw[i], label=f"t={plot_times[i]}days")
plt.legend()
plt.show()
The solid dark green line is the modelled concentration distribution at steady state, with a plume length somewhere around $160m$. To get a more exact plume length:
# If multiple minimum values, argmin gives the first index where array is minimum value
plume_length_index = np.argmin(results_mbt.cxyt[-1,len(results_mbt.y)//2,:])
print("Plume length is", results_mbt.x[plume_length_index], "m")
Plume length is 162 m
Thus, under these parameters, the modelled plume length is $162m$, slightly more than the on site measured $150m$ extent [4]. As a relatively simple model, the calculated plume length is quite accurate. Of course, given the parameter estimations and assumptions, the model is by no means conclusive.
[1] Essaid, H. I., Bekins, B. A., Herkelrath, W. N., & Delin, G. N. (2011). Crude Oil at the Bemidji Site: 25 Years of Monitoring, Modeling, and Understanding. Groundwater, 49(5), 706–726. https://doi.org/10.1111/j.1745-6584.2009.00654.x
[2] Essaid, H. I., Cozzarelli, I. M., Eganhouse, R. P., Herkelrath, W. N., Bekins, B. A., & Delin, G. N. (2003). Inverse modeling of BTEX dissolution and biodegradation at the Bemidji, MN crude-oil spill site. Journal of Contaminant Hydrology, 67(1), 269–299. https://doi.org/10.1016/S0169-7722(03)00034-2
[3] Zech, A., Attinger, S., Bellin, A., Cvetkovic, V., Dietrich, P., Fiori, A., Teutsch, G. & Da-gan, G. (2019). A critical analysis of transverse dispersivity field data. Groundwater 57 (4), 632–639. http://dx.doi.org/10.1111/gwat.12838
[4] Köhler, A. V., Craig, J. R., Yadav, P. K., & Liedl, R. (2026). An Analytic Element Method solution for simulating multiple steady-state groundwater contamination scenarios. Journal of Contaminant Hydrology, 276, 104733. https://doi.org/10.1016/j.jconhyd.2025.104733
[5] Bekins, B. A., Cozzarelli, I. M., Erickson, M. L., Steenson, R. A., & Thorn, K. A. (2016). Crude Oil Metabolites in Groundwater at Two Spill Sites. Groundwater, 54(5), 681–691. https://doi.org/10.1111/gwat.12419