Siyuan Wang | Research - Ch 11. Diffusion and Chemistry

Ch 11. Vertical Diffusion and Chemistry

OVERVIEW

- Okay in the previous chapter we learned how to setup a 1-D diffusion model, in which the vertical diffusion is solved by using a simple explicit scheme.
- In this chapter we'll spice things up by adding chemistry to it! There're many 1-D photochemistry model applications out there and I myself have a few publications too.
- To drive a 1-D model, you need a few inputs, such as temperature, pressure, turbulent diffusivity, etc. You can setup comprehensive physics to take care of all of those, but that's a bit beyond the scope of this tutorial, so I'll just grab these inputs from another model which is a very common thing to do. You can also make up fake (idealized) inputs too. The world is your oyster!
- One key input to drive a 1-D model, as you've seen in the previous chaper, is the height-dependent turbulent diffusivity (Kz). This is a fairly complex topic and and an active area of research as well. I'm not an expert on this, so I just took a fairly basic scheme from one of the most popular textbooks: Stull (2003). Here's the fun part: to calculate Kz, you need friction velocity (u*) and Obkuhov length (L) and their formulations kinda depend on each other. It's not entirely easy to solve it by hand and get the analytical solution. So we'll develop a Newton scheme to solve it iteratively. This is actually a quite useful skill!
- As you may have guessed, 1-D models are a lot more expensive than 0-D models, so here we'll use simpler chemical mechanisms: the Common Representative Intermediates Mechanism (CRI), which is also developed by the MCM folks. Link here. CRI is less complex than MCM. For instance, the mechanism we'll use in this chapter, containing C1-C3 akanes + isoprene + benzene + toluene, has ~400 reactions. The same MCM mechanism would have ~2000 reactions easily.
- Personally, this is one of the most exciting chapter of the tutorial, and I hope you enjoy it too!

1. Get all the files, load the mechanism, setup emissions/depositions

Go to the MCM/CRI website and mark a few species. We'll try to mimick an urban environment with a few simple VOCs, so let's mark CH4, C2H6, C3H8, isoprene (C5H10), benzene, and toluene. Then generate the KPP file. Don't forget to include inorganic reactions and generic rate coefficients! The KPP file is also uploaded to Github (cri_export_C1_C3_Isop_Aromatics.eqn).
Also download the Jupyter Notebook from Github (YAMCHA_tutorial_ch11.ipynb).
You also need a netCDF file with the inputs that I prepared (hrrr_info_for_1D_model.nc), which are from the HRRR model (High-Resolution Rapid Refresh), which is a NOAA-developed real-time convection-allowing, cloud-resolving weather model at 3-km resolution that is also updated hourly. HRRR outputs were obtained from Amazon Clouds (https://registry.opendata.aws/noaa-hrrr-pds/). You can certainly use other models, or use whatever fakeidealized inputs you'd like.
Run the first 3 cells which load a bunch of libraries and utility functions. You know the drill.
Drop the KPP file you just downloaded into the Google Colab storage (or your mounted Google Drive). Then in Cell #4, load the mechanism.

# ==========================================================

# Use this to load the chemical mechanism file in KPP format

# ==========================================================

ChemMech = load_mcm_kpp('/content/cri_export_C1_C3_Isop_Aromatics.eqn')

ChemMech.T

If loaded correctly, a summary will be printed. This mechanism has 425 reactions. Not too bad for a mechanism with these major VOCs!

Then in Cell #5, you'll setup emissions and depositions. I'm lazy so I'm not setting up emissions and depositions for all species. This is just to demo how emissions and depositions are setup.

ChemMech_EmisDep = pd.DataFrame({'E0': [' = C5H8', 'emission_rate_isoprene_molec_cm3_s', None],

'E1': [' = CO', 'emission_rate_CO_molec_cm3_s', None],

'E2': [' = NO', 'emission_rate_NO_molec_cm3_s', None],

'E3': [' = C2H6', 'emission_rate_C2H6_molec_cm3_s', None],

'E4': [' = C3H8', 'emission_rate_C3H8_molec_cm3_s', None],

'E5': [' = BENZENE', 'emission_rate_benzene_molec_cm3_s', None],

'E6': [' = TOLUENE', 'emission_rate_toluene_molec_cm3_s', None],

'D0': ['O3 = ', 'dep_rate_O3_sec', None],

'D1': ['NO2 = ', 'dep_rate_NO2_sec', None],

})

ChemMech_EmisDep.index = ['reaction', 'rate_coefficient', 'RO2']

ChemMech_EmisDep.T

Don't forget to merge this with the main mechanism!

ChemMech_merge = pd.concat([ChemMech,ChemMech_EmisDep],axis=1)

ChemMech_merge.T

2. Setup the vertical grid

Very important but fairly straightforward. We'll setup heights at layer centers and layer interfaces. Here the model top is 5km, with a vertical resolution (dz_m) of 100 m, so overall we'll have 51 vertical layers.

# --- height grid: layer center

model_top_m = 5000.

dz_m = 100

z_m = np.arange(dz_m/2, model_top_m+dz_m, dz_m)

nlev = len(z_m)

# --- height grid: layer interface

z_interface_m = np.zeros(len(z_m)+1)

for i in range(len(z_interface_m)):

if i==0: z_interface_m[i] = 0

else: z_interface_m[i] = 2.*(z_m[i-1]-z_interface_m[i-1])+z_interface_m[i-1]

3. Load the Input Data (netCDF)

Drop the netCDF file to the Google Colab storage (or in the mounted Google Drive), then load the file using xarray.

hrrr = xr.open_dataset('/content/hrrr_info_for_1D_model.nc')

hrrr

Familiarize yourself with this file and its contents. Should be mostly self-explanatory.

4. Process the Inputs

I grabbed the HRRR data every 3 hours (because I'm cheap), so these will be interpolated to the model timestamps. Also, these are on HRRR native vertical grids, and they need to be interpolated to the model vertical grids too. Note that time interpolation will be done at model runtime.
Firstly, let's interpolate the profiles (temperature, pressure, specific humidity) to the model vertical grid. I'll skip the source code here which are fairly basic. I'm using SciPy's built-in interp1d function, which is pretty nice.
The plots are here. I'm a careful person, and I do sanity checks often to make sure things are done right. The left plot is a sanity check, with blue lines representing HRRR profiles and orange lines representing the interpolated profiles. The right plot shows the curtain plots of pressure, temperature, and specific humidity interpolated to the model grids, which will be used in the model.

Then the next cell plots a few "surface" variables from HRRR that we'll use as well: planetary boundary layer height (pblh), 10-meter wind (U10), and sensible heat flux (hfx).
The figure is attached on the right. You can see PBLH reaches ~4 km at around noon time but remains low at night (a few hundred meters). Sensible heat flux also peaks at ~500 W/m2 at around noon, and is slightly negative at night. These are fairly typical for a good sunny day in the front range region. Basically the boundary layer is deep and unstable during daytime but shallow and stable at night.
Again these are not yet interpolated to the model timestamps, which will be done at the model runtime.

5. Calculate the Turbulent Diffusivity using Newton Approach

The formulations for friction velocity (ustar) and Obukhov length (L_Obukhov) are:

ustar = 0.4*uref/(np.log(zref/z0)-PSI(zol))

L_Obukhov = -1*(ustar**3.)*rho_air_kg_m3*Cp_J_kg_K*temp/0.4/9.8/HFX_W_m2

where uref is the wind speed measured at the reference height, zref, and z0 is the roughness. PSI is the integrated form of the stability correction function which is a function of the ratio between z and L_Obukhov (zol). Then L_Obukhov also depends on ustar, as well as a few other things like air density (rho_air_kg_m3), heat capcacity (Cp_J_kg_K), temperature (temp), and sensible heat flux (HFX).

The Newton approach is pretty interesting. Say you need to solve f(x) = 0. Firstly, you make an initial guess, x0. Then you draw a tangent line passing thru (x0, f(x0)). Next, you find x1 where the tangent line crossing y=0. Then you draw another tangent line passing thru (x1, f(x1))... rinse and repeat, until the result is "good enough": when f(x_n) <= tolerance, x_n is then the numerical approximation. Obviously I'm not doing a good job describing that, and this method has pros and cons, and you can find better documentations elsewhere. But since a figure is worth a thousand words, here I'm attaching three figures, with two unstable demos (z/L<0) and one stable demo (z/L>0). The accuracy is controlled by a tolerance here which is set to 1e-3. Mostly it converges <10 iterations.
A similar Newton solver was implemented in FORTRAN in a previous publication of mine (Wang et al. 2021).

Once this is done, we can move forward and calculate Kz!

Kz = 0.4*ustar*z_interface_m*((1.-z_interface_m/pblh_m)**1.5)/phi(zol)

where z_interface_m is the height at the layer interfaces (remember, Kz is represented at layer interfaces rather than layer centers), and phi is yet another stability correction function that depends on zol.

The results are shown on the right, with pblh also shown as a dashed line. You can see that Kz peaks at around noon but at ~2 km. The calculated Kz will later be used in the model to drive the vertical diffusion!

6. Setup A Few Diurnal Profiles

Now we'll create a few diurnal scaling factors which will be used to scale emissions (biogenic and anthropogenic) and solar zenith angle (SZA)!
First, the HRRR sensible heat flux is used as a "typical" diurnal profile shape, which remains low at night and peaks at around noon. This scaling factor will be used for biogenic emissions and SZA.
Second, in a typical urban setting, many anthropogenic emissions will follow traffic which has two rush hour peaks, one in the morning and the other in the afternoon. So I just made up a profile for that...
The plot is shown on the right!

7. Main Model

Finally we're here in the main 1-D model! Most of the stuff here in Cell #11 should look familiar to you. We have the functions that calculates the derivates and jacobian matrix for the chemistry solver, we also have the explicit scheme for diffusion (same as the previous chapter). Do make sure that you're loading the correct mechanism, i.e. the one with emissions/depositions! A few more things:
Initial conditions: because this is a 1-D model, the initial condition will be profiles too. Here we're setting constant mixing ratios across all layers. If you have profiles from other models, or if you want to set up profiles differently (e.g., different fix values in the boundary layer and free troposphere), be my guest!
Emissions and depositions: we discussed this before. For emissions, we'll provide maximum emission fluxes (in molec/cm2/s), and the diurnal profile shapes (anthropogenic and biogenic) will be applied later inside the main time loop. Depositions are provided as deposition velocities (cm/s).
Main time loop: this is where the time stamp is progressing. We'll also do time-interpolation on the inputs too. The HRRR inputs are at 3-hour resolution, but the model output frequency is set to be 900 seconds. So at each model time step, the closest two HRRR timestamps are grabbed, and HRRR variables are interpolated to the local model time step. Fairly straightforward.
Chemistry solver: this is called at every vertical level, and hence inside another loop. At each level, we grab the temperature, pressure, water vapor concentration first, and then the kinetics will be updated accordingly. Emissions and depositions are coupled with chemistry and done at only the lowest model layer.
Vertical diffusion: once chemistry is done in all vertical layers, we'll do vertical diffusion on all species. Note that in the previous chapter we discussed the (in)stability issue with the explicit scheme. Here because the chemistry output time step is fairly long (15 minutes here), and sometimes Kz values can be large, the chemistry time step may not work well for the diffusion! Therefore we dynamically determine the diffusion time step here:

dt_diffusion = 0.4* dz_m^2/max(kz_m2_s)

0.4 is an empirical factor which works in this demo and I've seen people using lower values such as 0.2. This is important to ensure the diffusion stability. When the diffusivity is large (and if you set higher vertical resolution, or smaller dz), you'll notice that dt_diffusion can be much smaller than the chemistry output time step. Fortunately the diffusion scheme is pretty fast!

Lastly, once the model is done, the results are packed into a xarray dataset.

8. Run the Model!

Let's run the model for 2 days, with 15 min output frequency. This mechanism has some 425 reactions (including emissions and depositions), and 51 vertical layers. With my free Google Colab account, it took me ~15 min to complete. Not too bad at all! If you have access to big clusters, it'll be a lot faster! For instance, on the NCAR super computer, this model runs roughly 3x faster!
However, if you'd like to run complex mechanisms (e.g. MCM), it can be slow.
The curtain plots and surface concentrations of several species are shown on the right, which are pretty interesting! Let's quickly go over it.
O3 was uniformly distributed in the model, initially. With the fresh NO emission in the early morning, a small fraction of O3 near the surface is quickly titrated, leading to the formation of NO2 (easier seen in the time series plot). Later in the day, the photolysis of NO2 produces NO and keeps NO2 low. O3 is also produced during daytime.
N2O5 is produced at night but mostly in the residual layer, because it's colder up there!
OH and HO2 peak during daytime, with OH approaching 1e+7 molec/cm3 and HO2 approaching 3-5e+8 molec/cm3. OH is a bit lower near the surface, because many OH scavengers (isoprene and other VOCs) are released at the surface!
Isoprene is elevated at the surface. Although its emission peaks at noontime, it's concentration is actually slightly lower at noon. We talked about this before, and this is because OH radical concentration is also very high at noontime.
HCHO is an oxidation produce from isoprene and other VOCs, and it's fairly well-mixed throughout the mixing layer, also remains high at the nighttime residual layer!

9. Summary

In this chapter, we integrated the vertical diffusion with the chemistry solver and put together a fun 1-D photochemistry model! You can expand and add additional components (e.g. aerosols, aqueous-phase chemistry, etc).
Compared to chemistry, diffusion is generally fast. But if you run the model with high vertical resolution (small dz), the diffusion timestep might need to be small which can slow down the model. You may do diffusion on "long-lived" species only and skip diffusion for radicals and very short-lived species, such as O1D, O3P, OH, HO2, etc. Some 3-D models do this to improve the computational efficiency. However, in this demo it doesn't matter -
You've probably noticed that I added a few timers in the main time loop to track how long would chemistry vs diffusion take. At night, all the chemistry steps take about 1-2 seconds, while the diffusion solver takes 0.02 seconds. During daytime, chemistry becomes more complicated which takes 5-6 seconds; diffusion also takes longer because of the dynamically adjusted diffusion timestep, which is up to 1-2 seconds. As you can see, chemistry is definitely the bottleneck here, which is often the case in other chemical transport models too. A few things that might reduce the computational cost. e.g. you can remove emissions/depositions and add them in the diffusion solver (I showed you this in the previous chapter), which could help a little. You may also relax the chemistry tolerance. You could also parallel computing because we're solving chemistry sequentially (in all vertical layers) and the chemistry in one layer does not depend on other layers.
We also covered a few technical stuff! A simple Newton method is implemented to iteratively solve friction velocity and Obukhov length. The Newton approach is a very useful numerical trick that I have used in quite a few occasions in my work in multiple models (WRF-Chem, CESM/CAM-Chem). My implementation here also has an option to make diagnostic plots which will help you better understand the process. When coupling chemistry with the diffusion, I also demonstrate how to dynamically adjust diffusion timestep to ensure the stability.

Google Sites

Report abuse