Skip to content

Near field: a short stack on a hospital roof

This page uses the near-field module for a typical cyclotron situation: an exhaust a few metres above a hospital roof, with the people of interest tens of metres away. It goes from the annual dispersion factor to an annual emission limit for F-18 that keeps the dose below 0.3 mSv. Every code block on this page is executed by the test-suite.

Not part of ENSI-G14

G-14 is built for tall stacks and receptors beyond about 200 m. The near-field module keeps the G-14 chain and only changes the spreading of the building-trapped plume, its behaviour at the source, the plume rise before 200 m, the path from the stack exit to nearby receptors, and optionally the wind speed. The equations, their sources (Ramsdell & Fosmire 1995 and the NRC code ARCON96 built on it; IAEA SRS-19) and the direction of each bias are documented in prototyping/near_field_dispersion.tex.

Source, building and weather

The stack is 3 m above a 20 m roof. The building that carries it is described once by its footprint; the module computes the face seen by the wind for every weather bin.

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

from ensi_g14_calculator.averaging.weather_averager import SourceParams, WeatherAverager, WeatherStats
from ensi_g14_calculator.near_field import Building, NearFieldWeatherAverager
from ensi_g14_calculator.orography import Orography
from ensi_g14_calculator.misc import Sv

source = SourceParams(
    H_a=23 * m,        # stack top above the ground at the stack base
    H_G=20 * m,        # building height: G-14 trapped fraction, and the building height of the wake model
    D_i=0.3 * m,
    D_e=0.35 * m,
    W_0=6.0 * m / s,
)
building = Building(
    length=60 * m,           # side along orientation_deg
    width=25 * m,
    orientation_deg=30.0,    # azimuth of the long side, clockwise from north
)                            # height=None: taken from source.H_G

orography = Orography(flat=True)
# meteorological convention: 0 deg is north, clockwise, direction the wind comes FROM
weather = WeatherStats.from_dataframe(pd.read_csv("data/weather_statistics.csv"))

Roofs

Here the ground is flat and receptors stand on it. Receptors on roofs are handled with building footprints on top of the orography (section Receptors on roofs below).

Annual dispersion close to the stack

Three averagers: plain G-14 (extrapolated below 200 m, which G-14 does not foresee), the near-field model, and the near-field model with all wind speeds reduced by 30 % for a built-up area. The maximum over azimuth is taken on rings around the stack, at 1.5 m above ground.

Receptors are evaluated with average_dispersion. For them, the near-field model measures the path of the building-trapped plume from the stack exit, \(S = \sqrt{x^2 + (z - H_a)^2}\), as ARCON96 does (NF-12 in the .tex; stretched_string=False turns it off). average_dispersion_at_points evaluates the plume as a cloud, for the point-kernel integral, and does not use it.

averagers = {
    "G-14": WeatherAverager(orography, source),
    "near field": NearFieldWeatherAverager(orography, source, building),
    "near field, f_U = 0.7": NearFieldWeatherAverager(orography, source, building, wind_speed_factor=0.7),
}

radii = np.array([20.0, 30.0, 50.0, 100.0, 200.0])
azimuths = np.deg2rad(np.arange(0.0, 360.0, 10.0))
R, AZ = np.meshgrid(radii, azimuths, indexing="ij")
X, Y = R * np.sin(AZ) * m, R * np.cos(AZ) * m    # east, north
Z = np.full(R.shape, 1.5) * m

chi_rings = {name: av.average_dispersion(X, Y, weather, height=Z)[0].to_value(s / m**3)
             for name, av in averagers.items()}
table = pd.DataFrame({name: chi.max(axis=1) for name, chi in chi_rings.items()},
                     index=pd.Index(radii, name="r (m)"))
print(table.to_string(float_format="{:.2e}".format))

          G-14  near field  near field, f_U = 0.7
r (m)
20.0  4.31e-04    2.83e-04               2.15e-04
30.0  2.34e-04    1.80e-04               1.38e-04
50.0  1.04e-04    9.07e-05               7.08e-05
100.0 3.88e-05    3.59e-05               3.29e-05
200.0 1.66e-05    1.58e-05               1.79e-05

(WeatherStats.from_dataframe scatters the wind directions within their bins, so the last digits change from run to run.)

Three things to notice:

  • With the broad Vogt parameters of G-14, the wake changes the annual average by a few percent only. The near-field module mainly makes the calculation below 200 m well-defined: it caps the plume at the source (where the G-14 power laws diverge) and evaluates the plume rise at the actual distance.
  • Within 50 m, the ground-level values are 10 to 35 % lower than with G-14: the building-trapped plume leaves the stack 23 m up, so a person on the street 20 m away is reached after about 29 m, not 20 m (NF-12).
  • A lower wind speed does not always raise the concentration. Through the G-14 rules, a smaller \(U\) means a larger \(W_0/U\): more plume rise and a smaller building-trapped fraction \(G_t\). Close to this stack that wins over the \(1/U\) dilution. The wind-speed factor is a sensitivity parameter, not a safety factor.

Receptor height matters as well, for example windows on a neighbouring façade 40 m from the stack, in the most exposed direction (the scenario runners take height too):

near_field = averagers["near field"]
k = int(np.argmax(chi_rings["near field"][1]))      # most exposed azimuth at 30 m
direction = azimuths[k]
heights = np.array([1.5, 5.0, 10.0, 15.0, 20.0, 25.0]) * m
x_w = np.full(heights.shape, 40.0) * np.sin(direction) * m
y_w = np.full(heights.shape, 40.0) * np.cos(direction) * m
chi_windows = near_field.average_dispersion(x_w, y_w, weather, height=heights)[0]
for h, c in zip(heights, chi_windows):
    print(f"z = {h.value:4.1f} m: chi = {c.to_value(s / m**3):.2e} s/m3")

z =  1.5 m: chi = 1.24e-04 s/m3
z =  5.0 m: chi = 1.30e-04 s/m3
z = 10.0 m: chi = 1.41e-04 s/m3
z = 15.0 m: chi = 1.60e-04 s/m3
z = 20.0 m: chi = 2.07e-04 s/m3
z = 25.0 m: chi = 3.98e-04 s/m3

The concentration grows with height: the building-trapped part is reached sooner the closer the window is to the stack exit, and at 25 m the window also sits in the elevated part of the plume. There it is three times the ground-level value at the same distance and higher than anywhere on the ground 20 m from the stack. Upper-floor windows of neighbouring buildings at stack height are candidate critical receptors and should be on the list. (With stretched_string=False, the trapped part is a ground-level plume a few metres thick, and a window 10 m up would get a fifth of the street value: too little inside the cavity of the building, which is mixed over about its height.)

Receptors on roofs

Observer heights place a receptor in the air above flat ground, which suits windows on façades. On a roof, the building-trapped plume flows over the roof itself, and the cloud seen by the point-kernel integral should not extend inside the building. Buildings can be added on top of the orography as rectangles (centre, length, width, orientation as for the wake Building, and the absolute elevation of the roof). Inside a footprint the terrain is the roof; the terrain grid is only interpolated outside the footprints. The reference elevation stays the bare ground at the stack, so H_a remains the stack height above the street, as G-14 uses it (wind speed at release height, trapped fraction, plume rise).

from ensi_g14_calculator.orography import BuildingFootprint

footprints = [
    # the hospital carrying the stack, 20 m high (centre at the stack here; latitude/longitude also work)
    BuildingFootprint(x=0 * m, y=0 * m, length=60 * m, width=25 * m, orientation_deg=30.0, roof_elevation=20 * m),
    # a taller neighbour 45 m east
    BuildingFootprint(x=45 * m, y=0 * m, length=20 * m, width=15 * m, orientation_deg=0.0, roof_elevation=26 * m),
]
with_roofs = NearFieldWeatherAverager(Orography(flat=True, buildings=footprints), source, building)

# 1.5 m above: the hospital roof 20 m along its long axis, the neighbour's roof, the street 60 m east
x_r = np.array([20.0 * np.sin(np.deg2rad(30.0)), 45.0, 60.0]) * m
y_r = np.array([20.0 * np.cos(np.deg2rad(30.0)), 0.0, 0.0]) * m
print("terrain:", with_roofs.orography.elevation_xy(x_r, y_r))
roofs = with_roofs.average_dispersion(x_r, y_r, weather, height=1.5 * m)[0]
heights_only = near_field.average_dispersion(x_r, y_r, weather, height=np.array([21.5, 27.5, 1.5]) * m)[0]
for label, a, b in zip(["hospital roof", "neighbour roof", "street"], roofs, heights_only):
    print(f"{label:15s} roofs as terrain {a.to_value(s / m**3):.2e}   heights only {b.to_value(s / m**3):.2e} s/m3")

terrain: [20. 26.  0.] m
hospital roof   roofs as terrain 5.11e-04   heights only 6.20e-04 s/m3
neighbour roof  roofs as terrain 1.90e-04   heights only 1.38e-04 s/m3
street          roofs as terrain 2.00e-05   heights only 2.00e-05 s/m3

For \(\chi_K\) the two agree to within a factor 1.4: the receptor is at the same height above the stack base either way, so the trapped part of the plume reaches it after the same slant distance (NF-12). What differs is the elevated part, which the G-14 orography correction lowers over a roof. The footprints matter most for the exact PKI: the trapped part of the cloud then lies on the roofs, and cloud points inside a footprint are below the terrain and are masked, so photons from inside the buildings are not counted. Steps in the terrain are a crude representation of the flow around buildings: the plume follows the roofs instead of flowing around them, and the G-14 orography correction (made for smooth hills) lowers the elevated plume over tall neighbours. Both are conservative on roofs.

Dose per becquerel and emission limit for F-18

For one receptor, the annual dose per emitted becquerel combines

  • inhalation from \(\chi_K\) at the receptor,
  • submersion from the exact photon point-kernel integral over the near-field cloud. The G-14 polynomial \(\chi_{K,S}\) is not valid below 200 m, and the fast PKI surrogate is not reliable within about 40 m of a short stack (next section). For comparison, the semi-infinite-cloud value with the local \(\chi_K\) is printed too.

Decay during transport

The G-14 dose functions decay the parent over a fixed flight time T_fz, 8.4 minutes for long-term scenarios. Close to a cyclotron stack the travel time is \(x/U\), tens of seconds: with the default, the dose of O-15 (\(T_{1/2}\) = 2 min) would be underestimated by a factor 17, that of N-13 by a factor 2. Below, decay is applied inside the plume instead: every weather bin and plume component is weighted by \(e^{-\lambda x/U}\) at its own travel time (time_weight), daughters by their in-growth, and the dose functions receive the result without T_fz.

from ensi_g14_calculator.dispersion.plume_decay import decayed_dispersion, parent_survival
from ensi_g14_calculator.dose import PointKernelIntegrator
from ensi_g14_calculator.dose.long_term_dose import LongTermDoseCalculator

integrator = PointKernelIntegrator()
dose_calculator = LongTermDoseCalculator()


def near_field_cloud(x, y, z):
    return near_field.average_dispersion_at_points(x, y, z, weather)


def annual_dose_per_bq(nuclide, x, y, z=1.5 * m):
    """Annual dose (Sv/Bq) at one receptor, for the most exposed age group, decay in the plume."""
    Q = 1.0 * Bq
    # chi of the parent and of its daughters, each decayed at its own travel time x / U
    chi = decayed_dispersion(
        lambda w: near_field.average_dispersion(x, y, weather, time_weight=w, height=z)[0],
        nuclide, dose_calculator,
    )
    # PKI over a cloud in which every point is decayed by its own travel time
    decay_constant, _ = dose_calculator.decay_chain(nuclide)
    decayed_cloud = near_field.physical_dispersion_field(weather, time_weight=parent_survival(decay_constant))
    pki_coefficient = integrator.dose_coefficient(nuclide, (x, y, z), decayed_cloud)

    inhalation = dose_calculator.inhalation(nuclide, Q, chi)
    submersion = dose_calculator.submersion_pki(nuclide, Q, pki_coefficient, T_fz=0 * year, k_s=1.0)
    semi_infinite = dose_calculator.submersion(nuclide, Q, chi, k_s=1.0)
    worst_group = max(inhalation, key=lambda g: inhalation[g] + submersion[g])
    return {
        "inhalation": inhalation[worst_group].to_value(Sv),
        "submersion (PKI)": submersion[worst_group].to_value(Sv),
        "submersion (semi-infinite)": semi_infinite[worst_group].to_value(Sv),
        "total": (inhalation[worst_group] + submersion[worst_group]).to_value(Sv),
    }


rows = {}
for i, r in enumerate(radii[:4]):
    k = int(np.argmax(chi_rings["near field"][i]))   # most exposed azimuth on this ring
    rows[r] = annual_dose_per_bq("F-18", X[i, k], Y[i, k])
doses = pd.DataFrame(rows).T.rename_axis("r (m)")
doses["limit (TBq/a)"] = 0.3e-3 / doses["total"] / 1e12
print(doses.to_string(float_format="{:.2e}".format))

       inhalation  submersion (PKI)  submersion (semi-infinite)    total  limit (TBq/a)
r (m)
20.0     5.61e-18          1.41e-18                    1.38e-17 7.02e-18       4.27e+01
30.0     3.57e-18          1.14e-18                    8.78e-18 4.72e-18       6.36e+01
50.0     1.80e-18          8.16e-19                    4.41e-18 2.61e-18       1.15e+02
100.0    7.10e-19          5.15e-19                    1.74e-18 1.22e-18       2.45e+02

The submersion dose from the finite cloud is 3 to 10 times lower than the semi-infinite-cloud value with the local concentration: close to a short stack the cloud is small compared with the photon mean free path in air (about 90 m at 511 keV). The receptor with the smallest limit sets the emission limit.

For the short-lived nuclides the treatment of decay matters. O-15 and N-13 have no inhalation coefficients in G-14 (they are counted by submersion only); here is the PKI submersion of O-15 at the most exposed point at 30 m, without decay, with decay in the plume, and with the G-14 flight time:

x30, y30 = X[1, int(np.argmax(chi_rings["near field"][1]))], Y[1, int(np.argmax(chi_rings["near field"][1]))]
lambda_o15, _ = dose_calculator.decay_chain("O-15")
no_decay = integrator.dose_coefficient("O-15", (x30, y30, 1.5 * m), near_field_cloud)
in_plume = integrator.dose_coefficient(
    "O-15", (x30, y30, 1.5 * m),
    near_field.physical_dispersion_field(weather, time_weight=parent_survival(lambda_o15)),
)
g14_flight_time = dose_calculator.submersion_pki("O-15", 1 * Bq, no_decay, k_s=1.0)["adult"]
print(f"no decay        {dose_calculator.submersion_pki('O-15', 1 * Bq, no_decay, T_fz=0 * year, k_s=1.0)['adult']:.2e}")
print(f"decay in plume  {dose_calculator.submersion_pki('O-15', 1 * Bq, in_plume, T_fz=0 * year, k_s=1.0)['adult']:.2e}")
print(f"G-14 T_fz       {g14_flight_time:.2e}")

no decay        1.18e-18 Sv
decay in plume  1.04e-18 Sv
G-14 T_fz       6.75e-20 Sv

Over the few tens of metres that matter, O-15 loses about 10 % of its activity; the G-14 flight time would remove 94 % of it.

Choices made on this page

k_s=1.0: no shielding by buildings for submersion (G-14 uses 0.4 for long-term scenarios). Occupancy is 100 %. The weather statistics are the package's test data, not a real site. Repeat the calculation for every nuclide of the source term (C-11, N-13, O-15, Ar-41 from activated air) and sum the fractions of the limit.

Submersion map with the PKI surrogate

The exact PKI takes about a second per receptor with a full weather statistics: fine for a list of candidate receptors, too slow for a map. The neural surrogate evaluates thousands of points at once. It sees a plume only through its widths at two distances, so the near-field plume is handed to it as the power law through those two widths (NF-9 in the .tex). Against the exact PKI this adds at most 5 % from 50 m on. Within about 40 m the surrogate itself is unreliable for short stacks (from 0.4 to 14 times the exact value, already for the plain G-14 plume): there, use the exact PKI as in the previous section.

near_field.submersion_model = "pki_surrogate"

map_radii = np.array([50.0, 75.0, 100.0, 150.0, 200.0, 300.0])
map_azimuths = np.deg2rad(np.arange(0.0, 360.0, 5.0))
RM, AM = np.meshgrid(map_radii, map_azimuths, indexing="ij")
XM, YM = RM * np.sin(AM) * m, RM * np.cos(AM) * m
coefficient_map = near_field.average_submersion_surrogate_dose_coefficient(XM, YM, weather, "F-18")

# check the surrogate against the exact PKI at the most exposed point of two rings
for i in (0, 2):
    k = int(np.argmax(coefficient_map[i]))
    exact = integrator.dose_coefficient("F-18", (XM[i, k], YM[i, k], 0 * m), near_field_cloud)
    print(f"r = {map_radii[i]:5.0f} m: surrogate {coefficient_map[i, k].to_value(Sv / Bq):.2e}, "
          f"exact {exact.to_value(Sv / Bq):.2e} Sv/Bq")

r =    50 m: surrogate 7.92e-19, exact 7.35e-19 Sv/Bq
r =   100 m: surrogate 5.53e-19, exact 5.26e-19 Sv/Bq

coefficient_map (Sv per emitted Bq) can be plotted directly, e.g. with matplotlib.pyplot.pcolormesh on XM, YM; multiply by the annual emission and the shielding factor to get annual doses. For decay during transport, pass time_weight=parent_survival(decay_constant) to average_submersion_surrogate_dose_coefficient (applied at the receptor, at \(x/U\)), then use submersion_pki(..., T_fz=0 * year).