Companion Python models for Mathematical Geoenergy: Discovery, Depletion, and Renewal by Pukite, Challou & Coyne (Wiley/IEEE Press, 2019), Chapter 15 (Latent Energy: Hydrological Cycle), Section 15.3 — Rainfall.
Entropic dispersion shapes the probability distribution of daily rainfall amounts, just as it shapes lake sizes and cloud areas elsewhere in Chapter 15. Two complementary models are derived.
Sigmoid / preferential-attachment model (Eqs. 15-12, 15-13). Clouds develop by preferential attachment — already-large clouds accumulate moisture faster. If the growth rate follows an exponential ramp \(g(x) = k\,(e^{ax}-1)\) then the critical-point PDF for reaching sufficient intensity takes the inverse-square (Lomax-like) form \(p(x) = \frac{r}{(r + g(x))^2}\) When integrated, this traces the descending (right-hand) side of a logistic sigmoid, matching North Carolina January rainfall histograms with parameters $r=1$, $a = 1/17$, $k = 2$ mm/h (GEM Fig. 15-5/15-6).
MaxEnt Bessel-K model (Eqs. 15-14 – 15-18). A cleaner, single-parameter derivation applies the super-statistics framework:
The corresponding PDF is the Bessel-K₀ distribution (also seen in terrain slopes, Ch. 16): \(p(E) = \frac{2}{\bar{E}}\,K_0\!\left(2\sqrt{\frac{E}{\bar{E}}}\right)\)
The sole free parameter $\bar{E}$ equals the sample mean directly from GEM Table 15-1 (combined mean = 2.29 mm/h), and the model outperforms the fractal/Hurst-Kolmogorov approach of Papalexiou et al. (2011) in matching empirical data.
Eq. 15-13 — Preferential-attachment growth rate
\[g(x) = k\!\left(e^{ax} - 1\right)\]The growth rate $g(x)$ encodes preferential attachment: its own derivative satisfies $dg/dx = a\,(g+k)$, meaning the growth rate is proportional to the current level plus a base $k$. Starting from zero at $x=0$, it rises exponentially with rainfall amount $x$.
Eq. 15-12 — Sigmoid rainfall PDF
\[p(x) = \frac{r}{\bigl(r + g(x)\bigr)^{2}}\]With $g(x)$ from Eq. 15-13, this is a monotone-decreasing density that begins at $p(0) = 1/r$ and decays faster than an exponential once the preferential-attachment ramp takes hold. When integrated it traces the descending slope of a logistic sigmoid (GEM Fig. 15-6), fitting the observed North Carolina data with $r=1$, $a=1/17$, $k=2$ mm/h.
Eq. 15-14 — Conditional energy CCDF (Boltzmann / MaxEnt)
\[P(E \mid E_i) = e^{-E/E_i}\]Given a storm with characteristic energy $E_i$ (proportional to its rainfall rate), the MaxEnt distribution with known mean $E_i$ is the exponential. This is the Boltzmann/Gibbs probability for a potential- energy activation process: holding moisture aloft against gravity costs energy proportional to the moisture mass.
Eq. 15-15 — Energy-rate proportionality
\[\text{Rate}_i \;\sim\; E_i\]The physical premise: the energy needed to sustain a storm is proportional to the moisture mass, hence to the rainfall rate. This identifies the conditional mean $E_i$ as the per-storm rainfall rate, so the prior over $E_i$ is simultaneously a prior over storm intensities.
Eq. 15-16 — MaxEnt prior on storm energy
\[p(E_i) = \alpha\, e^{-\alpha E_i}, \qquad \alpha = 1/\bar{E}\]Storm energies vary from event to event. Knowing only the long-run mean $\bar{E}$ (the single measurable quantity from GEM Table 15-1), MaxEnt selects the exponential distribution for $E_i$ with rate $\alpha=1/\bar{E}$.
Eq. 15-17 — Marginalisation over storm energy (mixture integral)
\[P(E) = \int_0^{\infty} P(E \mid E_i)\; p(E_i)\, dE_i\]The unconditional CCDF is found by integrating the conditional Boltzmann distribution against the MaxEnt prior. This mixture over a hidden variable (the per-storm energy $E_i$) is the super-statistics construction.
Eq. 15-18 — Bessel-K CCDF solution
\[P(E) = 2\sqrt{\frac{E}{\bar{E}}}\;K_1\!\!\left(2\sqrt{\frac{E}{\bar{E}}}\right)\]Result of evaluating Eq. 15-17 via the integral identity $\int_0^\infty \exp(-a/t - bt)\,dt = 2\sqrt{a/b}\,K_1(2\sqrt{ab})$ with $a=E$, $b=\alpha=1/\bar{E}$. Here $K_1$ is the modified Bessel function of the second kind of order 1. The CCDF satisfies $P(0)=1$ and decays sub-exponentially for large $E$.
PDF — Bessel-K₀ rainfall density
\[p(E) = -\frac{dP}{dE} = \frac{2}{\bar{E}}\,K_0\!\!\left(2\sqrt{\frac{E}{\bar{E}}}\right)\]Obtained by differentiating Eq. 15-18 using the Bessel recurrence $dP/dz = -z K_0(z)$ where $z = 2\sqrt{E/\bar{E}}$. The $K_0$ PDF normalises to 1 and has mean $\bar{E}$. The same $K_0$ form appears in the terrain slope PDF (GEM Ch. 16) and is the signature distribution of two nested MaxEnt exponentials.
| File | Purpose |
|---|---|
rainfall_symbolic.py |
SymPy derivation of Eqs. 15-12 through 15-18 with algebraic checks |
rainfall_numerical.py |
Numerical implementation, Monte-Carlo sampling, composite figure |
rainfall_model_output.png |
6-panel figure produced by rainfall_numerical.py |
# Install shared dependencies (from repository root)
pip install -r models/requirements.txt
# Symbolic derivation — verify all equations
python models/rainfall/rainfall_symbolic.py
# Numerical model — generate composite figure
MPLBACKEND=Agg python models/rainfall/rainfall_numerical.py
Two routes to the same data. The sigmoid model (Eqs. 15-12/13) is a three-parameter local fit, while the Bessel-K model (Eqs. 15-14 to 15-18) is a one-parameter global theory derived purely from the potential-energy argument. Both match the North Carolina data; the Bessel-K requires only the measured mean.
Super-statistics = two nested MaxEnt priors. Applying MaxEnt once (known mean $E_i$) gives an exponential for $E$. Applying it a second time (known mean $\bar{E}$ for the varying $E_i$) and marginalising yields the Bessel-K family — the same construction used for wind speed (Ch. 11) and terrain slopes (Ch. 16).
The tail is sub-exponential, not power-law. The Bessel-K CCDF decays as $\sqrt{z/\pi/2}\,e^{-z}$ with $z = 2\sqrt{E/\bar{E}}$, which is lighter than any power law. This is consistent with GEM Fig. 15-5 showing the rainfall tail stays close to exponential — not fractal — once sampling noise is accounted for.
One parameter. The sole free parameter $\bar{E}$ equals the arithmetic sample mean of the rainfall rate data. No fitting is required; the parameter can be read directly from GEM Table 15-1 (combined mean = 2.29 mm/h across 7 storm events, 29,536 observations).
Preferential-attachment dg/dx = a(g+k). The growth-rate identity shows that the cloud’s moisture-accumulation rate increases in proportion to its current size — the same feedback that drives Zipf’s law in city growth and citation networks.
Monte-Carlo verification. Drawing storm energies $E_i$ from the exponential prior and then drawing rainfall $E$ from the conditional exponential reproduces the Bessel-K distribution exactly, confirming the super-statistics construction with no approximation.
Pukite, P., Challou, D., & Coyne, C. (2019). Mathematical Geoenergy: Discovery, Depletion, and Renewal. Wiley/IEEE Press. ISBN 978-1-119-43429-0. Chapter 15, Section 15.3 — Equations 15-12 through 15-18.
Papalexiou, S. M., Koutsoyiannis, D., & Makropoulos, C. (2011). “How extreme is extreme? An assessment of daily rainfall distribution tails.” Hydrology and Earth System Sciences, 17(2), 851–862.
Gradshteyn, I. S., & Ryzhik, I. M. (2007). Table of Integrals, Series, and Products (7th ed.), Eq. 3.914.5. Academic Press.