Mathematical models for dispersive transport of solutes and contaminants through porous media derived from first principles in Mathematical Geoenergy (GEM), Chapter 20 (Dispersed Energy: Particulates and Transport in the Environment).
Dispersive transport through porous media is a prime example of a Maximum Entropy (MaxEnt) dispersion problem applied to the movement of fluids, contaminants, and particles in heterogeneous geological systems. The derivation in GEM Ch. 20 starts from the Fokker-Planck equation (FPE) and shows that, when the diffusion coefficient and/or drift velocity are themselves uncertain (disordered), the resulting concentration profiles are radically different from the standard Gaussian prediction.
Disorder in the diffusion coefficient — Applying MaxEnt to D (only its mean D₀ is known) and integrating the Gaussian FPE solution over all possible D values yields a Laplace (double-exponential) profile (Eq. 20-7) rather than a Gaussian.
Disorder in the drift velocity — Applying MaxEnt to the velocity ν (only the mean ν₀ is known) and using the delta-function FPE (no diffusion) yields a reciprocal exponential breakthrough curve (Eq. 20-11), or a Cauchy/Lorentzian profile when a dominant primary flow is present (Eq. 20-22).
Disorder in decomposition rates — Mixing materials with MaxEnt- distributed decay rates produces a hyperbolic decline (Eq. 20-30) rather than the standard single-rate exponential.
These results match breakthrough curve measurements from natural watersheds, laboratory porous-medium experiments, and contamination studies — all with only one or two free parameters.
Eq. 20-1 — MaxEnt exponential velocity distribution
\[p(\nu) = \frac{1}{\nu_0}\cdot e^{-\nu/\nu_0}, \quad \nu \geq 0\]Applying the Principle of Maximum Entropy with only a known mean ν₀ gives an exponential distribution for fluid velocity. This is the “disordered drift” (or disordered advection / convection) model.
Eq. 20-2 — Fokker-Planck equation
\[\frac{\partial n(x,t)}{\partial t} = -\mu F\frac{\partial n}{\partial x} + D\cdot\frac{\partial^2 n}{\partial x^2} - \frac{n}{\tau}\]General transport equation for solute concentration $n(x,t)$, with drift term (mobility $\mu$, force $F$), Fickian diffusion (coefficient $D$), and particle-loss term (time constant $\tau$).
Eq. 20-3 — Gaussian FPE solution (homogeneous medium)
\[n(x,t) = \frac{1}{\sqrt{4\pi D t}}\cdot e^{-\frac{(x-\mu F t)^2}{4Dt}}\cdot e^{-t/\tau}\]Standard Green’s function of the FPE for constant $D$ and $\mu F$. Replaced by the Laplace profile once disorder is introduced.
Eq. 20-4 — Integration over D distribution
\[n(x,t) = \int_0^\infty p(D)\cdot n(x,t\mid D)\\,dD\]Eq. 20-5 — MaxEnt diffusion coefficient distribution
\[p(D) = \frac{1}{D_0}\cdot e^{-D/D_0}\]In a disordered porous medium, local diffusion coefficients vary pore to pore. Without information beyond the mean $D_0$, MaxEnt gives an exponential distribution.
Eq. 20-6 — Integration identity
\[\int_0^\infty \frac{e^{-Ax - B/x}}{\sqrt{x}}\\,dx = \sqrt{\frac{\pi}{A}}\cdot e^{-2\sqrt{AB}}\]Standard Laplace-Bessel integral (Gradshteyn & Ryzhik §3.471), the key analytic tool used to evaluate the MaxEnt diffusion integral. Follows from the Bessel-K identity $\int_0^\infty x^{\nu-1}e^{-\beta x-\gamma/x}dx = 2(\gamma/\beta)^{\nu/2}K_\nu(2\sqrt{\beta\gamma})$ with $\nu=\tfrac{1}{2}$ and $K_{1/2}(z)=\sqrt{\pi/(2z)}\cdot e^{-z}$.
Eq. 20-7 — Dispersive transport formula (MaxEnt D)
\[n(x,t) = \frac{e^{-|x - \mu F t|/\sqrt{D_0 t}}}{2\sqrt{D_0 t}}\]Closed-form result of integrating Eq. 20-4 over the MaxEnt prior (Eq. 20-5) using identity Eq. 20-6 with $A = 1/D_0$, $B = (x-\mu F t)^2/(4t)$. The Gaussian profile of the homogeneous case is replaced by a Laplace (double-exponential) profile — the characteristic “fat tail” of dispersive transport.
Eq. 20-14 — Delta-function FPE (no diffusion)
\[n(x,t) = \int_0^\infty p(\nu)\cdot\delta(x - \nu t)\\,d\nu\]Eq. 20-16 — Closed-form delta result
\[n(x,t) = \frac{1}{t}\cdot p\left(\frac{x}{t}\right), \qquad J(x,t) = \frac{x}{t^2}\cdot p\left(\frac{x}{t}\right)\]Using the Jacobian of the delta function. This remarkably simple result converts any velocity PDF directly to a concentration profile.
Eq. 20-11 — No-diffusion reciprocal exponential
\[n(x,t) = \frac{e^{-x/(\nu_m t)}}{\nu_m t}\]Applying Eq. 20-16 with the MaxEnt velocity PDF (Eq. 20-1) gives a “reciprocal exponential” — the time parameter appears in the denominator of the exponent because we are dealing with rates, not times.
Eq. 20-17 — Primary flow with drag
\[\nu = \nu_0 - \nu_\text{drag}\]Eq. 20-20 — Lomax drag velocity distribution
\[p(\nu_\text{drag}) = \frac{\nu_d}{(\nu_d + \nu_\text{drag})^2}\]A scale-free (Lomax / Pareto-II) distribution for the drag velocity, derived from a ratio distribution of path increments.
Eq. 20-22 — Cauchy/Lorentzian breakthrough profile
\[n(x,t) = \frac{\nu_d\cdot t}{\bigl[(\nu_d + \nu_0)\cdot t - x\bigr]^2}\]Applying Eq. 20-16 with the Lomax velocity PDF (Eq. 20-21) gives a Cauchy/Lorentzian profile, which drops off as $1/t^2$ — the hallmark of heavy-tailed breakthrough curves with a dominant primary flow.
Eq. 20-24 — Laplace velocity distribution
\[p(\nu) = \frac{1}{\nu_m}\cdot e^{-|\nu - \nu_0|/\nu_m}\]MaxEnt distribution subject to a known mean $\nu_0$ and a known mean absolute deviation $\nu_m$. This is the Laplace (double-exponential) or bilateral exponential — it has higher entropy than the Gaussian and is the least-biased estimator for those two constraints.
Eq. 20-25 — Laplace breakthrough concentration
\[n(x,t) = \frac{1}{\nu_m t}\cdot \exp\left(-\left|\frac{x}{\nu_m t} - \frac{\nu_0}{\nu_m}\right|\right)\]Applying Eq. 20-16 with the Laplace velocity PDF (Eq. 20-24). On a log scale the breakthrough curve appears as a perfect isosceles triangle (symmetric decay both sides of the peak at $x = \nu_0 t$). The ratio $\nu_m/\nu_0 = 0.18$ gives an excellent fit to the Le Borgne [7] laboratory data spanning 4 orders of magnitude.
Eq. 20-23 — Particle loss by multiplicative factor
\[n(x,t\mid\tau) = n(x,t)\cdot e^{-t/\tau}\]The FPE loss term adds an exponential envelope to any concentration profile.
Eq. 20-26 — Waste half-life (compact form)
\[P(t) = \frac{1}{1 + k\cdot t}\]Eq. 20-27 — MaxEnt decay-energy distribution
\[p(E) = \frac{1}{E_0}\cdot e^{-E/E_0}\]Eq. 20-29 — Conditional survival probability
\[P(t\mid E) = e^{-r\cdot E\cdot t}\]Eq. 20-30 — MaxEnt mixture half-life (derived)
\[P(t) = \int_0^\infty P(t\mid E)\cdot p(E)\\,dE = \frac{1}{1 + r\cdot E_0\cdot t}\]Mixing materials (oil-spill compounds, radioactive waste) with MaxEnt-distributed decay rates yields a hyperbolic decline rather than a single exponential. For large $t$, $P(t)\sim 1/(k t)$ vs $e^{-kt}$ — a much fatter tail, consistent with observed long-term persistence of contaminants. The identical functional form appears as the hyperbolic decline of oil reservoirs (GEM Ch. 14).
Eq. 20-31 — Lomax / entroplet size distribution
\[p(x) = \frac{S}{(S + x)^2}\]In maximally disordered environments (e.g. ice-crystal growth in cirrus clouds), particle growth rates are MaxEnt distributed, yielding the Lomax (Pareto-II) size distribution with shape $\kappa = 1$. This is the same functional form as the lake-size distribution (GEM Ch. 15), demonstrating the universality of entropic dispersion.
CDF:
\[F(x) = \frac{x}{S + x}\]CCDF:
\[\bar{P}(x) = \frac{S}{S + x} \\;\to\\; \frac{S}{x} \quad (x \gg S)\]| File | Purpose |
|---|---|
flow_symbolic.py |
Symbolic derivation of all equations using SymPy |
flow_numerical.py |
Numerical implementation, validation, and composite figure |
flow_model_output.png |
Output figure (6 panels) generated by flow_numerical.py |
Install dependencies (from models/requirements.txt):
pip install -r ../requirements.txt
Run symbolic derivation (all assertions print ✓):
python flow_symbolic.py
Run numerical model and generate figure:
MPLBACKEND=Agg python flow_numerical.py
Disorder transforms Gaussian → Laplace: introducing MaxEnt uncertainty in the diffusion coefficient $D$ (Eq. 20-5) replaces the Gaussian concentration profile with the Laplace (double-exponential) profile (Eq. 20-7). The “fat tails” of breakthrough curves arise naturally from this disorder, not from anomalous diffusion models.
Simple delta-function rule: when diffusion is negligible, $n(x,t) = (1/t)\cdot p(x/t)$ (Eq. 20-16) converts any velocity distribution directly to a breakthrough curve. This provides a powerful one-parameter fitting tool.
Universal fat-tail character: MaxEnt velocity disorder always produces power-law or hyperbolic drop-offs in breakthrough curves. The Cauchy profile ($\sim 1/t^2$, Eq. 20-22) and the reciprocal exponential (Eq. 20-11) both emerge from the same delta-function rule with different velocity PDFs.
Hyperbolic decline extends to waste: the same $P(t) = 1/(1+kt)$ formula (Eq. 20-30) describes both oil-reservoir production decline (GEM Ch. 14) and the persistence of mixed-rate waste dumps or oil spills. Disorder in decay rates is as ubiquitous as disorder in reservoir permeabilities.
Lomax universality: the size distribution of ice crystals, lake sizes (GEM Ch. 15), cloud sizes, and oil-reservoir volumes all follow the same Lomax (Pareto-II) form. Disorder in growth rates → Lomax sizes is a recurring pattern wherever MaxEnt applies to dispersive aggregation.
Minimal parameterisation: the models require at most two free parameters (a mean velocity and a scale factor). Complex numerical models with many parameters are unnecessary — MaxEnt disorder smooths out the underlying complexity.