CMAQ output sanity check

[skunwar]

I just ran CMAQv5.3.2 for Jan 1 - March 1, 2014 using 2016 NEI emissions, and then used COMBINE and CALC_TMETRIC postprocessing tools to get monthly average of PM25_TOT and other PM25 components. I am getting very low average (of hourly) February PM2.5 for the Western US:

When I checked the inline emissions, windblown dust does not seem to show any emission during any hour of the day (28th Feb) in the CCTM_DUSTEMIS_* output file:

Also, COMBINE (or rather species definition file of cb6r3 mechanism) has this PM25_UNSPEC1 which is pretty high average hourly concentration:

Since I am very new to postprocessing CMAQ output, I am trying to figure if there maps make sense, so I have a few CMAQ quality check questions to see where things could have gone wrong:

1. What are some diagnostic tests that I can do in CMAQ itself to make sure things went right, especially with the windblown emissions? How can I make use of CCTM_PMDIAG_* files here?
2. How should I understand PM_UNSPEC1 from COMBINE tool as an entity?
3. I should probably use ground observations to compare these simulations - is there a one-stop-shop to get the ground observations like IMPROVE, AQS and CASTNET for the present decade to compare with? May I also ask if AMET is the right tool for me to compare CMAQ monthly average of hourly PM2.5 (components and total) with observations?

Thanks!

[hogrefe.christian]

Thanks for your post.

I checked one of our recent runs for February 2016, and while the general patterns of monthly average PM2.5 concentrations are quite similar to what you’re showing, the area with monthly average PM2.5 concentrations <2 ug/m3 in the Western U.S. isn’t quite as widespread in our runs as in yours. We do see 2-4 ug/m3 for much of Arizona and New Mexico and values > 8ug/m3 for the San Joaquin Valley and Southern California. In our simulations, wind-blown dust was not turned on, so that should bias our results on the low side.

Are you seeing any emissions in the CCTM_DUSTEMIS* output files on any day? If not, have you double checked that CTM_WB_DUST is set to TRUE in the run script?

PM25_UNSPEC1 as defined in the species definition file includes the non-carbon portion of OM (NCOM). Assuming a typical OM:OC ratio of 1.8-2.0 (though it varies in space and time), NCOM spatial patterns and magnitudes typically are very similar to those of OC, so for the Eastern U.S. where OC is often highest, NCOM tends to dominate UNSPEC1. UNSPEC2 excludes NCOM by subtracting it from UNSPEC1.

Yes, AMET is a good tool to perform model evaluation and further investigate where your simulations deviate from observations. You can obtain observations (incl. AQS, IMPROVE, and CASTNET) packaged for AMET from the CMAS data warehouse. For further information on AMET, see https://github.com/USEPA/AMET

PM2.5 concentrations especially during wintertime and at night can be very sensitive to PBL height and potentially even the thickness of your first model layer, so this is something you might want to look at.

[wyat.appel]

I agree with Christian. The pattern of PM2.5 looks correct, but I think the concentrations are on the low side, especially in the western U.S. I looked at some February 2014 observed PM2.5 data and the values in the western U.S. are 5-10 ugm-3 in most places, with higher values in the San Joaquin Valley and southern California. So I agree that your values do look a bit low.

AMET is a good analysis tool to use and I encourage you to try to get it setup. I do think it will reveal a large underprediction of PM2.5 compared to observations. We just need to figure out why.

Out of curiousity, have you looked at the ATOTIJ species in your combine file? That’s PM2.5 without the sharp cutoff applied. It should look very, very similar to the PM25_TOT species, but I just want to make sure that’s the case. And since you’re new to CMAQ post-processing, I’d make sure that everything worked correctly in both combine and the calc_metric scripts (check the log file for anything weird). Your plot looks more like a single day than a monthly average. I doubt that’s the problem either, but again good to rule these things out upfront before we dig deeper elsewhere.

Wyat

[skunwar]

Thanks a lot @hogrefe.christian and @wyat.appel for your quick and insightful responses. In response to your suggestions, I drew more plots (average of hourly Feb 2014) from the postprocessed CMAQ output:

ATOTIJ: @wyat.appel, this looks similar to PM25_UNSPEC1 above but not identical

PBLH: @hogrefe.christian, the PBL height plot below is the average for all of Feb, but you are probably interested in seeing nighttime and daytime PBLH separately? The bottom layer in CMAQ is between the pressure levels VGLVLS = 1.f, 0.995f. Is this bottom layer height OK with the PBLH shown below?

I looked at additional quantities just to see if they would give any clue about the low Western US PM.

O3: The pattern of higher O3 (average of Feb hourly O3) in the western Rocky Mountains than in the eastern US looks right. I am thinking I am missing the hotspots in eastern US cities because it is winter.

CO: I also took a quick check of CMAQ CO concentration with EPA website’s observed values, although the quantities are different. 1ppm (=1000ppb) 8 hour max CO (national average) in EPA charts below are probably comparable to monthly average of ~300ppb in the CMAQ simulation below.

[wyat.appel]

Thanks for these plots. They are very helpful. Your February average ozone looks pretty high relative to our 2016 simulation, even in the eastern U.S. After discussing with Christian, we agree, at least at this point, that given that you’re seeing low PM2.5 and high ozone, particularly in the west, perhaps there’s an issue with the meteorology you’re using. Specifically, your average February PBL heights are quite a bit higher than the February average PBL heights in our 2016 simulation. Also, it appears that your surface layer thickness is approximately 40m, which is twice the thickness that we use (20m), which could also be contributing to the low PM2.5 concentrations you’re see.

So, at this point we’re suspecting that there may be issues with the PBL and/or mixing in your meteorology data. Have you done any analysis or evaluation of your meteorology? Also, can you provide the details on what meteorology you’re using (e.g. WRF version, PBL scheme, LSM)? Those will be helpful to have as well.

If your meteorology is too aggressive with the vertical mixing, you could see low PM2.5 concentrations (from dilution) and high ozone values (mixing down higher ozone values from aloft). Anyway, it’s a good place to start looking given the performance issues you’re seeing.

Wyat

[skunwar]

Thanks for the tips! I think you just helped me spot an inconsistency- I checked the WRF output files, they use SF_SURFACE_PHYSICS = 2 which turns out to be Noah LSM; but the CMAQ runscript has set the environment variables such that setenv PX_VERSION Y (and setenv CLM_VERSION N; setenv NOAH_VERSION N). Could this inconsistent LSM setting in CMAQ cause all the issues in PM and O3 in the plots above?

[wyat.appel]

Others can probably weight in better than I about this specific issue, but having consistency between the WRF inputs and CMAQ is important. And those flags were put in there specifically to deal with the differences in the LSM models and maintain consistency within CMAQ. So while I’m not sure it’s the root cause of what you’re seeing, it could certainly contribute. Others who are more knowledgeable would be able to explain how. But you certainly should set the appropriate flags for your meteorology you’re using.

I would still encourage you to do some evaluation of your met fields to make sure there are not other potential issues. Hopefully others will chime in about the CMAQ settings you mention.

Wyat

[cgnolte]

I would not want to promise that setting PX_VERSION to N and NOAH_VERSION to Y will solve all your issues, but it might be a good start.

What is the horizontal grid spacing you used? What version of WRF and MCIP, and what vertical coordinate? How many vertical layers, and what is the depth of layer 1? (This is variable ZH in the METCRO3D file.)

[skunwar]

Thanks a lot for your quick help.

I had used downscaled meteorology from global model (2 deg by 2.5 deg lat-lon) to two nested domains -36km and 12km - in WRFv3.9.1.1 with 34 vertical layers. For the run that gave above plots, I am feeding that 12km WRF output to CMAQ. For context, I picked the following header information of one of the mcip output files:

LAY = 34 ;
VGTOP = 5000.f ;
VGLVLS = 1.f, 0.995f, 0.989f, 0.983f, 0.976f, 0.968f, 0.959f, 0.95f, 0.94f, 0.929f, 0.916f, 0.902f, 0.887f, 0.871f, 0.853f, 0.833f, 0.811f, 0.7
87f, 0.76f, 0.731f, 0.699f, 0.663f, 0.624f, 0.581f, 0.534f, 0.482f, 0.429f, 0.375f, 0.322f, 0.268f, 0.214f, 0.161f, 0.107f, 0.054f, 0.f ;
MCIP V5.1  FROZEN 11/21/2019
INPUT METEOROLOGY DATA FROM WRF ARW V3.9.1.1
CUMULUS PARAMETERIZATION:  MSKF with radiative feedback                                                                                                         
SHALLOW CONVECTION:  No shallow convection                                                                                                             
MICROPHYSICS:  WSM 6-Class                                                                                                             
LONGWAVE RADIATION:  RRTMg                                                                                                             
SHORTWAVE RADIATION:  RRTMg                                                                                                             
PBL SCHEME:  YSU                                                                                                             
SURFACE LAYER SCHEME:  Revised MM5 (Jimenez)                                                                                                             
LAND-SURFACE SCHEME:  NOAH Land-Surface Model                                                                                                             
URBAN MODEL:  No urban physics                                                                                                             
LAND USE CLASSIFICATION:  NLCD40                                                                                                             
3D ANALYSIS NUDGING:  SPECTRAL                                
WIND COEFF:   3.000E-04 s-1                                                     
TEMP COEFF:   3.000E-04 s-1                                
MOIS COEFF:   4.500E-05 s-1                                                     
GEOP COEFF:   3.000E-04 s-1                                                                                                             
SFC ANALYSIS NUDGING:  OFF                               
WIND COEFF:  not applicable                                                     
TEMP COEFF:  not applicable                                
MOIS COEFF:  not applicable                             
OBS NUDGING:  OFF                                                                  
WIND COEFF:  not applicable                                
TEMP COEFF:  not applicable                                                     
MOIS COEFF:  not applicable                                                                                 
EARTH RADIUS ASSUMED IN MCIP:  6370000.000 m

I also had this vertical coordinate setting in WRF namelist.input file:

eta_levels                           = 1.000, 0.995, 0.989, 0.983,
                                       0.976, 0.968, 0.959, 0.950,
                                       0.940, 0.929, 0.916, 0.902,
                                       0.887, 0.871, 0.853, 0.833,
                                       0.811, 0.787, 0.760, 0.731,
                                       0.699, 0.663, 0.624, 0.581,
                                       0.534, 0.482, 0.429, 0.375,
                                       0.322, 0.268, 0.214, 0.161,
                                       0.107, 0.054, 0.000

I took a snapshot of the ZH variable (mid-layer height in meters above the ground, hour 0 of Feb 1, bottom layer) from mcip output - I didn’t readily have the Feb month average for the ZH variable. Most of the US has value of ~19m except the Rocky Mountains with about 17m value for the timestamp:

I will try to use M3TPROC command to get ZH variable Feb average if the snapshot is not as useful.

[cgnolte]

It’s not necessary to compute the monthly-averaged layer height. As Wyat surmised earlier, you are using a height of around 38 m for layer 1, about twice as thick as has been typically used here the past few years. That could be contributing a little to your lower concentrations.

Your use of the Noah LSM could also be a factor. I am unfamiliar with the Jimenez surface layer scheme, has your group used that previously? And since this is downscaled meteorology from a global climate model, that could play a role as well.

[skunwar]

Thanks a lot for suggesting the different avenues for further exploration, namely (1) checking the parent global model, (2) changing the LSM model, (3) looking at how Jimenez surface layer scheme behaves compared to other schemes, (4) evaluating meteorological simulation with observations. I will look at these options and get back to you again, but I also wanted to ask few more things here before running CMAQ again:

1. It is recommended (in CMAQ user guide) that we inherit vertical coordinates from WRF. Since this is already done in WRF, I can rerun WRF for a short period (2 months only) with new vertical coordinate spacing that is denser at the bottom. Do you mind providing me the eta_levels that you have been using recently for WRF namelist.input in your own WRF simulations?

2. My windblown dust emissions are not showing any emissions in the CMAQ output file for windblown dust. Although this should not be the cause of the current low PM and high O3, I wanted to make sure it runs. I had used the setting:

   setenv CTM_WB_DUST Y
   setenv CTM_WBDUST_BELD UNKNOWN

I think BELD3 option for environment variable CTM_WBDUST_BELD did not work in the runscript, so I had thought above combination would use MCIP landcover database, but I may be missing something.

3. Does lightning NOx show as a variable in itself somewhere in the CMAQ output?

[hogrefe.christian]

1. We often use the following eta levels in our WRF simulations, the lowest layer resulting from this setup has a thickness around 19 m.
   eta_levels = 1.000, 0.9975, 0.995, 0.990, 0.985,
   0.980, 0.970, 0.960, 0.950,
   0.940, 0.930, 0.920, 0.910,
   0.900, 0.880, 0.860, 0.840,
   0.820, 0.800, 0.770, 0.740,
   0.700, 0.650, 0.600, 0.550,
   0.500, 0.450, 0.400, 0.350,
   0.300, 0.250, 0.200, 0.150,
   0.100, 0.050, 0.000

2. My understanding matches yours, I don’t have a good clue about what causes the absence of windblown dust emissions in your simulations.

3. You should be able to obtain a separate diagnostic file with lightning NOx emissions by setting LTNGDIAG to Y in your run script.

[skunwar]

I am following up with a more complete ‘6-month evolution’ of plots for average monthly PM2.5 total, O3 and PBLH for Jan-June 2014 (6 months), whereas I only had one month (Feb 2014) in above early posts from a while ago. Basically, I am looking to check (or rather get checked ) if my CONUS CMAQ-simulated monthly means of PM2.5/O3 and meteorological variable PBLH closely represent present day meteorology/concentration as per observation or other simulation for those months. January is basically spinup, and PM2.5 total had been underestimated for Feb 2014 especially in Western US, while O3 had been overestimated (high PBLH was suspected to cause so):

Mean monthly PM2.5 total first, in order from Jan to June, but in the same color scale:

Corresponding 6-month O3 evolution is below:

And the corresponding PBLH monthly means are:

[skunwar]

Hi @hogrefe.christian , I was curious what NEI emissions year+version was used in the reference February 2016 CMAQ simulation you referred to to check PM2.5_TOTAL monthly average?

[hogrefe.christian]

Unfortunately I do not recall which file I looked at back then. I am sure the NEI year would have been 2016, but I don’t know if the version was beta (2016ff) or 1 (2016fh).

[skunwar]

I was suspecting it might be 2011 or 2014, but thanks for letting me know!

[skunwar]

@wyat.appel

your average February PBL heights are quite a bit higher than the February average PBL heights in our 2016 simulation

Hi @wyat.appel, would it be possible to share the month-by-month or season-by-season average PBLH from your 2016 simulation?