Field site example: Bemidji oil pipeline burst¶

This notebooks outlines analysis of a contaminated field site near Bemidji (Minnesota, USA). It focusses on the analysis of the plume length under steady state conditions.

Authors: Jorrit Bakker, Alraune Zech

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$ of crude-oil remained in the soil.

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$^2$ between two wells ($1.6m$ apart) calculated a flow velocity of $0.06m/day$ and an estimated longitudinal dispersivity of $\alpha_l = 0.15m$. As the tracer test was performed over a short distance, this values is not necessarily representing macro-dispersivity, but the actual longitudinal macrodispersivity is expected to be higher. No information is given on transverse horizontal or transverse vertical dispersivities. Observed average porosity is $0.38$. The velocity determined by the tracer test is close to the range of $0.004 - 0.056 m/day$ which reflects calculated velocity based on measured hydraulic gradients and hydraulic conductivities.

Model setup¶

Goal is to setup a transient transport model for the contaminant transport at the Bemidji site using field site characteristics and reproduce the steady state plume distribution/length as done by Köhler et al., 2026$^4$.

Hydrogeological parameters¶

We chose hydraulic parameters following the site description. In this analysis, we 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$.

Attenuation and instant reaction model parameters¶

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 is reached, not the concentration distribution at steady state.

Following Köhler et al., 2026^$4$, we use only one electron acceptor, oxygen, with a concentration of $8 g/m^3$ and the default utilization factor of $3.14 g/g$

Source conditions¶

No source depletion is considered, since that is incompatible with steady state conditions. We use a source concentration of $6 g/m^3$ and a source dimensions of $1m$ following $^{1,4}$ Note that Bekins et al., 2016$^5$ reports initial source concentrations in the range of $4-7 g/m^3$ which is close to the values we use.

Model parameters¶

We model the plume within a distance of 200m and a total of $10$ years.

References¶

$^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.

$^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.

$^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.

$^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.

$^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.

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

Setup Data Input¶

In [2]:

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]
)

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] )

In [3]:

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
)

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 )

In [4]:

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]
)

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] )

In [5]:

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]
)

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 with instant reaction for (aerobic) biodegredation to calculate the concentration distribution.

In [6]:

mbt_bemidji = mbt.Mibitrans(hydro, att, source, model)
mbt_bemidji.instant_reaction(ea)
results_bemidji = mbt_bemidji.run()

mbt_bemidji = mbt.Mibitrans(hydro, att, source, model) mbt_bemidji.instant_reaction(ea) results_bemidji = mbt_bemidji.run()

Plotting¶

Concentration along plume center over time¶

In [7]:

cmap = plt.get_cmap('Dark2')
colors = cmap.colors   # list of RGB tuples

cmap = plt.get_cmap('Dark2') colors = cmap.colors # list of RGB tuples

In [8]:

for i,plot_time in enumerate([10, 7, 5, 3, 1]):
    results_bemidji.centerline(time=plot_time*365,
                               color=colors[i+1],lw = 2,
                               label=f"t={plot_time}y")
plt.legend()
plt.show()

for i,plot_time in enumerate([10, 7, 5, 3, 1]): results_bemidji.centerline(time=plot_time*365, color=colors[i+1],lw = 2, label=f"t={plot_time}y") plt.legend() plt.show()

Interpretation: The solid dark green line is the modelled concentration distribution at steady state, with a plume length somewhere around $160m$.

Calculate the exact plume length via:

In [9]:

# If multiple minimum values, argmin gives the first index where array is minimum value
plume_length_index = np.argmin(results_bemidji.cxyt[-1,len(results_bemidji.y)//2,:])
print("Plume length is", results_bemidji.x[plume_length_index], "m")

# If multiple minimum values, argmin gives the first index where array is minimum value plume_length_index = np.argmin(results_bemidji.cxyt[-1,len(results_bemidji.y)//2,:]) print("Plume length is", results_bemidji.x[plume_length_index], "m")

Plume length is 162 m

Steady state plume¶

In [10]:

results_bemidji.plume_2d(cmap="RdYlGn_r", shading="gouraud")
plt.ylim([-2,2])
plt.xlim([0,160])

results_bemidji.plume_2d(cmap="RdYlGn_r", shading="gouraud") plt.ylim([-2,2]) plt.xlim([0,160])

Out[10]:

(0.0, 160.0)

Interpretation: The modelled plume length is $162m$. This is only slightly more than the measured extend of $150m$, also modelled by Köhler et al, 2026$^4$.

Disclaimer: Given the parameter estimations, and model assumptions, the modelling results are by no means conclusive.