Updates to the CMAQ aqueous chemistry routine (aqch em)

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 1 0 6 + modeled hours, and

• Unexpected large differences in predicted sulfate m ass concentrations when comparing simulations employing a perturbation in p rimary 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 episo dic and variable.

After confirming with OAQPS that the emissions file s were not responsible, ORD focused on suspected numerical instabilities genera ted by the aqueous chemistry subroutine in CMAQ.

Initial testing by ORD focused on box modeling of g rid 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 an d 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 t he “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 becaus e the maximum number of iterations (i.e., 100) was reached in the “moves” c ase, this was causing too much sulfate production. During the last iteration in the oxida tion loop, the final time step is equal to the total time of the cloud minus the time spent al ready during previous iterations. For example, a one hour cloud – 20 minutes consumed dur ing the first 99 iterations results in a final 40 minute time step. Further box model tes ting demonstrated that increasing the max number of iterations from 100 to 10,000 solved the instability errors, however, model run time increased dramatically. ORD determi ned that the hydroxyl radical (OH) was controlling the aqchem oxidation timestep and t hat after approximately 30 iterations (in this test case), [OH] concentrations did not ch ange.

The aqueous chemistry subroutine in the multipollut ant 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 rea ction rate constants are very fast ~

10 10 (L/mol/sec) compared to other reactions in the CMA Q aqueous mechanism, (e.g.,

~10 4 for S(IV) oxidation by H 2 O 2 ). 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 insta bilities, as evidenced in the OAQPS
simulations.ORD determined that more precision was needed to ca pture the changes in OH

concentration and all OH related variables were cha nged to double precision (update 1). Testing with the change to double precision demonst rated that the aqueous solver only took 34 iterations for the "moves", same as base ca se 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 simul ation and for the simulation that used aqchem update 1. A difference plot of the two simu lations is shown in Figure 3.

ORD continued investigation of aqchem and incorpora ted additional updates and design changes. Testing demonstrated that mass was not be ing conserved for species produced via aqueous phase oxidation reactions during the la st 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 st ep = cloud life time – (sum of ∆ t for the previous 99 iterations for oxidation)) for the last (i.e., 100 th ) iteration. ORD implemented a design change so that now, after reaching 100 ite rations, 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 li mit for the oxidized amount. With this change, the remaining mass conservation problems we re 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 chan ge was implemented. Previously, OH reactions of glyoxal and methylglyoxal and Hg (f or the multipollutant model) controlled the time step of oxidation reactions. N ow only the sulfur reactions affect ∆ t. In CMAQ’s gas phase chemistry there is a pseudo-ste ady state assumption for OH and it is not a transported species. This is because hydr oxyl radical reactions tend to be catalytic (e.g., consumption and production). In t he 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 rep lenished via gas-to-cloud partitioning. However, due to operator splitting a queous and gas-phase chemistry are not solved simultaneously. To account for this and oth er uncertainties in predicted OH aqueous phase concentrations, (e.g., neglect of pro duction reactions (H 2 O 2 + h υ  2 OH) not currently implemented in aqchem).

ORD will be releasing an updated version of CMAQ (v 4.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 us ing the CMAQ multipollutant model. Note

maximum SO4 concentration of 283.4 ug/m3.

Figure 2. Simulation results for ASO4J using the C MAQ 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 -26 9.0 to 0.8 ug/m3.

Figure 4. Simulation results for ASO4J using the C MAQ 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 differ ences 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 ) SO 4 , HSO 4 , SO 3 , HSO 3 , CO 3 , HCO 3 , OH, NH 4 , HCO 2 , NO 3 , Cl

Start iteration and bisection (max. 3000 iterations )

Calc. final gas phase partial pressure of SO 2 , NH 3 , HNO 3 , HCOOH, CO 2

liquid conc. (mol L - 1 ) SO 4 , HSO 4 , SO 3 , HSO 3 , CO 3 , HCO 3 , OH, NH 4 , HCO 2 , NO 3 , Cl

Check for convergence

Compute ionic strength and activity coefficient (Da vies 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 fo r each species

TIME = TAUCLD (OR 100 iterations) Check for

100 max. iterations