Estimation of the Impact of Ozone on Four Economically Important Crops in the City Belt of Central Mexico

In this work, we report the economic impact of exposure to high ozone concentrations on four important crops in the area of influence of the Mexico City Megalopolis. Estimated yield losses were as follows: maize: 3%; oats: 26%; beans: 14%; sorghum: 15%. The information needed to estimate the impact of air pollution in Mexico is decidedly deficient. Regarding ozone, the coverage provided by the monitoring networks is strongly focused on urban monitoring and its consistency over time is highly irregular. Apart from the Mexico City Metropolitan Area (MCMA) and less than a handful of other cities, the quality of the data is poor. Ozone in rural areas can be estimated with air quality models. However, these models depend on a high-resolution emissions inventory, which has only been done through validation processes in the MCMA. With these limitations, we set out to estimate the economic impact of exposure to ozone in these crops with a varying degree of sensitivity to ozone in the city belt of Central Mexico. To this end, we developed a procedure that makes optimal use of the sparse information available for construction of AOT40 (accumulated exposure over the threshold of 40 ppb) exceedance maps for the 2011 growing season. We believe that, due to the way in which we dealt with the sparse information and the uncertainty regarding the available data, our findings lie on the safe side of having little knowledge such that they may be useful to decision-makers. We believe that this procedure can be extended to the rest of the country, and that it may be useful to developing countries with similar monitoring and modeling capacities. In addition, these impacts are not evenly distributed in the region and sometimes they were greater in municipalities that have a higher index of poverty. Air pollution arriving from urban areas increases the social inequalities to which these already vulnerable populations are exposed.


Introduction
Ozone is a delocalized free radical with high oxidant potential.In the troposphere it is formed by the photolysis of NO 2 and by lightning [1].Subsequent titration with NO yields a steady state driven by UV irradiance.In the presence of VOCs, the peroxy radicals (RO 2 • ), formed during their photo degradation, oxidize NO, opening a reaction path to NO 2 without consuming O 3 , thus allowing its accumulation.In urban areas VOC and NO x emissions enhance these processes, both formation and titration of O 3 .There, O 3 production is generally VOC-limited.In rural and periurban areas with lower NO x levels, O 3 production is NO x -limited [2].Increased ozone production occurs during times of high temperature and solar radiation, such as the stagnant high pressure systems in the spring and summer [3].The different spatial distributions of NO x and VOC production as well as NO destruction of ozone often result in higher ozone concentrations downwind of urban areas, rather than in the urban areas themselves, where they hinder plant growth [4,5].Ozone can also be brought into a region by long range transport or by stratosphere-troposphere exchange [6].
Plants exposed to O 3 show a reduction in growth [7,8], which in commercial crops means lesser yield and revenue.In the year 2000, economic losses attributed to the worldwide adverse effect of ground-level ozone on maize, soybean, and wheat were estimated at $11-18 billion USD [9].
One of the exposure indexes most closely associated with the injury observed in vegetation is the AOT40 cumulative exposure (accumulated exposure over the threshold of 40 ppb) index.This is calculated as the sum of the differences between the hourly mean ozone concentration and 40 ppb when the concentration exceeds 40 ppb during daylight hours.This value accumulates during periods of time that are defined based on the types of vegetation: usually three months for crops [10].The exposure-response functions for the AOT40 index were derived from studies conducted under experimental conditions in which environmental conditions favoring the development of the plant and the maximum absorption of the pollutant were ensured [11].
There are also dose-response functions, which indicate the response of the plant to the dose of ozone absorbed using the calculation of ozone flow rates from the air to the plant [12].By analogy with the exposure-response functions, a flow rate index for ozone accumulated during the plant's growth period is used.This flow rate is obtained by adding the hourly ozone flow rates over a given growth period that is defined depending on the type of vegetation.The calculation of the ozone flow rate takes into account the ambient ozone concentration and stomatal conductance [13].Dose-exposure functions are easier to use and are more widely used [2,14].
Geographic information systems (GISs) can be used in order to estimate the economic impact of ozone on crops [15,16].O 3 maps consistent with the chosen function are made and overlaid on crop maps in order to apply those functions [15][16][17].However, worldwide, there are much fewer rural ozone monitoring networks than urban ozone monitoring networks.Ozone monitoring networks are usually established in cities in order to protect the health of the population.This deficiency is even more visible in developing countries.As of 2011, in Mexico, there was not a single continuous air quality monitoring station specifically devoted to rural monitoring, and in the city belt of Central Mexico, there were three air quality monitoring systems (AQMSs) reporting ozone records from 36 monitoring stations to the "Sistema Nacional de Indicadores de Calidad del Aire" (SINAICA) (National Air Quality Information System) [18], not all of which have the same quality standards as the Mexico City Metropolitan Area (MCMA) network [19].The data thus obtained are not sufficient to create exceedance maps for the critical value AOT40 in the region.It is well established that ozone concentration increases downwind of urban areas due to the absence of fresh of NO x emissions, from which NO can titrate O 3 [20].Creating maps by interpolating only data from urban stations implies underestimation of the concentrations over all the rural areas in the city belt.Correlations between urban and rural sites [21] cannot be obtained in the complete absence of rural sites.The orographic barriers that separate the air basins in the center of the country impose further restrictions on the interpolation methods, when only using data from the available monitoring stations.Figure 1 shows how the monitoring stations (empty white circles) cluster in three groups.The largest group at the center corresponds to the MCMA, the cluster to the west corresponds to the Toluca Metropolitan Area, and the group located towards the east corresponds to the Puebla Metropolitan Area.Passive samplers are regarded as a practical choice for measuring ambient ozone concentrations for the purposes of atmospheric chemistry and ecological assessment [22,23].They are inexpensive and easy to deploy in the field, requiring no power supply or supporting equipment, and they also offer flexibility of placement.These features make them attractive for monitoring in rural areas and in difficult environments, such as within the forest canopy [24].However, they only provide continuous cumulative values without differentiating day from night.To obtain AOT40 estimates based on measurements with passive samplers, it is necessary to correct the total quantity of O3 provided by the passive sampler to the equivalent ozone levels for daylight hours above 40 ppb.The use of mathematical functions from the Weibull probability family applied to nearby continuous ozone measurements has been the tool used in order to approximate the distribution of hourly O3 concentrations [25].However, the resulting model must be verified against monitoring stations in the study area, which in the case of Central Mexico are available only in urban areas.Nevertheless, some authors have proposed the use of the concept of relative altitude.This is an empirical statistical approach based on using the difference between the altitude at each point at which ozone levels have been measured and the lowest altitude in each circular area around each measured point [26,27], but this approach may not apply in areas that are under the influence of nearby sources of ozone precursors; it is subject to considerable differences between each particular site in the region [26].Therefore, exceedance maps for the city belt of Central Mexico cannot be obtained using the statistical modeling nor the relative altitude approaches.
Air quality models have also been used to estimate the formation and transport of ozone over the areas covered by the vegetation of interest [15,28], but the results of the model need to be validated against observations.Additionally, if available, data from monitoring networks can be used to adjust the results of air quality models [29][30][31].Tarrasón et al. [29] created difference fields in the simulation domain by interpolating the differences observed between the results of a regional air quality model and the measurements taken.This technique allows for correction of the results of the model at any point in the simulation domain.In our case, the absence of rural monitoring stations precluded us from creating such difference fields outside of urban areas.Passive samplers are regarded as a practical choice for measuring ambient ozone concentrations for the purposes of atmospheric chemistry and ecological assessment [22,23].They are inexpensive and easy to deploy in the field, requiring no power supply or supporting equipment, and they also offer flexibility of placement.These features make them attractive for monitoring in rural areas and in difficult environments, such as within the forest canopy [24].However, they only provide continuous cumulative values without differentiating day from night.To obtain AOT40 estimates based on measurements with passive samplers, it is necessary to correct the total quantity of O 3 provided by the passive sampler to the equivalent ozone levels for daylight hours above 40 ppb.The use of mathematical functions from the Weibull probability family applied to nearby continuous ozone measurements has been the tool used in order to approximate the distribution of hourly O 3 concentrations [25].However, the resulting model must be verified against monitoring stations in the study area, which in the case of Central Mexico are available only in urban areas.Nevertheless, some authors have proposed the use of the concept of relative altitude.This is an empirical statistical approach based on using the difference between the altitude at each point at which ozone levels have been measured and the lowest altitude in each circular area around each measured point [26,27], but this approach may not apply in areas that are under the influence of nearby sources of ozone precursors; it is subject to considerable differences between each particular site in the region [26].Therefore, exceedance maps for the city belt of Central Mexico cannot be obtained using the statistical modeling nor the relative altitude approaches.
Air quality models have also been used to estimate the formation and transport of ozone over the areas covered by the vegetation of interest [15,28], but the results of the model need to be validated against observations.Additionally, if available, data from monitoring networks can be used to adjust the results of air quality models [29][30][31].Tarrasón et al. [29] created difference fields in the simulation domain by interpolating the differences observed between the results of a regional air quality model and the measurements taken.This technique allows for correction of the results of the model at any point in the simulation domain.In our case, the absence of rural monitoring stations precluded us from creating such difference fields outside of urban areas.
Damage to vegetation by ozone has been reported for Mexico.These reports have been of a qualitative nature [32,33], or with too coarse resolution in global modeled studies [17].The same lack of information about ozone in rural areas that hindered quantitative estimates in those early studies still prevails.In this work, we report for the first time in Mexico, the estimated cost of agricultural production loss due to the exposure of several crops of economic importance to high concentrations of ozone present in the city belt of Central Mexico, which includes the Mexico City megalopolis.We address the methodological challenges facing an emerging economy such as that of Mexico in the estimation of the economic impacts of crop ozone exposure and we describe the equity issues that such exposure exacerbates.

Study Area
The study region, located in the central region of Mexico, is demarcated by the coordinates 20.6200-18.2276• N and 100.2030-97.5693• W and has an area of 72,900 km 2 .It encompasses the states of Morelos, Mexico State, Tlaxcala, and Mexico City and includes parts of the states of Guerrero, Hidalgo, Queretaro, and Puebla.The total population of this region is 33.77 million inhabitants, which represents 31.2% of the total national population [34].It has an estimated surface area of rain-fed crops of 28,154 km 2 [35].

Passive Monitoring
The passive sampler (Ogawa & Co., USA, Inc., Pompano Beach, FL, USA) consists of a cylindrical polymer body (2 × 3 cm) and a plastic pin clip holder.There are two cavities on the ends of the cylinder, each of which holds one coated filter between two stainless steel screens.The passive sampler uses two nitrite-coated filters in the interior part.When exposed to ambient O 3 , nitrite ion is oxidized to nitrate ion.Both exposed filters are extracted and analyzed by ion chromatography to determine cumulative nitrate concentration for the two-week exposure period.This allows for a determination of average cumulative O 3 exposure in ppb [37].
For analytical determination, a Shimadzu HPLC model CDD-10AVP was used.It is equipped with an isocratic pump LC-20AD and an ionic conductivity detector CD-10A-VP.The column was Shodex SI-90 4E (4.0 × 250 mm, particle 9 µm).The mobile phase was 1.8 mM Na 2 CO 3 + a 1.7 mM NaHCO 3 solution, with a flow of 2.0 mL/min and a 200 mL injection volume.Analysis time was 15 min.
Our implementation of this method was tested by co-colocation of passive samplers with a standard continuous O 3 monitor.In two independent experiments, the Ogawa passive samplers always overestimated accumulated O 3 concentrations by no more than 10% in samples taken between 1 day and 3 weeks.The remaining nitrite in the passive sampler followed an exponential decay curve as the sampling period increased, staying in its linear-like part for up to two weeks.
The passive sampling was conducted in a smaller area, the mountain zone south of the MCMA (Figure 1, triangles with dots).

Continuous Monitoring
The MCMA continuous ozone monitoring data were obtained from the databases of the environment authority of Mexico City [38].Data from the Metropolitan Area of Puebla and from the Metropolitan Area of Toluca were obtained from SINAICA [18].All of them are represented by empty white circles in Figure 1.Short-term (up to five weeks) continuous data were obtained in several rural and periurban sites within the region using mobile units as part of studies funded for other purposes [39][40][41] (Figure 1, stars).

Maps
Maps were created using ArcGis software v.10.2 by ESRI, Redlands, CA, USA and consisted of five layers in shapefile format, with Lambert conformal conic projection (LCC), Datum International Terrestrial Reference Frame 92 (ITRF92) epoch 1988.0 (DATUM ITRF92).These layers are as follows: state, municipal, land use administrative boundaries [42], the study area limit (own elaboration), and the final layer, in which the spatial distribution of the ozone exposure index studied in this work, AOT40 for the year 2011, is calculated and represented, with our data.
The maps were created on a 1 ha scale.This exceeds the spatial resolution of the information available to us on agricultural production.It also exceeds the spatial resolution in our Mesoscale Climate and Chemistry Model (MCCM) runs (3 km × 3 km).However, we kept this resolution because we expect that the government agencies that have information at the parcel level will eventually make it available.We are developing the capability of running WRF-Chem, instead of MCCM at 1 km 2 for Central Mexico.For AOT40, the selected interpolation method was kriging, which is widely applied for interpolating air quality data [31,43,44].

Exposure-Response Functions
The exposure-response functions (ERFs) from the literature were derived from open-top chamber studies in which the effect on the growth of plants exposed to different concentrations of ozone was measured.The functions compiled by Mills et al. [12] were used; these functions give a linear exposure-response relationship as a function of the AOT40 metric for the crops studied (Table 1).For sorghum, oats and beans, in the absence of specific functions we choose to use the generic sensitive function suggested by Mills.The functions were modified to force RY = 1.0 at x = 0 and RY = 0.95 at the critical level [16,17].
Table 1.Exposure-response functions for the four crops of economic interest analyzed in this work in the city belt of Central Mexico [12].AOT40: accumulated exposure over the threshold of 40 ppb.
Production loss per unit area was calculated using the following equation: where CPL i is the crop production loss [t] per unit area i [ha], RYL i is the relative yield loss, and CP i is the actual crop yield [t/ha] [16].
Municipal crop production loss used in this study is the mean of the crop production loss in all the crop cells within each municipality.
Finally, the economic damage was calculated using the following equation: where EL i is the value of the economic loss [$/ha], and ARP i was the average rural price per municipality [$/t] in 2011.In Mexico, this is the highest resolution crop price available.Monetary units are current Mexican pesos in 2011.Final aggregate values are given in USD using the average 2011 exchange rate.

Statistical Information
The first tests were carried out using information from the database of "Oficina Estatal de Información para el Desarrollo Rural Sustentable" (OEIDRUS) [45], which provided the spatial distribution of the plots of the crops of interest, with a resolution of 1 ha.Subsequently, we realized that the information available to us regarding other states was very irregular.In order to standardize procedures, we decided to use statistical information available from "Secretaría de Agricultura, Ganadería, Desarrollo Rural, Pesca y Alimentación" (SAGARPA) [46], which provides the following municipal scale information by crop: hectares planted, hectares harvested, average yield per hectare, production value, and average production value.
In order to include a measure of environmental equity associated with externalities arising from air pollution we used the marginalization index, which is defined as a measurement summary that allows for differentiation between locations in the country according to the overall impact of the deficiencies suffered by the population as a result of their lack of access to education and assets, and the fact that many reside in inadequate housing, all of which is associated with a low income [47].

Air Quality Model
The air quality model (AQM) used was the MCCM [48] in its non-hydrostatic mode.The chemical model is consistent with the meteorological model, using both the same transport schemes and the same horizontal and vertical grid.The chemical mechanism was RADM2 [49].The meteorological information was obtained from the global forecast (GFS).The modeling domain has 73 × 73 cells of 3 × 3 km and 22 vertical levels.The parametrizations for the cloud radiation scheme were as follows: non-hydrostatic option, Grell's cumulus parameterization, the Burk-Thompson PBL scheme, warm-rain scheme, Reisner scheme for explicit moisture, and a five-layer soil model.
A comprehensive assessment of model performance may be done using a set of metrics.Among them, the index of agreement I c [50] allows for a quick and intuitive verification of model performance.
where p i is the model-predicted concentration at time interval i. o i is the hourly average of the observed concentration at time interval i.I c varies from 0 to 1.When I c = 1, there is a perfect match between modeled and observed concentration values.

Creation of the Exceedance Maps
The limiting factor for estimating crop loss due to ozone exposure is the availability of cumulative AT40 exceedance maps.The method developed for arriving at a reliable estimate, on the safe side of having sparse information, was the first result of this work.Figure 2 shows a conceptual model of the procedure followed in order to obtain the exceedance maps with the information available in Mexico.The following is a description of the contents of each of the model's nine boxes.

Box 1
The first layer of information necessary to obtain the AOT40 exceedance maps was the results from MCCM for January 2009 to July 2011 for the regular grid shown in Figure 1 (white circles with dots).The high spatial and temporal resolution emissions model used in the simulations (INEM-AR 2008, Spanish initials), was obtained from the National Emissions Inventory disaggregated to a high spatial resolution (3 km × 3 km) using by default temporal emissions profiles or local profiles if available.As is commonly necessary in the MCMA, NO x and VOC emissions were adjusted in order to optimize the adjustment of the O 3 profiles modeled using the observations of the local monitoring network [51,52].In any case, the first task when using the model in each specific application is to evaluate its performance by comparing the results of the model against observations of the MCMA monitoring network during the episode or scenario under study [53] (Figure 1, white empty circles).The O 3 comparison between the model and data from several MCMA monitoring stations has an index of agreement greater than 0.81 [54].The index value decreased to 0.5 < I c < 0.66 when compared with the rural and periurban sites monitored during two intensive campaigns [40,41].The aim of some tasks in the following boxes is to improve on this metric.
Atmosphere 2018, 9, x FOR PEER REVIEW 7 of 18 network [51,52].In any case, the first task when using the model in each specific application is to evaluate its performance by comparing the results of the model against observations of the MCMA monitoring network during the episode or scenario under study [53] (Figure 1, white empty circles).The O3 comparison between the model and data from several MCMA monitoring stations has an index of agreement greater than 0.81 [54].The index value decreased to 0.5 < < 0.66 when compared with the rural and periurban sites monitored during two intensive campaigns [40,41].The aim of some tasks in the following boxes is to improve on this metric.


Box 3 To prepare the results from the model for comparison with the measurements taken by the passive samplers, which are of lower temporal resolution, the hourly ozone concentrations reported by MCCM were averaged over the two months of sampling for each season, resulting in the MCCM profile shown in Figure 3.The averages over each point of the grid, 12 sites, in the part of the To prepare the results from the model for comparison with the measurements taken by the passive samplers, which are of lower temporal resolution, the hourly ozone concentrations reported by MCCM were averaged over the two months of sampling for each season, resulting in the MCCM profile shown in Figure 3.The averages over each point of the grid, 12 sites, in the part of the simulation domain that contains the passive monitor sampling region were interpolated in the GIS in order to create a field of modeled averages and thus extract the averages over the passive monitoring sites for comparison.
Atmosphere 2018, 9, x FOR PEER REVIEW 8 of 18 in order to create a field of modeled averages and thus extract the averages over the passive monitoring sites for comparison.

 Box 4
The average value of the modeled field extracted at the passive sampling sites (Cmi) was compared with the average of the results observed with the passive samplers during the two months at the site ( ) as follows: where Cm is the mean of the modeled concentrations from MCCM for the period of time similar to the passive monitoring time.Cp is the concentration predicted by the simulation model throughout the simulation period.n is the number of time intervals for the entire simulation period.Simultaneously, where Co is the mean of the observed concentrations obtained from the passive monitoring operations for the sampling period.Cj is the concentration obtained by a passive monitoring operation for a period of 15 days.l is the number of samples collected during a campaign, be it during the cold dry, warm dry, or wet seasons.
For each passive monitoring site, At this point, a field of adjustment quotients such as that used by Tarrasón et al.
[29] can be created, but a review of the exceedance map shows that the region where an interpolation can be applied is very small.Therefore, at this stage, we opted for an average adjustment factor for the MCCM results.An average adjustment factor for the results of the model in the entire region was obtained as follows:

Box 4
The average value of the modeled field extracted at the passive sampling sites (Cm i ) was compared with the average of the results observed with the passive samplers during the two months at the site (C o i ) as follows: where C m is the mean of the modeled concentrations from MCCM for the period of time similar to the passive monitoring time.C p is the concentration predicted by the simulation model throughout the simulation period.n is the number of time intervals for the entire simulation period.Simultaneously, where C o is the mean of the observed concentrations obtained from the passive monitoring operations for the sampling period.C j is the concentration obtained by a passive monitoring operation for a period of 15 days.l is the number of samples collected during a campaign, be it during the cold dry, warm dry, or wet seasons.
For each passive monitoring site, At this point, a field of adjustment quotients such as that used by Tarrasón et al.
[29] can be created, but a review of the exceedance map shows that the region where an interpolation can be applied is very small.Therefore, at this stage, we opted for an average adjustment factor for the MCCM results.An average adjustment factor for results of the model in the entire region was obtained as follows: Table 2 shows the correction factors obtained for the passive sampling sites.The wet season data could not be used because, in that season, at all those sites, saturated relative humidity conditions are persistent; this causes condensation in the passive sampler.There is no specific pattern.Tres Marías is a small community in the State of Morelos that is under the influence of the Mexico City-Cuernavaca toll and non-toll highways.However, Parres, in Mexico City, is also under the influence of the same highways.Only Desierto de los Leones site showed the opposite ratio.In 2011 and 2012, our group carried out three intensive monitoring campaigns at several sites in two regions of the city belt of Central Mexico.The sites selected were rural or periurban.From the 2011 campaign [40], the Amecameca site was chosen.This site, located in the State of Mexico at the Tenango del Aire mountain pass, connects the Valley of Mexico with the valleys of Cuautla and Cuernavaca.From the 2012 campaign [41], the Chipilo site was chosen, and they were used as described below.
With the adjustment factor obtained in Box 4, the results of the model were corrected at all grid points and at the coordinates of the mobile units in Amecameca and Chipilo for the periods of time during which the campaigns were carried out; these corrected results were compared with the hourly profiles observed at those sites.At the urban sites, the comparison showed that the model results corrected with the average adjustment value obtained from passive monitoring considerably overestimated the reported values; however, for Amecameca and Chipilo, the comparison showed that the daytime profiles of the corrected model are more similar to the observations, with an overestimate between 10 and 30%.For completeness, a small additional correction was made.The nighttime profiles are also overestimated, a possible sign of the absence of enough nocturnal NO emissions in the emissions inventory.During nighttime, excess modeled boundary layer height and/or wind intensity could also reduce the NO concentrations available to titrate O 3 .These results are also valuable in top-down validations of emissions inventories, but that is the object of another work to be reported elsewhere.With two sites, a nighttime adjustment factor of 0.47 was obtained, from 19:00 to 07:00 h.Its impact on the AOT40 indicator is small, producing an effect only at dawn and at dusk.Nevertheless, it was applied so as not to overestimate the calculation of the AOT40 indicator for those two hours.At this stage, comparison of modeled and observed concentrations in the rural and periurban available pairs of data yielded much better values of the index of agreement 0.62 < I c < 0.90.

Box 7
With the results from MCCM optimized for the MCMA, the hourly ozone profiles were obtained for the 2011 crop season.These profiles, for all sites in the regular grid, were corrected with the daytime and nighttime adjustment factors, and the AOT40 indicator was obtained for all points the regular grid, using the kriging interpolating method (Map a, Figure 4).

Box 8
The values of the AOT40 indicator were obtained for all monitoring sites with data from the local MCMA network (RAMA) and that of the Metropolitan Area of Puebla (REMA) and added to the interpolating data set yielding Map b in Figure 4.These allowed for the correction, in urban areas, of the overestimation in the interpolation introduced by the indiscriminate use of the adjustment factor applied to the results of the model.In the absence of sufficient data from a station, the exceedances for the growing period were estimated using a ratio proportional to the number of observation days [26].

Box 9
The AOT40 indicator calculated for the rural and periurban sites in the 2011 [40] and 2012 [56] campaigns were then appended to the interpolating data set, yielding Figure 4c.
At this point, we still have some from a few rural sites not yet used and must decide it is safe to use them.Some were from the campaign carried out in 2009 in the State of Morelos (CARIEM) [39].We also had data from some sites, obtained by the mobile unit of the air quality monitoring system of the State of Mexico between 2006 and 2009, two sites, Tula [57] and ININ [58] plus data for periurban sites during MILAGRO [59].The AOT40 indicator was estimated for all of these sites.We decided it was safe to use them after the following empirical analysis.
In total, we had 16 data sets from short time campaigns.Two (AME and CHI) were used to confirm the use of the single average adjustment factor.Two from the 2011 campaign, OZU and TEN, and three from the 2012 campaign, AMO, CA, and HUA, had already been used to obtain the map in Figure 4c.To produce the final map (Figure 4d), we added the estimated AOT40 exceedance of the remaining nine rural and periurban sites to the interpolating data and analyzed the resulting map.Altogether, when making Figure 4c,f, of the 16 sites, 11 of them warmed the colors to the already adjusted map.Five of them modulated an excess correction by the average adjustment factor In the absence of a region wide of several years of AOT40 exceedance data, Figure 4d offers the best representation possible of the AOT40 exceedance using the available information for July-September months.We call it a hybrid map, as it uses results of a mesocale air quality model fed with the 2008 INEM-AR emissions inventory and driven by 2011 meteorology.These estimates were corrected with passive sampling data collected during 2011 and 2012 and with short-term continuous monitoring data in rural and periurban sites, 16 sites, in four field campaigns in 2006, 2009, 2011, and 2012.Figure 4d is consistent with some known facts that help to explain the shapes and shades on the maps.The parts of the map with the greatest accumulation of exceedances, the hottest parts, are the high-altitude O 3 in the higher parts of the Sierra Nevada to the east of the MCMA, the Sierra Chichinautzin to the south, and the Sierra de las Cruces to the west of the MCMA; on the other hand, the red patch to the north of the Toluca Valley Metropolitan Area (TVMA) represents a net receptor area of air masses, rich in ozone and still photochemically active, arriving, depending on the dominant direction of surface winds, from the TVMA [60] to the south-southwest, from the MCMA [61,62] to the east, or from the Tula-Tepeji industrial zone [63] to the north.The impact of the Tula-Tepeji industrial zone on the MCMA is well documented [54,63]; however, due to the presence of the Sierra de las Cruces, a difference of only a few degrees can divert emissions from Tula-Tepeji towards the TVMA instead of MCMA.
The maps in Figure 4 show that the parts with the lower AOT40 exceedances systematically correspond to areas where there is no measurement of ozone, despite the upward correction with a constant factor of the MCCM results.Additionally, in Figure 4a-d, in most cases, as the scarce observational data from rural and periurban sites not previously used were added to the interpolation database, the map showed warmer colors meaning higher AOT40 exceedances.This supports the following hypothesis: It is highly likely that the cumulative estimate of injury to vegetation and economic costs due to ozone exposure in the city belt of Central Mexico will be an underestimate of the actual costs and damages.

Estimation of Economic Losses
Figure 5a-f shows the results of the application of Equations ( 1)-(3) to the corresponding layers of information in the GIS.Table 3 shows input values and results from the application of Equations ( 1)-(3) in the 20 municipalities most affected by the exposure of maize to ozone throughout the simulation domain.
The impacts are not homogeneous in the region; they depend on the production of the crop, on the planted area, the cumulative AOT40, and the price paid to the farmer, which is not the same in all municipalities.For example, in the case of maize, Table 3 shows that the farmers in the municipality of Huamantla in the State of Tlaxcala lost 3% of their production due to ozone exposure.The estimated value of their production is 156 million pesos MXN and the economic loss is 5 million pesos MXN.In the municipality of Villa de Allende in the State of Mexico farmers lost 7% of their production to ozone exposure.The estimated value of their production is 152 million pesos MXN, and the economic loss is million pesos MXN.In terms of equity, it can be observed that the degree of urban marginalization, used here as a euphemism for poverty, in Huamantla is low, and in Villa de Allende it is high; surely such a potential loss of income represents a high cost to already impoverished communities.The total production value for maize in these 20 municipalities was 1.590 billion pesos and, due to ozone exposure, 85 million pesos MXN (5%) were lost.In the entire simulation domain, the value of maize production in 2011 was 6.598 billion pesos MXN; 46,641 tons were lost to yield reduction, and the value of that unfulfilled production was 191 million pesos MXN.This represents 3% of the potential production.
The same procedure was performed to estimate the effect of ozone exposure in oats, beans, and sorghum, which are more sensitive to ozone than maize.Table 4 shows the economic loss for each of the four crops throughout the simulation domain.The total economic loss was 671 million pesos MXN, which at the 2011 exchange rate equals 50 million USD.For the other crops, the relative impact of the degree of marginalization was not uniform in the municipalities that make up the simulation domain.It is also clear that these costs may be underestimated due to the likely underestimation of the AOT40 exceedances in the areas of the map with colder colors.

Conclusions
This work is the first attempt to estimate the cost of air pollution in Mexico, in terms of its impact on agriculture.Maize, the most important cereal in the Mexican diet, showed low sensitivity to ozone and economic losses of up to 3%.However, if the other crops are as sensitive as their proxies, then farmers who grow beans may lose up to 14% in potential revenue.Those farming sorghum may lose 15%, and for oats the loss may be up to 26% of the yield under optimal growth conditions.The impact on other species such as vegetables, which are highly sensitive to ozone and are cultivated near the cities, was not estimated, nor was the impact on species that are native to Mexico and economically important in the region, such as nopal and agave, due to the lack of exposure-response functions, even by proxy.
This attempt comes more than 30 years after continuous ozone monitoring was first performed in Mexico.It has been done with almost non-existent information for the affected rural areas.Although it is incomplete, the evidence accumulated by the incorporation of data from rural and periurban sites seems to show that, in terms of AOT40 exceedance maps, the costs reported here may be an underestimation.
The other source of uncertainty or bias is the use of proxy species which points to the need for national or local determination of exposure-response functions for crops of economic importance not reported previously in the literature.
It is evident that the agricultural areas affected by the metropolitan areas of the city belt of Central Mexico are strongly influenced by the cumulative AOT40 exceedances.Some of the affected municipalities are ranked high on the marginality index, tightening the screw on already-impoverished farming communities.
It falls upon the decision-makers and the affected groups to decide if the information provided here is sufficient to initiate actions for the mitigation of the impacts that negatively affect the availability of food and reduce the sustainability of the agricultural activities in the region.The first of such actions should be the long-term rural monitoring of air pollutants nationwide.

Figure 1 .
Figure 1.Study area with orography and boundaries between federal states and Mexico City.Four data point categories are shown.

Figure 1 .
Figure 1.Study area with orography and boundaries between federal states and Mexico City.Four data point categories are shown.

Figure 2 .
Figure 2. Conceptual model of the process used to obtain the AOT40 exceedance maps.

Figure 2 .
Figure 2. Conceptual model of the process used to obtain the AOT40 exceedance maps.

Figure 3 .
Figure 3.The average passive monitoring values at any site (Co) and the Mesoscale Climate and Chemistry Model (MCCM) average (Cm), which is obtained by averaging the profile with the hourly means from MCCM calculated over the two-month passive sampling period.

Figure 3 .
Figure 3.The average passive monitoring values at any site (Co) and the Mesoscale Climate and Chemistry Model (MCCM) average (Cm), which is obtained by averaging the profile with the hourly means from MCCM calculated over the two-month passive sampling period.

Figure 4 .
Figure 4. Evolution of the hybrid exceedance maps as observational data were gradually added to correct the maps obtained only with the air quality model.(a) Using only modeled and corrected regular grid points.(b) Panel (a) + RAMA + REMA.(c) Panel (b) + García-Yee + Barrera-Huertas.(d) Panel (c) + MILAGRO + CARIEM + Tula.

Table 2 .
Empirical correction factors for x = co /cm showing that MCCM systematically underestimates O 3 in the rural sites sampled using passive monitors, with X = 2.06 ± 1.0.

Table 3 .
Results for 2011 in the 20 municipalities most in terms of maize production, showing the input values and the results of the application of Equations (1)-(3) in the geographic information system (GIS).

Table 4 .
Planted area, production value, and economic loss estimated for four crops of economic importance in the city belt of Central Mexico.