- Email: [email protected]
Mass conservation in the Community Multiscale Air Quality model
ARTICLE IN PRESS
Atmospheric Environment 40 (2006) 1199–1204 www.elsevier.com/locate/atmosenv
Mass conservation in the Community Multiscale Air Quality model Yongtao Hu, M. Talat Odman, Armistead G. Russell School of Civil and Environmental Engineering, Georgia Institute of Technology, Atlanta, GA 30332, USA Received 17 February 2005; received in revised form 15 October 2005; accepted 15 October 2005
Abstract Mass conservation errors were discovered when the Community Multiscale Air Quality (CMAQ) model (version 4.3) was coupled with the Fifth-Generation PSU/NCAR Mesoscale Model (MM5). These errors can be severe particularly over complex terrain. A recent correction made to the vertical advection module (in version 4.4) reduces the mass conservation errors, but does not completely eliminate them. This means that the method of renormalizing concentrations according to air density, which is currently applied to address the inconsistency problem, is not an ideal solution. Reconstruction of the wind fields by adjusting the vertical wind is proposed as a simple remedy to make CMAQ strictly mass conservative. r 2005 Elsevier Ltd. All rights reserved. Keywords: Regional air quality modeling; Inconsistency of meteorological fields; Mass conservation error; Source-receptor relationship; Emissions control strategy
1. Introduction Air quality models (AQMs) such as the Community Multiscale Air Quality model (CMAQ) (Byun and Ching, 1999) are expected to maintain a uniform mass mixing ratio field for an inert tracer after transport with the winds produced by meteorological models (MMs), such as the Fifth-Generation PSU/NCAR Mesoscale Model (MM5) (Grell et al., 1994). This expectation is realized only if the two models use the same discretization, i.e., grid, time step, and finite difference forms (Odman and Russell, 2000). In an integrated air quality modeling system, e.g. SMOG (Lu et al., 1997), the meteorolCorresponding author. Tel.: +1 404 894 1854; fax: +1 494 894 8266. E-mail address: [email protected] (Y. Hu).
1352-2310/$ - see front matter r 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.atmosenv.2005.10.038
ogy component and chemistry/transport (i.e., air quality) component naturally share the same grid structure. However, in several present-day air quality modeling applications, CMAQ and MM5 do not share the same grid structure and spatial interpolations are usually required. Also, when AQMs and MMs are off-line coupled, the outputs of the MMs are stored less frequently than the AQMs’ time step. For example, MM5 outputs are usually stored hourly while CMAQ has a maximum transport time step of 900 s. Input meteorological variables are reconstructed for AQMs’ time intervals, usually by linear interpolation between the stored values. Time step differences may also occur in on-line coupled integrated air quality modeling systems. For example, in SMOG, the dynamical equations are solved at a time step of 6 s while tracer transport is calculated every 300 s, accordingly,
ARTICLE IN PRESS 1200
Y. Hu et al. / Atmospheric Environment 40 (2006) 1199–1204
parameters predicted by the MM are averaged over the appropriate number of dynamical time steps to obtain values for the AQM (Lu et al., 1997). Finally, transport algorithms (e.g. finite difference forms) used for advection in AQMs are usually different from those in MMs. Because of all these differences in discretization, the wind components and the air density do not necessarily satisfy the discrete continuity equation in AQMs even though they may have done so in MMs. This is referred to as ‘‘the inconsistency of the meteorological fields’’ (Byun, 1999; Odman and Russell, 2000; Lee et al., 2004). As a result of this inconsistency, uniform mass mixing ratio fields are perturbed. The perturbation is usually more pronounced with data from non-hydrostatic MMs than with data from hydrostatic MMs for the same domains (Odman and Russell, 2000; Lee et al., 2004). To solve this inconsistency problem, some widely used AQMs, including CMAQ, adopted an approach where a uniform tracer is transported along with other species and the concentrations are renormalized using the deviation of the tracer from uniformity. In CMAQ, this uniform tracer is air density (Byun, 1999). As a result, a uniform mixing ratio field is maintained (Lee et al., 2004). However, the approach overlooks that non-linear advection schemes are used in AQMs for high accuracy. The non-linearities affect pollutant fields of different distribution differently. Most species of interest, especially the emitted ones, have spatial distributions far from being uniform. Therefore, renormalizing the concentration of a species based on the perturbation of a uniform field can artificially increase or decrease mass. Odman and Russell (2000) showed that this attempt to establish consistency and produce stable results in AQMs sacrifices mass conservation of pollutant species and found that the mass conservation errors introduced by this approach could be very large. Species mass conservation is not sacrificed if the consistency of meteorological fields is re-established through the continuity equation. To satisfy the discrete continuity equation in the AQMs, wind or density fields from MMs (or both) must be adjusted. Lu et al. (1997) prefer adjusting the density field, while Odman and Russell (2000) adjust either the vertical wind component or the vertical flux. Odman and Russell (2000) demonstrated that meteorological data could be made consistent in AQMs with species masses conserved at the same time. Also they found that the adjustments of wind fields have
insignificant impact on tracer transport. On the other hand, the renormalization approach could greatly increase or decrease the apparent emissions from single sources. As more applications use CMAQ for regional air quality modeling to study the relationships between receptors and precursor emissions at specific sources (Hu et al., 2004; Cohan et al., 2005), a mass conservative approach becomes more necessary. In this paper we examine the mass conservation errors in CMAQ by simulating tracer experiments, and evaluating the abilities of the two approaches (i.e., renormalization using air density and adjustment of vertical winds) in assuring mass conservation. 2. Approach An imaginary tracer experiment was simulated with CMAQ over southeastern United States using the meteorological data generated for the Fall-line Air Quality Study (FAQS) by Hu et al. (2004). In this experiment, inert tracers are released from four different locations in the Tennessee Valley (Fig. 1). These locations were chosen to represent major point sources of nitrogen oxides (NOx) emissions with greater potential of transport over complex terrain (Table 1). All releases start at 8:00 UTC and end at 14:00 UTC on 11 August 2000. The total
Fig. 1. Modeling domain (southeastern United States) with terrain heights shown as contour lines. Tracer release locations in Tennessee Valley are marked as dots. Also shown are the wind trajectories started at 8:00 UTC on 11 August 2000 from effective stack heights and calculated using the MM5 fields.
ARTICLE IN PRESS Y. Hu et al. / Atmospheric Environment 40 (2006) 1199–1204
1201
Table 1 Tracer release locations and effective stack heights Tracer name
Tracer-1 Tracer-2 Tracer-3 Tracer-4
Release location
Latitude
Longitude
36.0208 36.3767 35.8989 36.5063
84.1569 82.9639 84.5203 82.2538
Effective stack height (m)
Facility
NOx emissions in 2000 (Gg N year1)
807.9 377.7 1036.7 199.5
Plant Plant Plant Plant
4.3 3.0 7.0 2.9
A, Clinton, TN B, Rogersville, TN C, Harriman, TN D, Kingsport, TN
Also listed are the corresponding facilities and their NOx emissions in the year 2000.
amount of inert tracers released at each location is 24 kg. Starting two hours before the release of inert tracers and using the meteorological fields generated by MM5, a 42 h period is simulated using CMAQ following only advection and diffusion processes. Trajectories calculated by starting from effective stack heights for each tracer and using the MM5 fields are shown in Fig. 1. The CMAQ grid covers the horizontal domain shown in Fig. 1 with 75 66 12-km cells. It has 13 layers vertically, with 7 layers in the lowest kilometer and the surface layer having a thickness of 18 m. We follow the tracer concentrations and calculate the total mass of each tracer within the entire domain at every time step during the simulation. Note that while the release locations are close to the northern boundary of the domain, the winds are mostly northerly during the simulation period. This maximizes the residence time of the tracers within the domain. The same experiment was simulated with three different versions of CMAQ. First, we used version 4.3 (the ‘‘original’’ version). Note that CMAQ uses the piecewise parabolic method (PPM) (Collela and Woodward, 1984) both for horizontal and vertical advection. Second, we used version 4.4, which includes a correction to the implementation of PPM (the ‘‘corrected’’ version). In early versions, including 4.3, uniform grid spacing was assumed despite the non-uniform vertical layer spacing, leading to errors in vertical advection. Both of these versions of CMAQ applied the renormalization approach using air density to assure the consistency of the meteorological fields. In the third version, which we call the ‘‘mass conservative’’ version, we applied the inverse donor cell method (Odman and Russell, 2000) to satisfy the discrete continuity equation by adjusting the vertical winds. During the simulation with the ‘‘mass conservative’’ version, we
also tracked the adjustments to the vertical winds to check if there were significant deviations from the original winds. Both the renormalization method and the inverse donor cell method are briefly described below. In the renormalization method, air density is advected in addition to all the pollutant species and the ratio of the interpolated air density to the air density after advection is used to correct the advected species concentration (Byun, 1999): ccor ¼
rint T c . rT
(1)
Here ccor is the corrected species concentration, cT the species concentration after advection, rT the air density after advection and rint the air density linearly interpolated in time from the MM data. In CMAQ model’s advection modules, both c and r are coupled with the Jacobian of the coordinate transformation. The inverse donor cell method described in Odman and Russell (2000) requires advection of air density too, and, in addition, the use of the donor cell scheme for vertical advection. First, horizontal wind components are applied to the air density field, rn , to obtain an intermediate field r . Then, it is required that the air density rnþ1 be equal to rint after the vertical advection operator is applied to r . To satisfy this requirement the vertical wind must be adjusted, and the new vertical wind that would yield the desired air density is calculated as 8 B > > if BX0; > < r k (2) wkþ1=2 ¼ B > > if Bo0; > : r k1
ARTICLE IN PRESS Y. Hu et al. / Atmospheric Environment 40 (2006) 1199–1204
1202
where wk1=2 þ jwk1=2 j Dzk ðrk rint k Þþ 2 Dt wk1=2 jwk1=2 j rk . rk1 þ 2 Since w1=2 ¼ 0, Eq. (2) can be used starting from the surface layer ðk ¼ 1Þ and moving to the top, each time using the vertical wind computed for the previous layer. The winds at the top of the model will usually have a vertical component. The vertical velocities computed from Eq. (2) must be used in Eq. (3) below to advect all other pollutant species and obtain updated species concentrations. B¼
Dt ½ðwkþ1=2 jwkþ1=2 jÞckþ1 2Dzk ðwk1=2 jwk1=2 jÞck
cnþ1 ¼ ck k
þ ðwkþ1=2 þ jwkþ1=2 jÞck ðwk1=2 þ jwk1=2 jÞck1 .
ð3Þ
Note that the horizontal advection operator must always be applied before vertical advection. In CMAQ, the horizontal advection operator is split into two one-dimensional advection operators. When operator splitting is used, the density can change artificially even in non-divergent flows (Walcek and Aleksic, 1998). This may place unnecessary burden on the vertical wind adjustment resulting in a much larger deviation from the original velocity field. The deviation can be reduced by using a two-dimensional horizontal advection operator (e.g., Odman and Russell, 1991) or correcting the second one-dimensional advection operator for ‘‘cell depletion’’ that may occur after the application of the first one-dimensional operator (e.g., Walcek and Aleksic, 1998). Finally, note that since the donor cell scheme is subject to the Courant stability condition, the case when the Courant number exceeds unity requires special treatment. 3. Results and discussion If the model is mass conservative, the total mass of each inert tracer simulated in our test cases should remain constant at 24 kg after 14:00 UTC on August 11 until it starts to decrease when that tracer is transported out of the domain. However, the original version of CMAQ is artificially creating mass (Fig. 2). The amount of erroneous increase in tracer mass is large, especially for Tracer-4 whose total mass was artificially increased by about 5 folds
(Fig. 2(d)). Note that Tracer-4 encounters mountainous terrain right after its release (Fig. 1). The second largest mass conservation error is observed for Tracer-2 which also travels over very complex terrain but, after a short initial period, over more flat terrain. After correcting vertical advection, the mass conservation errors are reduced significantly in the simulation using the corrected version of CMAQ. However, mass is still artificially increased after an initial decrease (Figs. 2(a)–(c)). The largest artificial increase produced by the corrected version of CMAQ was 7% for Tracer-2 and the largest decrease was 8% for Tracer-3. These artificial fluctuations would misrepresent the NOx emissions shown in Table 1 for Plant B (Rogersville, TN) and Plant C (Harriman, TN), respectively, by +0.22 Gg N year1 and 0.58 Gg N year1 in the year 2000. Such large misrepresentation of emissions may have significant impacts on simulated responses to emission changes, i.e., for testing control strategies. In contrast to the fluctuations in the mass simulated using the corrected version of CMAQ, the mass conservative CMAQ with vertical wind adjustment kept total mass constant before the tracers were transported out of the domain either from the top or through lateral boundaries (Fig. 2). At the same time, the adjustments made to the vertical winds were small: the mean difference between adjusted and unadjusted vertical winds was 6.9 104 m s1 over the domain throughout the simulation, while the mean of unadjusted and adjusted vertical wind components are 2.4 103 and 3.0 103 m s1, respectively. Table 2 shows the distance between the trajectories calculated by using unadjusted and adjusted wind fields, respectively. Note that the trajectory end positions are in the same grid cells for all tracers after 12 h. After 24 h, the unadjusted and adjusted wind trajectories for Tracer-1 still end in the same grid cell but those for other tracers differ by one grid cell (diagonally in the case of Tracer-2). These results show that adjusting vertical winds introduce a small impact on tracer transport yet it assures strict mass conservation in CMAQ. 4. Conclusion Simulations with CMAQ prior to correcting the error in vertical advection may be subject to severe mass conservation errors over complex terrain. Using those versions of CMAQ (i.e., version 4.3 and earlier) to study source-receptor relationships
ARTICLE IN PRESS Y. Hu et al. / Atmospheric Environment 40 (2006) 1199–1204 30
40 35 Total Mass (Kg)
25 Total Mass (Kg)
1203
20 15 10
30 25 20 15 10
5
5 0 8/11/00 8/11/00 8/11/00 8/12/00 8/12/00 8/12/00 8/12/00 8/13/00
0 8/11/00 8/11/00 8/11/00 8/12/00 8/12/00 8/12/00 8/12/00 8/13/00 6:00
12:00
18:00
(a)
0:00
6:00
12:00
18:00
0:00
6:00
12:00
18:00
0:00
(b)
Time (UTC)
6:00
12:00
18:00
0:00
Time (UTC) 160
30
140 Total Mass (Kg)
Total Mass (Kg)
25 20 15 10
120 100 80 60 40
5
20
0 8/11/00 8/11/00 8/11/00 8/12/00 8/12/00 8/12/00 8/12/00 8/13/00 0:00 6:00 6:00 12:00 18:00 12:00 18:00 0:00 Time (UTC)
(c)
0 8/11/00 8/11/00 8/11/00 8/12/00 8/12/00 8/12/00 8/12/00 8/13/00 6:00
12:00
18:00
(d)
0:00
6:00
12:00
18:00
0:00
Time (UTC)
Fig. 2. Evolution of total mass of Tracer-1 (a), Tracer-2 (b), Tracer-3 (c) and Tracer-4 (d), as simulated using the original (in black), corrected (in gray), and mass conservative (in light gray) versions of CMAQ. Along the time series, the linear ramp-up between 8:00 and 14:00 UTC on 11 August 2000 is from increased mass during tracer release, while the gradual decrease later results from mass leaving the domain.
Table 2 Trajectory positions after 12 and 24 h, calculated using unadjusted and adjusted wind fields Starting point of trajectory
Trajectory age (h)
Trajectory position calculated using unadjusted winds
Trajectory position calculated using adjusted winds
Column
Row
Layer
Column
Row
Layer
Distance apart (m)
Release location of Tracer-1
12 24
33 46
35 23
6 9
33 46
35 23
6 9
202 4739
Release location of Tracer-2
12 24
39 48
46 42
8 9
39 47
46 43
8 9
474 13,724
Release location of Tracer-3
12 24
33 49
35 22
7 9
33 48
35 22
7 9
1871 7786
Release location of Tracer-4
12 24
46 45
48 36
4 9
46 44
48 36
4 9
2604 3350
All trajectories are started at 8:00 UTC on 11 August 2000.
ARTICLE IN PRESS 1204
Y. Hu et al. / Atmospheric Environment 40 (2006) 1199–1204
for the purpose of designing emission control strategies may lead to over/underestimates of impacts. The correction (i.e., version 4.4) reduces mass conservation errors but it does not completely eliminate them since the renormalization method, though it can preserve uniform mixing ratio fields, does not assure conservation of pollutant mass, suggesting continued issues with the application for source-specific analyses. Adjusting vertical winds proved to be an effective method to conserve mass in CMAQ. Acknowledgements This work was funded by the Visibility Improvement State and Tribal Association of the Southeast (VISTAS) under contract number V-2003-12, Environmental Protection Division of the Georgia Department of Natural Resources through the Fallline Air Quality Study, and the U.S. Environmental Protection Agency under Grants R-82897601, RD83096001, and RD83107601. We thank Dr. Peter Percell and Dr. Daewon W. Byun of the University of Houston for providing the corrected CMAQ code. References Byun, D.W., 1999. Dynamically consistent formulations in meteorological and air quality models for multiscale atmospheric studies. Part II: mass conservation issues. Journal of the Atmospheric Sciences 56, 3808–3820.
Byun, D.W., Ching, J.K.S., 1999. Science algorithms of the EPA Models-3 Community Multiscale Air Quality (CMAQ) Modeling System, U.S. EPA/600/R-99/030. Cohan, D.S., Hakami, A., Hu, Y., Russell, A.G., 2005. Nonlinear response of ozone to emissions: source apportionment and sensitivity analysis. Environmental Science and Technology 39, 6739–6748. Collela, P., Woodward, P.R., 1984. The piecewise parabolicmethod (PPM) for gas-dynamical simulations. Journal of Computational Physics 54, 174–201. Grell, G., Dudhia, J., Stauffer, D., 1994. A description of the Fifth-Generation Penn State/NCAR Mesoscale Model (MM5), NCAR Technical Note: NCAR/TN-398+STR. Hu, Y., Cohan, D., Odman, M.T., Russell, A.G., 2004. Air Quality Modeling of the August 11–20, 2000 Episode for the Fall Line Air Quality Study. Atlanta, GA, Prepared for Georgia Department of Natural Resources, Environmental Protection Division. Lee, S.M., Yoon, S.C., Byun, D.W., 2004. The effect of mass inconsistency of the meteorological field generated by a common meteorological model on air quality modeling. Atmospheric Environment 38, 2917–2926. Lu, R., Turco, R.P., Jacobson, M.Z., 1997. An integrated air pollution modeling system for urban and regional scales: 1. Structure and performance. Journal of Geophysical Research 102, 6063. Odman, M.T., Russell, A.G., 1991. A multiscale finite element pollutant transport scheme for urban and regional modeling. Atmospheric Environment 25A, 2385–2394. Odman, M.T., Russell, A.G., 2000. Mass conservative coupling of non-hydrostatic meteorological models with air quality models. In: Gryning, S.-E., Batchvarova, E. (Eds.), Air Pollution Modeling and its Application XIII. Kluwer/ Plenum, New York, pp. 651–660. Walcek, C.J., Aleksic, N.M., 1998. A simple but accurate mass conservative, peak-preserving, mixing ratio bounded advection algorithm with FORTRAN code. Atmospheric Environment 32 (22), 3863–3880.
Atmospheric Environment 40 (2006) 1199–1204 www.elsevier.com/locate/atmosenv
Mass conservation in the Community Multiscale Air Quality model Yongtao Hu, M. Talat Odman, Armistead G. Russell School of Civil and Environmental Engineering, Georgia Institute of Technology, Atlanta, GA 30332, USA Received 17 February 2005; received in revised form 15 October 2005; accepted 15 October 2005
Abstract Mass conservation errors were discovered when the Community Multiscale Air Quality (CMAQ) model (version 4.3) was coupled with the Fifth-Generation PSU/NCAR Mesoscale Model (MM5). These errors can be severe particularly over complex terrain. A recent correction made to the vertical advection module (in version 4.4) reduces the mass conservation errors, but does not completely eliminate them. This means that the method of renormalizing concentrations according to air density, which is currently applied to address the inconsistency problem, is not an ideal solution. Reconstruction of the wind fields by adjusting the vertical wind is proposed as a simple remedy to make CMAQ strictly mass conservative. r 2005 Elsevier Ltd. All rights reserved. Keywords: Regional air quality modeling; Inconsistency of meteorological fields; Mass conservation error; Source-receptor relationship; Emissions control strategy
1. Introduction Air quality models (AQMs) such as the Community Multiscale Air Quality model (CMAQ) (Byun and Ching, 1999) are expected to maintain a uniform mass mixing ratio field for an inert tracer after transport with the winds produced by meteorological models (MMs), such as the Fifth-Generation PSU/NCAR Mesoscale Model (MM5) (Grell et al., 1994). This expectation is realized only if the two models use the same discretization, i.e., grid, time step, and finite difference forms (Odman and Russell, 2000). In an integrated air quality modeling system, e.g. SMOG (Lu et al., 1997), the meteorolCorresponding author. Tel.: +1 404 894 1854; fax: +1 494 894 8266. E-mail address: [email protected] (Y. Hu).
1352-2310/$ - see front matter r 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.atmosenv.2005.10.038
ogy component and chemistry/transport (i.e., air quality) component naturally share the same grid structure. However, in several present-day air quality modeling applications, CMAQ and MM5 do not share the same grid structure and spatial interpolations are usually required. Also, when AQMs and MMs are off-line coupled, the outputs of the MMs are stored less frequently than the AQMs’ time step. For example, MM5 outputs are usually stored hourly while CMAQ has a maximum transport time step of 900 s. Input meteorological variables are reconstructed for AQMs’ time intervals, usually by linear interpolation between the stored values. Time step differences may also occur in on-line coupled integrated air quality modeling systems. For example, in SMOG, the dynamical equations are solved at a time step of 6 s while tracer transport is calculated every 300 s, accordingly,
ARTICLE IN PRESS 1200
Y. Hu et al. / Atmospheric Environment 40 (2006) 1199–1204
parameters predicted by the MM are averaged over the appropriate number of dynamical time steps to obtain values for the AQM (Lu et al., 1997). Finally, transport algorithms (e.g. finite difference forms) used for advection in AQMs are usually different from those in MMs. Because of all these differences in discretization, the wind components and the air density do not necessarily satisfy the discrete continuity equation in AQMs even though they may have done so in MMs. This is referred to as ‘‘the inconsistency of the meteorological fields’’ (Byun, 1999; Odman and Russell, 2000; Lee et al., 2004). As a result of this inconsistency, uniform mass mixing ratio fields are perturbed. The perturbation is usually more pronounced with data from non-hydrostatic MMs than with data from hydrostatic MMs for the same domains (Odman and Russell, 2000; Lee et al., 2004). To solve this inconsistency problem, some widely used AQMs, including CMAQ, adopted an approach where a uniform tracer is transported along with other species and the concentrations are renormalized using the deviation of the tracer from uniformity. In CMAQ, this uniform tracer is air density (Byun, 1999). As a result, a uniform mixing ratio field is maintained (Lee et al., 2004). However, the approach overlooks that non-linear advection schemes are used in AQMs for high accuracy. The non-linearities affect pollutant fields of different distribution differently. Most species of interest, especially the emitted ones, have spatial distributions far from being uniform. Therefore, renormalizing the concentration of a species based on the perturbation of a uniform field can artificially increase or decrease mass. Odman and Russell (2000) showed that this attempt to establish consistency and produce stable results in AQMs sacrifices mass conservation of pollutant species and found that the mass conservation errors introduced by this approach could be very large. Species mass conservation is not sacrificed if the consistency of meteorological fields is re-established through the continuity equation. To satisfy the discrete continuity equation in the AQMs, wind or density fields from MMs (or both) must be adjusted. Lu et al. (1997) prefer adjusting the density field, while Odman and Russell (2000) adjust either the vertical wind component or the vertical flux. Odman and Russell (2000) demonstrated that meteorological data could be made consistent in AQMs with species masses conserved at the same time. Also they found that the adjustments of wind fields have
insignificant impact on tracer transport. On the other hand, the renormalization approach could greatly increase or decrease the apparent emissions from single sources. As more applications use CMAQ for regional air quality modeling to study the relationships between receptors and precursor emissions at specific sources (Hu et al., 2004; Cohan et al., 2005), a mass conservative approach becomes more necessary. In this paper we examine the mass conservation errors in CMAQ by simulating tracer experiments, and evaluating the abilities of the two approaches (i.e., renormalization using air density and adjustment of vertical winds) in assuring mass conservation. 2. Approach An imaginary tracer experiment was simulated with CMAQ over southeastern United States using the meteorological data generated for the Fall-line Air Quality Study (FAQS) by Hu et al. (2004). In this experiment, inert tracers are released from four different locations in the Tennessee Valley (Fig. 1). These locations were chosen to represent major point sources of nitrogen oxides (NOx) emissions with greater potential of transport over complex terrain (Table 1). All releases start at 8:00 UTC and end at 14:00 UTC on 11 August 2000. The total
Fig. 1. Modeling domain (southeastern United States) with terrain heights shown as contour lines. Tracer release locations in Tennessee Valley are marked as dots. Also shown are the wind trajectories started at 8:00 UTC on 11 August 2000 from effective stack heights and calculated using the MM5 fields.
ARTICLE IN PRESS Y. Hu et al. / Atmospheric Environment 40 (2006) 1199–1204
1201
Table 1 Tracer release locations and effective stack heights Tracer name
Tracer-1 Tracer-2 Tracer-3 Tracer-4
Release location
Latitude
Longitude
36.0208 36.3767 35.8989 36.5063
84.1569 82.9639 84.5203 82.2538
Effective stack height (m)
Facility
NOx emissions in 2000 (Gg N year1)
807.9 377.7 1036.7 199.5
Plant Plant Plant Plant
4.3 3.0 7.0 2.9
A, Clinton, TN B, Rogersville, TN C, Harriman, TN D, Kingsport, TN
Also listed are the corresponding facilities and their NOx emissions in the year 2000.
amount of inert tracers released at each location is 24 kg. Starting two hours before the release of inert tracers and using the meteorological fields generated by MM5, a 42 h period is simulated using CMAQ following only advection and diffusion processes. Trajectories calculated by starting from effective stack heights for each tracer and using the MM5 fields are shown in Fig. 1. The CMAQ grid covers the horizontal domain shown in Fig. 1 with 75 66 12-km cells. It has 13 layers vertically, with 7 layers in the lowest kilometer and the surface layer having a thickness of 18 m. We follow the tracer concentrations and calculate the total mass of each tracer within the entire domain at every time step during the simulation. Note that while the release locations are close to the northern boundary of the domain, the winds are mostly northerly during the simulation period. This maximizes the residence time of the tracers within the domain. The same experiment was simulated with three different versions of CMAQ. First, we used version 4.3 (the ‘‘original’’ version). Note that CMAQ uses the piecewise parabolic method (PPM) (Collela and Woodward, 1984) both for horizontal and vertical advection. Second, we used version 4.4, which includes a correction to the implementation of PPM (the ‘‘corrected’’ version). In early versions, including 4.3, uniform grid spacing was assumed despite the non-uniform vertical layer spacing, leading to errors in vertical advection. Both of these versions of CMAQ applied the renormalization approach using air density to assure the consistency of the meteorological fields. In the third version, which we call the ‘‘mass conservative’’ version, we applied the inverse donor cell method (Odman and Russell, 2000) to satisfy the discrete continuity equation by adjusting the vertical winds. During the simulation with the ‘‘mass conservative’’ version, we
also tracked the adjustments to the vertical winds to check if there were significant deviations from the original winds. Both the renormalization method and the inverse donor cell method are briefly described below. In the renormalization method, air density is advected in addition to all the pollutant species and the ratio of the interpolated air density to the air density after advection is used to correct the advected species concentration (Byun, 1999): ccor ¼
rint T c . rT
(1)
Here ccor is the corrected species concentration, cT the species concentration after advection, rT the air density after advection and rint the air density linearly interpolated in time from the MM data. In CMAQ model’s advection modules, both c and r are coupled with the Jacobian of the coordinate transformation. The inverse donor cell method described in Odman and Russell (2000) requires advection of air density too, and, in addition, the use of the donor cell scheme for vertical advection. First, horizontal wind components are applied to the air density field, rn , to obtain an intermediate field r . Then, it is required that the air density rnþ1 be equal to rint after the vertical advection operator is applied to r . To satisfy this requirement the vertical wind must be adjusted, and the new vertical wind that would yield the desired air density is calculated as 8 B > > if BX0; > < r k (2) wkþ1=2 ¼ B > > if Bo0; > : r k1
ARTICLE IN PRESS Y. Hu et al. / Atmospheric Environment 40 (2006) 1199–1204
1202
where wk1=2 þ jwk1=2 j Dzk ðrk rint k Þþ 2 Dt wk1=2 jwk1=2 j rk . rk1 þ 2 Since w1=2 ¼ 0, Eq. (2) can be used starting from the surface layer ðk ¼ 1Þ and moving to the top, each time using the vertical wind computed for the previous layer. The winds at the top of the model will usually have a vertical component. The vertical velocities computed from Eq. (2) must be used in Eq. (3) below to advect all other pollutant species and obtain updated species concentrations. B¼
Dt ½ðwkþ1=2 jwkþ1=2 jÞckþ1 2Dzk ðwk1=2 jwk1=2 jÞck
cnþ1 ¼ ck k
þ ðwkþ1=2 þ jwkþ1=2 jÞck ðwk1=2 þ jwk1=2 jÞck1 .
ð3Þ
Note that the horizontal advection operator must always be applied before vertical advection. In CMAQ, the horizontal advection operator is split into two one-dimensional advection operators. When operator splitting is used, the density can change artificially even in non-divergent flows (Walcek and Aleksic, 1998). This may place unnecessary burden on the vertical wind adjustment resulting in a much larger deviation from the original velocity field. The deviation can be reduced by using a two-dimensional horizontal advection operator (e.g., Odman and Russell, 1991) or correcting the second one-dimensional advection operator for ‘‘cell depletion’’ that may occur after the application of the first one-dimensional operator (e.g., Walcek and Aleksic, 1998). Finally, note that since the donor cell scheme is subject to the Courant stability condition, the case when the Courant number exceeds unity requires special treatment. 3. Results and discussion If the model is mass conservative, the total mass of each inert tracer simulated in our test cases should remain constant at 24 kg after 14:00 UTC on August 11 until it starts to decrease when that tracer is transported out of the domain. However, the original version of CMAQ is artificially creating mass (Fig. 2). The amount of erroneous increase in tracer mass is large, especially for Tracer-4 whose total mass was artificially increased by about 5 folds
(Fig. 2(d)). Note that Tracer-4 encounters mountainous terrain right after its release (Fig. 1). The second largest mass conservation error is observed for Tracer-2 which also travels over very complex terrain but, after a short initial period, over more flat terrain. After correcting vertical advection, the mass conservation errors are reduced significantly in the simulation using the corrected version of CMAQ. However, mass is still artificially increased after an initial decrease (Figs. 2(a)–(c)). The largest artificial increase produced by the corrected version of CMAQ was 7% for Tracer-2 and the largest decrease was 8% for Tracer-3. These artificial fluctuations would misrepresent the NOx emissions shown in Table 1 for Plant B (Rogersville, TN) and Plant C (Harriman, TN), respectively, by +0.22 Gg N year1 and 0.58 Gg N year1 in the year 2000. Such large misrepresentation of emissions may have significant impacts on simulated responses to emission changes, i.e., for testing control strategies. In contrast to the fluctuations in the mass simulated using the corrected version of CMAQ, the mass conservative CMAQ with vertical wind adjustment kept total mass constant before the tracers were transported out of the domain either from the top or through lateral boundaries (Fig. 2). At the same time, the adjustments made to the vertical winds were small: the mean difference between adjusted and unadjusted vertical winds was 6.9 104 m s1 over the domain throughout the simulation, while the mean of unadjusted and adjusted vertical wind components are 2.4 103 and 3.0 103 m s1, respectively. Table 2 shows the distance between the trajectories calculated by using unadjusted and adjusted wind fields, respectively. Note that the trajectory end positions are in the same grid cells for all tracers after 12 h. After 24 h, the unadjusted and adjusted wind trajectories for Tracer-1 still end in the same grid cell but those for other tracers differ by one grid cell (diagonally in the case of Tracer-2). These results show that adjusting vertical winds introduce a small impact on tracer transport yet it assures strict mass conservation in CMAQ. 4. Conclusion Simulations with CMAQ prior to correcting the error in vertical advection may be subject to severe mass conservation errors over complex terrain. Using those versions of CMAQ (i.e., version 4.3 and earlier) to study source-receptor relationships
ARTICLE IN PRESS Y. Hu et al. / Atmospheric Environment 40 (2006) 1199–1204 30
40 35 Total Mass (Kg)
25 Total Mass (Kg)
1203
20 15 10
30 25 20 15 10
5
5 0 8/11/00 8/11/00 8/11/00 8/12/00 8/12/00 8/12/00 8/12/00 8/13/00
0 8/11/00 8/11/00 8/11/00 8/12/00 8/12/00 8/12/00 8/12/00 8/13/00 6:00
12:00
18:00
(a)
0:00
6:00
12:00
18:00
0:00
6:00
12:00
18:00
0:00
(b)
Time (UTC)
6:00
12:00
18:00
0:00
Time (UTC) 160
30
140 Total Mass (Kg)
Total Mass (Kg)
25 20 15 10
120 100 80 60 40
5
20
0 8/11/00 8/11/00 8/11/00 8/12/00 8/12/00 8/12/00 8/12/00 8/13/00 0:00 6:00 6:00 12:00 18:00 12:00 18:00 0:00 Time (UTC)
(c)
0 8/11/00 8/11/00 8/11/00 8/12/00 8/12/00 8/12/00 8/12/00 8/13/00 6:00
12:00
18:00
(d)
0:00
6:00
12:00
18:00
0:00
Time (UTC)
Fig. 2. Evolution of total mass of Tracer-1 (a), Tracer-2 (b), Tracer-3 (c) and Tracer-4 (d), as simulated using the original (in black), corrected (in gray), and mass conservative (in light gray) versions of CMAQ. Along the time series, the linear ramp-up between 8:00 and 14:00 UTC on 11 August 2000 is from increased mass during tracer release, while the gradual decrease later results from mass leaving the domain.
Table 2 Trajectory positions after 12 and 24 h, calculated using unadjusted and adjusted wind fields Starting point of trajectory
Trajectory age (h)
Trajectory position calculated using unadjusted winds
Trajectory position calculated using adjusted winds
Column
Row
Layer
Column
Row
Layer
Distance apart (m)
Release location of Tracer-1
12 24
33 46
35 23
6 9
33 46
35 23
6 9
202 4739
Release location of Tracer-2
12 24
39 48
46 42
8 9
39 47
46 43
8 9
474 13,724
Release location of Tracer-3
12 24
33 49
35 22
7 9
33 48
35 22
7 9
1871 7786
Release location of Tracer-4
12 24
46 45
48 36
4 9
46 44
48 36
4 9
2604 3350
All trajectories are started at 8:00 UTC on 11 August 2000.
ARTICLE IN PRESS 1204
Y. Hu et al. / Atmospheric Environment 40 (2006) 1199–1204
for the purpose of designing emission control strategies may lead to over/underestimates of impacts. The correction (i.e., version 4.4) reduces mass conservation errors but it does not completely eliminate them since the renormalization method, though it can preserve uniform mixing ratio fields, does not assure conservation of pollutant mass, suggesting continued issues with the application for source-specific analyses. Adjusting vertical winds proved to be an effective method to conserve mass in CMAQ. Acknowledgements This work was funded by the Visibility Improvement State and Tribal Association of the Southeast (VISTAS) under contract number V-2003-12, Environmental Protection Division of the Georgia Department of Natural Resources through the Fallline Air Quality Study, and the U.S. Environmental Protection Agency under Grants R-82897601, RD83096001, and RD83107601. We thank Dr. Peter Percell and Dr. Daewon W. Byun of the University of Houston for providing the corrected CMAQ code. References Byun, D.W., 1999. Dynamically consistent formulations in meteorological and air quality models for multiscale atmospheric studies. Part II: mass conservation issues. Journal of the Atmospheric Sciences 56, 3808–3820.
Byun, D.W., Ching, J.K.S., 1999. Science algorithms of the EPA Models-3 Community Multiscale Air Quality (CMAQ) Modeling System, U.S. EPA/600/R-99/030. Cohan, D.S., Hakami, A., Hu, Y., Russell, A.G., 2005. Nonlinear response of ozone to emissions: source apportionment and sensitivity analysis. Environmental Science and Technology 39, 6739–6748. Collela, P., Woodward, P.R., 1984. The piecewise parabolicmethod (PPM) for gas-dynamical simulations. Journal of Computational Physics 54, 174–201. Grell, G., Dudhia, J., Stauffer, D., 1994. A description of the Fifth-Generation Penn State/NCAR Mesoscale Model (MM5), NCAR Technical Note: NCAR/TN-398+STR. Hu, Y., Cohan, D., Odman, M.T., Russell, A.G., 2004. Air Quality Modeling of the August 11–20, 2000 Episode for the Fall Line Air Quality Study. Atlanta, GA, Prepared for Georgia Department of Natural Resources, Environmental Protection Division. Lee, S.M., Yoon, S.C., Byun, D.W., 2004. The effect of mass inconsistency of the meteorological field generated by a common meteorological model on air quality modeling. Atmospheric Environment 38, 2917–2926. Lu, R., Turco, R.P., Jacobson, M.Z., 1997. An integrated air pollution modeling system for urban and regional scales: 1. Structure and performance. Journal of Geophysical Research 102, 6063. Odman, M.T., Russell, A.G., 1991. A multiscale finite element pollutant transport scheme for urban and regional modeling. Atmospheric Environment 25A, 2385–2394. Odman, M.T., Russell, A.G., 2000. Mass conservative coupling of non-hydrostatic meteorological models with air quality models. In: Gryning, S.-E., Batchvarova, E. (Eds.), Air Pollution Modeling and its Application XIII. Kluwer/ Plenum, New York, pp. 651–660. Walcek, C.J., Aleksic, N.M., 1998. A simple but accurate mass conservative, peak-preserving, mixing ratio bounded advection algorithm with FORTRAN code. Atmospheric Environment 32 (22), 3863–3880.