CMAQ 4.3 Release Notes

Description of Model Changes for the 2003 Release of the CMAQ Model

1. Introduction

The Community Multiscale Air Quality Model (CMAQ) has been extensively revised in 2003. Changes include: updated science, corrected implementations, efficiency enhancements, and bug fixes. There is a preliminary evaluation of the new release of the CMAQ model on the EPA Models-3 website comparing model results to observations of gas and aerosol species. There are also model-to-model comparisons between the 2003 release and the 2002 release. The biggest changes involve aerosol modeling, particularly nitrate aerosols and secondary organic aerosols (SOA). Nitrate modeling was updated so it is consistent with the most recent literature and the SOA implementation was corrected to allow for reversible semi-volatility. These changes resulted in substantially lower concentrations of both aerosol nitrates and SOA. Minor changes have also been made to aqueous processes and dry deposition. The details of changes to chemistry/transport modeling are presented in Section 2.

There have been major modifications to improve model efficiency. A new fast gas phase chemistry solver, known as the Euler Backward Iterative (EBI) scheme, has been developed for the CB4 mechanism. Also, some of the fastest reacting species have been dropped from the transport processors. The time step for operator splitting has been revised to allow different advective time steps by vertical layer. Details of these changes are described in Section 3.

Note that other components of the CMAQ system such as SMOKE and MCIP have also been revised recently. Revisions to MCIP are especially important because serious errors in the processing of the horizontal and vertical wind fields have been corrected. Thus, one needs to use version 2.2 of MCIP (see MCIP release notes for details) for CMAQ modeling.

2. Chemistry/Transport Modeling Changes

2.1 Secondary Organic Aerosols (SOA)

Background: The formation of secondary organic aerosols (SOA) in the atmosphere occurs via the oxidation of organic compounds to form semi-volatile gases that can subsequently condense on particles. Because of their volatility, these compounds are also capable of evaporating from the particle-phase back to the gas-phase as temperatures increases. The SOA gas-particle partitioning model used in CMAQ is based on the Secondary Organic Aerosol Model (SORGAM) described by Schell et al. ( 2001). The CMAQ SOA module accounts for SOA formation from the oxidation of six primary organic groups: alkanes, alkenes, cresol, high-yield aromatics, low-yield aromatics, and monoterpenes. Ten semi-volatile, gaseous SOA precursors are produced via these reactions - two each for olefins, monoterpenes, and aromatics, and one each for alkanes and cresol. The gas-phase reaction yields and the volatility properties of the individual semi-volatile et al., 1997; Kalberer et al., 2000; Strader et al., 1999; Griffin et al.,1999).

Updates: Two major problems with the SOA module contained in the 2002 release of CMAQ were identified. First, the manner in which the gas-particle partitioning model was implemented allowed for condensation of the semi-volatile gases to the particle phase, but not for the reverse process of evaporation. Second, the amount of semi-volatile SOA precursors formed from the atmospheric oxidation of some of the organic compounds of anthropogenic origin could be over-estimated in some cases. The latest release of CMAQ contains a number of changes to correct these problems. In addition, some minor revisions to the treatment of SOA from biogenic sources have also been made. The net effect of these changes is to reduce CMAQ predictions of both anthropogenic and biogenic SOA under most conditions. These changes are all internal to the model, and, hence, do not require any specific user action. Each is described individually in more detail below.

Reversibility. As noted above, the current release of CMAQ includes modifications to the SOA module to make the gas-particle partitioning of semi-volatile compounds fully reversible (i.e., allow for both evaporation and condensation). Accounting for reversibility necessitated making some changes to the treatment of semi-volatiles in the model. In the previous version of CMAQ, the gas-phase components of the semi-volatile SOA precursors were represented by ten model species, and the corresponding particle-phase components of the semi-volatiles were lumped into the two generic anthropogenic and biogenic SOA species. The names of the gas-phase species were tagged with the prefix VAP_ to indicate they represented species in the gas-phase. In the new reversible version of the SOA model, the ten gas-phase species have been replaced by ten species that represent the sum of the gas-phase and the particle-phase concentrations. (Note that only the total gas- plus particle-phase concentration is needed for partitioning, not the concentrations in each phase. For computational efficiency, only total concentrations are tracked in the revised model.) Hence, the prefixes of the ten semi-volatile species have been changed from VAP_ to SGTOT_ to indicate that these species represent the sum of the gas- and particle-phase mechanism. The units for these species in the CMAQ output files are in ppmV however. The representation of the anthropogenic and biogenic components of the SOA remains unchanged (i.e., SOA is represented by the model species named AORGxy, where x is either A or B for anthropogenic or biogenic, respectively, and y is either I or J for Aitken mode or accumulation mode, respectively). The corresponding output units for these species are still ug/m3. Finally, the replacement of the set of gas-phase species with the set of species representing the sum of the gas and particle phases necessitated changing the dry deposition computations for these species. In the previous release, a gas-phase deposition velocity was used for the ten semi-volatile species since they were treated as gases. In the new approach, a mass-based, weighted average of gas-phase and particle-phase deposition velocities is used for these species since compounds in the different phases should dry deposit at significantly different rates under most conditions. Note that all of these changes are internal to the model.

SOA production from alkenes. As noted in the background section, the atmospheric reactions of alkenes are assumed to form two semi-volatile SOA precursor species. In the previous version of CMAQ, these semi-volatiles were formed via the reactions of the SAPRC99 and RADM2 mechanism species OLE2 and OLI, respectively. (SOA production from alkenes in the CB4 mechanism was not included in the previous release.) The OLE2 and OLI species were assumed to be of anthropogenic origin only, and, hence, the SOA produced from the semi-volatile species were included in the anthropogenic SOA category. A review of emission inventory data shows that not all compounds from anthropogenic sources that are lumped into the OLE2 or OLI mechanism species form SOA precursors. In addition, some compounds emitted from biogenic sources that do not produce SOA precursors are also lumped into the SAPRC99 OLE2 species. As a result, production of SOA precursors from alkenes was likely overestimated in the previous version of CMAQ when the SAPRC99 or RADM2 mechanisms were used. In addition to the emissions inconsistencies, the vapor pressures of the semi-volatile gases formed via this pathway were set too low in the previous release. Tests conducted with a special version of CMAQ using corrected data have shown that the contribution to SOA formation from alkenes is very small compared to other chemical pathways. Note that changes to the emissions processing method and the addition of special "bookkeeping" type changes to the CMAQ model would be required to properly apportion the production of SOA from alkenes into their biogenic and anthropogenic components. To avoid introducing this additional level of complexity for essentially very little change in model results, the production of SOA precursors via the reactions of alkenes has been set to zero in the new release. Place holders are maintained in the model code should this pathway be included in future releases. Consequently, the CMAQ output files in the new release will include the semi-volatile, SOA precursor species produced by alkene reactions (SGTOT_OLI_1 and SGTOT_OLI_2), but their concentrations will be zero.

SOA production from alkanes. In previous versions of CMAQ, the semi-volatile SOA precursor produced by the oxidation of alkanes was based on the reaction of the ALK5 species in the SAPRC99 mechanism and the HC8 species in the RADM2. In that implementation, all emitted compound slumped into ALK5 and HC8 were assumed to produce a semi-volatile SOA precursor gas. A review of emissions data shows that many compounds lumped into these model species do not form SOA precursors. Therefore, in the new release, the molar yields of the semi-volatiles produced by reactions of ALK5 and HC8 are multiplied by the average mole fraction of SOA producers that are lumped into these species. The average mole fractions (0.57 for ALK5 and 0.56 for HC8) were derived from an annual U.S. inventory of anthropogenic emissions, and are applied uniformly in time and space. Note that production of SOA from alkanes with the CB4 mechanism was omitted in the previous version of CMAQ, and that remains the case with this release. Applying the same approach to CB4 is complicated by the fact that the CB4 species into which alkanes are lumped, PAR, also includes contributions from many non-SOA producing biogenic and anthropogenic compounds. Methods to account for SOA production via this pathway in CB4 will be examined for future releases.

SOA production from aromatics. Similar to alkanes, the production rate of semi-volatile gases from high-yield aromatics (ARO1 in SAPRC99 and TOL in RADM2) has been adjusted downward to account for the fact that benzene, a non-SOA producer, is lumped into these mechanism species. In the new release, the molar yields of semi-volatile SOA precursors from ARO1 and TOL are multiplied by 0.93 and 0.94, respectively, to account for the inclusion of benzene in these species. These fractions were derived using an annual U.S. inventory of anthropogenic emissions, and are applied uniformly in time and space. No adjustment factor is needed for CB4 since benzene is not lumped into either of the aromatic compounds in that mechanism.

SOA production from monoterpenes. The molecular weights of the two semi-volatile SOA precursors produced by the oxidation of monoterpenes have been changed from 184.0 to 177.0 to be consistent with the molecular weight of the corresponding particle-phase SOA. In addition, changes have been made to the rate constants used in the CB4 and RADM2 mechanisms for the reactions of monoterpenes. Since monoterpenes were not treated explicitly in either CB4 or RADM2, special tracking reactions were added to the CMAQ versions of two mechanisms to obtain the amounts of semi-volatile SOA precursor produced by the reaction of monoterpenes (Gipson and Young, 1999). In previous versions of CMAQ, the rate constants used for the monoterpene reactions in both mechanisms were based on the rate constants for the reactions of the RADM2 species OLI, the species into which monoterpenes are lumped. In the new release, the monoterpene rate constants in both CB4 and RADM2 have been changed to the values used in SAPRC99 for the species TRP1, which is the SAPRC99 explicit species for monoterpenes. Hence, the three mechanisms are now consistent in this regard. It should be noted that the choice of rate constants for monoterpenes in CB4 and RADM2 does not affect the gas-phase chemistry in either mechanism since the oxidants reacting with monoterpene in these reactions are regenerated.

2.2 Nitrate Modeling

Background: A significant pathway for the formation of nitrate aerosol involves the heterogeneous conversion of N2O5 to HNO3. In the pre-2002 versions of CMAQ, heterogeneous conversion was not explicitly included. However, the earlier versions of the RADM2 and CB4 gas-phase mechanisms included a homogeneous hydrolysis reaction converting N2O5 to HNO3 at a rate that may have implicitly included some effects of heterogeneous conversion (e.g., Stockwell et al., 1990). In the 2002 CMAQ release, the heterogeneous pathway was explicitly added using the model described by Dentener and Crutzen (1993). In that approach, the heterogeneous conversion of N2O5 to HNO3 is modeled as a pseudo first-order loss process (i.e., N2O5--> 2HNO3), where the reaction rate constant is a function of the aerosol surface area, the N2O5 gas-phase diffusion coefficient, particle radius, and reaction probability. The amount of HNO3 retained in the aerosol phase is determined by equilibrium considerations. With the implementation of the heterogenous pathway, the rate constants for the corresponding gas-phase hydrolysis reactions in CB4 and RADM2 were lowered to the value used in SAPRC99. The SAPRC99 rate constant is almost an order of magnitude lower and reflected a more recent estimate of the upper limit on the homogeneous, gas-phase conversion rate.

Updates: The same approach is incorporated in the latest release to CMAQ, but three major changes have been made. First, the amount of HNO3 produced by the pseudo first-order reaction in the 2002 CMAQ release was a factor of 1.7 too high due to a units conversion error. This problem has been corrected in the latest release. The second change pertains to the reaction probability used in calculating the pseudo first-order rate constant for converting N2O5 to HNO3. In the 2002 CMAQ release, the reaction probability was set to 0.1 after Dentener and Crutzen (1993). For the 2002 CMAQ release, the fixed reaction probability has been replaced by the method of Riemer et al. (2003) in which the reaction probability is computed as a function of aerosol composition. It is based on experimental results suggesting that the heterogeneous N2O5conversion rate is inhibited by the presence of nitrates in the aerosol phase (Mentel, 2000). In the Riemer method, the reaction probability is computed as function of the mass ratio of aerosol sulfate to the sum of aerosol sulfate and nitrate. If this ratio is 1.0 (i.e., all sulfate) the reaction probability is set to an upper limit of 0.02. If the ratio is zero (i.e., all nitrate), the reaction probability is set to a lower limit of 0.002. Reaction probabilities for intermediate ratios are assumed to scale linearly between these two limits. Thus, using this method always results in a lower reaction probability than used in the 2002 release. The third change affects the homogeneous, gas-phase conversion of N2O5 to HNO3. Since evidence exists to suggest that the conversion of N2O5 to HNO3 proceeds in hydrated aerosols as opposed to the gas-phase (e.g., Jacob, 2000), the rate constants for the gas-phase hydrolysis reactions of N2O5have now been set to zero in all of the gas-phase mechanisms with linkages to the latest CMAQ aero3 module. The original mechanism rate constants are used in all other situations. The net effect of all these changes is to significantly lower CMAQ predictions of aerosol nitrates under most conditions.

2.3 Clouds: Chemistry, Scavenging and Wet Deposition

Several modifications were made to the CMAQ cloud module. The most significant changes were made to the routine which calculates the Henry's law constants. This routine is now capable of calculating the effective Henry's law value using dissociation equilibrium constants from Seinfeld and Pandis (1998). This change will affect species that do not participate in cloud chemistry (e.g. HCHO, HNO2, HO2) by allowing them to dissociate, which should increase their scavenging and wet deposition. However, for other species (e.g. SO2, HNO3, NH4, H2O2, O3, MPH, HCOOH), the aqueous chemistry routine already calculated their dissociation to compute the ionic strength of the cloud water, so this change will not affect those species. A few other changes were made to the Henry's law routine, including updates to the coefficients of the table using data from Sanders (1999), and the addition of MEK, MVK, and METHACROLEIN to the table.

Within the aqueous chemistry module, the section of the code that calculates the incremental time-step for the cloud was changed to be more stable than in the previous version. This change was necessary because at very low oxidant concentrations (near the machine precision), the program crashed because of division by zero. The new algorithm is more streamlined and eliminates some of the needless divisions.

Changes were made to the washout of coarse aerosol number concentrations (NUMCOR). Previously, the NUMCOR washout was scaled against the SO4 removal rate---but SO4 is also produced within clouds. Linking NUMCOR in this way created a problem for the diameter of coarse-mode particles, since the mass of coarse particles was washed out depending on the rainfall rate. To solve this problem, NUMCOR is now washed out in the same manner as the coarse mass. This should preserve the aerosol diameter during the washout process.

2.4 Dry deposition

here are two dry deposition models available in MCIP, M3dry, which includes direct linkages to the Pleim-Xiu land surface model (PX LSM) in MM5, and RADMdry, which was originally developed for the RADM model. Version 2.2 of MCIP includes updates to both dry deposition schemes with an emphasis on winter simulations. In particular, the M3dry routine has significant updates for MM5 runs that are run with or without the PX LSM. The effects of snow cover are now considered by M3dry as well as RADMdry.

There are three new dry deposition species: N2O5, NO3, and generic aldehyde. The nitrogen species, N2O5and NO3, were added as part of the aerosol nitrate review and revisions. In addition, there were several species for which dry deposition velocities were previously computed in MCIP but were not used in CMAQ for dry deposition: PAN, HONO, CO, and Methanol. The dry deposition surrogates for all chemical mechanisms have been updated to
include these species.

2.5 Corrections for vertical diffusion of aerosols

Previously aerosol species were treated as mass density units in the vertical diffusion processing. This is incorrect because vertical diffusion was developed initially for gas-phase concentrations in molar mixing ratio. We converted the aerosols to molar mixing ratio (ppmV) units, including emissions sources and dry deposition sinks. The vertical profiles are now more realistic, typically with higher concentrations at lower levels.

2.6 Re-ordered processes for operator splitting

Previously vertical diffusion was done after advection, mass-adjustment and horizontal diffusion in the time-splitting configuration for CMAQ. This was then followed by the cloud, gas-phase chemistry and aerosol processes. It is now felt that a better characterization of the time-splitting would be effected by doing the vertical distribution of concentrations before advection, primarily to get emissions injected into the proper layers before advection takes place (emissions are injected in the vertical diffusion processing).

In addition, the symmetric operator splitting processing option was eliminated. This option has never been fully tested, and it is felt that it is less accurate when recording the outputs. The chemistry process should be the last before writing the concentrations to the output file, since then the concentrations would be more in chemical balance than otherwise.

The result of these changes should show improvement in the characterization of time-splitting, and preliminary tests indicate this. Because of these changes, the couple/decouple processes are now called only in sciproc, and the decoupling of the average concentration array, AGRID, is no longer required. These reductions reduce the CPU time somewhat, but not significantly.

3. Modifications to Achieve Improved Efficiencies

3.1 New Gas Phase Chemistry Solver

Background: The previous releases of the CMAQ model have included a non-generalized but fast chemistry solver - the Modified Euler Backward Iterative (MEBI) method as originally developed by Hertel et al. (1993) and modified by Huang and Chang (2000). As noted in previous releases, the MEBI solver has proven to be significantly faster than the other chemistry solvers in CMAQ and quite accurate when bench marked against the high accuracy Gear method. The only major shortcoming is that it cannot be completely generalized to handle any chemical mechanism, but rather must be individually customized for each mechanism. Nevertheless, it has become the solver of choice in most model applications because of its speed and accuracy.

Update: The MEBI solver included in previous releases is based on the methods described by Hertel et al (1993), except that the analytical solutions derived for the NO/NO2/O3/O3P and OH/HO2/HNO4/HONO groups of species were replaced by numerical solutions as formulated by Huang and Chang (2000). Use of the numerical solution enhances generalization, but may increase computational costs. This release includes a new EBI solver in which the numerical methods used for the two groups were replaced by the original analytical formulations of Hertel et al (1993) as applied to the CMAQ version of CB4. Note, however, that the basic structure of the EBI solver is similar to MEBI's. Tests conducted with EBI show it to be about 2 times faster than the MEBI solver. Differences in model predictions between the two solvers were found to be minimal (e.g., maximum differences on the order of 1 ppb for ozone). Hence, it is expected that the EBI solver will replace the MEBI solver even though the latter is also included in this release.

Like MEBI, the CMAQ EBI solver has certain limitations. First, it can only be used with the CB4 family of mechanisms. Extensions to other mechanisms are planned, but are not yet available. Second, the EBI solver cannot be used to conduct integrated reaction rate (IRR) analysis because of the fairly large time step used in EBI. Integrated process rate (IPR) analysis can be invoked however. Finally, the EBI solver includes error control parameters that affect both speed and accuracy. A default set of error control parameters has been incorporated in the solver to enhance model speed while maintaining accuracy. If the error control parameters were to be changed for a particular application, it is highly recommended that testing and comparison with a Gear type solver be performed first to ensure that the desired degree of accuracy is achieved.

Usage: The use of the EBI solver is completely analogous to using the MEBI solver module. Normally, the chemical solver module is selected in the configuration file or run script that is used to build the executable CMAQ. To use the EBI solver, the module ebi_cb4 should be selected. No other CMAQ processors are affected. The only current restriction is that the ebi_cb4 module must be used in conjunction with one of the following CB4 mechanisms: cb4, cb4_aq, cb4_ae2, cb4_ae2_aq, cb4_ae3, or cb4_ae3_aq. At present, there is no consistency check when the model is built, but execution of the model will be stopped with an error message at run time if an inconsistency is detected. Similarly, if an attempt is made to invoke IRR analysis, the model simulation will abort with an error message. Finally, it was noted above that the EBI solver was roughly 2 times faster than MEBI. Reductions in model run times will vary, depending on the model configuration. Tests conducted with the CMAQ model configuration that included aerosols gave overall model runtime reductions of about 15% when MEBI was replaced with EBI.

3.2 Elimination of transport for fast-reacting CB4 species

In previous releases of CMAQ, all of the species included in the gas-phase chemistry mechanism have been subjected to transport (i.e., advection and diffusion in the horizontal and vertical). Some of these species are at very low concentrations and react very rapidly however. If the fast reacting species concentrations are governed primarily by the concentrations of slower reacting species, it might be possible to forego transport of the fast reacting species, thereby reducing computational costs. A series of tests was conducted with CMAQ using the CB4 to test this hypothesis. The results of these tests indicate little loss in accuracy when transport is eliminated for the following CB4 species: O, O1D, OH, HO2, C2O3, ROR, CRO, TO2, XO2, and XO2N. For example, maximum differences in ozone were less than 1 ppb over a multi-day simulation period. In fact, the differences were sufficiently small to negate the need to invoke steady-state approximations, thereby maximizing computational savings. As a result, the CB4 mechanism files have been configured to eliminate the transport of these 10 species in the latest release. The run-time savings obtained with this change will vary depending on the model configuration. When aerosols are included, model run time reductions should be on the order of about 5%. Thus, when this configuration is combined with the new EBI solver, run time reductions on the order of 20% can be realized as compared to using the MEBI solver and transporting all species.

3.3 Changes in Time-stepping

Previously, the time-step for operator splitting was determined by the shortest Courant-condition limited advective time-step anywhere in the 3-d domain. Thus, the synchronization time-step tends to be dictated by high winds aloft, leading to unnecessarily short time-steps for advection in the lower layers where wind speeds are much less. This can be mitigated to some extent by using a longer user-specified minimum synchronization time-step, which then may require multiple horizontal advection steps per sync step. This approach, however, can waste computation time with no accuracy benefit because all layers are still required to have the same minimum Courant-condition limited advective time-step. Since the region of greatest modeling interest is the planetary boundary layer (PBL), the algorithm was modified to try to get reasonable sync steps near the surface with equal advection steps and allow multiple advection steps aloft, if needed.

There is a new control parameter which specifies the model layers which determine the synchronization time-step. The new input parameter (SIGMA_SYNC_TOP) is the sigma value of the highest model layer to be used in determining the sync time-step. All model layers from the ground up to this sigma level will have the same advective time-step which is the minimum Courant-condition limitation in these layers. The model synchronization time-step is also set to this value such that these model layers will execute horizontal advection once per sync time-step. The default value for this parameter is 0.7 which reflects a height above ground of about 2.5 km. This is intended to include most of the daytime PBL in most areas. For model simulations where the maximum PBL height is expected to be very different from this height (e.g. dry climates and/or high elevations) the user can specify a different value for this parameter. The ZF parameter on the METCRO3D file from MCIP can provide guidance on the heights of the sigma layers.

There also has been a modification of model defaults for maximum and minimum synchronization time-step values. On the basis of model testing with 32 km grid resolution compared to benchmark runs with 6 minute synchronization time-steps, the maximum recommended (default) sync time-step is now set to 12 minutes rather than the previous 15 minutes. Also, based on testing at 2 km grid resolution, the recommended (default) minimum sync time-step is now 1 minute rather than the previous 5 minutes. The following table shows an example of what the new algorithm does for part of a 32 km horizontal resolution run with 21 layers. For this model integration step the synchronization time-step was determined to be 12 minutes on the basis of the minimum Courant-condition time-step in the lowest 14 layers (up to sigma = 0.7):

Layer

Advection Step
HHMMSS

Per Sync Step

21

000600

2

20

000600

2

19

000400

3

18

000600

2

17

000600

2

16

000600

2

15

001200

1

14

001200

1

13

001200

1

12

001200

1

11

001200

1

10

001200

1

9

001200

1

8

001200

1

7

001200

1

6

001200

1

5

001200

1

4

001200

1

3

001200

1

2

001200

1

1

001200

1

Initial tests compared with fixed, short synchronization steps show an improvement in the accuracy in the lower layers. Also, there is some decrease in CPU time for coarser resolutions, but the opposite is sometimes true for the finer resolutions.

3.4 Parallel Processing Changes

When the code runs in parallel in a so-called MPICH cluster, synchronization between the node that launches the run and the other participating nodes can be a problem, evidently due to network file system (NFS) latency. A spin-wait code was employed to attempt to deal with this problem by delaying execution for 60 sec. after processor 0 issued an IO/API call to open a new CONC file. However, the spin-wait solution allows incorrect IO/API header data to be written to the output CONC and ACONC files. Specifically, the netCDF TFLAG variable, representing the timestamp record), can have invalid or no data at the first date/time for a number of file variables. A correction obtained from Zion Wang of CERT-UC-Riverside has been applied, which involves simplifying the opening and closing of the output concentration file for the case that the file does not previously exist. The problematic spin-wait code is no longer used.