Skip to content

Retrospective example: the Chernobyl release

A demonstration of the API, not a dose reconstruction

This page shows how to use the lower-level components of the package for a retrospective case that is not the directive's 24 h scenario: a ten-day release with a day-by-day source term, wind and release height, and populations that are evacuated at different times. The ENSI-G14 atmospheric model was designed for Swiss short-range scenarios, not for a lofted, burning-graphite release, and the inputs below were collected from web summaries of the literature (several primary sources could not be opened), with assumptions where nothing was found. The results are order-of-magnitude, partly calibrated numbers with large uncertainty, see Limitations. Do not quote them as doses received by anyone.

We compute, for four settlements close to the reactor, the dose to a resident (adult, 10-year-old child and 1-year-old toddler) from inhalation, submersion in the cloud and ground-shine, from the explosion on 26 April 1986 until the settlement was evacuated.

Inputs

Source term

Released activities (as quoted for the reactor accident, see the sources at the end of the page), for the three nuclides modelled here:

Nuclide Released activity
Cs-137 85 PBq
I-131 1760 PBq
Xe-133 6500 PBq

The time profile is taken from the published day-by-day I-131 release (26 April to 5 May, in PBq): 704, 204, 150, 102, 69, 62, 102, 107, 130, 130. We assume that all nuclides follow the same profile. The release height is reported as above 1 km on 26 April and about 600 m on the following days.

Settlements

Coordinates are taken from the settlements' Wikipedia pages, evacuation times from the literature: Pripyat and Yaniv on 27 April (afternoon), Chernobyl town on 5 May. The evacuation of Kopachi is only known to be "in the beginning of May"; we assume 3 May.

Weather

Only qualitative information was found: on 26 April the plume went west to north-west in a stable night-time atmosphere (about 8 m/s), on 27 April towards the north-east (2-3 m/s), and from 29 April the wind pushed it south towards Kyiv. We assume neutral conditions (category D) and 2.5 m/s from 27 April on, and no information on the direction after 2 May (all directions equally likely). Nothing was found on the rain near the reactor: we model it as a single free parameter, the share of the time during which it rains at 1 mm/h, calibrated below on the measured fallout.

Setting up

import warnings
warnings.filterwarnings('ignore')   # silence astropy / pandas deprecation noise

import numpy as np
import pandas as pd
from astropy.units import Bq, m, s, year, day

from ensi_g14_calculator.misc import Sv
from ensi_g14_calculator.orography import Orography
from ensi_g14_calculator.averaging.weather_averager import SourceParams, WeatherAverager, WeatherStats
from ensi_g14_calculator.dose.short_term_dose import ShortTermDoseCalculator
from ensi_g14_calculator.deposition.deposition import DepositionModel

# the reactor is the origin of the local (x east, y north) frame
REACTOR_LAT, REACTOR_LON = 51.38972, 30.09917
orography = Orography(flat=True, latitude0=REACTOR_LAT, longitude0=REACTOR_LON)

# time is counted in days after 26 April 00:00; the explosion was at 01:23
t_explosion = 1.39 / 24

villages = pd.DataFrame([
    # name,       lat,       lon,       evacuated (days after 26 Apr 00:00)
    ('Pripyat',   51.40611,  30.05722,  1.65),   # 27 April, 15:30
    ('Yaniv',     51.3942,   30.0564,   1.65),   # 27 April, 15:30
    ('Kopachi',   51.357278, 30.121972, 7.5),    # ASSUMED: 3 May, midday
    ('Chernobyl', 51.27222,  30.22417,  9.5),    # 5 May, midday
], columns=['name', 'lat', 'lon', 't_evacuation'])

x, y = orography.latlon_to_xy(villages.lat.values, villages.lon.values)
villages['distance_km'] = np.hypot(x.value, y.value) / 1e3
villages['bearing_deg'] = np.degrees(np.arctan2(x.value, y.value)) % 360   # from the reactor, clockwise from north
print(villages.round(2).to_string(index=False))

which prints (version 0.5.0):

     name   lat   lon  t_evacuation  distance_km  bearing_deg
  Pripyat 51.41 30.06          1.65         3.43       302.05
    Yaniv 51.39 30.06          1.65         3.01       279.53
  Kopachi 51.36 30.12          7.50         3.94       156.32
Chernobyl 51.27 30.22          9.50        15.68       146.42

# released activity in Bq, and the daily profile (share of the total released each day, from 26 April)
total_release = {'Cs-137': 85e15, 'I-131': 1760e15, 'Xe-133': 6500e15}
i131_per_day = np.array([704, 204, 150, 102, 69, 62, 102, 107, 130, 130.])
daily_share = i131_per_day / i131_per_day.sum()
day_start = t_explosion + np.arange(10)     # each 24 h slice starts at the time of the explosion

# one entry per day: weather category, wind speed (m/s, at release height), release height (m),
# and the documented direction the plume travels TO as (bearing, +- spread), or None if unknown.
days = [
    ('E', 8.0, 1000., (300, 30)),                       # 26 Apr: west / north-west, stable night, > 1 km
    ('D', 2.5,  600., (45, 30)),                        # 27 Apr: north-east
    ('D', 2.5,  600., (45, 30)),                        # 28 Apr: north-east (ASSUMED, as the 27th)
    ('D', 2.5,  600., (160, 30)), ('D', 2.5, 600., (160, 30)),   # 29 Apr - 2 May: south towards Kyiv
    ('D', 2.5,  600., (160, 30)), ('D', 2.5, 600., (160, 30)),
    ('D', 2.5,  600., None), ('D', 2.5, 600., None), ('D', 2.5, 600., None),   # 3 - 5 May: unknown
]

dose_calculator = ShortTermDoseCalculator()
deposition_model = DepositionModel()

The calculation

For each day and each settlement, the dispersion factor of that day's weather is computed with the WeatherAverager (the "bins" are a spread of wind directions, all with the same probability), then converted to doses with the same dose calculator the directive scenarios use:

  • inhalation and submersion from ShortTermDoseCalculator.inhalation and .submersion (the release of each day is integrated over the passage of the plume; only the part released before the evacuation counts),
  • ground-shine, from the activity deposited by dry fallout (\(F = \chi v_g\)) and by washout \(W\), and the soil weathering model of the directive, integrated from the deposition until the evacuation. Note that, unlike dry fallout, washout scavenges the whole plume column: it does not depend on the release height and does not need the plume to touch the ground.

The plume argument selects where the plume goes:

  • 'documented': the directions from the table above (all directions for the days without information),
  • 'early_aligned': as above, except that on 26 and 27 April, which carry more than half of the release, the plume is blowing straight at the settlement (a worst case for the early release; it is documented that the plume enveloped Pripyat on 27 April),
  • 'uniform': all directions equally likely on all days, an expectation when nothing is known about the wind.

height_scale multiplies all release heights, rain_fraction is the share of the time during which it rains at 1 mm/h (washout), k_s is the shielding factor (None: the directive's 0.4, 1.0: outdoors).

RAIN_MM_H = 1.0   # intensity of the rain, when it rains


def wind_statistics(category, speed, from_directions, rain_fraction=None):
    """
    One weather bin per direction, with equal probabilities. Without `rain_fraction` the bins sum to a probability
    of 1 (what the dispersion needs); with it they sum to `rain_fraction` and carry a rainfall intensity
    (what the washout needs: it is not normalised, "it does not rain all the time").
    """
    probability = (1.0 if rain_fraction is None else rain_fraction) / len(from_directions)
    table = pd.DataFrame({
        'weather_category': category,
        'wind_speed': speed,
        'wind_direction': np.asarray(from_directions) % 360,   # direction the wind comes FROM
        'probability': probability,
    })
    if rain_fraction is not None:
        table['rainfall_intensity'] = RAIN_MM_H
    return WeatherStats.from_dataframe(table, check_deg=False, randomize_wind_direction_within_bins=False)


def chernobyl_dose(plume='documented', height_scale=1.0, rain_fraction=0.0, k_s=None):
    """Returns a long dataframe: one row per village, day, nuclide, pathway and age group."""
    rows = []
    for v, village in villages.iterrows():
        xv, yv = x[v:v + 1], y[v:v + 1]
        for d, (category, speed, height, documented) in enumerate(days):
            # share of this day's release that happens before the evacuation
            released_before_evacuation = float(np.clip(village.t_evacuation - day_start[d], 0, 1))
            if released_before_evacuation == 0.0:
                continue

            # directions the wind comes from (= direction the plume travels to + 180 degrees)
            if plume == 'early_aligned' and d < 2:
                from_directions = np.array([village.bearing_deg + 180.])
            elif plume == 'uniform' or documented is None:
                from_directions = np.arange(5, 360, 10.)
            else:
                to_bearing, spread = documented
                from_directions = to_bearing + 180. + np.linspace(-spread, spread, 7)

            # a source without plume rise (W_0 = 0, no buildings): the effective height is the given height
            source = SourceParams(H_a=height * height_scale * m, H_G=0 * m, D_i=1 * m, D_e=1 * m, W_0=0 * m / s)
            averager = WeatherAverager(orography=orography, source_params=source)
            averager.correct_wind_speed_depending_on_height = False   # the speeds are given at release height
            chi_k, chi_ks = averager.average_dispersion(xv, yv, wind_statistics(category, speed, from_directions))
            rain = wind_statistics(category, speed, from_directions, rain_fraction) if rain_fraction > 0 else None
            time_of_flight = village.distance_km * 1e3 / speed * s

            for nuclide, total in total_release.items():
                Q = total * daily_share[d] * released_before_evacuation * Bq
                inhalation = dose_calculator.inhalation(nuclide, Q, chi_k, T_fz=time_of_flight)
                submersion = dose_calculator.submersion(nuclide, Q, chi_ks, T_fz=time_of_flight, k_s=k_s)

                # ground-shine, from the dry deposition until the evacuation
                coefficients = dose_calculator._get_nuclide_row_table41(nuclide)
                decay = coefficients['lambda_a-1'] / year
                e_bs = coefficients['external_ground_eBS_(Sv/a)/(Bq/m2)'] * (Sv / year) / (Bq / m**2)
                deposited = Q * deposition_model.fallout_factor(chi_k, nuclide=nuclide)         # dry, Bq / m2
                if rain is not None:    # plus washout, which scavenges the whole column: it does not need ground contact
                    deposited = deposited + Q * averager.average_washout(xv, yv, rain, nuclide)
                residence = max(village.t_evacuation - (day_start[d] + 0.5), 0) * day   # deposition at mid-day
                shielding = dose_calculator.K_S_BS if k_s is None else k_s
                ground_shine = deposited * shielding * e_bs * (
                    dose_calculator.f_f * dose_calculator.integral_exp(decay + dose_calculator.lam_f, residence)
                    + dose_calculator.f_s * dose_calculator.integral_exp(decay + dose_calculator.lam_s, residence)
                )

                for age in dose_calculator.age_groups:
                    rows += [
                        (village['name'], d, nuclide, 'inhalation', age, np.nan_to_num(inhalation[age].to_value(Sv))[0]),
                        (village['name'], d, nuclide, 'submersion', age, submersion[age].to_value(Sv)[0]),
                        (village['name'], d, nuclide, 'ground_shine', age, ground_shine.to_value(Sv)[0]),
                    ]
                rows.append((village['name'], d, nuclide, 'deposition_Bq_m2', 'n/a', deposited.to_value(Bq / m**2)[0]))
    return pd.DataFrame(rows, columns=['village', 'day', 'nuclide', 'pathway', 'age_group', 'value'])


def cs137_deposition_kBq_m2(df):
    """Total modelled Cs-137 deposition at each village until its evacuation."""
    cs = df[(df.pathway == 'deposition_Bq_m2') & (df.nuclide == 'Cs-137')]
    return cs.groupby('village').value.sum().reindex(villages.name) / 1e3

Step 1: calibrating the rain on the measured fallout

The Cs-137 ground deposition measured in 1986-87 exceeded about 1500 kBq/m² within the 30 km zone (up to about 5000 kBq/m² in some villages), and was 600-650 kBq/m² at 15 km from the reactor. Let us compare with the model, with the reported release heights and the documented plume directions, for several shares of rainy time:

shares = [0.0, 0.002, 0.01, 0.05]
check = pd.DataFrame({
    f"rain {share:.1%} of the time": cs137_deposition_kBq_m2(chernobyl_dose('documented', rain_fraction=share))
    for share in shares
})
check.index.name = 'Cs-137 deposition [kBq/m2]'
print(check.round(0).to_string())

which prints (version 0.5.0):

                            rain 0.0% of the time  rain 0.2% of the time  rain 1.0% of the time  rain 5.0% of the time
Cs-137 deposition [kBq/m2]                                                                                            
Pripyat                                       0.0                  140.0                  702.0                 3511.0
Yaniv                                         0.0                  154.0                  769.0                 3844.0
Kopachi                                       0.0                  193.0                  961.0                 4802.0
Chernobyl                                   335.0                  390.0                  608.0                 1700.0

Without rain, the plume of the model does not touch the ground within a few km (1000 m, then 600 m high): there is almost no fallout close to the reactor, whereas the measurements show a heavy one. Rain changes this completely, because washout scavenges the whole plume column and is independent of the release height: with about 1 % of the time at 1 mm/h (a quarter of an hour per day), the modelled deposition is of the order of the measured values everywhere. We use this "calibrated" rain from now on. It is only one of several things that can produce this deposition (a lower release, fallout of heavy particles near the reactor, ...): the data used do not identify the cause. For Pripyat and Yaniv, which sit at the edge of the documented plume, the position of the plume matters as well, see the plumes below.

Step 2: the doses

rain_fraction = 0.01
results = {
    'documented': chernobyl_dose('documented', rain_fraction=rain_fraction),
    'early_aligned': chernobyl_dose('early_aligned', rain_fraction=rain_fraction),
    'uniform': chernobyl_dose('uniform', rain_fraction=rain_fraction),
}

def dose_table(df, age_group):
    """Dose in mSv, one row per village, one column per pathway, plus the total."""
    d = df[(df.age_group == age_group)]
    table = d.pivot_table(index='village', columns='pathway', values='value', aggfunc='sum').reindex(villages.name)
    table['total'] = table.sum(axis=1)
    return table[['inhalation', 'submersion', 'ground_shine', 'total']] * 1e3

for plume, df in results.items():
    print(f"\n--- adult, plume: {plume} (mSv)")
    print(dose_table(df, 'adult').round(3).to_string())

which prints (version 0.5.0):

--- adult, plume: documented (mSv)
pathway    inhalation  submersion  ground_shine   total
name                                                   
Pripyat         0.000       0.303         0.152   0.455
Yaniv           0.000       0.317         0.167   0.484
Kopachi         0.038       0.665         0.424   1.126
Chernobyl      25.549       0.464         1.193  27.205

--- adult, plume: early_aligned (mSv)
pathway    inhalation  submersion  ground_shine    total
name                                                    
Pripyat         0.010       2.557         0.517    3.083
Yaniv           0.001       2.760         0.582    3.344
Kopachi         0.162       3.852         5.551    9.565
Chernobyl     106.487       3.120         8.614  118.222

--- adult, plume: uniform (mSv)
pathway    inhalation  submersion  ground_shine   total
name                                                   
Pripyat         0.000       0.122         0.031   0.153
Yaniv           0.000       0.127         0.035   0.162
Kopachi         0.016       0.383         0.395   0.794
Chernobyl      12.807       0.582         0.691  14.081

The same, summed over pathways, for the three age groups:

totals = pd.DataFrame({
    (plume, age): dose_table(df, age)['total']
    for plume, df in results.items() for age in ['adult', 'child_10_year', 'toddler_1_year']
})
print((totals).round(2).to_string())

which prints (version 0.5.0):

          documented                              early_aligned                              uniform                             
               adult child_10_year toddler_1_year         adult child_10_year toddler_1_year   adult child_10_year toddler_1_year
name                                                                                                                             
Pripyat         0.45          0.45           0.45          3.08          3.09           3.09    0.15          0.15           0.15
Yaniv           0.48          0.48           0.48          3.34          3.34           3.34    0.16          0.16           0.16
Kopachi         1.13          1.15           1.17          9.56          9.68           9.74    0.79          0.80           0.81
Chernobyl      27.21         44.85          54.71        118.22        191.75         232.86   14.08         22.92          27.87

Which nuclide matters, for the worst-case plume (adult, mSv)?

by_nuclide = results['early_aligned'].query("age_group == 'adult'").pivot_table(
    index=['village', 'nuclide'], columns='pathway', values='value', aggfunc='sum') * 1e3
print(by_nuclide.round(3).to_string())

which prints (version 0.5.0):

pathway            ground_shine  inhalation  submersion
village   nuclide                                      
Chernobyl Cs-137          0.444       1.071       0.147
          I-131           8.170     105.416       2.184
          Xe-133          0.000       0.000       0.789
Kopachi   Cs-137          0.660       0.002       0.181
          I-131           4.891       0.160       2.695
          Xe-133          0.000       0.000       0.976
Pripyat   Cs-137          0.051       0.000       0.120
          I-131           0.466       0.010       1.789
          Xe-133          0.000       0.000       0.648
Yaniv     Cs-137          0.058       0.000       0.130
          I-131           0.525       0.001       1.931
          Xe-133          0.000       0.000       0.699

Sensitivity to the shielding

The directive's shielding factor (0.4, a person spending most of the time indoors) applies to submersion and ground-shine. People outdoors are not shielded (k_s = 1):

outdoors = chernobyl_dose('early_aligned', rain_fraction=rain_fraction, k_s=1.0)
print(dose_table(outdoors, 'adult').round(3).to_string())

which prints (version 0.5.0):

pathway    inhalation  submersion  ground_shine    total
name                                                    
Pripyat         0.010       6.392         1.292    7.694
Yaniv           0.001       6.901         1.456    8.358
Kopachi         0.162       9.630        13.877   23.669
Chernobyl     106.487       7.801        21.534  135.823

Limitations

  • The plume position is a major uncertainty. The three plumes differ by a factor 8 to 20 for the same village (Pripyat, adult: 0.45, 3.1 and 0.15 mSv). The 1986 plume was lofted, moved, and its trajectory over the first days is not captured by a constant wind direction per day.
  • The rain is a calibrated parameter, not a measurement: nothing was found on the rain near the reactor, and the 1 % of rainy time was tuned on loosely sourced values of the measured Cs-137 deposition (which a lower release or heavy particles could also explain). Ground-shine scales roughly linearly with it: try other values of rain_fraction in Step 1.
  • Only three nuclides are modelled. Te-132/I-132, Ba/La-140, Zr/Nb-95, Ru and others, which dominated the early ground-shine, are missing, so ground-shine is a lower bound. Their activities are not given above and they would need their own release and profile.
  • The temporal profile of I-131 is assumed for all nuclides, the I-131 chemical form is the package default.
  • Constant weather during each day, transport time not resolved, no depletion of the plume by the washout and the fallout, no sheltering of the inhaled air.
  • Doses are only accumulated until evacuation, in the settlement itself. Ingestion is not computed (nobody is assumed to eat local food), neither are thyroid doses, which dominated the reported health impact: the package computes effective doses.
  • Inhalation is the weakest pathway. It comes from I-131 in the air at ground level. For an elevated release the model has almost none within a few km (Pripyat: 0.01 mSv even with the plume on the village), while the plume is reported to have enveloped Pripyat and the reported thyroid doses of its evacuees are large (a mean of about 0.17 Gy). At 15 km (Chernobyl town) the plume has reached the ground and inhalation dominates. The inhalation doses close to the reactor are therefore likely underestimated, and this is not fixed by the rain, which only acts on the deposition.
  • For comparison, the reported mean external dose of the evacuees of Pripyat and the 30 km zone is about 17 mSv. The model's submersion plus ground-shine is 0.5 to 3 mSv for Pripyat (adult): 5 to 30 times lower, as expected from the missing nuclides and the uncertain plume.

Sources

Values were gathered from web search summaries of the following references (some of the original pages could not be opened); verify against the originals before reuse.