Environ. Sci. Technol. 1996, 30, 1163-1175
Air Quality Predictions of the Urban Airshed Model Containing Improved Advection and Chemistry Algorithms SANDRA L. WINKLER AND DAVID P. CHOCK* Ford Research Laboratory, Ford Motor Company, P.O. Box 2053, Mail Drop 3083, Dearborn, Michigan 48121-2053
We have carried out an extensive study on the effect of replacing the advection and chemistry algorithms of the Urban Airshed Model (UAM) with more accurate algorithms, using the SCAQS episode of August 26-28, 1987, and input data supplied by the California Air Resources Board. Replacing the UAM chemistry solver by the implicit-explicit hybrid solver (IEH) tends to lower the afternoon ozone concentrations somewhat, the afternoon PAN concentrations more noticeably, and the H2O2 concentrations significantly. IEH is also much less likely to have a convergence failure than the UAM solver. Replacing the Smolarkiewicz scheme by the forward Euler Taylor Galerkin scheme (FETG) yields a greater impact. The replacement increases the local peak NO concentrations and lowers the local peak ozone and CO concentrations significantly. Replacing both chemistry and advection algorithms lowers the peak ozone, PAN, and H2O2 concentrations quite noticeably. For example, the predicted peak ozone for August 28 is reduced from 234 to 200 ppb. The complicated wind field for this episode generates two high-ozone subdomains in the modelsone in the east and one in the northwest. Reducing the NOx emission by 50% shifts the predicted domainwide peak ozone from the east to the northwest. If we disregard this shift, then the use of more accurate modules in the UAM leads to a 50% increase in the domainwide peak ozone while the original UAM yields a 27% increase. Both the original and improved UAM versions give a 20% increase in peak hourly ozone and a 30% increase in peak 8-h ozone in the eastern subdomain, but much greater and vastly different increases for each of the peak 1-h and peak 8-h ozone concentrations in the northwestern subdomain. The best available numerical methods should be used in modeling a complex system like air quality so that potential hidden errors may be identified and the design of air quality control may be more reliable.
0013-936X/96/0930-1163$12.00/0
1996 American Chemical Society
Introduction The 1990 Clean Air Act Amendments require the use of a photochemical grid model for ozone air quality control planning in ozone nonattainment areas that are classified as serious, severe, or extreme. Several photochemical grid models exist for the study of air pollution episodes in urban areas: the Urban Airshed Model (UAM, developed by Systems Applications International), the CIT Photochemical Airshed Model (CIT, developed by California Institute of Technology and Carnegie Mellon University), CALGRID (developed by Sigma Research Corporation), etc. Presently, only the UAM has been approved by the U.S. Environmental Protection Agency for use in control planning applications. Because of the enormous impact of the UAM predictions, it is important that they be based on the best possible science and input data and that the numerical methods for the model solutions be accurate. The grid models describe the transformation and transport of air pollutants. Transformation typically entails several tens of chemical species in some 80 to more than 100 chemical reactions with a very wide range of reaction rate constants. Accurate chemistry solvers exist but require too much computing time to be practical for air quality simulations; less accurate but fast chemistry solvers have been used in existing grid models. The present version of the UAM (1) uses the Carbon Bond IV chemical mechanism, which is solved by a hard-wired UAM solver. This solver was found (2) to be fairly accurate. But it also has an execution time that depends sensitively on some reaction rate constants and on the initial concentrations of each time step and has some tendency not to converge readily. Transport refers to the advection and diffusion of air pollutants. Solving the advection-dominated transport equation has been a major numerical hurdle in air quality modeling. This is because spatial discretization can alter both the amplitudes and the phase speeds of the Fourier components of the concentration profile, leading to artificial diffusion and dispersion of the solution. The present version of the UAM uses Smolarkiewicz’s scheme (3) (SMOL, denoted STS3 in ref 4) to solve the horizontal advection step in the model. This scheme preserves the sign of the concentration but is highly diffusive when the concentration background is zero and yet can yield an oscillatory solution that may even enhance the solution peak near the region of a steep concentration gradient when the background is not zero (4). Because of the perceived diffusive nature of this scheme, the horizontal diffusivity of the UAM may have been made small to avoid excessive horizontal diffusion. Uncertainties in basic assumptions and input data aside, nonlinear chemistry in the model can magnify inaccuracies introduced by either the chemistry solver or the transport algorithm (5). Use of deficient solvers or algorithms can also create compensating errors, mask or obscure errors in model assumptions and input, and lead to failure in control strategies. Since more accurate chemistry solvers and advection algorithms have recently become available, it behooves us to test and implement them in the UAM in hope of improving our understanding of the atmosphere and the model performance. In the UAM, the top of the second layer generally corresponds to the mixed-layer height or a minimum height
VOL. 30, NO. 4, 1996 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
1163
of 50 m, below which the layers are of equal depth. This layer depth can vary from cell to cell in the horizontal direction. Most other models use vertical layers of fixed depth. In a separate paper (6), the impact on air quality prediction using a fixed vertical layer structure in the UAM will be explored. In this paper, we investigate the potential impact on air quality predictions and emissions control of replacing the present UAM chemistry solver and advection scheme with other viable candidates. It is important to note that improving the numerical methods in an air quality model should improve the model performance when all the relevant model assumptions and input data are sufficiently accurate. Otherwise, the model performance may become worse because some hidden errors may be exposed. (As a case in point, if the physical horizontal diffusivity is made small to accommodate the perceived diffusive nature of the Smolarkiewicz scheme, then this error by design will express itself when a more accurate advection scheme is used in the model.) The following sections describe our selection of the numerical methods to replace the corresponding modules in the UAM; the assumptions made in the modified UAM models and the simulation conditions; the results of the simulations and the conclusion.
Selection of Numerical Methods Chemistry Solvers. Chemistry is the most time-consuming portion of the modeling process. Gear-based codes such as LSODE (7) accurately solve stiff systems of chemical equations but with excessively long computing times. Fast chemical kinetic solvers used in current grid models such as the hybrid method (8), the quasi-steady-state approximation (9), and the UAM chemistry solver reduce computation time but may also reduce accuracy (2). Recently, Sun et al. (10) developed a highly accurate implicit-explicit hybrid (IEH) fast solver. This new solver has been favorably compared with the above fast solvers (2). It is very accurate, robust, and fast under a wide range of initial conditions. It also has a greater tendency of reaching a convergent solution than the UAM solver. Accordingly, the IEH solver will be chosen to substitute the UAM solver in our UAM modification. Advection Algorithms. Advection algorithms are continuously being invented, modified, and in many cases, abandoned. Numerous algorithms have been tested (see, e.g., refs 4, 11, and 12). A few rather promising algorithms for solving the advection equation have been identified recently. They include the forward Euler Taylor Galerkin (FETG, denoted FCFTG/FE in ref 4) method (13) and the accurate space derivative (ASD, denoted FASD/FE in ref 4) method (14). Both schemes, especially the latter, retain the advected concentration peaks well. But FETG is still subject to numerical dispersion resulting in oscillations in the solution. ASD, on the other hand, encounters the Gibbs phenomenon when discontinuity is present in the concentration profile and/or its low-order derivatives, leading to oscillations in the solution around the point of discontinuity. In either case, negative concentrations often appear. The resulting concentration oscillations may be suppressed by the use of the Forester filter (15) so that the occurrence and magnitude of negative concentrations may be reduced. In our applications, all remaining negative concentrations will have to be reset to zero to assure stability in the chemistry. The spurious gain in mass is generally away from the concentration peaks and is generally small. There are also flux-limiting schemes designed to prevent
1164
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 30, NO. 4, 1996
the occurrence of negative concentrations. One of them is the fully multi-dimensional flux-corrected transport algorithm of Zalesak (16). This scheme, to be denoted ZAL, was also studied previously (denoted MFCT/PC in ref 11) and was found to perform moderately well but to also have a tendency to suppress concentration peaks. More recently, Bott’s scheme (17), also flux limiting and monotonic, has gained some popularity in some air quality modeling applications. The scheme uses time splitting and assumes the concentration of a grid cell to be represented by an area-preserving integral (over the grid cell) of a polynomial function of space. Thus, the area below the polynomial function determines the outflow or inflow advective flux. However, although the polynomial function assigned to each grid is determined by the concentrations of the grid and the surrounding grids, the functions for neighboring grids are generally not continuous across the grid boundary. Moreover, the polynomial function within a grid can have negative values. As a result, the flux calculated by this method can be quite inaccurate, especially when the concentration profile is complex and the degree of the polynomial functions is high. Consequently, unexpected distortion of the concentration profile may occur in the application of Bott’s scheme. To avoid negative concentrations, monotone flux limitations are imposed, but this can cause a distortion of the profile when time splitting is used because the flow field is not necessarily divergence free in each direction. Bott has proposed a remedy to reduce the inaccuracy due to time splitting (18). ZAL, on the other hand, corrects the multi-dimensional flux without time splitting. Based on the above information, we selected FETG, ZAL, and ASD as our advection schemes to replace Smolarkiewicz’s scheme in the UAM.
Methodology Our UAM modification process has three parts. Part 1 examines the effect of changing only the chemistry solver used in the UAM. Part 2 evaluates the effect of changing the advection scheme. Part 3 combines changes for both the chemistry solver and the advection algorithm and evaluates their combined effect. Chemistry. We began evaluating the effect of changing chemistry solvers by replacing the UAM solver in the original UAM with a low-tolerance version (relative tolerance ) 10-7; absolute tolerance ) 10-11 ppm) of the IEH solver. This version was found to yield solutions within 10-6 of the practically “exact” LSODE01 (LSODE with a relative tolerance of 10-7, an absolute tolerance of 10-11 ppm, and a Jacobian update at every integration step) solutions (10) at eight times the speed. The UAM solution with this version of IEH will be denoted UAM/Ex. The solution of UAM coupled with the high-tolerance version of the IEH solver will be denoted UAM/I. For this model version, the relative and absolute tolerances of 10-2 and 10-6 ppm, respectively (IEH26 in ref 2), are used for the daytime simulations, and the relative and absolute tolerances of 10-3 and 10-7 ppm, respectively (IEH37 in ref 2) are used for the nighttime simulations. UAM/Ex provides a standard against which both UAM and UAM/I will be compared. Advection. Comparison of different advection algorithms for horizontal transport in the UAM is somewhat complicated by the fact that the vertical wind field in the UAM is internally calculated. The UAM does not use input vertical wind but instead calculates it by requiring that a
nonreactive calibration gas concentration be constant and equal before and after the advective transport in each time step. If the originally uniformly-constant calibration gas concentration deviates from the constant after a horizontal advection, then the vertical wind field must be adjusted such that the calibration gas concentration is restored to the constant after the application of the donor-cell differencing scheme for vertical advection. This is an ingeneous scheme to remove any local mass imbalance (of the calibration gas) regardless of its source. However, since the source includes the concentration distortion created by the horizontal advection algorithm, errors generated by the algorithm can therefore become concealed by the calculated vertical wind field. How seriously these errors influence the vertical wind field is difficult to ascertain. We have considered the case of the UAM with a fixed vertical layer structure and an externally-supplied vertical wind field elsewhere (6). In this paper, replacing the Smolarkiewicz scheme in the UAM by FETG, ZAL, and ASD yields solutions denoted by UAM/FETG, UAM/ZAL, and UAM/ASD, respectively. Chemistry and Advection. Having evaluated the effect of chemistry solvers and advection algorithms separately, we introduced both types of modifications simultaneously and examined their combined effect. Replacing the chemistry solver in the UAM by IEH and the horizontal advection scheme by FETG, ZAL, and ASD leads to results denoted by UAM/IFETG, UAM/IZAL, and UAM/IASD, respectively.
Model Assumptions and Simulation Conditions All models are based on the Urban Airshed Model-CBM IV code version 6.21 and dated 93287 (1). The code contains a correction for an error in the UAM chemistry solver of version 6.20 that underestimated the XO2 loss rate by a factor of 2 in the XO2 + XO2 reaction (19). The code was modified to accommodate different advection algorithms and/or chemistry solver as the case might be. Each of the chemistry solvers requires different assumptions. The UAM solver has a single relative tolerance input as the convergence criterion, while the IEH solver uses a relative and an absolute tolerance to determine convergence. For the UAM solver, we used a tolerance of 0.02, as recommended in Volume II of the UAM User’s Manual (20). For the IEH solver, we used IEH26 for the daytime and IEH37 for the nighttime simulations as mentioned earlier. In some initial test runs, we did use IEH26 for nighttime simulations as well and obtained results considerably more accurate than those of UAM but not as accurate as those with IEH37 for some fast species concentrations. We chose to retain IEH37 for nighttime runs to preserve high accuracy in the solutions. A default value of 10-16 ppm (2) was assumed for the initial concentrations of the nonadvected fast (UAM’s steady-state) species in each time step. It turns out that with the present version of CBM-IV, the concentrations of XO2N and, to a considerably lesser degree, those of NO3 and N2O5 are dependent on the initial concentrations, especially for nighttime. The use of a default initial value of 10-16 ppm for XO2N tends to suppress its subsequent concentration and leads to a slight exaggeration of the concentrations of O3. To ameliorate this problem without having to transport the species, we assigned 10-6 ppm as the initial conditions for XO2N, NO3, and N2O5 and used their concentrations at the end of each time step as the initial concentrations for
the next time step in the same grid cell. (Note: In a more recent correction of the XO2N mechanism in CBM-IV, the default value of 10-16 ppm for the initial concentrations of all nonadvected fast species does not lead to any significant deterioration of the solutions for all species.) The FETG and ASD algorithms were substituted into the code and coupled with the diffusion algorithm (unchanged) in the same manner as in the original UAM. FETG and Smolarkiewicz’s scheme as applied in the UAM are based on time splitting while ZAL and ASD are fully twodimensional. In UAM, concentrations that are small after advection are set to a lower bound. For FETG and ASD coupled with the IEH chemistry solver, negative concentrations due to the advection scheme are set to zero. At each time step, two consecutive iterations (for both x and y) of the Forester scheme were applied to FETG and ASD to smooth out ripples represented by 2∆x waves (12). The choice of the smoothing coefficient, µ, in the Forester scheme requires some experimentation. Too high a value will suppress the oscillation drastically and may also suppress the peak concentrations significantly. Too low a value may leave essentially untouched the concentration profile generated by the advection scheme; with a resulting concentration profile that may be checkerboard-like or stripy but with well-maintained peaks. Obviously, concentration profiles containing many steep and flat gradients (like the primary pollutants) or having poorly matched boundary conditions tend to generate oscillations and are highly sensitive to the choice of µ while species with a smoother concentration distribution would not be. For example, in our 1-day model trial runs with UAM/IFETG, compared to a µ of 0.001, a µ of 0.01 led to a suppression of the domainwide daily maximum concentration of
