Skip to content

Long term computation example

Here we compute the dose experienced yearly by groups situated near a cluster of chimneys emitting activity.

Parameters

Let us first go over the required input parameters. For each set of parameter, the package routine we will see below accepts either a path to a CSV file containing the values, or a pandas dataframe with the same information.

Position and characteristics of the chimneys

Prepare a CSV file following the template:

name,h,W,D_i,D_e,h_building,lat,long
bag1,20,10.2,1.75,1.8,15,46.93026329420960,7.420436181732923
bag2,18,7.2,1.65,1.7,15,46.92953305334195,7.421908702163402

That is, for each chimney:

  • its name
  • its height (h) in meters
  • its exhaust velocity (W) in meters / second
  • its interior and exterior diameters (D_i and D_e) in meters.
  • a representative height characterising the surrounding buildings (h_building) in meters.
  • its latitude and longitude (WGS84) in degrees.

Position of the exposed groups

Please format these in another CSV file, here is an example:

name,lat,long
mcdo,46.92626527964869,7.41639267991673
park,46.92834332443608,7.41955510310637

Each group is labelled by its name, latitude and longitude.

An optional height column (metres above local ground, outside G-14) places a group above the ground, e.g. a window on a façade. Inhalation and submersion (also the photon point-kernel observer) are then evaluated at that height; deposition and ingestion stay at ground level. Missing values mean ground level. The PKI surrogate does not support elevated observers.

Emissions

In yet another CSV file, prepare the emissions in this way: 

nuclide,bag1,bag2
H-3_HTO,1e6,1e7
Be-7,1e3,3e2
C-11,1,1
C-14,1,1
Ar-41,1,1
Sr-90,1,1
That is, a 2D table of activities (in becquerels) indexed by

  • the names of the nuclide (rows),
  • the names of the chimneys (columns).

Weather statistics

This is the most complex part, as there are many perfectly valid ways of parametrising weather characteristics.

Each weather situation is labelled by

  • its turbulence category, which we will call weather_category in the package,
  • its intensity of rainfall, rainfall_intensity (given in mm/h),
  • its wind speed, wind_speed (given in m/s),
  • the direction FROM which the wind is coming, wind_direction (given in degrees, north = 0°, east = 90°).

Moreover, each weather situation (let us call it a weather bin) is weighted by a probability, probability.

So, for the 4-parameter weather statistic, an input CSV file would look like this: 

weather_category,rainfall_intensity,wind_speed,wind_direction,probability
F,0.5,1.5,25,0.0208920964316648
D,1.5,3.5,205,0.0172028049725891
F,1.,1.5,15,0.0153606844254326
F,0.1,1.5,215,0.0130703440055958
A,0.2,1.5,205,0.0128506531384229
F,0.2,0.5,5,0.0127546962079336
D,0.5,4.5,205,0.0123064763352533
  (...)

The 3-parameter weather statistic described in the ENSI-G14 would look similar, simply minus the rainfall_intensity column.

Note: 3-parameter statistic does not mean dry weather. We use it internally to compute the dispersion \(\chi\).

On the other hand, 4-parameter statistic does not mean wet weather only. However we use this one internally to compute the washout factor \(W\). So, weather bins with 0 rainfall will indeed contribute nothing in our case. Thus, you can remove all the weather bins with no rainfall from your 4-parameter weather statistic. See the weather clustering tool in averaging/weather_measurement_clustering.py to cluster "raw" weather measurements.

In summary: I would prepare a 3-parameter weather statistic containing all the weather bins, simply omitting the rainfall_intensity information for clarity. Next, I would select all the weather bins that do contain rainfall to constitute the 4-parameter weather statistic.

Orography

Use the Shuttle Radar Topographic Mission (download from opentopography) data to get an elevation map of the zone of interest around your chimneys. Do not make it too big, but make sure to include all your chimneys and exposed groups. Download as an ASC file.

If an exposed group finds itself outside the footprint of your downloaded topographic file, the interpolator will throw a cryptic error: ValueError: One of the requested xi is out of bounds in dimension 1 -- Just so you know. Download a bigger topographic footprint to fix it.

Running the calculator

As always, use your favourite Python editor or a Jupyter notebook to script these commands.

Let us do the data loading first. We assume having prepared all the required data described above in a subdirectory called data.

The chimneys

import numpy as np
import pandas as pd

chimneys = pd.read_csv('data/chimneys.csv')

# by the way, say we are missing the D_e info: per the guideline we can set D_e = D_i
# chimneys['D_e'] = chimneys['D_i']

The exposed groups and emissions

exposed_groups = pd.read_csv('data/exposed_groups.csv')
emissions = pd.read_csv('data/emission.csv')

The weather

weather = pd.read_csv('data/weather_statistics.csv')

# say we have one file with all weather bins, split it in two
# as described above
weather_3_parameter = weather.drop('rainfall_intensity', axis=1)
weather_4_parameter = weather[weather['rainfall_intensity'] > 0.0]

The actual run

from ensi_g14_calculator.long_term_scenarios import LongTermEvaluator
from astropy.units import m

# instantiate the calculator, give it the path to the
# the database that will contain the results 
# (the said database will be created),
# and the path to your orography ASC file.
LE = LongTermEvaluator('dose_database.sqlite3', 'data/orography.asc') 

LE.run(chimneys, exposed_groups, 
       weather_3_parameter,weather_4_parameter, 
       emissions, summer_total_precipitation=0.3*m)
This last parameter, summer_total_precipitation, is only required if your emissions contain tritium. Here I set it to 300mm.

Optional photon point-kernel integral (PKI) for submersion

This package can also evaluate submersion dose through point-kernel integrals, which is useful for photon submersion doses when the observer is less than 200m away from the stack or in general far from the plume with respect to its size. Replace the LE.run call above with this extended call to compute PKI alongside the regular submersion for a given nuclide, here Ar-41:

```python notest from ensi_g14_calculator.dose import PointKernelConfig

LE.run( chimneys, exposed_groups, weather_3_parameter, weather_4_parameter, emissions, summer_total_precipitation=0.3 * m, pki_nuclides=["Ar-41"], pki_config=PointKernelConfig( # this is optional, default works good enough. radial_order=12, polar_order=12, azimuthal_order=16, ), ) 

!!! warning "Elevation map extent"
    PKI samples the plume out to several photon mean free paths around each
    observer (up to a few km for high-energy lines such as Ar-41). Every
    sampled point needs elevation data, so the `ASC` file must extend well
    beyond your exposed groups, otherwise you get the `out of bounds` error
    described in the orography section above.

PKI results are stored separately as
`model="pki_photon"` and `is_reference=False`. Query them explicitly when
comparing against the ENSI-G14 reference results:

```python
from ensi_g14_calculator.dose.dose_database import DoseDatabase

with DoseDatabase("dose_database.sqlite3") as db:
        ensi_total = db.get_dose_breakdown(
            group='mcdo', age_group='adult', chimney='bag', nuclide='Cs-137', pathway='submersion', model='ensi_g14'
        )
        pki_total = db.get_dose_breakdown(
           group='mcdo', age_group='adult', chimney='bag', nuclide='Cs-137',pathway='submersion', model='pki_photon', 
           include_non_reference=True  # pass this else the PKI dose gets ignored. 
        )

The long-term runner evaluates PKI at the exposed-group coordinates and uses local terrain for observer height. For arbitrary 3-D points that are not exposed groups, use the direct PKI example.

Note finally one why PKI: the ENSI model was fitted on PKI simulations by CERN, and attempts at correcting the semi-infinite cloud submersion model by introducing a cloud-shape and observer-position dependent correction factor. It works well in most situations, but breaks down severely below 200m, and whenever the observer is far enough from the plume (vertically or horizontally).

PKI should the be preferred, e.g. in urban settings where the distances are small, and the geometry generally very different from what the ENSI model was thought for. Final note, while we can evaluate the ENSI model on a grid and get an idea of what the submersion dose looks like on the 2D plane, PKI is more expensive and is evaluated only at given points.

Visualisation of the results

Dose inspection

Of course, you can conduct further analyses using the produced sqlite3 database with the tool of your choice.

We also provide a webapp that loads the produced sqlite3 file and shows a dynamic breakdown of the dose. See the files and the readme at https://github.com/duxfrederic/ensi_g14_calculator/tree/main/ensi_g14_dashboard.

Mapping dispersion and elevation

You can get a quick overview of the calculated dispersion with the LongTermEvaluator class, after having run LongTermEvaluator.run at least once:

LE.plot_dispersion('bag1')  # 'bag1' is here a chimney name
(example output) png

There you can inspect the order of magnitudes of your dispersion, as well as its coherence with respect to the topography and to the wind distribution.

Likewise, use the deposition plot to inspect the balance between washout and fallout: 

LE.plot_deposition('bag1', 'H-3_HTO')
(example output) png

LE.plot_deposition('bag1', 'Cs-137')
(example output) png

Miscellaneous

The wind measurement height is taken to be 10m by default, as it is fairly standard. You can change it after having initialised your LongTermEvaluator:

from astropy.units import m
LE = LongTermEvaluator('dose_database.sqlite3', 'data/orography.asc') 
LE.wind_measurement_height_zm = 20 * m

Then run the calculation as above.