Individual Risk Contours

Technical methodology for calculating Individual Risk (IR) contour maps using TekRiskGen2 — grid-based IR aggregation, Marching Squares contouring, and GeoJSON output for QRA visualization

1. Executive Summary

Individual Risk (IR) contours are iso-risk lines drawn on a map that connect all geographic points sharing the same annual probability of fatality due to industrial hazards. They are the cornerstone of Quantitative Risk Assessment (QRA) in the chemical and oil & gas industries.

TekRiskGen2 computes IR contours by:

Aggregating all defined risk scenarios (fires, explosions, toxic clouds) at a facility

Evaluating the fatality probability at thousands of grid points around the site

Extracting contour lines at standard risk levels (e.g., 1×1051 \times 10^{-5}, 1×1061 \times 10^{-6} per year)

Projecting results onto a geographic map as GeoJSON polygons

Map visualization

The result is a set of nested contour polygons — each representing a different risk level — overlaid on a Mapbox GL map, allowing engineers to visually assess whether vulnerable receptors (schools, hospitals, residential areas) fall within unacceptable risk zones.


2. Key Concepts

Individual Risk (IR)

Annual probability that a hypothetical person, continuously present and unprotected at a specific location, will be killed due to an industrial accident. Expressed as a dimensionless number per year (e.g., IR=1×105IR = 1 \times 10^{-5}/year = 1 in 100,000 chance).

Scenario Frequency (f)

How often a particular accident scenario is expected to occur, in events/year. Derived from historical failure rate databases (OREDA, OGP/IOGP) combined with event tree analysis.

Fatality Probability (Pf)

Probability of death at a specific location given a scenario occurs. Derived from probit models (thermal/blast) or point-in-polygon tests (flash fire — Pf=1.0P_f = 1.0 inside the LEL envelope, 00 outside).

Contour Levels

Specific IR values at which iso-risk lines are drawn. TekRiskGen2 uses seven standard levels: 10210^{-2} through 10810^{-8} per year.

2.1 Core Formula

The Individual Risk at any point (x,y)(x, y) is the sum of contributions from all scenarios:

IR(x,y)=ifi×Pf,i(x,y)IR(x, y) = \sum_{i} f_i \times P_{f,i}(x, y)
SymbolDescriptionUnits
IR(x,y)IR(x, y)Individual risk at location (x,y)(x, y)per year
fif_iFrequency of scenario iievents/year
Pf,i(x,y)P_{f,i}(x, y)Probability of fatality at (x,y)(x, y) given scenario ii occursdimensionless [0,1][0, 1]
iiScenario index (each risk model is a scenario)

3. Supported Risk Models

TekRiskGen2 supports five consequence model types, each contributing to the overall IR calculation:

ModelHazard TypeSymmetryPfP_f MethodWind Dependence
FireballThermal radiationRadialDistance-based interpolation on fatality profileNone — omni-directional
Pool FireThermal radiationRadialDistance-based interpolation on fatality profileNone — radially symmetric
Jet FireThermal radiationRadialDistance-based interpolation on fatality profileNone — radially symmetric
VCEBlast overpressureRadialDistance-based interpolation on fatality profileNone — radially symmetric
Flash FireFlame engulfmentDirectionalPoint-in-polygon on rotated LEL envelopeCritical — uses full 16-direction wind rose

Fireball, Pool Fire, Jet Fire, VCE — These models assume the hazard effect decays with distance from the source in all directions equally. The fatality probability at any point depends only on the Euclidean distance to the source, following a pre-computed fatality profile (fatality percentage vs. distance).


4. Step-by-Step Methodology

Cargando diagrama...

Step 1: Scenario Definition

Each enabled risk model defines a scenario characterized by:

ParameterDescriptionSource
Frequency (ff)Annual occurrence rateRisk model configuration
Source locationGeographic coordinates (lat, lon)Source definition
Model typefireball, poolfire, jetfire, vce, flashfireRisk model type
Fatality profilePfP_f vs. distance curve (radial) or LEL polygon (flash fire)Prior consequence calculation results

Fatality profiles are extracted from previously completed consequence calculations stored in RiskCalculation.results.fatalidades. Each profile consists of distance-fatality percentage pairs:

Distance (m)Fatality (%)
0100
2595
5070
10030
2005
3500

Flash Fire profiles

For flash fire models, the fatality "profile" is replaced by a LEL polygon — the geographic outline of the flammable cloud at the Lower Explosive Limit concentration for a single wind direction.


Step 2: Wind Rose Integration (Flash Fire)

Flash fire scenarios require integration with a Pasquill-Gifford wind rose to account for the directional nature of vapor cloud propagation.

Wind Rose Structure

The wind rose divides the compass into 16 equi-spaced directions at 22.5° intervals. Each direction carries a probability (fraction of time the wind blows from that direction), split into day and night periods:

PeriodHoursStability Classes
Day06:00 – 18:00A, B, C, D
Night18:00 – 06:00D, E, F

Polygon Rotation

For each of the 16 wind directions, the base LEL polygon is rotated around the source point:

x=xcos(θ)ysin(θ)x' = x \cos(\theta) - y \sin(\theta) y=xsin(θ)+ycos(θ)y' = x \sin(\theta) + y \cos(\theta)

Where (x,y)(x, y) is the original polygon vertex relative to source, θ\theta is the rotation angle, and (x,y)(x', y') is the rotated vertex.

The result is 16 rotated LEL polygons, each associated with the wind probability for its direction.

Flash Fire PfP_f Formula

Pf,ff(x,y)=min ⁣(1.0,  jp(θj)1{(x,y)Polygon(θj)})P_{f,\text{ff}}(x, y) = \min\!\left(1.0,\; \sum_{j} p(\theta_j) \cdot \mathbb{1}\{(x,y) \in \text{Polygon}(\theta_j)\}\right)

Where:

SymbolDescription
θj\theta_jWind direction j{0°,22.5°,,337.5°}j \in \{0°, 22.5°, \ldots, 337.5°\}
p(θj)p(\theta_j)Probability of wind from direction θj\theta_j (0(0 to 1)1)
Polygon(θj)\text{Polygon}(\theta_j)LEL polygon rotated to direction θj\theta_j
1{}\mathbb{1}\{\cdot\}Indicator function: 1 if point is inside polygon, 0 otherwise

Probability cap

The sum is capped at 1.0 because the maximum fatality probability cannot exceed certainty.


Step 3: Grid Generation

A uniform 2D rectangular grid is generated to evaluate IR across the study area.

ParameterDescriptionDefault
CenterCentroid of all source coordinatesAuto-calculated
ResolutionSpacing between adjacent grid points25 m
ExtentHalf-size of the grid from centerAuto-calculated

Available resolutions: 1, 5, 10, 25, 50, 100 meters.

Auto-Calculation of Extent

Find the maximum fatality profile distance across all scenarios

Multiply by 1.3 (30% padding factor)

Round up to the nearest 100 m

Grid Dimensions

Points per side=2×extentresolution+1\text{Points per side} = \frac{2 \times \text{extent}}{\text{resolution}} + 1 Total grid points=(Points per side)2\text{Total grid points} = (\text{Points per side})^2

Example: With extent = 2,000 m and resolution = 25 m → 161 points per side → 25,921 total points.

Coordinate System

The grid uses a local Cartesian coordinate system in meters:

  • Origin: Grid center (centroid of sources)
  • X-axis: West → East (positive eastward)
  • Y-axis: South → North (positive northward)

Flat-earth approximation

This approximation is valid for distances up to ~50 km from the center — adequate for typical industrial sites.


Step 4: IR Calculation at Each Grid Point

For every point in the grid, the system accumulates IR contributions from all enabled scenarios.

Fireball, Pool Fire, Jet Fire, VCE — For each grid point (xi,yj)(x_i, y_j) and each radial scenario ss:

Step 1 — Euclidean distance: d=(xixs)2+(yjys)2d = \sqrt{(x_i - x_s)^2 + (y_j - y_s)^2}

Step 2 — Quick reject: If d>dmaxd > d_{\max} → skip (Pf=0P_f = 0)

Step 3 — Interpolate fatality profile: If dd0d \leq d_0: Pf=P0/100P_f = P_0 / 100. Otherwise find interval [dk,dk+1][d_k, d_{k+1}] containing dd:

t=ddkdk+1dk,Pf=Pk+t(Pk+1Pk)100t = \frac{d - d_k}{d_{k+1} - d_k}, \qquad P_f = \frac{P_k + t \cdot (P_{k+1} - P_k)}{100}

Step 4 — Accumulate: IR(xi,yj)+=fs×PfIR(x_i, y_j) \mathrel{+}= f_s \times P_f


Step 5: Contour Extraction (Marching Squares)

Once the IR grid is fully populated, contour lines at each target level are extracted using the Marching Squares algorithm.


Step 6: Geographic Conversion

Contour polygons are converted from local meter coordinates back to geographic coordinates for map display.

Conversion Formulas

lat=latcenter+y111,320\text{lat} = \text{lat}_{\text{center}} + \frac{y}{111{,}320} lon=loncenter+x111,320×cos(latcenter×π/180)\text{lon} = \text{lon}_{\text{center}} + \frac{x}{111{,}320 \times \cos(\text{lat}_{\text{center}} \times \pi / 180)}

Where 111,320111{,}320 is the approximate meters per degree of latitude (WGS84).

GeoJSON Output

Each contour polygon is packaged as a GeoJSON Feature:

{
  "type": "Feature",
  "properties": {
    "level": 1e-5,
    "levelFormatted": "1e-5",
    "color": "#EA580C",
    "opacity": 0.3,
    "type": "ir_contour"
  },
  "geometry": {
    "type": "Polygon",
    "coordinates": [[[lon1, lat1], [lon2, lat2], "...", [lon1, lat1]]]
  }
}

Area Calculation

The area enclosed by each contour is computed using the shoelace formula on local meter coordinates:

A=12i(xiyi+1xi+1yi)[m2]A = \frac{1}{2} \left| \sum_{i} (x_i \cdot y_{i+1} - x_{i+1} \cdot y_i) \right| \quad [\text{m}^2]
Area RangeDisplay Unit
<10,000< 10{,}000
10,0001,000,00010{,}000 – 1{,}000{,}000hectares (ha)
1,000,000\geq 1{,}000{,}000km²

5. Risk Tolerability Criteria

The following table summarizes international standards for individual risk tolerability:

StandardRegionIntolerable LimitBroadly AcceptableReference
UK HSE R2P2United Kingdom1×1041 \times 10^{-4} (public)1×1061 \times 10^{-6}HSE (2001)
ASEA / NTE-002Mexico1×1031 \times 10^{-3}1×1061 \times 10^{-6}NTE-002-CGPC
RIVMNetherlands1×1051 \times 10^{-5}1×1081 \times 10^{-8}RIVM guidelines
Hong KongHong Kong1×1051 \times 10^{-5}1×1061 \times 10^{-6}HK risk guidelines
HIPAP No. 4Australia1×1051 \times 10^{-5}1×1061 \times 10^{-6}NSW HIPAP No. 4
EPA / OSHAUnited States1×1041 \times 10^{-4}1×1061 \times 10^{-6}EPA RMP / OSHA PSM

ALARP Principle

ALARP — As Low As Reasonably Practicable

Between the intolerable and broadly acceptable limits lies the ALARP region. Within this zone, risk is tolerated only if further reduction is impracticable or disproportionately costly.

ZoneRisk LevelRequired Action
Intolerable1×103\geq 1 \times 10^{-3} to 1×1041 \times 10^{-4} (varies by standard)Risk cannot be justified except in extraordinary circumstances
ALARPBetween intolerable and acceptable thresholdsTolerable only if reduction is impracticable
Broadly Acceptable1×106\leq 1 \times 10^{-6} to 1×1081 \times 10^{-8} (varies by standard)No further action needed

6. Contour Color Code

TekRiskGen2 uses a seven-level color scheme to distinguish risk zones on the map:

Level (per year)ColorHex CodeRisk Interpretation
1×1021 \times 10^{-2}Indigo#4B0082Extreme — immediate action required
1×1031 \times 10^{-3}Dark Red#8B0000Intolerable — exceeds all standards
1×1041 \times 10^{-4}Red#DC2626Intolerable for public (UK HSE, US EPA)
1×1051 \times 10^{-5}Orange#EA580CUpper ALARP bound (many jurisdictions)
1×1061 \times 10^{-6}Yellow#EAB308Broadly acceptable threshold
1×1071 \times 10^{-7}Light Green#84CC16Low risk
1×1081 \times 10^{-8}Green#22C55ENegligible risk

Non-standard levels

Contour levels not in this standard set default to gray (#6B7280).


7. Limitations


8. Bibliographic References