Note

This notebook is located in the ./examples directory of the gwtransport repository.

Macrodispersion and Microdispersion in Solute Transport#

This notebook explains the distinction between macrodispersion (spreading from aquifer-scale heterogeneity, captured by the APVD), microdispersion (spreading from pore-scale velocity variations, characterized by \(\alpha_L\)), and molecular diffusion (\(D_m\)), and provides formulas to convert between them. See dispersion scales and advection-dominated assumption for background.

Dispersion as Scale-Dependent Heterogeneity#

All spreading in groundwater transport arises from velocity heterogeneity. What we call “dispersion” vs “heterogeneity” is a matter of observation scale:

Scale	Mechanism	Dispersion type	gwtransport representation
Molecular	Brownian motion	Molecular diffusion	\(D_m\) [m²/day]
Pore (~mm-cm)	Velocity variations in pores	Microdispersion	\(\alpha_L\) [m] (longitudinal dispersivity)
Aquifer (~m-km)	Different streamline paths	Macrodispersion	\(\sigma_V\) [m³] (std of pore volume distribution)

Key properties:

Microdispersion (\(\alpha_L\)) is an aquifer property — independent of hydrological boundary conditions
Macrodispersion depends on both aquifer properties and boundary conditions
Molecular diffusion depends on the compound and temperature

Standard approach: The diffusion_fast module is the recommended transport module. It combines advection (via APVD) with Gaussian smoothing for microdispersion and molecular diffusion. The diffusion module provides exact solutions but is much slower. The advection module ignores microdispersion and molecular diffusion entirely.

Prerequisite: The APVD is only well-defined under the steady-streamlines assumption.

Spreading Components#

Breakthrough curve spreading arises from three processes, expressed as equivalent standard deviation in volume units [m³]:

Component	Symbol	Physical cause	Formula [m³]
Aquifer heterogeneity	\(\sigma_{V,apv}\)	Different streamline lengths	\(\sigma_{V,apv}\) (given)
Molecular diffusion	\(\sigma_{V,diff}\)	Brownian motion	\(\frac{\bar{V}}{L} \sqrt{2 D_m \tau}\)
Mechanical dispersion	\(\sigma_{V,disp}\)	Pore-scale velocity variations	\(\bar{V} \sqrt{\frac{2 \alpha_L}{L}}\)

where \(\bar{V}\) is mean pore volume, \(L\) is streamline length, \(\tau = R \bar{V} / Q\) is mean solute residence time.

Important: \(\sigma_{V,disp}\) is constant regardless of flow velocity. \(\sigma_{V,diff}\) decreases with higher \(Q\) (less time for diffusion).

Unit conversions#

Unit	Symbol	Conversion from spatial \(\sigma_x\)
Spatial [m]	\(\sigma_x\)	—
Volume [m³]	\(\sigma_V\)	\(\sigma_V = \sigma_x \cdot \bar{V}/L\)
Time [day]	\(\sigma_t\)	\(\sigma_t = \sigma_x / v = \sigma_x \cdot \tau / L\)

Choosing the Right Approach#

By APVD source#

APVD source	Recommendation
Calibration from measurements	`diffusion_fast.gamma_` (standard). Or `advection.gamma_` if diffusion is negligible — the fitted std then absorbs all mixing.
Gamma from streamline analysis	`diffusion_fast.gamma_` (standard). Or `advection.gamma_` with combined std (variance addition, see below).
Discrete pore volumes	`diffusion_fast.infiltration_to_extraction` (standard) or `diffusion.infiltration_to_extraction` (exact). Variance addition cannot be used with discrete volumes.

Cross-compound calibration: When the APVD was calibrated using one compound (e.g., temperature with \(D_m \approx 0.1\) m²/day) and is used to predict another (e.g., a dissolved solute with \(D_m \approx 10^{-4}\) m²/day), the molecular diffusion contribution in the calibrated std is too large for the solute. When calibrating with diffusion_fast or diffusion, the three components are separated automatically.

Modeling approaches#

Method	Function	Dispersion included	Speed	When to use
Gamma + fast diffusion	`diffusion_fast.gamma_*`	Macro + micro + molecular	Fast	Standard approach — recommended for most cases
Discrete + fast diffusion	`diffusion_fast.infiltration_to_extraction`	Macro + micro + molecular	Fast	Discrete APVD from streamline analysis
Gamma APVD only	`advection.gamma_*`	Macro only (or combined via variance addition)	Fast	When diffusion is negligible, or combined std from variance addition
Discrete + exact diffusion	`diffusion.infiltration_to_extraction`	Macro + micro + molecular	Slow	Validation, highly variable flow, or when high precision is needed

Variance addition (alternative fast approach)#

When using the advection module, microdispersion and molecular diffusion can be included by increasing the gamma std parameter: \(\sigma_{total} = \sqrt{\sigma_{apv}^2 + \sigma_{diff}^2 + \sigma_{disp}^2}\). This is only valid for gamma-parameterized (continuous) distributions from streamline analysis, not from calibration (where all mixing is already included) and not for discrete APVD.

[1]:

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

from gwtransport import advection, diffusion, diffusion_fast

np.random.seed(42)
plt.style.use("seaborn-v0_8-whitegrid")

1. System Parameters#

The input parameters are organized in three groups:

Aquifer geometry: streamline length \(L\) and aquifer pore volume distribution (mean \(\bar{V}\), std \(\sigma_{apv}\))
Flow: mean discharge \(Q\)
Transport: retardation factor \(R\), molecular diffusivity \(D_m\), and longitudinal dispersivity \(\alpha_L\)

[2]:

# Aquifer geometry
streamline_length = 100.0  # L [m]
mean_apv = 10000.0  # V_mean [m³] Mean aquifer pore volume
std_apv = 800.0  # sigma_apv [m³] Std from aquifer heterogeneity

# Flow
mean_flow = 120.0  # Q [m³/day]

# Transport
retardation = 2.0  # R [-] Retardation factor
diffusivity_molecular = 1e-4  # D_m [m²/day] Molecular diffusion (~1e-9 m²/s)
dispersivity = 1.0  # alpha_L [m] Longitudinal dispersivity

2. Converting Microdispersion and Molecular Diffusion to Equivalent APVD Spreading#

To compare microdispersion and molecular diffusion with macrodispersion (APVD heterogeneity), we convert spatial spreading to equivalent standard deviation in volume units [m³]. This allows direct variance addition.

Conversion from spatial to volume units#

Spatial spreading is \(\sigma_x\) [m]. To convert to volume units [m³]:

\[\sigma_V = \sigma_x \cdot \frac{\bar{V}}{L} \quad \text{[m³]}\]

This uses the relationship: pore volume \(V\) = cross-sectional area × length, so \(\bar{V}/L\) [m²] converts length to volume.

Formulas for equivalent spreading in volume units#

After substitution (see derivation below), the formulas simplify to expressions using only input parameters:

\[\sigma_{V,diff} = \frac{\bar{V}}{L} \sqrt{2 D_m \tau} = \frac{\bar{V}}{L} \sqrt{\frac{2 D_m R \bar{V}}{Q}} \quad \text{[m³]}\]

\[\sigma_{V,disp} = \frac{\bar{V}}{L} \sqrt{2 \alpha_L L} = \bar{V} \sqrt{\frac{2 \alpha_L}{L}} \quad \text{[m³]}\]

\[\sigma_{V,diff+disp} = \sqrt{\sigma_{V,diff}^2 + \sigma_{V,disp}^2} = \bar{V} \sqrt{\frac{2 D_m R}{L Q} + \frac{2 \alpha_L}{L}} \quad \text{[m³]}\]

Flow-dependence:

\(\sigma_{V,disp}\) (microdispersion) is independent of flow \(Q\): the plume traverses the same pore-scale heterogeneity regardless of velocity — this is a fixed aquifer property
\(\sigma_{V,diff}\) (molecular diffusion) decreases with higher \(Q\) (proportional to \(1/\sqrt{Q}\)): faster flow means less time for molecular diffusion

Derivation

Starting from \(\sigma_x = \sqrt{2 D \tau}\) [m] and converting to volume units via \(\sigma_V = \sigma_x \cdot \bar{V} / L\):

Molecular diffusion (\(D = D_m\), constant):

\[\sigma_{x,diff} = \sqrt{2 D_m \tau} = \sqrt{2 D_m R \bar{V} / Q} \quad \text{[m]}\]

\[\sigma_{V,diff} = \frac{\bar{V}}{L} \sqrt{\frac{2 D_m R \bar{V}}{Q}} \quad \text{[m³]}\]

Microdispersion (\(D_{disp} = \alpha_L v = \alpha_L L / \tau\)):

\[\sigma_{x,disp}^2 = 2 D_{disp} \tau = 2 \alpha_L \frac{L}{\tau} \tau = 2 \alpha_L L \quad \text{[m²]}\]

\[\sigma_{x,disp} = \sqrt{2 \alpha_L L} \quad \text{[m]}\]

\[\sigma_{V,disp} = \frac{\bar{V}}{L} \sqrt{2 \alpha_L L} = \bar{V} \sqrt{\frac{2 \alpha_L}{L}} \quad \text{[m³]}\]

Note: The retardation factor \(R\) cancels out for microdispersion because \(D_{disp} \propto v \propto 1/\tau\) and the product \(D_{disp} \cdot \tau\) is independent of \(\tau\). Physically: a retarded compound moves slower through the same pore structure, experiencing the same heterogeneity over the same distance.

For molecular diffusion, \(R\) does not cancel: a retarded compound spends more time in the aquifer, giving more time for diffusion to act.

Combined (variances add since processes are independent):

\[\sigma_{V,diff+disp}^2 = \sigma_{V,diff}^2 + \sigma_{V,disp}^2\]

[3]:

# Equivalent spreading in volume units from molecular diffusion
sigma_v_diff = (mean_apv / streamline_length) * np.sqrt(2 * diffusivity_molecular * retardation * mean_apv / mean_flow)

# Equivalent spreading in volume units from mechanical dispersion
# Note: R cancels out because D_disp = alpha_L * v = alpha_L * L / tau, and D_disp * tau = alpha_L * L
sigma_v_disp = mean_apv * np.sqrt(2 * dispersivity / streamline_length)

# Combined spreading (direct formula)
sigma_v_diff_disp = mean_apv * np.sqrt(
    2 * diffusivity_molecular * retardation / (streamline_length * mean_flow) + 2 * dispersivity / streamline_length
)

print("Equivalent spreading contributions in volume units [m³]:")
print(f"  sigma_v_diff (molecular diffusion):   {sigma_v_diff:6.0f} m³")
print(f"  sigma_v_disp (mechanical dispersion): {sigma_v_disp:6.0f} m³")
print(f"  sigma_v_diff_disp (combined):         {sigma_v_diff_disp:6.0f} m³")

Equivalent spreading contributions in volume units [m³]:
  sigma_v_diff (molecular diffusion):       18 m³
  sigma_v_disp (mechanical dispersion):   1414 m³
  sigma_v_diff_disp (combined):           1414 m³

3. Combining Macrodispersion and Microdispersion#

Section 2 computed \(\sigma_{diff,disp}\) [m³]—the spreading from microdispersion and molecular diffusion only.

Now we combine this with the original APVD heterogeneity (macrodispersion) \(\sigma_{apv}\) [m³] to get the total effective standard deviation:

\[\sigma_{total} = \sqrt{\sigma_{apv}^2 + \sigma_{diff}^2 + \sigma_{disp}^2} = \sqrt{\sigma_{apv}^2 + \sigma_{diff,disp}^2} \quad \text{[m³]}\]

This is the std parameter to use in advection.gamma_infiltration_to_extraction when you want to include microdispersion and molecular diffusion effects in the fast advection module.

When is this valid?

Only when APVD was obtained from streamline analysis (explicit geometry), not from calibration on measurements (where all mixing is already included in the calibrated std). See dispersion scales.
Only for gamma-parameterized (continuous) distributions where the std parameter can be adjusted. For discrete streamline volumes from a flow model, variance addition is not meaningful—use the diffusion module instead.

[4]:

# Total std in volume units combining all sources
sigma_v_total = np.sqrt(std_apv**2 + sigma_v_diff_disp**2)

# Variance fractions
var_total = sigma_v_total**2
frac_apv = std_apv**2 / var_total * 100
frac_diff = sigma_v_diff**2 / var_total * 100
frac_disp = sigma_v_disp**2 / var_total * 100

print("Variance contributions to total spreading:")
print(f"  sigma_v_apv  = {std_apv:7.0f} m³  ({frac_apv:5.1f}% of variance) - aquifer heterogeneity")
print(f"  sigma_v_diff = {sigma_v_diff:7.0f} m³  ({frac_diff:5.1f}% of variance) - molecular diffusion")
print(f"  sigma_v_disp = {sigma_v_disp:7.0f} m³  ({frac_disp:5.1f}% of variance) - mechanical dispersion")
print("  " + "-" * 60)
print(f"  sigma_v_total= {sigma_v_total:7.0f} m³  (100.0% of variance)")

print("\n" + "=" * 62)
if frac_diff + frac_disp < 5:
    print("RECOMMENDATION: Use APVD only (sigma_v_apv)")
    print("  -> advection.gamma_* with std=sigma_v_apv")
elif frac_diff < 1:
    print("RECOMMENDATION: Use APVD + dispersion")
    print("  -> advection.gamma_* with std=sqrt(sigma_v_apv^2 + sigma_v_disp^2)")
elif frac_diff + frac_disp < 25:
    print("RECOMMENDATION: Use APVD + diff + disp (constant flow assumed)")
    print("  -> advection.gamma_* with std=sigma_v_total")
else:
    print("RECOMMENDATION: Use full diffusion module for accuracy")
    print("  -> diffusion.infiltration_to_extraction")
print("=" * 62)

Variance contributions to total spreading:
  sigma_v_apv  =     800 m³  ( 24.2% of variance) - aquifer heterogeneity
  sigma_v_diff =      18 m³  (  0.0% of variance) - molecular diffusion
  sigma_v_disp =    1414 m³  ( 75.8% of variance) - mechanical dispersion
  ------------------------------------------------------------
  sigma_v_total=    1625 m³  (100.0% of variance)

==============================================================
RECOMMENDATION: Use APVD + dispersion
  -> advection.gamma_* with std=sqrt(sigma_v_apv^2 + sigma_v_disp^2)
==============================================================

Including Microdispersion via Variance Addition (Alternative Fast Approach)#

When using the advection module (instead of the standard diffusion_fast), microdispersion and molecular diffusion can be included by adjusting the gamma std parameter. This works because macrodispersion and microdispersion are independent processes whose variances can be added. See dispersion scales.

# Combined std for advection module:
std_total = np.sqrt(sigma_apv**2 + sigma_diff**2 + sigma_disp**2)
cout = advection.gamma_infiltration_to_extraction(..., std=std_total, ...)

Validity conditions:

APVD from streamline analysis (not calibration — where all mixing is already included)
Gamma-parameterized distribution (not discrete volumes — use diffusion_fast or diffusion instead)
Flow relatively constant (for \(\sigma_{diff}\) calculation)

4. Numerical Validation#

[5]:

n_days = 350
_tedges = pd.date_range("2019-12-31", periods=n_days + 1, freq="D")

_flow = np.full(n_days, mean_flow)
_cin = np.zeros(n_days)
_cin[50] = 100.0

cin_itedges = np.flatnonzero(np.diff(_cin, prepend=1.0, append=1.0))
flow_itedges = np.flatnonzero(np.diff(_flow, prepend=1.0, append=1.0))
itedges = np.unique(np.concatenate([cin_itedges, flow_itedges]))
tedges = _tedges[itedges]

cin = _cin[itedges[:-1]]
flow = _flow[itedges[:-1]]

cout_tedges = _tedges.copy()

nbins = 5000
streamline_lengths = np.full(nbins, streamline_length)

[6]:

# 1. APVD only (no diffusion/dispersion) - baseline showing aquifer heterogeneity
cout_apvd = advection.gamma_infiltration_to_extraction(
    cin=cin,
    flow=flow,
    tedges=tedges,
    cout_tedges=cout_tedges,
    mean=mean_apv,
    std=std_apv,
    n_bins=nbins,
    retardation_factor=retardation,
)

# 2. Full diffusion solution with a SINGLE pore volume (no APVD heterogeneity)
# This is the reference for validating the equivalent std approximation
cout_diffusion_single = diffusion.infiltration_to_extraction(
    cin=cin,
    flow=flow,
    tedges=tedges,
    cout_tedges=cout_tedges,
    aquifer_pore_volumes=np.array([mean_apv]),
    streamline_length=np.array([streamline_length]),
    molecular_diffusivity=diffusivity_molecular,
    longitudinal_dispersivity=dispersivity,
    retardation_factor=retardation,
    asymptotic_cutoff_sigma=4.0,
)

# 3. APVD approximation using sigma_v_diff_disp (diffusion + dispersion spreading only)
# This should match cout_diffusion_single since both have no APVD heterogeneity
cout_apvd_approx = advection.gamma_infiltration_to_extraction(
    cin=cin,
    flow=flow,
    tedges=tedges,
    cout_tedges=cout_tedges,
    mean=mean_apv,
    std=sigma_v_diff_disp,  # Equivalent std from diffusion + dispersion
    n_bins=nbins,
    retardation_factor=retardation,
)

# 4. Full diffusion solution with MULTIPLE pore volumes (full APVD + diffusion/dispersion)
# This is the slow but physically rigorous approach
from gwtransport import gamma

gamma_bins = gamma.bins(mean=mean_apv, std=std_apv, n_bins=nbins)
pore_volumes = gamma_bins["expected_values"]
cout_diffusion_multi = diffusion.infiltration_to_extraction(
    cin=cin,
    flow=flow,
    tedges=tedges,
    cout_tedges=cout_tedges,
    aquifer_pore_volumes=pore_volumes,
    streamline_length=streamline_lengths,
    molecular_diffusivity=diffusivity_molecular,
    longitudinal_dispersivity=dispersivity,
    retardation_factor=retardation,
    asymptotic_cutoff_sigma=4.0,
    suppress_dispersion_warning=True,  # We know what we're doing here
)

# 5. APVD approximation with combined std (APVD + diffusion + dispersion)
# This is the fast approximation - should match cout_diffusion_multi
cout_apvd_total = advection.gamma_infiltration_to_extraction(
    cin=cin,
    flow=flow,
    tedges=tedges,
    cout_tedges=cout_tedges,
    mean=mean_apv,
    std=sigma_v_total,  # Combined std: sqrt(sigma_apv^2 + sigma_diff^2 + sigma_disp^2)
    n_bins=nbins,
    retardation_factor=retardation,
)

# 6. diffusion_fast with MULTIPLE pore volumes (standard approach)
# flow_out is needed because tedges has non-uniform bins (from signal compression above)
flow_out = np.full(len(cout_tedges) - 1, mean_flow)
cout_diffusion_fast_multi = diffusion_fast.gamma_infiltration_to_extraction(
    cin=cin,
    flow=flow,
    tedges=tedges,
    cout_tedges=cout_tedges,
    mean=mean_apv,
    std=std_apv,
    n_bins=nbins,
    mean_streamline_length=streamline_length,
    mean_molecular_diffusivity=diffusivity_molecular,
    mean_longitudinal_dispersivity=dispersivity,
    retardation_factor=retardation,
    suppress_dispersion_warning=True,
    flow_out=flow_out,
)

[7]:

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# Plot around breakthrough
plot_start, plot_end = 100, 349
t_plot = cout_tedges[plot_start:plot_end]

# Left panel: Single pore volume validation
ax1.plot(
    t_plot[:-1],
    cout_apvd[plot_start : plot_end - 1],
    label=f"APVD only (std={std_apv:.0f} m³)",
    linewidth=2,
    alpha=0.6,
    color="gray",
)
ax1.plot(
    t_plot[:-1],
    cout_diffusion_single[plot_start : plot_end - 1],
    label="Diffusion (single pore volume)",
    linewidth=2.5,
    linestyle="-",
    color="black",
)
ax1.plot(
    t_plot[:-1],
    cout_apvd_approx[plot_start : plot_end - 1],
    label=f"Approx: std={sigma_v_diff_disp:.0f} m³",
    linewidth=2,
    linestyle="--",
    color="C0",
)
ax1.set_xlabel("Date")
ax1.set_ylabel("Concentration")
ax1.set_title("Single Pore Volume: Diffusion vs Equivalent Std")
ax1.legend()
ax1.grid(True, alpha=0.3)

# Right panel: Multiple pore volumes validation
ax2.plot(
    t_plot[:-1],
    cout_apvd[plot_start : plot_end - 1],
    label=f"APVD only (std={std_apv:.0f} m³)",
    linewidth=2,
    alpha=0.6,
    color="gray",
)
ax2.plot(
    t_plot[:-1],
    cout_diffusion_multi[plot_start : plot_end - 1],
    label="Diffusion exact (multi PV)",
    linewidth=2.5,
    linestyle="-",
    color="black",
)
ax2.plot(
    t_plot[:-1],
    cout_diffusion_fast_multi[plot_start : plot_end - 1],
    label="diffusion_fast (standard)",
    linewidth=2,
    linestyle="-",
    color="C2",
)
ax2.plot(
    t_plot[:-1],
    cout_apvd_total[plot_start : plot_end - 1],
    label=f"Approx: std={sigma_v_total:.0f} m³",
    linewidth=2,
    linestyle="--",
    color="C1",
)
ax2.set_xlabel("Date")
ax2.set_ylabel("Concentration")
ax2.set_title("Multi Pore Volume: Exact vs diffusion_fast vs Combined Std")
ax2.legend()
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

../_images/examples_05_Diffusion_Dispersion_13_0.png

Summary#

Formulas for equivalent APVD spreading [m³]#

\[\sigma_{diff} = \frac{\bar{V}}{L} \sqrt{\frac{2 D_m R \bar{V}}{Q}} \quad \text{(molecular diffusion — decreases with higher } Q \text{)}\]

\[\sigma_{disp} = \bar{V} \sqrt{\frac{2 \alpha_L}{L}} \quad \text{(microdispersion — independent of } Q \text{; } R \text{ cancels)}\]

\[\sigma_{total} = \sqrt{\sigma_{apv}^2 + \sigma_{diff}^2 + \sigma_{disp}^2} \quad \text{(total spreading)}\]

Key takeaways#

Standard approach: Use diffusion_fast for transport modeling — it handles macrodispersion, microdispersion, and molecular diffusion
Macrodispersion (APVD) depends on aquifer properties and boundary conditions; microdispersion (\(\alpha_L\)) is an aquifer property; molecular diffusion (\(D_m\)) depends on compound
Variance combination (\(\sigma_{total}\)) is an alternative fast approach but only valid for gamma-parameterized APVD from streamline analysis — not for calibrated APVD or discrete volumes
When APVD is calibrated from measurements using advection, all mixing is already included — do not add on top. When calibrating with diffusion_fast or diffusion, the components are separated automatically

Macrodispersion and Microdispersion in Solute Transport

Contents

Macrodispersion and Microdispersion in Solute Transport#

Dispersion as Scale-Dependent Heterogeneity#

Spreading Components#

Unit conversions#

Choosing the Right Approach#

By APVD source#

Modeling approaches#

Variance addition (alternative fast approach)#

1. System Parameters#

2. Converting Microdispersion and Molecular Diffusion to Equivalent APVD Spreading#

Conversion from spatial to volume units#

Formulas for equivalent spreading in volume units#

3. Combining Macrodispersion and Microdispersion#

Including Microdispersion via Variance Addition (Alternative Fast Approach)#

4. Numerical Validation#

Summary#

Formulas for equivalent APVD spreading [m³]#

Key takeaways#

References#