pykpp Tutorial

goals: Predict ozone at a monitor using HYSPLIT and pykpp authors: Barron Henderson, Robert Nedbor-Gross, Qian Shu

# Overview

  1. Create an idealized environment (Knote et al., 2015)

  2. Create a simple trajectory model

  • Chemistry - CB05

  • Emissions from NEI

  • Deposition from Knote et al., 2015

  • Dilution following Knote et al., 2015

Assumes pykpp is installed

Knote et al. 2015

This example application is based on the Knote et al. (2015) simulation of a typical summer day with CB05. The example here uses the “Summer” PBL and temperature profile. At present, the example uses the a relatively “clean” initial environment (NOx(t=0) = 1ppb).

Trajectory model

We’ll use a KPP language to implement the model. Each piece will be saved as a string and then combined to produce an active model.

  • prepare an air pracel and model environment,

  • select gas-phase chemical reactions,

  • add emissions as chemical reactions,

  • add deposition as chemical reactions

import time
import numpy as np
from pykpp.mech import Mech
from io import StringIO
import pandas as pd
import matplotlib.pyplot as plt
from pykpp.updaters import func_updater


infile = StringIO("""
#INLINE PY_INIT
TEMP = 298.15
P = 99600.
t = TSTART = 0 * 3600.
TEND = TSTART + 3600. * 24
DT = 60.
MONITOR_DT = 3600
StartDate = 'datetime(2010, 7, 14)'
Latitude_Degrees = 40.
Longitude_Degrees = 0.00E+00
KDIL = 0.
PBLH_OLD = 250.
vd_CO = 0.03
vd_SO2 = 1.8
vd_NO = 0.016
vd_NO2 = 0.1
vd_NO3 = 0.1
vd_N2O5 = 4
vd_HONO = 2.2
vd_HNO3 = 4
vd_HNO4 = 4
vd_O3 = 0.4
vd_H2O2 = 0.5
vd_CH3OOH = 0.1
molps_to_molpcm2ps = Avogadro / 1200000**2

add_world_updater(func_updater(Update_M, incr = 360., verbose=False))
add_world_updater(func_updater(Update_THETA, incr = 360., verbose=False))
add_world_updater(interpolated_from_csv(
    'summerenv.tsv', 'time', incr = 360., delimiter = '\\t', verbose=False
))
add_world_updater(interpolated_from_csv(
    'mean_emis.csv', 'time', incr = 360, verbose=False
))
add_world_updater(func_updater(
    Monitor, incr = 7200., allowforce = False, verbose=False
))
#ENDINLINE

#INITVALUES
CFACTOR = P * Avogadro / R / TEMP * centi **3 * nano {ppb-to-molecules/cm3}
ALL_SPEC=1e-32*CFACTOR;
TOTALNOx=1.
M=1e9
O2=.21*M
N2=.79*M
H2O=0.01*M
CH4=1750
CO=100
O3=30.
NO = 0.1 * TOTALNOx
NO2 = 0.9 * TOTALNOx
SO2 = 1
N2O = 320
{B = 210.; what to do about B?}

#MONITOR t/3600.; THETA; PBLH; TEMP; O3;
#LOOKAT THETA; t/3600; PBLH; TEMP; O1D; OH; HO2; O3; NO; NO2; eNOx; TUV_J(6, THETA); emis_NO; emis_NO * molps_to_molpcm2ps / PBLH / 100 / CFACTOR * 3600;

#include cb05cl.eqn
<EMISALD2> EMISSION = ALD2 : emis_ALD2 * molps_to_molpcm2ps / PBLH / 100 ;
<EMISALDX> EMISSION = ALDX : emis_ALDX * molps_to_molpcm2ps / PBLH / 100 ;
<EMISCH4> EMISSION = CH4 : emis_CH4 * molps_to_molpcm2ps / PBLH / 100 ;
<EMISCL2> EMISSION = CL2 : emis_CL2 * molps_to_molpcm2ps / PBLH / 100 ;
<EMISCO> EMISSION = CO : emis_CO * molps_to_molpcm2ps / PBLH / 100 ;
<EMISETH> EMISSION = ETH : emis_ETH * molps_to_molpcm2ps / PBLH / 100 ;
<EMISETHA> EMISSION = ETHA : emis_ETHA * molps_to_molpcm2ps / PBLH / 100 ;
<EMISETOH> EMISSION = ETOH : emis_ETOH * molps_to_molpcm2ps / PBLH / 100 ;
<EMISFORM> EMISSION = FORM : emis_FORM * molps_to_molpcm2ps / PBLH / 100 ;
<EMISHONO> EMISSION = HONO : emis_HONO * molps_to_molpcm2ps / PBLH / 100 ;
<EMISIOLE> EMISSION = IOLE : emis_IOLE * molps_to_molpcm2ps / PBLH / 100 ;
<EMISISOP> EMISSION = ISOP : emis_ISOP * molps_to_molpcm2ps / PBLH / 100 ;
<EMISMEOH> EMISSION = MEOH : emis_MEOH * molps_to_molpcm2ps / PBLH / 100 ;
<EMISNH3> EMISSION = NH3 : emis_NH3 * molps_to_molpcm2ps / PBLH / 100 ;
<EMISNO> EMISSION = NO + eNOx : emis_NO * molps_to_molpcm2ps / PBLH / 100 ;
<EMISNO2> EMISSION = NO2 + eNOx : emis_NO2 * molps_to_molpcm2ps / PBLH / 100 ;
<EMISOLE> EMISSION = OLE : emis_OLE * molps_to_molpcm2ps / PBLH / 100 ;
<EMISPAR> EMISSION = PAR : emis_PAR * molps_to_molpcm2ps / PBLH / 100 ;
<EMISSESQ> EMISSION = SESQ : emis_SESQ * molps_to_molpcm2ps / PBLH / 100 ;
<EMISSO2> EMISSION = SO2 : emis_SO2 * molps_to_molpcm2ps / PBLH / 100 ;
<EMISSULF> EMISSION = SULF : emis_SULF * molps_to_molpcm2ps / PBLH / 100 ;
<EMISTERP> EMISSION = TERP : emis_TERP * molps_to_molpcm2ps / PBLH / 100 ;
<EMISTOL> EMISSION = TOL : emis_TOL * molps_to_molpcm2ps / PBLH / 100 ;
<EMISXYL> EMISSION = XYL : emis_XYL * molps_to_molpcm2ps / PBLH / 100 ;
<DDEPCO> CO = DUMMY : vd_CO / PBLH / 100 ;
<DDEPSO2> SO2 = DUMMY : vd_SO2 / PBLH / 100 ;
<DDEPNO> NO = DUMMY : vd_NO / PBLH / 100 ;
<DDEPNO2> NO2 = DUMMY : vd_NO2 / PBLH / 100 ;
<DDEPNO3> NO3 = DUMMY : vd_NO3 / PBLH / 100 ;
<DDEPN2O5> N2O5 = DUMMY : vd_N2O5 / PBLH / 100 ;
<DDEPHONO> HONO = DUMMY : vd_HONO / PBLH / 100 ;
<DDEPHNO3> HNO3 = DUMMY : vd_HNO3 / PBLH / 100 ;
<DDEPHNO4> PNA = DUMMY : vd_HNO4 / PBLH / 100 ;
<DDEPO3> O3 = DUMMY : vd_O3 / PBLH / 100 ;
<DDEPH2O2> H2O2 = DUMMY : vd_H2O2 / PBLH / 100 ;
<DDEPCH3OOH> MEPX = DUMMY : vd_CH3OOH / PBLH / 100 ;

""")

keywords = [
    'DUMMY', 'EMISSION', 'BNZHRXN', 'BNZNRXN', 'ISOPRXN', 'SESQRXN', 'SULRXN',
    'TOLHRXN', 'TOLNRXN', 'TRPRXN', 'XYLHRXN', 'XYLNRXN'
]

mech = Mech(
    infile, mechname='summercb05', incr=360, add_default_funcs=False,
    keywords=keywords
)

# Create time start/stop matrix
nhour = 24.
start_end_ts = np.linspace(
    0, nhour, int(nhour) * 12 + 1
).repeat(2, 0)[1:-1].reshape(-1, 2) * 3600

# Capture initial values as "background"
# concentrations
mech.world['ybkg'] = mech.get_y(mech.parsed_world)

# Create specific tolerance values for radical species
# great care should be taken with O1D.
# note that CMAQ EBI solver has no absolute tolerance and
# uses 1 rtol for radicals and rxn counters.
#
atol = np.array([
    10 if (spc in ('O', 'O1D', 'HCO3', 'NTR') or spc[-3:] == 'RXN') else 1e-3
    for spc in mech.allspcs
])
rtol = np.array([
    0.1 if (spc in ('O', 'O1D', 'HCO3', 'NTR') or spc[-3:] == 'RXN') else 1e-5
    for spc in mech.allspcs
])

# Start timing the process
runstart = time.time()

# Update world for all interpolated values
mech.world['t'] = start_end_ts.min()
mech.update_y_from_world()
mech.Update_World(forceupdate=True)
Included /home/runner/work/pykpp/pykpp/pykpp/models/cb05cl.eqn
/home/runner/work/pykpp/pykpp/pykpp/parse.py:391: UserWarning: Ignoring DUMMY
  warn('Ignoring %s' % spc)
/home/runner/work/pykpp/pykpp/pykpp/parse.py:391: UserWarning: Ignoring ISOPRXN
  warn('Ignoring %s' % spc)
/home/runner/work/pykpp/pykpp/pykpp/parse.py:391: UserWarning: Ignoring TRPRXN
  warn('Ignoring %s' % spc)
/home/runner/work/pykpp/pykpp/pykpp/parse.py:391: UserWarning: Ignoring SULRXN
  warn('Ignoring %s' % spc)
/home/runner/work/pykpp/pykpp/pykpp/parse.py:391: UserWarning: Ignoring TOLNRXN
  warn('Ignoring %s' % spc)
/home/runner/work/pykpp/pykpp/pykpp/parse.py:391: UserWarning: Ignoring TOLHRXN
  warn('Ignoring %s' % spc)
/home/runner/work/pykpp/pykpp/pykpp/parse.py:391: UserWarning: Ignoring XYLNRXN
  warn('Ignoring %s' % spc)
/home/runner/work/pykpp/pykpp/pykpp/parse.py:391: UserWarning: Ignoring XYLHRXN
  warn('Ignoring %s' % spc)
/home/runner/work/pykpp/pykpp/pykpp/parse.py:391: UserWarning: Ignoring BNZNRXN
  warn('Ignoring %s' % spc)
/home/runner/work/pykpp/pykpp/pykpp/parse.py:391: UserWarning: Ignoring BNZHRXN
  warn('Ignoring %s' % spc)
/home/runner/work/pykpp/pykpp/pykpp/parse.py:391: UserWarning: Ignoring SESQRXN
  warn('Ignoring %s' % spc)
/home/runner/work/pykpp/pykpp/pykpp/parse.py:385: UserWarning: Ignoring EMISSION
  warn('Ignoring %s' % spc)
Species: 65
Reactions: 223
{t:0,t/3600.:0,THETA:1.2E+02,PBLH:2.5E+02,TEMP:2.9E+02,O3:30}
{t:0,t/3600.:0,THETA:1.2E+02,PBLH:2.5E+02,TEMP:2.9E+02,O3:30}

Define a function to dilute with PBL

  • In a trajectory model, we need a PBL dilution

  • It should only dilute when PBL goes up

  • When PBL collapses, it should not concentrate

  • This will be called by the mechanism in between integrations

def UpdatePBL(mech, world):
    """
    Defining updater for PBL rise
    """
    if 'y' not in world:
        return
    PBLH_OLD = world['PBLH_OLD']
    ybkg = world['ybkg']
    y = world['y']
    PBLH = mech.world['PBLH']
    if PBLH != PBLH_OLD:
        KDIL = np.minimum(1, PBLH_OLD / PBLH)
        mech.world['KDIL'] = KDIL
        world['PBLH_OLD'] = PBLH
        ydil = (KDIL * y) + (1 - KDIL) * ybkg
        nonzero = ybkg != 2.4195878568940516e-22
        y[nonzero] = ydil[nonzero]
    else:
        mech.world['KDIL'] = 1
    mech.update_world_from_y(y)

# Add an updater that depends on the y vector
mech.add_world_updater(func_updater(UpdatePBL, incr=360., verbose=False))
mech.world['PBLH_OLD'] = mech.world['PBLH']

Run Model

# Archive initial values
mech.archive()

# For each start/stop combination, update the world and integrate
for t0, t1 in start_end_ts:
    # set current time to t0
    t = t0
    # Get y vector from world state
    y = mech.update_y_from_world()
    # integrate from t to t1 with initial state of y
    # ts, Y = mech.integrate(
    #     t, t1, y0=y, solver='odeint', mxords=3, mxordn=3, atol=atol,
    #     rtol=rtol, mxstep = 1000, hmax = 300, verbose=False
    # )
    ts, Y = mech.integrate(
        t, t1, y0=y, solver='lsoda', atol=atol, rtol=rtol, nsteps=1000,
        max_step=300, max_order_ns=3, max_order_s=3, verbose=False
    )
    # Update world to match new time
    mech.world['y'] = Y[-1]
    mech.update_world_from_y()
    mech.Update_World(forceupdate=True)
    # Update y vector for mixing of PBL
    # Set new time to last time of integration
    t = ts[-1]
    # Update world from y vector
    mech.update_world_from_y()
    mech.archive()


# Optionally save to disk
# mech.output()

runend = time.time()
print((runend - runstart), 'seconds')
{t:3625,t/3600.:1,THETA:1.2E+02,PBLH:2.5E+02,TEMP:2.9E+02,O3:28}
{t:7210,t/3600.:2,THETA:1.1E+02,PBLH:2.5E+02,TEMP:2.9E+02,O3:26}
{t:7228,t/3600.:2,THETA:1.1E+02,PBLH:2.5E+02,TEMP:2.9E+02,O3:26}
{t:10830,t/3600.:3,THETA:1.1E+02,PBLH:2.5E+02,TEMP:2.9E+02,O3:25}
{t:14412,t/3600.:4,THETA:97,PBLH:2.5E+02,TEMP:2.9E+02,O3:23}
{t:14434,t/3600.:4,THETA:97,PBLH:2.5E+02,TEMP:2.9E+02,O3:23}
{t:18035,t/3600.:5,THETA:87,PBLH:3E+02,TEMP:2.9E+02,O3:23}
{t:21612,t/3600.:6,THETA:76,PBLH:5.5E+02,TEMP:2.9E+02,O3:26}
{t:21636,t/3600.:6,THETA:76,PBLH:5.5E+02,TEMP:2.9E+02,O3:26}
{t:25238,t/3600.:7,THETA:65,PBLH:8.5E+02,TEMP:2.9E+02,O3:27}
{t:28812,t/3600.:8,THETA:54,PBLH:1.1E+03,TEMP:3E+02,O3:28}
{t:28839,t/3600.:8,THETA:54,PBLH:1.1E+03,TEMP:3E+02,O3:28}
{t:32442,t/3600.:9,THETA:42,PBLH:1.3E+03,TEMP:3E+02,O3:28}
{t:36032,t/3600.:10,THETA:32,PBLH:1.4E+03,TEMP:3E+02,O3:29}
{t:36043,t/3600.:10,THETA:31,PBLH:1.4E+03,TEMP:3E+02,O3:29}
{t:39644,t/3600.:11,THETA:22,PBLH:1.5E+03,TEMP:3E+02,O3:30}
{t:43233,t/3600.:12,THETA:18,PBLH:1.5E+03,TEMP:3E+02,O3:33}
{t:43246,t/3600.:12,THETA:18,PBLH:1.5E+03,TEMP:3E+02,O3:33}
{t:46847,t/3600.:13,THETA:22,PBLH:1.5E+03,TEMP:3E+02,O3:35}
{t:50435,t/3600.:14,THETA:31,PBLH:1.5E+03,TEMP:3E+02,O3:39}
{t:50448,t/3600.:14,THETA:31,PBLH:1.5E+03,TEMP:3E+02,O3:39}
{t:54052,t/3600.:15,THETA:42,PBLH:1.5E+03,TEMP:3E+02,O3:41}
{t:57635,t/3600.:16,THETA:54,PBLH:1.5E+03,TEMP:3E+02,O3:44}
{t:57655,t/3600.:16,THETA:54,PBLH:1.5E+03,TEMP:3E+02,O3:44}
{t:61258,t/3600.:17,THETA:65,PBLH:1.5E+03,TEMP:3E+02,O3:45}
{t:64836,t/3600.:18,THETA:76,PBLH:1.5E+03,TEMP:3E+02,O3:45}
{t:64860,t/3600.:18,THETA:76,PBLH:1.5E+03,TEMP:3E+02,O3:45}
{t:68465,t/3600.:19,THETA:87,PBLH:1.5E+03,TEMP:2.9E+02,O3:45}
{t:72036,t/3600.:20,THETA:97,PBLH:1.5E+03,TEMP:2.9E+02,O3:45}
{t:72066,t/3600.:20,THETA:97,PBLH:1.5E+03,TEMP:2.9E+02,O3:45}
{t:75667,t/3600.:21,THETA:1.1E+02,PBLH:2.5E+02,TEMP:2.9E+02,O3:43}
{t:79241,t/3600.:22,THETA:1.1E+02,PBLH:2.5E+02,TEMP:2.9E+02,O3:41}
{t:79270,t/3600.:22,THETA:1.1E+02,PBLH:2.5E+02,TEMP:2.9E+02,O3:41}
{t:82872,t/3600.:23,THETA:1.2E+02,PBLH:2.5E+02,TEMP:2.9E+02,O3:38}
15.545097827911377 seconds

Make Plots

data = mech.get_output()
data['h_LST'] = data['t'] / 3600
data['NOx'] = data['NO'] + data['NO2']
data['jNO2*100'] = data['TUV_J(6,THETA)'] * 100
data['PBLH [km]'] = data['PBLH'] / 1e3
data['TEMP/25 [C]'] = (data['TEMP'] - 273.15) / 25.
data['OH/10 [ppqv]'] = data['OH'] * 1e5
data['HO2 [pptv]'] = data['HO2'] * 1e3

fig, ax = plt.subplots()
data.set_index('h_LST')[['emis_NO', 'jNO2*100', 'PBLH [km]', 'TEMP/25 [C]']].plot(ax=ax)
fig.savefig('knote_physical.png')

fig, ax = plt.subplots()
data.set_index('h_LST')[['NOx', 'NO', 'NO2']].plot(ax=ax)
ax.set(xlabel='hour [LST]', ylabel='ppb', ylim=(.1, 5), yscale='log')
fig.savefig('chemical_nox.png')

fig, ax = plt.subplots()
data.set_index('h_LST')[['O3', 'OH/10 [ppqv]', 'HO2 [pptv]']].plot(ax=ax)
ax.set(xlabel='hour (LST)')
fig.savefig('chemical_ozone.png')
  • plot knoteae2015
  • plot knoteae2015
  • plot knoteae2015

Total running time of the script: ( 0 minutes 16.635 seconds)

Gallery generated by Sphinx-Gallery