ABSTRACT (original article): We have analyzed 39 major ammonium nitrate (AN) accidents (1916–2022) that reveals systemic catastrophic risks have not diminished over the past century. This challenges the prevailing narrative of improved process safety of AN on account of iterative guidelines and technology advances. Applying extreme value theory (EVT), we model temporal occurrence as a homogeneous Poisson process (maximum likelihood estimate # # 0.23 events/year) and casualty severity as heavy-tailed (generalized Pareto distribution, shape parameter # # 1.13). These findings indicate persistent systemic patterns, with conventional scenario-based risk assessments underestimating tail risks realized. This is shown by the implied “shadow mean” exceeds the sample mean by 1.5–19-fold, and societal risk profiles occupy the “intolerable” region per frequency-number (F/N) curve criteria. Such approaches, reliant on expert-defined failure pathways, fail to capture emergent stochastic dynamics evident in historical ensembles, such as governance lapses contributing to the 218 fatalities in Beirut (2020). We propose an EVT-based systemic ensemble framework as a complementary method for chemical risk governance. We need to monitor realized systemic risks over mechanism-driven simulations or narrative expectations. CITATION (original article): Gürkan Sin (2025), Persistent Catastrophic Risks in Ammonium Nitrate: Extreme Value Theory Challenges Reductionist Safety Paradigms, ChemRxiv. https://doi.org/10.26434/chemrxiv-2025-nj4hf-v2

Truncated GPD - Exact closed-form solution

ClearAll["Global`*"];
Parameters (kept symbolic)
ξ = ξ; σ = σ; μ = μ; U = U;
Strong assumptions for perfect simplification
$Assumptions = ξ > 0 && σ > 0 && U > μ && ξ ≠ 1;​​​​Print["=== Step 1: PDF of the Truncated GPD ==="];​​fX[x_] := (1/σ) (1 + ξ (x - μ)/σ)^(-(1 + 1/ξ))​​FX[x_] := 1 - (1 + ξ (x - μ)/σ)^(-1/ξ)​​​​fX[x] // TraditionalForm​​FX[x] // TraditionalForm
=== Step 1: PDF of the Truncated GPD ===
Out[]//TraditionalForm=
-
1
ξ
-1
ξ(x-μ)
σ
+1
σ
Out[]//TraditionalForm=
1-
-1/ξ
ξ(x-μ)
σ
+1
denom = FX[U] - FX[μ];​​Print["Normalizing factor:"];​​denom // TraditionalForm
Normalizing factor:
Out[]//TraditionalForm=
1-
-1/ξ
ξ(U-μ)
σ
+1
Step 2: Expected Value Definition
Print["\n=== Step 2: Expected Value (Original Integral) ==="];​​Etrunc = 1/denom * Integrate[x fX[x], {x, μ, U}];​​Etrunc // TraditionalForm
=== Step 2: Expected Value (Original Integral) ===
Out[]//TraditionalForm=
1
ξ
σ
-1/ξ
(σ(-μξ+σ+ξU))

1
ξ
σ
(-μξ+σ+U)-(μ(-ξ)+μ+σ)
1
ξ
(-μξ+σ+ξU)

(ξ-1)1-
-1/ξ
ξ(U-μ)
σ
+1
if μξ<σ∧μ>0
Step 3: Substitution
Print["\n=== Step 3: Substitution y = 1 + ξ(x-μ)/σ ==="];​​yU = 1 + ξ (U - μ)/σ;​​xSub = μ + (σ/ξ)(y - 1);​​integrand = xSub * (1/ξ) y^(-1 - 1/ξ);​​Print["Integrand after substitution:"];​​integrand // TraditionalForm
=== Step 3: Substitution y = 1 + ξ(x-μ)/σ ===
Integrand after substitution:
Out[]//TraditionalForm=
-
1
ξ
-1
y
μ+
σ(y-1)
ξ
ξ
Step 4-6: Integration and Simplification
Print["\n=== Steps 4-6: Solving the Integrals ==="];​​numer = Integrate[x fX[x], {x, μ, U}]; (* full numerator *)​​finalE = Simplify[numer / denom];​​​​Print["Final expression after dividing by normalizing factor:"];​​finalE // TraditionalForm
=== Steps 4-6: Solving the Integrals ===
Final expression after dividing by normalizing factor:
Out[]//TraditionalForm=
-
-1/ξ
(-μξ+σ+ξU)

1
ξ
σ
(-μξ+σ+U)-(μ(-ξ)+μ+σ)
1
ξ
(-μξ+σ+ξU)

(ξ-1)
1
ξ

σ
-μξ+σ+ξU

-1
if μξ<σ∧μ>0
Step 7: Clean closed-form with y_U
Print["\n=== Step 7: Final Closed-Form Expression (y_U form) ==="];​​yUform = μ - σ/(ξ - 1) + σ (yU - 1)/(ξ (ξ - 1) (yU^(1/ξ) - 1));​​yUform // TraditionalForm​​​​Print"\n✅ PROOF: This y_U closed form is mathematically identical to the original integral";​​Print["(verified by Simplify[Etrunc - yUform] === 0)"];​​Simplify[finalE - yUform] // TraditionalForm​​​​Print["\nThe expression above is the exact analytical solution of the truncated GPD mean."];
=== Steps 4-6: Solving the Integrals ===
Final expression after dividing by normalizing factor:
Out[]//TraditionalForm=
-
-1/ξ
(-μξ+σ+ξU)

1
ξ
σ
(-μξ+σ+U)-(μ(-ξ)+μ+σ)
1
ξ
(-μξ+σ+ξU)

(ξ-1)
1
ξ

σ
-μξ+σ+ξU

-1
if μξ<σ∧μ>0
=== Step 7: Final Closed-Form Expression (y_U form) ===
Out[]//TraditionalForm=
μ-
σ
ξ-1
+
U-μ
(ξ-1)
1
ξ
ξ(U-μ)
σ
+1
-1
✅ PROOF: This y_U closed form is mathematically identical to the original integral
(verified by Simplify[Etrunc - yUform] === 0)
Out[]//TraditionalForm=
0 if μξ<σ∧μ>0
The expression above is the exact analytical solution of the truncated GPD mean.

Verification - approach 1: Shadow Mean of Truncated GPD

Exact symbolic integration (Mathematica) vs Paper’s Numerical/MC
ClearAll["Global`*"];
Parameters exactly as reported in the paper (Table 1)
ξ = 1.13;​​σ = 24.09;​​μ = 0; (* GPD location parameter = 0 (Pareto Type II / Lomax) *)
Upper truncation bounds from Table 2 of the paper
U1 = 10000; (* Bhopal-scale bound *)​​U2 = 8.19 * 10^9; (* World population bound *)
GPD PDF and CDF (standard form used in the paper)
gpdPDF[x_] := (1/σ) (1 + ξ (x - μ)/σ)^(-(1 + 1/ξ)) /; x >= μ;​​gpdCDF[U_] := 1 - (1 + ξ (U - μ)/σ)^(-1/ξ);
1. ANALYTICAL SOLUTION – exact symbolic integration
This follows the derivation PDF step-by-step (substitution y = 1 + ξ (x-μ)/σ)
analyticShadowMean[U_] := Module[{int},​​ int = Integrate[x * gpdPDF[x], {x, μ, U}, Assumptions -> U > μ && ξ > 0];​​ int / gpdCDF[U]​​];
2. NUMERICAL SOLUTION – high-precision Monte Carlo style
This replicates the paper’s numerical/MC method for the truncated mean
3. Comparison Table (exactly matches Table 2 structure)

Verification - approach 2: yU Closed-Form Shadow Mean

Matches paper Table 2 exactly (within rounding)
yU definition
Simplified closed-form (the one you showed)
Comparison with paper values

CITE THIS NOTEBOOK

Persistent catastrophic risks in ammonium nitrate: extreme value theory challenges reductionist safety paradigms​
by Gürkan Sin​
Wolfram Community, STAFF PICKS, March 25, 2026
https://community.wolfram.com/groups/-/m/t/3669326