Short term scenario example

In this example, we will compute the dose deposited onto a grid around the chimney for a (prospective) short term scenario. We will use a flat orography.

Let us start by importing the useful components from the calculator:

from ensi_g14_calculator.orography import Orography
from ensi_g14_calculator.averaging.weather_averager import SourceParams
from ensi_g14_calculator.short_term_scenarios.runner import ScenarioEvaluator, get_max_report
from ensi_g14_calculator.dose import PointKernelConfig
from ensi_g14_calculator.misc.plots import plot_dose_summary

We will also need some units from astropy, as all the quantities inside the calculators must be given with units. Additionally, we import numpy to define the grid later on.

from astropy.units import Bq, m, s
import numpy as np

Parameters

Let us define the parameters of our short term scenario. We first need a vector of nuclides being released, define them in this way:

nuclides = {
    'H-3_HTO': 1e10 * Bq, 
    'Cs-137': 1e7 * Bq, 
    'Ar-41': 1e10 * Bq
}

Next, we need to define the parameters of the chimney from which the activity is emitted:

source = SourceParams(
    H_a = 75 * m,
    H_G = 15 * m,
    D_i = 1.2 * m,
    D_e = 1.2 * m,
    W_0 = 10.0 * m / s,
)

Where

H_a is the height of the chimney (total of the stack),
H_G is the height of the tallest adjacent buildings or structures, to calculate potential ground effects,
D_i and D_e are the internal and external diameters of the chimney, respectively,
W_0 is the exhaust speed of the chimney,
and the optional z_m is the reference height at which the wind speed (set internally by the calculator to 1m/s as per the directive) was measured. It defaults to 10 * m.

We will also instantiate an orography object, which takes care of the geometry:

orography = Orography(flat=True)

Please look at the Orography API documentation should you need add an elevation map.

Defining the grid on which we will compute the dispersion and dose

For this example we will evaluate the dose on a square grid.

extent = 3000.0 * m
x = np.linspace(-extent, extent, 151)
y = np.linspace(-extent, extent, 151)
x, y = np.meshgrid(x, y)  # a square grid going from -extent to +extent

Note that, if you wanted to, you could instead evaluate the dose on a line, or at a single position, e.g. x = 300.0 * m, y = 0.0 * m.

Using the scenario evaluator

The ScenarioEvaluator class handles the definition of the scenario: short term, in two weather phases, with the activity spread over 24 hours.

evaluator = ScenarioEvaluator(source, orography)

Now we need to try each weather category to find the maximum dose. The evaluator will compute dispersion, deposition and ultimately dose maps for each. For the first 8 hours, a 2mm/h rain is assumed; then between 8 and 24h, 1mm/hour (per the directive, section 4.1.1.2.b.1).

Here, the direction of the wind can be chosen arbitrarily. But in a real scenario, with the actual position of the exposed group being known, it will have to be chosen in the way that maximises dose for the exposed group. (Usually, wind blowing in their direction from the release point.)

results = evaluator.evaluate_short_term_dose(nuclides_bq=nuclides,
                                             x=x, y=y, wind_direction=-90)

Optional photon point-kernel comparison

PKI was implemented for dealing with submersion dose in distances shorter than 200m. It accounts for photon dose only. For nuclear power plants, the ENSI-G14 model remains the reference calculation and the default dose totals exclude PKI records. Add the following optional arguments to the call above to evaluate Cs-137 at the direct and residence hotspots selected by G-14:

results_with_pki = evaluator.evaluate_short_term_dose(
    nuclides_bq=nuclides,
    x=x,
    y=y,
    wind_direction=-90,
    pki_nuclides=["Cs-137"],
    pki_config=PointKernelConfig(
        radial_order=12,
        polar_order=12,
        azimuthal_order=16,
    ),
    pki_hotspot_age_group="adult",
    min_radius_direct=200 * m,
    min_radius_residence=500 * m,
)

The two radius values are the G-14 spatial exclusions used to select the hotspots; they are not PKI integration settings. PKI runs after that selection and does not change the ENSI-G14 grids or hotspot locations. Compare the photon result at the selected direct point with:

weather = results_with_pki.worst_weather
direct_pki_dose = results_with_pki.hotspot_report[weather][
    "submersion_comparison_direct_0_8h"
]["Cs-137"]["pki_photon"]["adult"]

Point records for both models are also available with results_with_pki.get_point_doses(). The PKI records have model="pki_photon" and is_reference=False; use reference_only=False or an explicit model filter when comparing them. For arbitrary 3-D coordinates rather than G-14 hotspots, call evaluator.evaluate_short_term_pki(...); see the direct PKI example.

The above evaluate_short_term_dose admits a few optional parameters that you can fine tune if needed:

T_exp (time, with units), the exposure time after release, either 1 * year (default) or 50 * year.
release_date (string), the date of the release. Default 2 months before the end of harvest season to be conservative: '08-01'.

The resulting results object holds dispersion grids and dose grids, organised by nuclide, weather category, interval, pathway and age group (query it with sum_doses, below).

Inspecting the results

You can get the dose grid calculated for any pathway, nuclide, age group, weather category, and time interval, for example:

dose_inhalation_hto_first8h = results.sum_doses(weather='A', nuclide='H-3_HTO', interval='0-8h', 
                                                pathway='inhalation', age_group='adult')
print(dose_inhalation_hto_first8h.max())  # or you can plot it, or save it to a file and analyse it in your favourite tool ...

The procedure described in the directive is a bit complicated, so there is a helper for this: get_max_report. It performs the following:

Find the maximum dose for the combination of inhalation and submersion in the first 8h (interval "0-8h"), no closer than 200m to the release point. This is also done in the absence of shielding (\(k_s = 0\)), as it depicts the representative person witnessing the event outdoors and close to the release point.
Find the maximum dose for the combination of all ingestion pathways (and ground shine) in the medium term (interval "medium_term"), no closer than 500m to the release point. This location defines the living area of the representative person.
Compute the submersion and inhalation at the location of the living area of the representative person, for the interval "8-16h", this time with shielding (\(k_s\) according to directive).

Here is how to use the get_max_report utility:

worst_weather, reports = get_max_report(results, 'adult', x, y, 200.0*m, 500.0*m)

direct_dose = reports[worst_weather]['dose_direct_0_24h']  # inhalation + submersion
medium_term_dose = reports[worst_weather]['dose_medium_term']  # food + ground shine
direct_pos = reports[worst_weather]['pos_direct']
residence_pos = reports[worst_weather]['pos_residence']

print(f"The worst short term dose occurs in weather {worst_weather} "
      f"with a combined dose of {direct_dose:.2e} within the first 24h; then {medium_term_dose:.2e} from food/ground shine.\n"
      f"Evaluated at hotspot (first 8h) at y={y[direct_pos]}, x={x[direct_pos]}; "
      f"then, hours 8h->T_exp evaluated at y={y[residence_pos]}, x={x[residence_pos]}.")

The output looks like:

The worst short term dose occurs in weather D with a combined dose of 2.19e-07 Sv within the first 24h; then 1.26e-04 Sv from food/ground shine.
Evaluated at hotspot (first 8h) at y=0.0 m, x=660.0 m; then, hours 8h->T_exp evaluated at y=0.0 m, x=500.0 m.

We can also plot the dose map to inspect the contribution of the various pathways:

plot_dose_summary(results, worst_weather, 'adult', '0-8h', x, y, 
                  min_radius_direct=200*m, min_radius_residence=500*m)

png

Only the pathways corresponding to the given interval (0-8h) are displayed.

We can query other intervals:

plot_dose_summary(results, worst_weather, 'adult', 'medium_term', x, y, 
                  min_radius_direct=200*m, min_radius_residence=500*m)

png

The pathways are updated accordingly. We can also pass nuclide= to inspect the dose attributable to a specific nuclide.

Keep in mind that all "grid of values" carry units. Sv for dose grids, s/m³ for dispersion grids. You can strip the unit with .value, e.g.,

grid = results.sum_doses(weather='A', interval='0-8h', age_group='adult').value