blpd.chem

  1"""
  2Chemistry and emissions calculations
  3on LPD (Lagrangian particle dispersion) model output.
  4
  5>>> ds0 = blpd.Model().run().to_xr()  # LPD output
  6>>> ds = blpd.chem.calc_relative_levels_fixed_oxidants(ds0)
  7"""
  8
  9import numpy as np
 10import xarray as xr
 11from scipy import stats
 12
 13from .utils import auto_bins_xy
 14from .utils import cache
 15from .utils import calc_t_out
 16from .utils import load_p
 17
 18__all__ = (
 19    "calc_areal_emission_rates_canola",
 20    "calc_gridded_conc_canola",
 21    "calc_relative_levels_fixed_oxidants",
 22    "load_canola_species_data",
 23    "load_canola_flower_basal_emissions",
 24    "load_default_chemical_species_data",
 25)
 26
 27
 28@cache
 29def load_default_chemical_species_data():
 30    """Chemical species data used by default in `calc_relative_levels_fixed_oxidants`.
 31    Returns
 32    -------
 33    dict
 34        *keys*: ASCII chemical species names;
 35        *values*: another dict with *values*:
 36        `'display_name'`, `'kO3'`, `'kOH'`, `'kNO3'`
 37    """
 38    # key, display name string, kO3, kOH, kNO3
 39    # for now must have space after `,` !
 40    s = """
 41apinene, "α-pinene", 8.1e-17, 5.3e-11, 6.2e-12
 42bocimene, "β-ocimene", 5.4e-16, 2.52e-10, 2.2e-11
 43bpinene, "β-pinene", 2.4e-17, 7.8e-11, 2.5e-12
 44carene, "3-carene", 3.8e-17, 8.7e-11, 9.1e-12
 45cineole, "1,8-cineole", 6.0e-20, 1.0e-11, 1.7e-16
 46farnesene, "α-farnesene", 1.0e-15, 3.2e-10, 5.5e-11
 47limonene, "d-limonene", 2.5e-16, 1.6e-10, 1.3e-11
 48linalool, "linalool", 4.3e-16, 1.6e-10, 1.1e-11
 49myrcene, "β-myrcene", 4.7e-16, 2.1e-10, 1.3e-11
 50sabinene, "sabinene", 9.0e-17, 1.2e-10, 1.0e-11
 51thujene, "α-thujene", 4.4e-16, 8.7e-11, 1.1e-11
 52    """.strip()
 53    # TODO: store in separate file that makes it easier to add more (and other data). YAML?
 54
 55    d = {}
 56    for line in s.split("\n"):
 57        parts = [p.strip() for p in line.split(", ")]
 58        dl = {
 59            "display_name": eval(parts[1]),
 60            "kO3": float(parts[2]),
 61            "kOH": float(parts[3]),
 62            "kNO3": float(parts[4]),
 63        }
 64        key = parts[0]
 65        d[key] = dl
 66
 67    return d
 68
 69
 70# TODO: for relative levels, should always be 0--1 (or 0--100% optionally), so should remove `fv_0_default` and coordinate with gridded calc
 71def calc_relative_levels_fixed_oxidants(
 72    ds,
 73    species=None,
 74    *,
 75    t_out=None,
 76    p_overrides=None,
 77    fv_0_default: float = 100.0  # percentage mode by default
 78):
 79    """Calculate relative levels (rt. values at particle release)
 80    using fixed (over space and time) oxidant levels.
 81
 82    The fixed oxidant levels greatly simplify the problem—fractional chemical destruction only depends on time-since-release.
 83    Since we release particles in a known non-random way (fixed release rate),
 84    we can derive time-since-release after the LPD run has finished,
 85    without even needing the trajectory data.
 86    (This assumes the source strengths remain constant as well.)
 87
 88    Parameters
 89    ----------
 90    ds : xr.Dataset
 91        LPD model output generated by `blpd.model.Model.to_xarray`.
 92    species : dict
 93        *keys*: the species to calculate levels for;
 94        *values*: a `dict` for each species that must at least have the rate constants
 95        ``'kO3'``, ``'kOH'``, ``'kNO3'``
 96        but can also have ``'fv_0'``, the in-particle concentration at release.
 97        By default, the standard dataset `CHEMICAL_SPECIES_DATA` is used.
 98    fv_0_default
 99        Set to `1.0` instead for fractional mode (default is percentage mode)
100        or pick a different value.
101        Used for species in `species` that don't provide `'fv_0'` in their dict.
102
103    Returns
104    -------
105    xr.Dataset
106    """
107    if species is None:
108        species = load_default_chemical_species_data()
109
110    p = load_p(ds)
111
112    if p_overrides is not None:
113        p.update(p_overrides)
114
115    # unpack needed model options/params
116    Np_tot = p["Np_tot"]
117
118    if t_out is None:
119        t_out = calc_t_out(p)
120
121    # (molec cm^-3)
122    conc_O3 = p["conc_oxidants"]["O3"]
123    conc_OH = p["conc_oxidants"]["OH"]
124    conc_NO3 = p["conc_oxidants"]["NO3"]
125
126    # Save position dataset to merge with
127    ds0 = ds.copy()
128
129    ip = np.arange(Np_tot)  # for now
130    spcs = sorted(species.keys())  # species keys
131    n_spc = len(spcs)
132
133    def d0():
134        return np.empty((ip.size, n_spc))  # to be filled; or create d0 and use .copy()
135
136    ds = xr.Dataset(
137        coords={
138            "ip": ("ip", ip, {"long_name": "Lagrangian particle index"}),
139            "spc": ("spc", spcs, {"long_name": "chemical species"}),
140        },
141        data_vars={
142            "f_r": (("ip", "spc"), d0(), {"long_name": "fraction remaining"}),
143            "f_d_o3": (("ip", "spc"), d0(), {"long_name": "fraction destroyed by O3"}),
144            "f_d_oh": (("ip", "spc"), d0(), {"long_name": "fraction destroyed by OH"}),
145            "f_d_no3": (("ip", "spc"), d0(), {"long_name": "fraction destroyed by NO3"}),
146        },
147        attrs={
148            "conc_O3": conc_O3,  # molec cm-3
149            "conc_OH": conc_OH,
150            "conc_NO3": conc_NO3,
151            "p_json": ds0.attrs["p_json"],
152        },
153    )
154    for spc, d_spc in species.items():
155        conc0_val = p["fv_0"].get(spc, fv_0_default)
156
157        kO3 = d_spc["kO3"]
158        kOH = d_spc["kOH"]
159        kNO3 = d_spc["kNO3"]
160
161        # k_sum = kO3 + kOH + kNO3
162        # k_inv_sum = 1/kO3 + 1/kOH + 1/kNO3
163        kc_sum = conc_O3 * kO3 + conc_OH * kOH + conc_NO3 * kNO3
164        kc_O3 = kO3 * conc_O3
165        kc_OH = kOH * conc_OH
166        kc_NO3 = kNO3 * conc_NO3
167
168        # fraction remaining
169        f_remaining = (
170            1.0 * np.exp(-kc_O3 * t_out) * np.exp(-kc_OH * t_out) * np.exp(-kc_NO3 * t_out)
171        )
172
173        # fraction destroyed
174        # breakdown by oxidant depends on the rate const and oxidant levels
175        f_destroyed = 1 - f_remaining  # total
176        f_d_O3 = kc_O3 / kc_sum * f_destroyed
177        f_d_OH = kc_OH / kc_sum * f_destroyed
178        f_d_NO3 = kc_NO3 / kc_sum * f_destroyed
179        # note: transforming from f_r to these is linear
180        # so could be done after binning, but we might want to plot these particle-wise
181
182        f_r_spc = np.full((Np_tot,), conc0_val) * f_remaining
183
184        # update dataset values
185        ds["f_r"].loc[:, spc] = f_r_spc
186        ds["f_d_o3"].loc[:, spc] = f_d_O3
187        ds["f_d_oh"].loc[:, spc] = f_d_OH
188        ds["f_d_no3"].loc[:, spc] = f_d_NO3
189
190    ds_ret = xr.merge([ds, ds0[["x", "y", "z"]]])
191    ds_ret.attrs.update(ds.attrs)
192
193    # Add display names
194    display_names = [d["display_name"] for d in species.values()]
195    ds_ret["display_name"] = (
196        "spc",
197        display_names,
198        {"long_name": "better chemical species names (UTF-8)"},
199    )
200
201    return ds_ret
202
203
204@cache
205def load_canola_species_data():
206    """Canola species data as `dict` of `dict`s.
207
208    Returns
209    -------
210    dict
211        *keys*: ASCII chemical species names;
212        *values*: another dict with *values*:
213        `'MW'`, `'kO3'`, `'kOH'`, `'kNO3'`, `'OH_yield'`, `'HCHO_yield'`, `'display_name'`
214    """
215
216    def try_float(v):
217        try:
218            return float(v)
219        except ValueError:
220            return v
221
222    fields = ["key", "MW", "kO3", "kOH", "kNO3", "OH_yield", "HCHO_yield", "display_name"]
223    # ,mw,ko3,koh,kno3,oh_yield,hcho_yield,display_name
224    s = """
225apinene,136,8.09e-17,5.33e-11,6.16e-12,0.83,0.19,α-pinene
226bpinene,136,2.35e-17,7.81e-11,2.51e-12,0.35,0.45,β-pinene
227carene,136,3.8e-17,8.7e-11,9.1e-12,0.86,0.21,3-carene
228cineole,154,6e-20,1e-11,1.7e-16,0.01,0,"1,8-cineole"
229farnesene,204,1e-15,3.2e-10,5.5e-11,0.6,0.15,α-farnesene
230limonene,136,2.5e-16,1.61e-10,1.3e-11,0.67,0.19,d-limonene
231linalool,154,4.3e-16,1.6e-10,1.1e-11,0.72,0.36,linalool
232myrcene,154,4.7e-16,2.1e-10,1.27e-11,0.63,0.74,β-myrcene
233sabinene,136,9e-17,1.17e-10,1e-11,0.33,0.5,sabinene
234thujene,136,4.4e-16,8.7e-11,1.1e-11,0.69,0.36,α-thujene
235    """.strip()
236
237    import csv
238
239    lines = s.splitlines()
240    d = {}
241    for row in csv.DictReader(lines, fieldnames=fields):
242        dl = {k: try_float(v) for k, v in row.items() if k != "key"}
243        d[row["key"]] = dl
244
245    return d
246
247
248@cache
249def load_canola_flower_basal_emissions():
250    """Speciated basal emission rate data for canola plant."""
251
252    # basal emissions in ng per floret-hour (ng floret^-1 hr^-1)
253    # from Table 1 of Jakobsen et al. (1994) (the values for light conditions)
254    # doi: https://doi.org/10.1016/S0031-9422(00)90341-8
255    emiss_ng = {
256        "myrcene": 12.3,
257        "limonene": 7.2,
258        "sabinene": 6.5,
259        "farnesene": 5.0,
260        "apinene": 4.2,
261        "linalool": 3.0,
262        "carene": 1.9,
263        "bpinene": 2.0,
264        "thujene": 1.9,
265        "cineole": 1.7,
266    }
267
268    # convert to ug/s units
269    emiss_ug_s = {}
270    for spc, e_ng in emiss_ng.items():
271        # calculate emission in ug (micrograms) per second
272        e_ug_s = e_ng / 1000 / (60 * 60)
273        emiss_ug_s[spc] = e_ug_s
274
275    # reference temperature
276    # near the end of the paper they say that the
277    # Sampling was performed at 15 deg
278    T_S_C = 15.0  # deg. C
279    T_S_K = T_S_C + 273.15  # Kelvin
280
281    return {
282        "emiss_ug_per_floret-s": emiss_ug_s,
283        "T_S_K": T_S_K,
284    }
285
286
287def calc_areal_emission_rates_canola(
288    plant_density: float = 30.0,
289    num_flowers_avg: float = 25.0,
290    T: float = 298.0,
291    beta: float = 0.06,
292):
293    """Calculate areal emission rates (m$^{-2}$ s$^{-1}$) for canola floral volatiles.
294
295    Parameters
296    ----------
297    plant_density
298        Average number of plants per m$^2$.
299    num_flowers_avg
300        Average number of flowers in the plants.
301    T
302        Temperature (K).
303
304        Temperature sensitivity parameter for emissions (K$^{-1}$).
305
306    Returns
307    -------
308    dict
309        *values*: species ASCII; *keys*: dict of areal emissions in different units
310    """
311
312    # load basal data
313    basal = load_canola_flower_basal_emissions()
314    E_si_ug = basal["emiss_ug_per_floret-s"]
315    T_S = basal["T_S_K"]
316
317    # other needed
318    N_A = 6.02214076e23  # also in scipy.constants
319
320    # load canola species additional data
321    data = load_canola_species_data()
322
323    # loop through species
324    E_i = {}
325    for spc, e_si in E_si_ug.items():
326        # ug per floret-s at temperature T
327        e_i = e_si * np.exp(beta * (T - T_S))
328
329        # compute areal emissions
330
331        # ug per s per m^2
332        e_i_areal = e_i * num_flowers_avg * plant_density
333
334        # mol " "
335        mw = data[spc]["MW"]  # molecular weight (g/mol)
336        e_i_areal_mol = e_i_areal / 1e6 / mw  # note ug -> g
337
338        # molec " "
339        e_i_areal_molec = e_i_areal_mol * N_A
340
341        E_i[spc] = {
342            "emiss_areal_ug": e_i_areal,
343            "emiss_areal_mol": e_i_areal_mol,
344            "emiss_areal_molec": e_i_areal_molec,
345        }
346
347    return E_i
348
349
350def calc_gridded_conc_canola(ds):
351    """Some calculations for canola floral volatiles
352    using the fixed oxidants chem.
353
354    Parameters
355    ----------
356    ds : xr.Dataset
357        LPD model output generated by `blpd.model.Model.to_xarray`.
358    """
359    # ! 2-D for now (not real 3-D/volume conc.)
360    p = load_p(ds)
361
362    # 1. Determine floral volatile levels per particle
363    emiss_rates = calc_areal_emission_rates_canola()
364    dt = p["dt"]
365    dNp_per_dt_per_source = p["dNp_per_dt_per_source"]
366    Np_per_sec_per_source = dNp_per_dt_per_source / dt
367    fv_0 = {}  # floral volatile initial amounts in particle
368    for spc, d in emiss_rates.items():
369        # assume 1 m^2 source area
370        # TODO: source's effective area for emissions calculation can be an input somewhere?
371        # and that it does not overlap with other sources
372        molec_per_p = d["emiss_areal_molec"] / Np_per_sec_per_source
373        fv_0[spc] = molec_per_p
374
375    # 2. Compute particle-relative levels
376    #    decreased due to chemical destruction by oxidation
377    canola = load_canola_species_data()
378    ds_fo = calc_relative_levels_fixed_oxidants(ds, species=canola, fv_0_default=1.0)
379    # ^ can also pass fv_0
380    spcs_fo = ds_fo["spc"].values
381
382    # 3. grid
383    X, Y = ds.x.values, ds.y.values
384    dx = X[1] - X[0]
385    dy = Y[1] - Y[0]
386    pos = (X, Y)
387    bins = auto_bins_xy(pos)
388
389    # 4. Calculate concentrations from particle concentration field
390    #    and speciated fv_0
391    res = {}  # store results here
392    ret = None  # for first run we don't have a result yet
393    stats_to_calc = ["sum", "count", "mean", "median", "std"]
394    for spc in spcs_fo:
395        f_r_spc = ds_fo["f_r"].sel(spc=spc).values  # fraction remaining, by particle
396        f_d_o3_spc = ds_fo["f_d_o3"].sel(spc=spc).values  # fraction destroyed by O3
397        # Calculate each of the stats
398        rets = {}
399        variables = {"f_r_spc": f_r_spc, "f_d_o3": f_d_o3_spc}
400        for vn, values in variables.items():
401            rets_vn = {}
402            for stat in stats_to_calc:
403                ret = stats.binned_statistic_dd(
404                    np.column_stack(pos),  # binning by
405                    values,  # calculating the statistic on
406                    statistic=stat,
407                    bins=bins,
408                    binned_statistic_result=ret,
409                )
410                rets_vn[stat] = ret
411            rets[vn] = rets_vn
412
413        # Gridded particle count
414        # should be same for all, don't really need to re-calc this way...
415        n_p_g = rets["f_r_spc"]["count"].statistic
416
417        # Calculate conc. by summing the relative concs. of particles
418        # and scaling by the particle initial # of molec
419        f_r_sum = rets["f_r_spc"]["sum"].statistic
420        fv_spc_g = f_r_sum * fv_0[spc]
421
422        # Link molec destroyed by O3 to OH and HCHO yields
423        # via the ozonolysis oxidation pathway
424        f_d_o3_sum = rets["f_d_o3"]["sum"].statistic
425        oh_yield_oz_spc = canola[spc]["OH_yield"]
426        hcho_yield_oz_spc = canola[spc]["HCHO_yield"]
427        oh_y_tot = f_d_o3_sum * fv_0[spc] * fv_0[spc] * oh_yield_oz_spc
428        hcho_y_tot = f_d_o3_sum * fv_0[spc] * fv_0[spc] * hcho_yield_oz_spc
429
430        # xr.Datset per species then combine later??
431
432        # collect
433        res[spc] = {
434            "particle_count": n_p_g,
435            "molec_count": fv_spc_g,
436            "conc_molec_areal": fv_spc_g / (dx * dy),
437            "oh_yield_count": oh_y_tot,
438            "oh_yield_areal": oh_y_tot / (dx * dy),
439            "hcho_yield_count": hcho_y_tot,
440            "hcho_yield_areal": hcho_y_tot / (dx * dy),
441        }
442
443    # 5. Specify grid centers
444    x, y = ret.bin_edges[0], ret.bin_edges[1]
445    xc = x[:-1] + 0.5 * np.diff(x)
446    yc = y[:-1] + 0.5 * np.diff(y)
447
448    # Create xr.Dataset
449
450    # Specify units and long names for the above-calculated variables
451    units_long_names = {
452        "particle_count": ("", "Lagrangian particle count"),
453        "molec_count": ("molec", "molecule count"),
454        "conc_molec_areal": ("molec m^-2", "areal conc. (molec.)"),
455        "oh_yield_count": ("radical", "OH yield from ozonolysis"),
456        "oh_yield_areal": ("radical m^-2", "areal OH yield from ozonolysis"),
457        "hcho_yield_count": ("molec", "HCHO yield from ozonolysis"),
458        "hcho_yield_areal": ("molec m^-2", "areal HCHO yield from ozonolysis"),
459    }
460    units = {k: v[0] for k, v in units_long_names.items()}
461    long_names = {k: v[1] for k, v in units_long_names.items()}
462
463    # Aggregate: species as a third dim
464    spc_sorted = sorted(spc for spc in res.keys())
465    n_spc = len(spc_sorted)
466    varnames = res[spc_sorted[0]].keys()
467    res_agg = {}
468    for vn in varnames:
469        data = np.zeros((n_spc, xc.size, yc.size))
470        for i, (spc, data_spc) in enumerate(res.items()):
471            data[i, :, :] = data_spc[vn]
472        res_agg[vn] = data
473
474    ds = xr.Dataset(
475        coords={
476            "xe": ("xe", x, {"units": "m", "long_name": "x (bin edge)"}),
477            "ye": ("ye", y, {"units": "m", "long_name": "y (bin edge)"}),
478            "x": ("x", xc, {"units": "m", "long_name": "x (bin center)"}),
479            "y": ("y", yc, {"units": "m", "long_name": "y (bin center)"}),
480            "spc": ("spc", spc_sorted, {"long_name": "chemical species"}),
481        },
482        data_vars={
483            vn: (
484                ("spc", "x", "y"),
485                values,
486                {
487                    "units": units[vn],
488                    "long_name": long_names[vn],
489                },
490            )
491            for vn, values in res_agg.items()
492        },
493        attrs={
494            "t_tot_s": p["t_tot"],
495        },
496    )
497
498    # Fix particle count (doesn't vary with species)
499    ds["particle_count"] = ds["particle_count"].sel(spc="apinene")
500
501    # Add display names
502    display_names = [canola[spc]["display_name"] for spc in spc_sorted]
503    ds["display_name"] = (
504        "spc",
505        display_names,
506        {"long_name": "better chemical species names (UTF-8)"},
507    )
508
509    return ds
510
511
512# Define chem calculation options for the main model class
513chem_calc_options = {"fixed_oxidants": calc_relative_levels_fixed_oxidants}