An investigation of the sub-grid variability of trace gases and aerosols for global climate modeling

One fundamental property and limitation of grid based models is their inability to identify spatial details smaller than the grid cell size. While decades of work have gone into developing sub-grid treatments for clouds and land surface processes in climate models, the quantitative understanding of sub-grid processes and variability for aerosols and their precursors is much poorer. In this study, WRFChem is used to simulate the trace gases and aerosols over central Mexico during the 2006 MILAGRO field campaign, with multiple spatial resolutions and emission/terrain scenarios. Our analysis focuses on quantifying the sub-grid variability (SGV) of trace gases and aerosols within a typical global climate model grid cell, i.e. 75 ×75 km2. Our results suggest that a simulation with 3-km horizontal grid spacing adequately reproduces the overall transport and mixing of trace gases and aerosols downwind of Mexico City, while 75-km horizontal grid spacing is insufficient to represent local emission and terrain-induced flows along the mountain ridge, subsequently affecting the transport and mixing of plumes from nearby sources. Therefore, the coarse model grid cell average may not correctly represent aerosol properties measured over polluted areas. Probability density functions (PDFs) for trace gases and aerosols show that secondary trace gases and aerosols, such as O 3, sulfate, ammonium, and nitrate, are more likely to have a relatively uniform probability distribution (i.e. smaller SGV) over a narrow range of concentration values. Mostly inert and longlived trace gases and aerosols, such as CO and BC, are more likely to have broad and skewed distributions (i.e. larger SGV) over polluted regions. Over remote areas, all trace gases and aerosols are more uniformly distributed compared to polluted areas. Both CO and O 3 SGV vertical profiles are nearly constant within the PBL during daytime, indicating Correspondence to: Y. Qian (yun.qian@pnl.gov) that trace gases are very efficiently transported and mixed vertically by turbulence. But, simulated horizontal variability indicates that trace gases and aerosols are not well mixed horizontally in the PBL. During nighttime the SGV for trace gases is maximum at the surface, and quickly decreases with height. Unlike the trace gases, the SGV of BC and secondary aerosols reaches a maximum at the PBL top during the day. The SGV decreases with distance away from the polluted urban area, has a more rapid decrease for long-lived trace gases and aerosols than for secondary ones, and is greater during daytime than nighttime. The SGV of trace gases and aerosols is generally larger than for meteorological quantities. Emissions can account for up to 50% of the SGV over urban areas such as Mexico City during daytime for less-reactive trace gases and aerosols, such as CO and BC. The impact of emission spatial variability on SGV decays with altitude in the PBL and is insignificant in the free troposphere. The emission variability affects SGV more significantly during daytime (rather than nighttime) and over urban (rather than rural or remote) areas. The terrain, through its impact on meteorological fields such as wind and the PBL structure, affects dispersion and transport of trace gases and aerosols and their SGV.


Introduction
One fundamental property and limitation of all Eulerian models is their inability to identify spatial details smaller than the grid cell size, known as sub-grid variability (SGV).SGV is present for meteorological variables as well as trace gases and aerosols, even when very small grid spacings are employed (Haywood et al., 1997;Karamchandani et al., 2002;Ching et al., 2006).For weather and climate models, with grid spacing ranging from a few to hundreds of kilometers, the SGV of meteorological variables arises from withingrid variability of terrain, land surface type and properties, Published by Copernicus Publications on behalf of the European Geosciences Union.Cloud-induced variability 1 km to 100 km, e.g.shallow cumulus to stratus and organized synoptic systems and clouds.For example, Leung and Qian (2003) found that spatial resolution has a significant impact on simulation results over mountainous areas in the western US.Higher spatial resolution improves the simulation, especially for orographic precipitation and snowpack, due to better reproduction of the temperature gradient by resolving the complex terrain and mesoscale forcings.In the case of complex terrain, one cannot simply apply a bias correction to account for errors in the model terrain caused by smoothing or the effects of the sub-grid terrain variability itself.
For chemistry models, the SGV of trace gases and aerosols results from both the traditional sub-grid processes affecting meteorology and specific chemistry processes, such as emissions, chemical transformation, and removal.In most chemistry models, primary emissions are usually instantaneously and uniformly diluted over an entire grid cell volume.In the real world, the spatial distribution of anthropogenic emissions (e.g.SO 2 , NO x , CO, black carbon, organic matter) from point, area, and mobile sources are quite inhomogeneous and dissimilar within a model grid cell.Models that currently employ grid spacings of 1 to 100 km cannot resolve all the small-scale variations in anthropogenic, biomass burning, and biogenic emissions.The point-like nature of many emission sources makes them particularly difficult.In addition, the subsequent dispersion and mixing of trace gas and aerosol plumes in the horizontal and vertical dimensions occurs at highly variable rates.Different spatial resolutions may result in different predictions of secondary products such as ozone and sulfates, since many atmospheric chemistry processes are nonlinear and frequently diffusionlimited.Table 1 lists some processes that contribute to subgrid variability of trace gases and aerosols.Depending on the resolution of the model grid employed some or all of these, and other, processes can introduce model uncertainty and error if they are neglected.
Climate models in the Intergovernmental Panel on Climate Change Fourth Assessment Report typically employ a horizontal grid spacing on the order of 100 km.This is much larger than the width of pollution plumes near source regions, such as power plant stacks (e.g.Chapman et al., 2009), and also larger than the spatial scales of many of the other processes causing aerosol spatial variability, such as aqueous generation and scavenging of aerosol induced by cloud, and transport and mixing influenced by terrain features.While the reduction of grid spacing is foreseen in future modeling activities, the question of how much resolution is needed to accurately reproduce aerosol impacts on climate is not known.There are significant nonlinearities in the chemistry that arise from changing grid cell mean concentrations of gases and aerosols as a given amount of material is spread throughout a given grid cell.This issue has been known for years, but has not been adequately addressed, although there have been different attempts (e.g., Calbo et al., 1998;Mayer et al., 2000).
Another confounding factor when evaluating climate or chemistry models is the comparison of point measurements with model grid cell averages.This is considered a "change of support" problem in which inferences are made about differences between point-based measurements to modelpredicted values that represent volume average concentration (Gelfand et al., 2001).The quantitative comparisons of modeled concentrations of trace gases or aerosols with observations will change merely due to a different choice in the size of the grid cell chosen for the simulation.Moreover, any observation reflects an instantaneous event out of a population, while model predictions represent an average of the population during a time step.When SGV is significant, the comparison of grid model outputs against one or more point measurements can result in comparisons seeming worse than they really are (Ching et al., 2006;Touma et al., 2006).
Decades of work have gone into developing sub-grid treatments for clouds or land surface process in climate models (Slingo, 1980;Randall et al., 2003;Avissar and Pielke, 1989;Seth et al., 1994).It has come to be accepted that high resolution is needed to improve the handling of clouds in climate models (e.g.Shukla et al., 2009).The quantitative understanding for sub-grid processes and variability of aerosols and their precursors are much poorer.For example, is the same resolution needed for aerosols as for clouds to resolve the SGV?If aerosols could be modeled at a coarser resolution than the clouds, a significant savings could be had since many more variables must be stored and advected for the aerosols and the associated trace gas chemistry than are needed for the traditional meteorological fields.Quantification of subgrid variability of aerosols and better understanding of the extent to which different sub-grid-scale processes contribute to uncertainty in aerosols and their radiative forcing within climate models are needed to assist in developing parameterization schemes that could account for at least a portion of the neglected sub-grid aerosol processes in these models.
While a systematic approach has not yet been employed to document the impact of neglected subgrid aerosol variability, some work has been done to parameterize it.The most prominent of these attempts is the plume-in-grid concept, which was initially developed for the handling of ozone plumes within a grid cell (Karamchandani et al., 2002).The concept has later been extended to other constituents such as mercury (Vijayaraghavan et al., 2008) and particulate matter (Karamchandani et al., 2006).Another approach is the Explicit-Cloud Parameterized-Pollutant parameterization which is specifically designed for use in multiscale modeling framework (MMF) models (Gustafson et al., 2008).This latter technique uses statistics from high resolution cloud models embedded within each coarse GCM column to improve the treatment of vertical mixing and cloud processing of aerosols.
In this study, the chemistry version of the Weather Research and Forecasting model (WRF-Chem) is applied to simulate the aerosols and other trace gases over the vicinity of Mexico City during the 2006 Megacity Initiative: Local and Global Research Observations (MILAGRO) using multiple spatial resolutions and scenarios that examine the effect of SGV of emissions and terrain.Our analysis in this study focuses on quantifying the sub-grid variability of trace gases and aerosols within a typical global climate model (GCM) grid cell, i.e. 75 × 75 km.Nested domains with grid spacing representative of mesoscale models (i.e. 15 km), and cloud-system resolving models (i.e. 3 km) are used to identify how the simulated aerosol characteristics change with spatial scale.The paper is organized as follows: in Sect. 2 we briefly introduce the WRF-Chem model configuration, experiment design, and observational data.In Sect. 3 we evaluate WRF-Chem simulations at various spatial resolutions against observations, with the objective of quantifying the uncertainty caused by spatial variability for trace gases and aerosols when comparing point measurements to grid cell volumes at different spatial resolutions.In Sect. 4 we present the basic characteristics of the SGV of trace gases and aerosols, including their spatial pattern, diurnal and vertical variation, and dependences on the distance to emission sources and on the spatial resolution.In Sect.5, based on a series of sensitivity experiments, we analyze the factors affecting sub-grid scale processes, including local forcings such as topography and emission variability.Discussions and conclusions are presented in Sects.6 and 7, respectively.The results of this study improve our understanding of subgrid process of trace gases and aerosols and provide useful information guiding parameterization development designed to reduce the uncertainty in estimating the aerosol forcing of climate.

Model description
The non-hydrostatic Weather Research and Forecasting (WRF) community model includes various options for dynamic cores and physical parameterizations so that it can be used to simulate atmospheric processes over a wide range of spatial and temporal scales (Skamarock et al., 2008).WRF-Chem, the chemistry version of the WRF model (Grell et al., 2005), simulates trace gases and particulates interactively with the meteorological fields.WRF-Chem contains several treatments for photochemistry and aerosols developed by the user community.
The modules in WRF-Chem version 3 used in this study are: the CBM-Z gas-phase chemistry mechanism (Zaveri and Peters, 1999), the MOSAIC aerosol model that employs a sectional approach for the aerosol size distribution (Zaveri et al., 2008), and the Fast-J photolysis scheme (Wild et al., 2000).The aerosol direct effect is coupled to the Goddard shortwave scheme (Fast et al., 2006).The interactions between aerosols and clouds, such as the first and second indirect effects, activation/resuspension, wet scavenging, and aqueous chemistry (Gustafson et al., 2007;Chapman et al., 2009), are not turned on.Aerosol-cloud interactions were probably negligible prior to the cold surge on 23 March when mostly sunny conditions were observed over the central Mexican plateau (Fast et al., 2007).Prognostic species in MO-SAIC include sulfate, nitrate, ammonium, chloride, sodium, other (unspecified) inorganics, organic matter (OM), black carbon (BC), aerosol water, and aerosol number.Eight size bins are used for each aerosol specie.Aerosols are assumed to be internally mixed and volume-averaging is used to compute optical properties that influence radiation.It should be noted that no secondary organic aerosol (SOA) treatment is included in the version of MOSAIC used for this paper.The configuration of WRF-Chem is similar to Fast et al. (2007Fast et al. ( , 2009)), except that our study does not include observational data assimilation.While simulated plume locations would be in better agreement with the observations by including data assimilation, the forcing associated with data assimilation would vary among domains with different resolutions so that the results would not be directly comparable.We therefore chose to neglect data assimilation because it would confound the interpretation of SGV of trace gases and aerosols.
The following meteorological physics options were employed: the Rapid Radiative Transfer Model (RRTM) for longwave (Mlawer et al., 1997), the Goddard shortwave scheme (Chou and Suarez, 1994), the Noah land surface model (Chen and Dudhia, 2001) for land surface processes, the Kain-Fritsch cumulus and shallow convection scheme (Kain, 2004) (for domains with grid spacing greater than 10 km), the Yonsei University nonlocal boundary layer turbulence transfer scheme (Hong et al., 2006), and the Lin mixed phase cloud microphysics scheme.Advection included the www.atmos-chem-phys.net/10/6917/2010/Atmos.Chem.Phys., 10, 6917-6946, 2010 positive definite limiter (Skamarock, 2006) for both the water and chemistry species as was found necessary by Chapman et al. (2009) to prevent spurious mass production.

Experiment design
The period of the simulations, from 06:00 UTC (00:00 LT, Local Time) 1 March to 06:00 UTC 30 March, coincides with most of the airborne and surface measurements during MILAGRO.Only the results from 00:00 UTC 5 March to 00:00 UTC 30 March are averaged and used in the analysis shown later in this paper.Two computational domains are employed.The larger domain, which encompasses Mexico, Southern Texas, and a portion of Central America, has 75km grid spacing.A smaller domain, encompassing central Mexico, a portion of the Gulf of Mexico and includes a large fraction of the aircraft flight paths is used with smaller grid spacings.Simulations for the smaller domain use either 3-km or 15-km grid spacing.The analysis is conducted for a series of four locations that lie along the dominant synoptic flow pattern from Mexico City towards the Gulf of Mexico, with each station farther from the large source of emissions over Mexico City.The sites are referred to as T1, representing an area of large urban emissions, T2, which is a downwind with significantly lower local emissions than T1, T3, which is over the Mexican Plateau border, and T4, which represents a remote region influenced by marine processes.The 3-km and 15-km domains are setup with identical corner locations that coincide with the corners of cells on the 75-km grid.Each grid cell in the 75-km simulation consists of a 5 × 5 set of cells in the 15-km grid, and each grid cell in 15-km simulation consists of a 5×5 set of cells in 3-km simulation (see Fig. 1).This allows us to easily compare values over equivalent regions between grids.The statistics comparing the high-and low-resolution grids are all for a 75 km by 75 km square equivalent to the grid cell from the 75-km grid cell containing the site.We define this area as GC75.For the 75-, 15-, and 3-km grids, the resulting area contains 1, 25, and 625 cells, respectively, that can contribute to variability within the region.
To ensure identical boundary forcings to the central region of the 75-km grid where all three grids cover identical areas, the initial and boundary conditions are handled differently between the 75-km and the other two grids.For the 75-km grid, initial and boundary conditions are provided at 6 h intervals for the meteorological variables from the National Center for Environmental Prediction's Global Forecast System (GFS) model on a 1 by 1-degree grid; the initial and boundary conditions for the trace gas and aerosol species are provided every 6 h by the MOZART-4 global chemistry model (Pfister et al., 2008) as in Fast et al. (2009).Initial and boundary conditions for meteorological, trace gas, and aerosol variables for the 15-km and 3-km grids are derived once per hour from the 75-km grid using one-way nesting.This procedure ensures that the large-scale forcing for the region of comparison is identical, allowing us to attribute differences between the simulations to local impacts and SGV derived from differences between the grid resolutions.
Emissions for this study are identical to those used by Fast et al. (2009) with the exception that they have been regridded to the domains used in this study.Emissions of anthropogenic trace gases and particulates were obtained from two inventories: the 2002 Mexico City Metropolitan Area (MCMA) inventory, developed by the Comisión Ambiental Metropolitana (CAM, 2004, Lei et al., 2007), and the 1999 National Emissions Inventory (NEI), developed by Mexico's Secretariat of the Environment and National Resources, the US Environmental Protection Agency, and several other groups (http://www.epa.gov/ttn/chief/net/mexico. html).Emissions of CO, NO x , SO 2 , volatile organic compounds (VOCs), NH 3 , PM 2.5 , and PM 10 are available for point, area, and mobile sources.The right column of Fig. 2 illustrates the area emission of CO for the three grid spacings.Biomass burning emissions are included and are based on the MODIS thermal anomalies product (Wiedinmyer et al., 2006).Biogenic emissions were calculated using the MEGAN v2.04 model.Dust emissions were calculated interactively during the model simulation based on grid cell wind speed, moisture, and other relevant conditions using the dust module in WRF-Chem based on Shaw et al. (2008) (Zhao et al., 2010).Note that because the dust and biogenic emissions are determined interactively during the model simulation, the total mass of emitted dust differs between domains, whereas the mass of other emissions is consistent between the three grids for equivalent areas.
Table 2 summarizes the experiments done in this study.C75, C15, and C3 are control simulations, in which the model configuration is identical except for the grid spacings of 75 km, 15 km, and 3 km, respectively.E15 is the same as C15 except that the anthropogenic and biomass burning emissions are on the 75-km grid instead of the 15-km grid, i.e. with a uniform emission over each GC75.E15 is a

Observational data
The Mexico City metropolitan area (MCMA), with a population of ∼20 million, is the largest metropolitan area in North America and is located within a basin on the central Mexican plateau at an elevation of 2200 m above sea level.Mountain ranges that are ∼1000 m higher than the basin floor border the west, south, and east sides of the city, affect the local and regional circulations.Remarkably large pollutant emission and the surrounding mountains make the MCMA and the surrounding region an opportune place to study the SGV of trace gases and aerosols.Air pollution in MCMA has been studied for many years (Raga et al., 2001;Salcedo, 2006;Molina et al., 2007).MILAGRO, composed of several collaborative field experiments supported by various Mexican institutions and the US National Science Foundation (NSF) and Department of Energy (DOE), is the largest of a series of campaigns in and around the MCMA (Molina et al., 2008).The month of March was selected for the field campaign period because of the dry, mostly sunny conditions observed over Mexico at this time of the year.A comprehensive set of meteorological, trace gas, and aerosol measurements was obtained at the surface and aloft over a wide range of spatial scales.Extensive surface chemistry and meteorological profiling measurements were made at three "supersites" denoted by T0, T1, and T2, of which the latter 2 are shown in Fig. 2 (e.g.Doran et al., 2007;Shaw et al., 2007).A detailed list of instruments and instrument platforms is given in Molina et al. (2010) and research findings derived from the measurements have been in a special section of "Atmospheric Chemistry and Physics" at http://www.atmos-chem-phys.org/specialissue83.html.
The observations used in this study were collected by the many scientists who participated in MILAGRO, and were ported into the Aerosol Modeling Testbed (AMT) developed at the Pacific Northwest National Laboratory (http://www.pnl.gov/atmospheric/research/aci/amt/).The AMT (Fast et al., 2010) collects all the MILAGRO measurements into a central location and reformats the data into a single format, significantly reducing the time needed for analysis and graphing (Rishel et al., 2009).Table 3 highlights the primary data used in this study, which consists of surface observations at the T0 and T1 MILAGRO supersites and measurements from the DOE G-1 aircraft.

Evaluation of WRF-Chem simulations at various spatial resolutions against observations
We first compare the performance of the model at various spatial resolutions.Fast et al. (2007Fast et al. ( , 2009)), Tie et al. (2009), andZhang et al. (2009) have comprehensively evaluated WRF-Chem simulations at 3 and 6 km grid spacing against observations of meteorology, trace gases, and aerosols during MILAGRO, with the simulations by Fast most closely resembling those in this study.Here we focus on comparing the present simulations against observations, with the objective of investigating uncertainty that arises from comparing point measurements to model grid cell estimates at different grid spacing.

Meteorological fields
High-pressure systems, weak synoptic forcing in the subtropics and horizontal temperature gradients over the central Mexican plateau are favorable for the development of local and regional thermally-driven flows (Fast et al., 2007).Several studies have evaluated simulations of near-surface winds and PBL structure over the MCMA (e.g.Fast and Zhong, 1998;de Foy et al., 2006).It remains a challenging task to simulate the details of near-surface winds at specific locations and times over areas with complex terrain, whereas most mesoscale models can capture the primary thermallydriven circulations and their interaction (Zhang et al., 2009).
As summarized in Fast et al. (2007)   aerosols also modify boundary layer properties (Jauregui, 1997) and subsequently near-surface transport and mixing of pollutants.

Wind speed
The left column of Fig. 3 shows the near-surface wind fields (10-m height) in C75, C15, and C3, averaged from 5 to 30 March.Generally, the wind speed is less than 4 m s −1 over land.The differences for wind speed and direction among the three simulations are smaller over the ocean and coastal plains where the surface is flat.In contrast, C75 is incapable of capturing many local wind features associated with complex terrain.The average wind speed simulated in C75 is very weak over the MCMA as shown in Fig. 3a, which would lead to an underprediction of pollutant transport.C3 shows much larger spatial variability in wind speed and direction over central and southern Mexico since it better captures the local thermally-driven down-slope and up-slope flows due its higher spatial resolution.Small-scale heating and terrain geometry associated with the mountain ranges leads to localscale circulations.While different wind patterns associated with various synoptic conditions (e.g.cold surge events) occurred during March (Fast et al., 2007), southerly winds can be found in the vicinity of Mexico City when averaged from 5 to 29 March.Indeed, Fast et al. (2007) suggested that the Mexico City pollutant plume is transported northeastward 20-30% of the time during March.
The left column of Fig. 4 compares the observed and simulated wind speed at T1. Associated with the change of PBL structure, the near-surface wind speed over Mexico City exhibits a strong diurnal cycle, with minimum wind speed during early morning and a maximum during late afternoon.The observed wind is between 1-5 m s −1 most of time and the maximum wind speed is less than 10 m s −1 .The magnitude and variability of wind speed at T2 is similar to T1.While all three simulations capture the diurnal cycle of wind speed, C75 significantly underestimates the variability and diurnal range of wind speed.C75 underpredicts the maximum wind speed by 40-50% during late afternoon and overpredicts the minimum wind speed during morning.As shown in Fig. 4, C3 and C15 capture the peak wind speed most days.It should be noted that it remains a difficult task to quantitatively compare the simulated wind speed at a specific location and time with the observation, especially when wind speed is low (Zhang et al., 2009).Surface wind measurements in an area with a complex underlying surface are not likely to be representative of a larger area.

PBL height
The PBL height (PBLH) is often used to describe the depth of the vertical mixing that affects the dispersion of pollutants.Mixing heights during the MCMA-2003 campaign reached around 3000 m and vigorous vertical mixing implied pollutants were well mixed in the PBL during the daytime (de Foy et al., 2006).Our simulations show that PBLH in MC and downward is usually lower than 0.5 km during nighttime, and it starts growing after sunrise (i.e.07:00 LT), after which it grows rapidly to 1.0-2.0km between 11:00 to 13:00 LT and reaches a peak of 2.5 km around 15:00 LT.The simulated variation and magnitude of PBLH are very similar over T1 and T2, which is consistent with the measurements of Doran et al. (2007).It has been noted that the YSU PBL scheme used in WRF has a tendency to collapse the afternoon PBL too quickly (Fast et al., 2009).
The middle and right columns of Fig. 3 compare the simulated PBLH for the three spatial resolutions for night and daytime hours.We can see the PBLH is usually less than 400 m over night and larger than 600 m over daytime.Since PBLH is strongly influenced by the terrain, the detailed features of PBLH over mountainous areas are not captured in C75.Generally, PBLH is higher over the southwestern portion of the grids and gradually decreases to the northeastward as shown in all three simulations, with the lowest values over the Gulf of Mexico.C75 overpredicts the PBLH over the eastern coastal plains and fails to capture the minimum PBLH along 19 • N near MC.

Relative humidity
Relative humidity (RH) is an important meteorological variable because it directly affects uptake and evaporation of aerosol water, thus significantly affecting aerosol optical properties.The right column of Fig. 4 compares the observed and simulated RH at T1.The RH is usually lower than 60% before 15 March, and afternoon minimum RH often drops to below 10%.The daily maximum RH rises to above 80% around 15 March because of an El Norte event, which transports moisture to the plateau from the Gulf of Mexico.The averaged RH, especially maximum RH during later nighttime, is higher during the second half of the month than during the first half of month.
All three simulations reproduce the diurnal cycle of RH, with the maximum value associated with a lower temperature around sunrise and the minimum value associated with a higher temperature during afternoon.However, C75 overpredicts maximum RH by 20-30% and minimum RH by 5-15%, respectively over T1.C15 slightly overpredicts the maximum RH in the first half of the month and underpredicts the maximum RH in the second half of month, while it captures the minimum RH most days.C3 captures both maximum and minimum RH well for the majority of days.Generally RH is higher over the Gulf of Mexico and gradually decreases to southwestward as shown in all three simulations, with lowest RH over the central Mexican plateau.

CO
The chemical lifetime of carbon monoxide (CO) is about 2 months, thus over the few days relevant here it can be considered as a passive tracer that is emitted from the surface, mixed in the PBL, and transported by the prevailing winds (Tie et al., 2009).We first examine the predictions of CO to evaluate the impact of transport on the SGV of trace gases.Fast et al. (2009), Tie et al. (2009) and Zhang et al. (2009) all evaluated their CO simulations against the observations from the RAMA operational monitors in Mexico City.Here, we focus on the comparison of simulated CO at various spatial resolutions.Generally, the three simulations reproduce the diurnal cycle of CO reasonably well.The surface CO concentration reaches a peak during 07:00-09:00 LT, because of the morning rush-hour traffic combined with accumulation during nighttime from shallow PBL depth and lower wind speed.As suggested by Tie et al. (2007), the diurnal variation of surface CO concentrations is mainly controlled by the daily variability of PBL height and emission of CO.As morning progresses, the PBL height increases allowing rapid dilution of CO concentrations.Overall, the discrepancies for surface CO concentration between model and observation, for both the mean and percentiles, is smaller than 20% for C3 over Mexico City.The consistency of the observed and simulated CO suggests that the overall emission estimates of CO are reasonable over the city.
Bias of simulated CO outside the city is larger than in the city, but the model qualitatively captures the magnitude and temporal variation of CO.This is probably related to uncertainties in the emission inventories outside the city.Rapid changes in urban growth at the edge of the city and traffic along the highway just to the south of T1 during the morning rush hour period may not be represented well.Over the suburban (i.e.T2) site, CO concentration does not show the same diurnal variation as in the city and the simulated peak value in the morning is much lower than observed.
Predictions of CO further downwind are also evaluated using aircraft measurements.The 10th, 25th, 75th, and 90th percentiles of CO concentration show that C3 overpredicts the range of observed CO on some days and underpredicts the range on others.When averaged among all the aircraft flights, the percentiles are slightly larger than the measurements and the median value is somewhat lower than observed.The percentiles and mean value of CO concentration is 10-20% lower in C75 than measured when all aircrafts data are averaged.The results suggest that C3 adequately reproduces the overall transport and mixing of CO downward of Mexico City, although there are errors in space and time for the exact position and magnitude of plumes.The spatial distribution of bottom model level CO for the three simulations are shown in the left column of Fig. 5, in which we find that surface CO concentration simulated in C75 is a factor of 3-4 lower than in C15 and C3 over Mexico City, whereas the overall pattern of plume transport is similar among the simulations.Indeed, the percentiles and median value of CO in C75 near the surface is a factor of 4-5 lower than the measurements, while C3 captures the median and extreme values of CO well.Poor performance of C75 at the surface indicates that 75-km horizontal grid spacing is insufficient to represent   local emissions and terrain-induced flows along the mountain ridge, subsequently affecting the transport and mixing of plumes from nearby sources.

O 3
In contrast to CO, ozone (O 3 ) and nitrogen oxides (NO x ) are more reactive trace gases and their transport and mixing are influenced by more factors than CO. Figure 6 shows variability that is broadly consistent between the model and observations for the O 3 (left) and NO x (right) concentrations along the G-1 flight path on two example days, 20 March and 9 March, respectively.On 20 March the observations have four major peaks when the aircraft passed through plumes from the city during the 3-h flight, with the peak values twice as large as the average O 3 concentration.While both C15 and C3 capture these peaks, C15 underpredicts the peak values by 20-30% and C3 only slightly underpredicts the peak values.Because of its coarse spatial resolution, C75 fails to reproduce any peaks over the entire flight path.These results indicate that the simulation with coarse spatial resolution is not able to capture the spatial variability of O 3 .

NO x
Aircraft measurement for NO x on 9 March show multiple peaks with various magnitudes during the 4-h flight path, with the peak values 5-10 times larger than the average NO x concentration.C3 captures the variation and magnitudes well, including the maximum NO x mixing ratio around 17:30 UTC.C75 simulates higher NO x along polluted portions of the flight track but fails to capture any peaks.The performance of C15 is between C3 and C75. Figure 6 shows that C15 captures the multiple peaks but underpredicts the magnitude for each peak.

Black carbon
The pollutants over the region are mainly from man-made emissions in the vicinity of Mexico City and biomass burning.The right column of Fig. 5 shows the surface BC concentration for the three control simulations.Similar to CO, the spatial distribution of BC is similar between C3 and C15, with maximum mixing ratios over MC, Puebla, and Orizaba, but C3 provides more detail and spatial variability.The BC concentration of C75, however, is 50-100% less than in C3 over the central Mexico Plateau and eastern Mexico.The center of maximum BC mixing ratio is shifted to the northeast of Mexico City, with a maximum value 0.7-0.8µg kg −1 in C75.Over Mexico City (e.g.T0, T1), observed BC exhibits a strong diurnal variation (Yu et al., 2009), with the concentration increasing during the night and the largest values occurring in the morning hours around 08:00 LT.Because of the diurnal variation of PBL depth and wind, the pollutants are trapped in the city overnight and during early morning hours in a shallow surface layer before the rapid mixed layer growth commences in the morning.Comparing with observations over Mexico City (e.g.T0, Fig. 7), C3 captures the variation and maximum values of BC very well.C15 generally captures the observed variation of BC but significantly underestimates the magnitude during the morning (06:00-10:00 LT).C75 fails to capture any peaks during entire week.
Although the behavior of BC at T2 is less regular between days, the averaged diurnal variation of BC is very similar downwind of the city, with the largest concentration around 08:00 LT, since the PBL structure and evolution at T2 is similar with that at T1. Overall, the simulated BC concentration at T2 is 20-30% lower than at T1.Those features are consistent with the observations at T1 and T2 (Yu et al., 2009).

Organic matter
The spatial distribution (not shown) and temporal variation (Fig. 7, right panel) of OM are very similar with those for BC.This is not surprising since BC and OM (excluding SOA in this study) share many common emission sources and their transport and mixing are often linked together (Hodzic et al., 2009).When comparing the three simulations against observations, the performance for OM is similar to BC (see Fig. 7 left panel).Transport from Mexico City to T2 appears to account for a substantial fraction of the BC and OM at T2.Since SOA is not included in the WRF-Chem simulations, it is not surprising that simulated organic aerosol mass from all three simulations are lower than observed.However, on the positive side, primary organic aerosols in similar simulations have matched well to hydrocarbon-like organic aerosols derived from Aerosol Mass Spectrometer measurements (Fast et al., 2009;Hodzic et al., 2009).

Aerosol optical depth
Figure 8 (left) shows scatter plots of AOD at 500 nm for the observations at T1 versus the three simulations.The observed AOD ranges from 0.1 to 0.7, however, AOD simulated in C75 is below 0.3 for most cases, which is significantly underpredicted.C15 performs better than C75, with AOD ranging from 0 to 0.45.AOD simulated by C3 ranges from 0.05 to 0.95 and compares even better with observations.Recently, AERONET data have been widely used to evaluate or constrain AOD for global aerosol modeling and data assimilation (e.g.Dubovik et al., 2002).This study suggests that point measured AOD may reasonably represent the model grid cell average if grid spacing is 3 km or smaller.When the grid resolution becomes larger, e.g.larger than 75 km, the model grid cell average may not correctly estimate AOD measured over polluted areas.When comparing the simulations against observations over the downwind T2 site, however, performance is closer among the three simulations (not shown).This indicates that the point measured AOD over remote area could be used to reasonably represent the average of model grid cell with larger grid spacing.

Aerosol number concentration
The aerosol number concentration, N a , is critical information to investigate the indirect affect of aerosol on cloud and precipitation.The right column of Fig. 8 shows percentiles comparing observed (with a PCASP) and simulated N a for MOSAIC size bin 3 (0.15625 µm to 0.3125 µm dry particle diameter) along the G-1 flight paths.C75 predicted N a , including median, percentile values and ranges, are more than 2-3 times lower than observed for almost all the flights during the field campaign.C15 performs better at simulating N a , especially during the first half of the month.Comparing C75 and C15 with C3 shows that C3 reasonably captures the median and percentile values of Na for most flights.All three simulations underpredict N a , possibly for many reasons.For example, SOA is not included which might impact nucleation and there is uncertainty in the homogeneous nucleation parameterization.Also, uncertainties exist in the size distribution of emitted particles and some emission sources are  missing, e.g.small biomass burning events, are not included in the emission inventories used for this study.The underprediction of N a is much more serious near the surface than aloft as measured by aircraft, especially for the simulation with coarser resolution (not shown).

Sub-grid variability (SGV) of trace gases and aerosols
In this section we first present the probability density function (PDF) of trace gases and aerosols, and analyze the characteristics of their SGV, including the spatial distribution, diurnal and vertical variations, and dependences on the distance to emission sources and the spatial resolution.

Definition of SGV
Each grid cell in C75 covers an area of 75 × 75 km 2 (GC75) and contains a set of 5 × 5 C15 grid cells and 25 × 25 C3 grid cells (see Fig. 1).In this study, SGV is defined as a normalized standard deviation (SD) within a 75 km×75 km grid cell at a particular hour, which is then averaged across all the hours in the simulation.
Here N = 25 for C15 and N = 625 for C3.T represents the number of hours in the analysis.Except where noted, this the analysis period is from 5-30 March 2008, so T = 600.
x i,t refers to the value of a given species (e.g.trace gas or particulate) for a C15 or C3 grid cell at a particular time t.xt is an average for the high resolution cells residing within the given C75 cell, GC75 at time t.This average consists of 25 C15 or 625 C3 grid cells for each GC75: Atmos.Chem.Phys., 10, 6917-6946, 2010 www.atmos-chem-phys.net/10/6917/2010/Spatial variability is normalized by the average within GC75 so that the magnitude of SGV is not dependent on the species types or values.For a few variables, such as elevation, we also use time averaged standard deviation (SD) to describe their spatial variability.

PDF
Figure 9 shows the frequency of occurrence of C3 simulated trace gases and aerosols at the surface, distributed as a function of concentration, over T1 GC75, which covers an area of 625 C3 grid cells.Simulated CO concentration for C3 is most frequently on the lower end of the distribution with approximately 50% of the cells having values between 0.15 and 0.45 ppm over Mexico City.The frequency of occurrence of higher CO concentrations dramatically decreases, and maximum CO concentrations extends beyond 2.1 ppm only over regions close to emission sources.The maximum CO mixing ratio is 10 times larger than the minimum one, which is consistent with the observations and simulations in Tie et al. (2009).
O 3 is distributed more evenly over a smaller range of values, i.e. the probabilities are more similar across the O 3 concentration values and the maximum concentration is around twice as large as the minimum one.The mixing ratio of O 3 shows an approximately normal distribution.In contrast to CO, O 3 is a secondary product with more diverse sources and sinks and has a shorter lifetime.Similar to CO, the frequency of occurrences of BC dramatically decreases with the increase of concentration and the maximum BC concentration is 9 times larger than the minimum one.Contrasting with BC are the properties of SO 4 , NO 3 , and NH 4 (sulfate + nitrate + ammonium = SNN), of which SO 4 is the dominant contributor by mass.SNN exhibit similar SGV properties with PDFs that are more evenly and narrowly distributed, and the maximum mixing ratio is around twice as large as the minimum one, although the frequency has an overall decreasing trend with the increase of SNN mixing ratio.
In summary for the urban location, secondary trace gases and aerosols, such as O 3 and SNN, are more likely to have a normal distribution over a small range of values.For primary trace gases and aerosols, such as CO and BC, it is more likely to have a large range of values over polluted source regions, with the frequency decreasing with increasing concentration.
For periods of southwesterly winds, T1, T2, T3, to T4 show the PDFs of the Mexico City plume as it is transported downwind.For other periods, the PDFs at these sites characterize the variability of trace gases and aerosols associated with high (T1) to low (T4) local emissions.For all periods, the PDFs from T1, T2, T3, and T4 gradually become more evenly distributed for CO, BC, O 3 and SNN.Over remote areas in general, all trace gases and aerosols are more evenly distributed compared to polluted areas.
Figure 9 also shows the frequency of occurrence of wind speed, RH, PBLH, and elevation over the 625 C3 grid cells surrounding the T1 GC75.Wind speed typically falls between 2.5 and 5.0 m s −1 with a negative skew.RH has a positive skew with values typically between 40% and 50%.Maximum PBLH is two times higher than the minimum, which indicates a strong spatial variability of PBL structure over MC.The elevation varies from 2200 m to 3600 m, with around 30% of the area higher than 2500 m and a much larger amount of area with lower elevations.

Spatial pattern
Figure 10 shows the C3 spatial distribution of SGV for surface CO, BC, and SNN, and the SD for wind speed, PBLH, and terrain height.Maximum SGV for CO is centered over the MCMA, with a maximum value of 0.6-0.8.The spatial distribution of SGV for CO coincides with the maximum CO concentration (Fig. 5) and emission rates (Fig. 2), even though SGV is normalized by the mean CO mixing ratio.This makes sense given that the strongest gradients in concentration typically exist where pollutants are emitted.High values exist in a grid cell containing emissions while nonemitting neighboring upwind cells can have very low concentrations in the most extreme case.Farther downwind pollutants become more spatially mixed, and therefore more uniform with smaller SGV.SGV for CO is smaller than 0.2 in other regions.Overall, the SGV for BC is larger than for CO, with a maximum value greater than 1.0 over Puebla (southeast of Mexico City), attributed to the higher emission and complex terrain along the southern border of the Mexican Plateau.The SGV of BC is larger than 0.3 over other remote regions, because of more widely spread emission sources (e.g.biomass burning) and larger variability in surface concentrations of BC (Fig. 5).
SGV of SNN is larger than 0.2 over the majority of land areas, with a maximum value around 0.5 over Puebla.The SGV of SNN, with a range of 0.4 to 0.5, is much smaller than for BC and CO over the MCMA, which is consistent with the more evenly distributed PDF for SNN as shown in Fig. 9.This implies that the overall features of the PDF for SNN at T1 are representative over a large portion of central Mexico.Conversely, the PDFs of CO and BC at T1 are not representative over as large an area as seen in Fig. 10.
The SD for wind speed over land is around 1.0 to 1.6 m s −1 (Fig. 10), which is half of the mean wind speed at 10-m height.The SGV for wind speed is between 0.4 and 0.55 over the majority of the ocean (not shown).The SD for PBLH is between 140 and 280 m over land, with SGV around 0.4 to 0.65 (not shown).Maximum SD for elevation is along the southern, eastern and northeastern borders of the Mexican Plateau, where T3 is located.

Vertical profile and diurnal variation of SGV
Figures 11 and 12 show the vertical profiles of SGV over the four GC75 for 6 variables during daytime and nighttime, respectively.During daytime (e.g.15:00 LT as shown in Fig. 11), the SGV for CO over T1 is around 0.6 within the PBL, i.e. below 2.75 km, but substantially drops to around 0.1 above the PBL top.The near vertically constant SGV within the PBL indicates that simulated CO is very efficiently transported and mixed vertically by turbulence but not well mixed horizontally, even in upper layers of the PBL.Indeed, the horizontal wind speed decreases only slightly with height and vertical gradient of CO concentration is very small within the PBL (not shown).The nighttime vertical profile of CO, as shown at 03:00 LT in Fig. 12, is significantly different than during daytime.The SGV for CO over T1 reaches 0.8 at the surface, and quickly decreases with height.The shallower PBL and lower wind speed during nighttime does not facilitate the dispersion of pollutants vertically and horizontally, so the CO mixing ratio is much larger at the surface during nighttime than during daytime.The CO concentration dramatically decreases in the free troposphere over T1 because of efficient horizontal transport, resulting in a significant decrease of CO SGV with height.The vertical variations of O 3 SGV are similar to CO during both daytime and nighttime, except for a maximum SGV for O 3 observed in the lower stratosphere.
Although the concentration of BC averaged over T1 GC75 does not vary significantly vertically within the PBL during daytime, which is similar to CO, the vertical variation of SGV for BC is different than for CO.The vertical variability of SGV for BC is also smaller below 1.5 km, which is similar as for CO.However, SGV for BC dramatically increases with height above 1.5 km and reaches a maximum value around the PBL top (around 3 km MSL), which is not observed for CO.The vertical variability of SGV of CO is relatively straightforward within the PBL because CO concentration is mainly controlled by its emission rate and meteorological conditions.The concentration of BC is also affected by deposition and interactions with other types of particles.An example of this type of interaction is the mixing of two different types of particulate species that have different spatial patterns for emissions.If only one of these species is present, there would be a given SGV, but when both species Y. Qian et al.: Sub-grid variability of trace gases and aerosols for GCMs Daytime (15:00 LT) SGV Profiles are present the SGV is modified, and typically increased, if the second species has a different emission pattern.Because MOSAIC uses an internal mixture for representing the particles, the net effect is an increased SGV for particles compared to gases.
The vertical structure of SGV for SNN is very similar to BC during daytime over all four GC75, with a maximum SGV at the PBL top.During nighttime, however, the SGV is almost constant vertically with no large second peak as observed during daytime.The different vertical variations of SGV for aerosols between daytime and nighttime within the PBL, where most of particulates are suspended, are important for estimating the SGV of direct radiative forcing of aerosol by reflecting and/or absorbing solar radiation during only daytime.
SGV of RH also exhibits double peaks in its vertical profile during daytime over T1 and T2 GC75, one at the surface and one at the PBL top.The vertical variability of SGV for RH is small and steady during nighttime over all four GC75.After being normalized by mean wind speed, the SGV for wind speed gradually decreases with height during nighttime.The SGV for wind speed during daytime is smaller at the surface and varies slightly vertically within the PBL.
Understanding the SGV peak at the top of the PBL is related to the combination of strong vertical gradients of concentration and horizontal gradients of the PBL top.For variables with a strong vertical gradient at the top of the PBL, undulations of the PBL top lead to strong horizontal gradients within the GC75.This is particularly true for RH and the particulate species.Animations of all the RH profiles contributing to a given GC75 (not shown) reveal that the RH typically is highest at the top of the PBL and then strongly decays above the PBL.Because the PBL height is not constant within a given GC75, small changes in the PBL height cause a strong gradient in RH leading to large SGV.At night when the PBL decays, the RH typically peaks near the surface and the peak aloft goes away.
It should be noted that there is a strong relationship between the SGV vertical profile for RH and the particulates.Close examination of Figs.11 and 12 reveal that the existence of the SGV peak at the top of the PBL typically correspond between RH, BC, and SNN.All of these tend to form or not form peaks consistently at the same time of day and locations.For example, in the nighttime profile for T3 the RH maintains the SGV peak similar to during the daytime and the aerosol SGV also maintains the peak aloft at night.There are mechanisms that connect RH to aerosol processes such as particle coagulation and chemical reaction rates.So, the RH SGV likely directly impacts the particulate SGV.However, since the particulates also have a strong gradient at the top of the PBL, the undulating PBL most likely plays a larger role in establishing the SGV peak at the PBL top than the RH interactions.
Figure 13 shows the diurnal cycle of SGV at the surface.SGV of CO exhibits two peaks over T1.The first peak with maximum SGV of 1.2 occurs at 07:00 LT when the CO concentration peaks and the PBL starts growing.The second SGV peak of about 1.0 occurs around midnight.The diurnal variation pattern of SGV for CO shown in Fig. 12 is consistent with that for CO concentration (not shown), but it should be remembered that SGV is normalized by mean CO concentration.While CO concentration during nighttime is larger than during daytime, it exhibits a noticeable minimum around 03:00 LT.In theory, the CO concentration should continuously increase from 18:00 to 06:00 LT with the accumulation of pollutants within the shallow nighttime PBL.The minimum CO concentration results from the much smaller CO emission rates late at night in Mexico City.In fact, surface emissions between 00:00 and 06:00 LT for CO as well as for EC and SO 2 all are 5-10 times lower than during the daytime.O 3 , which has higher concentrations during daytime and lower concentrations during nighttime, exhibits an opposite diurnal cycle compared to CO.However, SGV for O 3 shows a similar diurnal cycle to CO, except for a smaller magnitude, with one peak around 07:00 LT and another around midnight over T1 GC75.SGV for BC is larger during nighttime than during daytime, with one SGV peak of 1.1 around 07:00 LT.The SGV for BC exhibits a different diurnal variation compared to CO, although the concentration of BC has a very similar diurnal cycle as CO, with a low values around 03 affected by lower BC emission between 00:00 and 06:00 LT.The maximum SGV for SNN occurs in early afternoon with a magnitude around 0.75, smaller than CO, O 3 , or BC.The SGV for RH is around 0.2 and a maximum value occurs later in the afternoon that is less than 0.4, which is smaller than SGV for trace gases and aerosols.The SGV for wind speed reaches a maximum (above 0.7 for T1) around early morning over all GC75, which partly contributes to the maximum SGV for trace gases and aerosols around 07:00 LT.

Dependences of SGV on the distance to the polluted sources
Figures 11-13 also can be interpreted as vertical profiles and diurnal variations of SGV for increasing distance from the large urban emission source as once progresses from T1 to T2, T3, and T4.The vertical and diurnal variations over T1 GC75 are described in Sect.4.4.Here the discussion focuses on the differences of SGV among the four areas.
The vertical variation of CO and O 3 over T2 is similar with that over T1 during both daytime and nighttime, but with a smaller magnitude of SGV within the PBL over T2, especially at the surface during nighttime.For example, within the PBL, the SGV for both CO and O 3 over T2 is nearly constant vertically during daytime but at night decreases quickly with distance above the surface, which is consistent with T1.But the SGV for CO and O 3 at T2 is 30-40% as large as at T1 during daytime, resulting from smaller emission variability over T2.The SGV for CO and O 3 over T3 and T4 is less than 0.1 and the vertical variability is also small during both daytime and nighttime.
The vertical profile of SGV for BC and SNN over T2 is similar as over T1 and a SGV peak also appears at the PBL top during daytime, except with a smaller magnitude over T2.A SGV peak also appears at the PBL top over T3 in both daytime and nighttime, but vertical variability of SGV is smaller over T4.Generally, the SGV for BC and SNN are larger than for CO and O 3 because of larger and more variable emissions of BC and SNN precursors over rural or remote areas.For RH and wind speed, SGV over T2 exhibits similar vertical pattern as over T1, but with a smaller magnitude at T2 during daytime.Overall SGV for RH and wind speed does not differ as much over the four GC75 as for trace gases and aerosols because the spatial variability of meteorological variables does not depend on local emission rates.
The vertical profiles show that major differences in SGV among the four GC75 occur within the PBL, especially near the surface.Figure 13 compares the daily mean SGV over the four GC75 at the lowest model level, which reflects the state of the entire PBL for CO and O 3 .Surface SGV of CO at T1 is about two times larger than at T2 and 10 times larger than at T3 and T4.Surface SGV of O 3 at T1 is about two and three times as large as at T2 and at T3 and T4, respectively.SGV of BC and SNN over T2, T3, and T4 are smaller during daytime because of more efficient ventilation, are larger during nighttime, and do not have the two diurnal peaks.This is different than at T1 where two peaks exist at midnight and 07:00 LT.The SGV for BC and SNN are larger at T1 than at the other three GC75, but the difference between urban and downwind remote regions is much smaller than for CO and O 3 .For BC and SNN, the difference of SGV is very small among T2, T3, and T4.
As summarized in Fig. 14, the SGV for trace gases and aerosols generally decreases with the distance away from the urban area with large emission sources, even with the normalization by the mean concentrations.As described in Sect.2.2, the model lateral boundaries provide time dependent inflow conditions for pollutants from other portions of the globe, which is important for species with longer chemical life times, e.g.CO and O 3 .Aerosols and short-lived trace gases over T3 and T4 are not likely to be affected by the boundary conditions or contribute to their SGV over central Mexico.
Our simulations show that the decreasing rate of SGV with the distance is more significant for trace gases than for aerosols.Among the trace gases, the decrease of SGV with distance for CO is more dramatic than for O 3 , while among the aerosols the decrease of SGV with distance for BC is more dramatic than for SNN.This implies that the SGV primary trace gases and aerosols with a longer lifetime decrease with distance more quickly than for secondary trace gases and aerosols.This is probably related to the emission source, whether it is more localized or more widespread spatially, and the interaction and deposition processes of aerosols species.In addition, the rate at which SGV decreases with the distance away from the polluted urban area is more significant during daytime than at nighttime.

Dependences of SGV on the spatial resolution
Figure 14 also compares SGV for trace gases and aerosols over the four GC75 between two simulations with 3-km and 15-km grid spacing (i.e.C3 and C15), respectively.The SGV for CO, BC and SNN are 60-100% larger for C3 than for C15 over the T1 urban site, even though the SGV of emissions is only 25-35% larger at C3 compared to C15 over the same grid cell (Fig. 15).Over the other GC75 (i.e.T2, T3, and T4), the SGV for trace gases and aerosols are 30-60% larger for C3 than for C15, except for CO, which has a much lower SGV over these rural or remote areas.For meteorological variables, the SGV is 20-30% larger in C3 than in C15 for RH and PBLH and 20-60% larger for wind speed over the four GC75.The increase of SGV from C15 to C3, ranging from 150% over T1 to 60% over T3 (Fig. 15), is much larger for cloud optical depth than for other meteorological variables.Overall, the SGV for C3 is larger than for C15 for all variables, which numerically implies that the model solution has not converged at the 3-km grid spacing.However, it is unreasonable to expect full convergence, i.e. the result does not change with a further increase in resolution, at any affordable model resolution.Until convergence is reached, it appears that the magnitude of SGV generally increases with the spatial resolution of the model.This is due to the increased detail of the emission sources and additional small-scale forcings from clouds, topography, etc. that get introduced as the resolution increases.Most importantly, the increase of SGV for trace gases and aerosols is stronger than for most meteorological variables.This sensitivity is most likely due to the increased SGV of emissions at higher resolutions.

Impacts of emission and terrain on the SGV of trace gases and aerosols
In Sect. 4 we discussed the spatial and temporal variations of SGV for trace gases, aerosols, and meteorological variables.However, what factors affect the subgrid processing and variability of trace gases and aerosols, and how significant are each of those factors?In this section we discuss the contributions of emissions and orography on the SGV based on the results of sensitivity experiments.

Emissions
It would be expected that the spatial variability of emissions has a great impact on the SGV of trace gases and aerosols over urban areas.Figure 16 shows the vertical profiles of 24-h mean SGV over T1 for three simulations with 15-km grid spacing.The settings for E15 are exactly same as for C15, except that the emission rates in E15 are averaged to match the emissions from the 75-km grid.In effect, this makes the emission values constant for each 5 × 5 set of grid cells in E15, corresponding to the single 75-km grid cell from C75.So, the difference between output from C15 and E15 reflects the contribution of emission spatial variability between global model and regional model grid spacings on the SGV of trace gases and aerosols.
As shown in Fig. 16, the vertical profile of SGV for C15 is similar to the profile for C3 over T1 when we combine the day and nighttime profiles in Figs.11 and 12, except for a smaller magnitude in C15 as discussed in Sect.4.6.With a uniform emission rate used in T1 GC75 (i.e.E15), the SGV decreases from 0.38 to 0.21 for CO and from 0.34 to 0.18 for BC at the surface during daytime (not shown).The differences of SGV (15-25%) between C15 and E15, however, are much smaller during nighttime over T1, partly because of the minimum emissions rate after midnight (see Sect. 4.4).As shown in Fig. 16, the daily averaged SGV difference between E15 and C15 is around 37% for CO and BC at the surface, and decreases with height in the PBL.
Differences in SGV between C15 and E15 for O 3 and SNN decrease with height within the PBL during daytime (not shown) and are smaller than for CO and BC.When uniform emissions are used in each GC75, the SGV drops by 25% for O 3 and by 15% for SNN during daytime.The changes in SGV are minor for O 3 and SNN during the night.The daily averaged SGV changes are 10-20% for O 3 and SNN at the surface (Fig. 16).The differences of SGV between C15 and E15 are near zero in the free troposphere for trace gases, aerosols, and meteorological variables.The SGV for Y. Qian et al.: Sub-grid variability of trace gases and aerosols for GCMs meteorological variables, including RH and wind, are not affected by the emission rates in the simulations since there is minimal feedback between aerosols and meteorology in the chosen model configuration, only the aerosol direct effect.
The impact of emissions on SGV of trace gases and aerosols is much less significant over rural or remote areas.As shown in Fig. 17, the SGV differences are almost indistinguishable over T3, except for near the surface where SGV actually increases in E15 for CO, BC, and SNN.However, both the emission amount and the SGV are small over T3.At T3, which has complex and varied terrain, the SGV for trace gases and aerosols is primarily determined by terrain rather than by emissions.
In summary, the spatial variability of emissions can account for up to 50% of the SGV during daytime for longlived trace gases and aerosols, such as CO and BC, over urban areas like MC.The impact of emissions on secondary trace gases and aerosols, such as O 3 and SNN, are less than for CO and BC.The impact of emissions on the SGV of trace gases and aerosols decays with altitude in the PBL and has almost no effect in the free troposphere.Also, emission variability affects the SGV more significantly during daytime rather than during nighttime.This last point is interesting in that it represents a balance between processes operating on differing spatial and time scales that impact the SGV.For a given region, there is a diurnal cycle and spatial distribution associated with the emissions.The emissions are typically lower during nighttime.But simultaneously, PBL mixing is also minimum at night.This leads to the spatial structure imposed by the small-scale emission sources being maintained, and thus the impact on SGV of the emissions.Alternatively during the day, the emissions are typically higher, which implies a greater contribution to the SGV from the emissions.However, mixing within the PBL is also higher during the day, which would work to smooth out the spatial gradients, and thus reduce the SGV.Because the simulated SGV is actually stronger during the day, this implies that the increased mixing within the PBL is insufficient to counteract the higher emission rates.Whether or not this is a universal finding, or is specific to the Mexico City area, is unknown.

Terrain
The PBL evolution and regional flows, which are strongly affected by topography, have a significant impact on pollutant dispersion.Importantly, this impact is nonlinear and not easily generalized.For example, just because there is a lot of variability in the terrain height within a given grid cell, one cannot know a priori what the bias will be on the flow and PBL structure.The direction of mountain ridges or valleys within the cell, in combination with how these connect to features in neighboring grid cells, and the current meteorological conditions, will alter the terrain induced SGV for a given cell.
Figure 16 compares the vertical profiles of SGV between the simulations with different treatments of terrain (i.e.C15 vs. T15).The settings for T15 are exactly same as in C15, except that the terrain in T15 is identical to the terrain in C75.The terrain is constant for each set of 5 × 5 cells in T15 corresponding to a single cell in C75.So, the difference between C15 and T15 reflects the contribution of the terrain variability on the SGV of trace gases and aerosols, at least for variability differences between GCM-like and regional model-like grid spacings.
With the uniform terrain used in T1 GC75 (i.e.T15), the SGV significantly decreases for trace gases and aerosols over both urban and rural areas.In contrast to the emissions effect, the impact of terrain on SGV of CO and BC is minor at the surface over T1, but gradually increases with altitude in the PBL and above.Figure 16 shows that the daily mean SGV decreases 10-30% for CO, O 3 , and BC in the upper PBL and lower free troposphere.The terrain-induced reduction of SGV is 25-35% for SNN in the PBL, which is more significant than for trace gases over T1, especially at the surface.
The impacts of terrain variability on SGV of trace gases and aerosols are more significant over areas with more variable terrain, i.e.T3, which is located near the Mexican plateau border.The SD of terrain for T3 is maximum among the four GC75 (Fig. 15).As shown in Fig. 17, the trace gas and aerosol SGV over T3 decreases 30-40% in T15 compared to C15.The differences of SGV between C15 and T15 are similar from the surface through the free troposphere for CO and O 3 .The SGV decreases more for BC and SNN than for CO and O 3 , especially near the PBL top where the reduction of SGV approaches 50% for BC and SNN.
The terrain, by modifying the meteorological fields, such as wind and PBL structure, affects the dispersion and trans-port of trace gases and aerosols and their SGV.Figure 17 shows that the coarse terrain reduces the SGV of meteorological variables, including RH and wind speed.It is interesting to note that the impact of terrain on the SGV of RH and wind speed is maximum at the surface, and decreases with altitude up to the top of atmosphere over T1.Over T3, however, the difference of SGV for RH and wind speed between C15 and T15 is very small at the surface.Note that the T15 simulation, in which terrain is flat within each GC75 region consisting of 5 × 5 15-km grid cells, is just a sensitivity experiment to test how much the SGV of trace gases and aerosols could be changed due to the terrain.Since the terrain affects the meteorological fields nonlinearly, and thus the impacts on trace gases and aerosols, the quantitative change of SGV obtained in T15 could change if a different method were used to coarsen the terrain, e.g.filtering out short wavelengths from the C15 terrain instead of using the flat, stair-step method of substituting terrain from C75.

Discussion
With the expected improvement of high performance computational resources, the use of high spatial resolution is perceived as a solution to partially address problems with climate simulations, including for the aerosols and their radiative forcing.How much one gains by going to very high spatial resolution modeling, however, will remain an important question.Leung and Qian (2003) suggest that increasing spatial resolution does not appear to lead to uniform improvements in precipitation and snowpack simulations over complex terrain.They found an overprediction of precipitation along windward slopes of the Cascades as model goes to higher resolution.Meanwhile their results show that errors in the snow simulation are not simply explained by elevation bias; there is a tendency for the model to grossly underpredict snow.In numerical weather forecasts for the same region, Colle et al. (1999) obtained similar findings with the resolution increasing.By analogy, these findings imply that better results cannot be guaranteed in air quality and aerosol modeling as just increasing the model spatial resolution.This is why effort needs to go into understanding the SGV.
There is also the issue of higher resolution models leading to more stringent comparisons for certain statistical comparisons.Just because the statistics look worse for higher resolution does not always mean that the model behaves worse.This is essentially the opposite of the problem raised in the preceding paragraph.For example, contingency table based metrics, such as equitable threat scores (Gandin and Murphy, 1992), that rely on "hits" and "misses" of forecasted values can have more misses for fine-scale features when simulated at high resolution.With coarser simulations, the results are smoothed out over larger areas leading to higher chances of a hit, whereas a finer grid might be closer to reality but have a feature located incorrectly.Discerning improved model behavior is therefore very difficult when considering both the issue of true error introduced by changing the grid from error based on the analysis methodology.
Another point, that was noted in the Introduction, the SGV of trace gases and aerosols results from both the traditional subgrid processes affecting meteorology (e.g.clouds) and specific chemistry processes (e.g.emissions).Many of these processes are correlated due to mutual interactions.For example, the land use variability could potentially affect the SGV of aerosols and their precursors by changing the biogenic emissions.The vegetation can emit climate-sensitive biogenic VOCs that are oxidized in the atmosphere to form organic aerosols, or SOA.Therefore, any SGV of vegetation is related to the SGV of emissions.
A number of factors associated with our modeling study should be taken into account when assessing the results presented in this paper.The primary assumption is that the full effect of clouds on SGV of trace gases and aerosols is not included.The feedbacks between aerosols and clouds, specifically activation/resuspension, wet scavenging, aqueous chemistry, and the two indirect effects of aerosols in WRF-Chem are turned off.However, the impact of clouds on wind, temperature, and humidity on the local meteorological environment do affect the aerosol life cycle.The feedbacks were turned off to minimize the changing behavior of clouds at different resolutions since clouds behave very differently at different grid spacings.For example, the 75-km grid requires a cumulus parameterization but the 3-km grid does not.Since this is one of the first studies of subgrid variability for aerosols, we felt it appropriate to isolate the contributors to variability, as much as possible, to those that we could control more easily for sensitivity studies, such as emissions and terrain spatial variability.Clouds were not prevalent over the Mexican plateau prior to the cold surge on 23 March, but they did occur more frequently over the ocean and along the coast.Figure 15 shows that SGV of cloud optical depth, which ranges from 2.5 to 5.3 at 3-km grid spacing, almost one order of magnitude higher than for other meteorological variables.This is partially due to the threshold-like nature of cloud development, and also the fact that clouds are not a continuous field like other variables such as wind speed.Qian et al. (2001) shows that aqueous chemistry and wet removal are the primary factors regulating aerosols such as sulfate in all-sky conditions over East Asia.It can be expected that the cloud variability will significantly affect the SGV of trace gases and aerosols, especially for soluble aerosols.Therefore, the conclusions obtained in this study only represent the conditions under cloud-free skies, and should be interpreted as conservative values of SGV.
The second factor also relates to clouds.All settings (including lateral boundary conditions for both chemistry and meteorology) and parameters in WRF-Chem have been kept identical in the three simulations (i.e.C3, C15 and C75) except for grid spacing and the convection cloud parameterization.Convective cloud parameterization is used in C15 and C75 but not in C3, since it is assumed that clouds could be resolved in C3.Although the cloud-aerosol interactions are turned off in this study, clouds still influence other meteorological variables including wind and PBLH, thus affecting trace gases and aerosols.This is an inherent problem for comparing simulations under cloud-resolving and noncloud-resolving spatial resolutions.
The third factor is the poor performance of the model in simulating dust.Dust in northern Mexico accounts for the majority of mass in the simulated internally mixed aerosols particles, especially for the large size bins.The uncertainty in dust simulation affects the estimated AOD as well as the radiative forcing and deposition of aerosols.While we do not have direct measurements of crustal materials, the total PM 2.5 mass of the model is overpredicted (not shown), which we attribute to excessive dust mass.

Summary and conclusion
One fundamental property and limitation of grid based models is their inability to resolve spatial gradients smaller than twice the grid cell size due to aliasing, and for some processes four or more times the grid spacing is required due to numerical diffusion (Pielke, 1991).Sub-grid variability (SGV), as illustrated in this study for meteorology, trace gases, and aerosols is an inherent problem of all grid models.For air quality or chemistry models, the SGV is affected not only by the traditional sub-grid processes affecting meteorology, such as the within-grid variability of terrain, land surface properties, and clouds, but also by specific chemical processes, such as emissions.While decades of work have gone into developing sub-grid treatments for clouds or land surface process in climate models, the quantitative understanding of sub-grid processes and variability for aerosols and their precursors is much poorer.In this study, WRF-Chem is applied to simulate the aerosols and other trace gases over the vicinity of MC during the 2006 MILAGRO field campaign using multiple spatial resolutions and scenarios that examine SGV of emissions and terrain.Our analysis focuses on quantifying the SGV of trace gases and aerosols within a typical GCM grid cell, i.e. 75 × 75 km 2 .
We first compared the model performance at three grid spacings (i.e. 3 km, 15 km, and 75 km) in simulating meteorological variables and trace gases and aerosols.C75 significantly underestimates the diurnal variability of wind speed, PBLH, and RH and did not capture many local features of these meteorological variables associated with complex terrain in central Mexico, while C3 captured their spatial and diurnal variability reasonably well.The surface CO concentration, which reached a peak during 07:00-09:00 LT, was associated with morning rush-hour traffic, low wind speed, and weak vertical mixing within the shallow PBL.Overall, the bias of C3 simulated surface CO concentration, including mean and percentiles, was smaller than 20% over the MCMA.The C3 simulation also captured multiple peaks of BC measured at the surface and of O 3 and NO x in the atmosphere observed by aircraft flights, which indicates that C3 adequately reproduced the overall transport and mixing of trace gases downwind of the MCMA.Conversely, C75 failed to capture any peaks for BC and OC at the surface or for trace gases aloft compared to aircraft observations.The poor performance of C75 at the surface indicates that 75km horizontal grid spacing was insufficient to represent local emission patterns and terrain-induced flows along the mountain ridge and subsequently affected the transport and mixing of plumes from nearby sources.All three simulations underpredicted the aerosol number concentration, partly because SOA is not included and some emission sources are missed in the WRF-Chem simulations.The underprediction of aerosol number concentration is much more serious near the surface than aloft, as measured by aircraft, especially for the simulation with coarse resolution.The evaluation of simulated AOD at various spatial resolutions against measurements suggests that point measured AOD may reasonably represent the model grid cell average if grid spacing is around or smaller than 3 km.When the grid spatial resolution becomes larger, e.g.larger than 75 km, the model grid cell average may not correctly represent AOD measured over areas with high emission rates.
PDFs for trace gases and aerosols show that more reactive and better mixed trace gases or aerosols, such as O 3 and SNN, were more likely to have evenly distributed and narrow distributions of concentrations (i.e.smaller SGV).Primary and less reactive trace gases and aerosols, such as CO and BC, were more likely to have broad distributions (i.e.larger SGV) over polluted regions.Over remote areas, all trace gases and aerosols were more evenly and narrowly distributed compared to polluted areas.Along the path from T1, T2, T3, to T4, the distributions became more evenly and narrowly distributed (i.e.smaller SGV) for CO, BC, O 3 , and SNN as the distance from high emission sources increases.Both CO and O 3 SGV vertical profiles were nearly constant within the PBL during daytime, indicating that trace gases were very efficiently transported and mixed vertically by turbulence but not well mixed horizontally in the PBL.During nighttime, the SGV for trace gases reached a maximum at the surface and quickly decreased with height.The shallower PBL and lower wind speed during nighttime did not facilitate the dispersion of pollutants vertically and horizontally.Unlike the trace gases, the SGV of BC and SNN reached a maximum at the PBL top during the day.The SGV of trace gases and aerosols generally decreased downwind of the MCMA, which is a large emission source, and the rate of SGV decrease is greater for trace gases than for aerosols.The SGV decreased with downwind distance away from the polluted urban area, had a quicker decrease for less-reactive trace gas or aerosol than for reactive ones, and was greater during daytime than nighttime.The SGV for C3 was larger than for C15 for trace gases and aerosols (by 60-100%) and for me-teorological variables (by 20-60%), which indicates that the magnitude of SGV generally increased with the spatial resolution of model.
Emissions can account for up to 50% of the SGV over urban areas such as the MCMA during daytime for lessreactive trace gases and aerosols, such as CO and BC.The impact of emission spatial variability on SGV decayed with altitude in the PBL and was insignificant in the free troposphere.The emission variability affects SGV more significantly during daytime (rather than nighttime) and over urban (rather than rural or remote areas).The terrain, through its impact on meteorological fields such as wind and the PBL structure, affected the transport and mixing of trace gases and aerosols and their SGV.The impacts of terrain spatial variability on SGV were more significant over areas with more variable terrain, i.e.T3.The SGV decreased more for BC and SNN than for CO and O 3 , especially within the PBL where the reduction of SGV can be up to 50% for BC and SNN at the upper PBL layers.
The above results that quantify the basic characteristics of SGV and their causes will improve our understanding of sub-grid processes that affect trace gases and aerosols.The results will also provide useful information for parameterization developers who need to take into account sub-grid scale variability for aerosols and their precursors so that they can reduce uncertainty in estimating aerosol radiative forcing on climate at larger scales.

Fig. 1 .
Fig. 1.Conceptual subdivision of a coarse grid cell from the 75-km grid into smaller grid points for the finer grids.Each 75-km grid cell is subdivided into 25 15-km grid cells.Each of the 15-km grid cells are further refined into 25 3-km grid cells for a total of 625 3-km grid cells per 75-km cell.
and de Foy et al. (2008), clear skies, low humidity, and weak winds aloft associated with high-pressure systems are usually observed over Mexico during March.The near surface winds over the central Mexican plateau are influenced by the thermally-driven circulation associated with terrain and large-scale synoptic flow.The thermal and dynamic effects of urbanization and

Y.Fig. 8 .
Fig.8.Scatter plots of observed versus the closest model grid cell simulated AOD at 500 nm at T1 (left column) for the 75-km (top), 15-km (middle), and 3-km (bottom) simulations, respectively.Comparison of observed (blue) and simulated (red) aerosol number concentration at for MOSAIC size bin 3 (particle dry diameters between 0.15625 and 0.3125 µm, unit: # cm −3 ) along the G-1 flight paths, where horizontal lines denote the median, boxes denote 25th and 75th percentiles, and vertical lines denote 10th and 90th percentiles (right column).

Fig. 9 .
Fig. 9. Probability density functions (PDF) of values from the 3-km simulation that lie within the 75-km host grid cell for the T1 site for the following variables: CO, O 3 , BC, SNN (SO 4 + NO 3 + NH 4 ), wind speed, relative humidity, PBL height, and terrain height.The PDFs encompass the time period 5-30 March 2006.

Fig. 10 .
Fig. 10.The spatial distribution of SGV for the 3-km simulation within each 75-km grid cell from the coarser domain for bottom model level CO (ppb), BC (µg kg −1 ), and SNN (µg kg −1 ) (left column) and of Standard Deviation (SD) for wind speed (m s −1 ), PBL height (m), and terrain elevation (m) (right column).The SGV here is averaged over the time period 5-30 March 2006.

Fig. 11 .
Fig. 11.Vertical profiles of SGV for CO, O 3 , BC, SNN, RH, and wind speed over T1, T2, T3, and T4 at 15:00 LT.The vertical axis is the model η level, which makes comparison between the points easier, since each point is at a different elevation.The SGV is calculated from the single hour each day for the time period 5-30 March 2006.

Fig. 15 .
Fig. 15.The SGV in C3 and C15 simulations for surface emission (CO, EC and SO 2 ), cloud optical depth, and PBL height, and SD for elevation.The SGV and SD are calculated over the time period 5-30 March 2006.

Fig. 16 .
Fig. 16.Comparison of vertical profiles of daily mean SGV simulated in C15, E15, and T15 simulations for CO, O 3 , BC, SNN, RH, and wind speed over the T1 site.The SGV is calculated over the time period 5-30 March 2006.

Table 1 .
Example influences on aerosol subgrid variability.

Table 3 .
Primary observations used in this study.