Updates to the CMAQ aqueous chemistry routine (aqchem)

August 5, 2009 

Shawn Roselle and Ann Marie Carlton 

During application of the multipollutant version of CMAQv4.7, OAQPS encountered 

several issues, including: 

• Abnormally high sulfate concentrations (>400 ug/m3), 

• Large temporal discontinuities in predicted sulfate aerosol mass (e.g., ∆ > 40 

(ug/m3) in one hour), encountered in a few of the 10

6

+ modeled hours, and 

• Unexpected large differences in predicted sulfate mass concentrations when 

comparing simulations employing a perturbation in primary PM emissions, below 

referred to at the “MOVES” case. In these modeling scenarios there were no 

changes in sulfur emissions. 

In general, these observed abnormalities were episodic and variable. 

After confirming with OAQPS that the emissions files were not responsible, ORD 

focused on suspected numerical instabilities generated by the aqueous chemistry 

subroutine in CMAQ. 

Initial testing by ORD focused on box modeling of grid cells in Alberta, Canada where 

the problem seemed to be the most pronounced. Box model testing showed the aqueous 

chemistry solver took 34 iterations to solve for and complete the oxidation reactions (see 

Appendix 1 for a map of aqchem solver). In aqchem, there is a maximum of 100 

iterations for oxidation chemistry. When testing the “moves” emissions scenario, the 100 

iteration limit was removed and ORD discovered that 1148 iterations was required for the 

"moves" emissions case. ORD determined that because the maximum number of 

iterations (i.e., 100) was reached in the “moves” case, this was causing too much sulfate 

production. During the last iteration in the oxidation loop, the final time step is equal to 

the total time of the cloud minus the time spent already during previous iterations. For 

example, a one hour cloud – 20 minutes consumed during the first 99 iterations results in 

a final 40 minute time step. Further box model testing demonstrated that increasing the 

max number of iterations from 100 to 10,000 solved the instability errors, however, 

model run time increased dramatically. ORD determined that the hydroxyl radical (OH) 

was controlling the aqchem oxidation timestep and that after approximately 30 iterations 

(in this test case), [OH] concentrations did not change. 

The aqueous chemistry subroutine in the multipollutant version of CMAQv4.7 includes 

new reactions with the hydroxyl radical to convert mercury (Hg) covalence states and to 

make “cloud SOA” (Carlton et al., (2008)). The reaction rate constants are very fast ~ 

1010 (L/mol/sec) compared to other reactions in the CMAQ aqueous mechanism, (e.g., 

~104 for S(IV) oxidation by H2O2). The coupling of faster and slower processes in a 

forward Euler solver (i.e., the solver employed by aqchem.F for oxidation reactions) can 

create a stiff system that leads to numerical instabilities, as evidenced in the OAQPS 

simulations. 
ORD determined that more precision was needed to capture the changes in OH 

concentration and all OH related variables were changed to double precision (update 1). 

Testing with the change to double precision demonstrated that the aqueous solver only 

took 34 iterations for the "moves", same as base case when double precision was 

employed. OAQPS conducted simulations of the same time period with the double 

precision version of aqchem and observed that this update corrected roughly 70% of 

observed problems. Figures 1 and 2 show the ASO4J concentrations at hour 7 on March 

4, 2005 for the original multipollutant model simulation and for the simulation that used 

aqchem update 1. A difference plot of the two simulations is shown in Figure 3. 

ORD continued investigation of aqchem and incorporated additional updates and design 

changes. Testing demonstrated that mass was not being conserved for species produced 

via aqueous phase oxidation reactions during the last cloud timestep. This issue is only 

important episodically, when the maximum iterations is reached (100). During the 

aqueous chemistry solving routine, timestep for the later iteration is expanded to the 

remainder of the cloud lifetime (i.e., last time step = cloud life time – (sum of ∆t for the 

previous 99 iterations for oxidation)) for the last (i.e., 100th) iteration. ORD implemented 

a design change so that now, after reaching 100 iterations, a minimum 1 second timestep 

for the remainder of the cloud lifetime is imposed. Mass is conserved in the later part of 

the cloud lifetime, as total mass available is a limit for the oxidized amount. With this 

change, the remaining mass conservation problems were eliminated. Figures 4 and 5 

compare SO4 for the update 2 change in aqchem. 

Considering that the aqueous chemistry routine was designed to predict sulfate and not 

stiff OH-reacting mechanisms, a further design change was implemented. Previously, 

OH reactions of glyoxal and methylglyoxal and Hg (for the multipollutant model) 

controlled the time step of oxidation reactions. Now only the sulfur reactions affect ∆t. 

In CMAQ’s gas phase chemistry there is a pseudo-steady state assumption for OH and it 

is not a transported species. This is because hydroxyl radical reactions tend to be 

catalytic (e.g., consumption and production). In the aqueous phase chemistry, OH is 

absorbed by cloud water and in an open cloud model (i.e., the design approach currently 

employed), absorbed species (e.g., OH) would be replenished via gas-to-cloud 

partitioning. However, due to operator splitting aqueous and gas-phase chemistry are not 

solved simultaneously. To account for this and other uncertainties in predicted OH 

aqueous phase concentrations, (e.g., neglect of production reactions (H2O2 + hυ  2 OH) 

not currently implemented in aqchem). 

ORD will be releasing an updated version of CMAQ (v4.7.1) in early Fall 2009, that will 

include the aqueous chemistry updates. In addition, we are developing a more robust 

solver for aqueous chemistry (ROS3) that we plan to include in the next major CMAQ 

release (2011). 
Figure 1. Original simulation results for ASO4J using the CMAQ multipollutant model. Note 

maximum SO4 concentration of 283.4 ug/m3. 

Figure 2. Simulation results for ASO4J using the CMAQ multipollutant model with aqchem update 

1. Note maximum SO4 concentration of 24.3 ug/m3. 
Figure 3. Difference in ASO4J between simulations using the CMAQ multipollutant model produced 

by aqchem update 1. SO4 differences range from -269.0 to 0.8 ug/m3. 

Figure 4. Simulation results for ASO4J using the CMAQ multipollutant model with aqchem update 

2. Maximum SO4 concentration of 24.3 ug/m3. 
Figure 5. Difference in ASO4J between simulations using the CMAQ multipollutant model w/ 

aqchem update 2 compared with update 1. SO4 differences range from -2.8 to 1.4 ug/m3. 
Appendix 1. 

Map of aqchem.F 

Molar conc. = initial amt. – amt. deposited (mol L-1) 

bisection for pH, initial guesses between 0.01 – 10

liquid conc. (mol L-1) SO4, HSO4, SO3, HSO3, CO3, HCO3, OH, NH4, HCO2, NO3, Cl

Start iteration and bisection (max. 3000 iterations)

Calc. final gas phase partial pressure of SO2, NH3, HNO3, HCOOH, CO2 

liquid conc. (mol L-1) SO4, HSO4, SO3, HSO3, CO3, HCO3, OH, NH4, HCO2, NO3, Cl

Check for convergence

Compute ionic strength and activity coefficient (Davies Eqn.)

Calculate liquid concentrations and final gas phase concs. of oxdidants 

Kinetic calcs Cal. Min time step – check for large time step 

Don’t let DT > TAUCLD 

Compute wet depositions and phase concentrations for each species 

TIME = TAUCLD (OR 100 iterations) 

Check for 

100 max. iterations