*Article* **Land Use Changes in a Peri-Urban Area and Consequences on the Urban Heat Island**

**Marianna Nardino 1,\* and Nicola Laruccia <sup>2</sup>**


Received: 17 October 2019; Accepted: 12 November 2019; Published: 13 November 2019

**Abstract:** The effect of urbanization on microclimatic conditions is known as "urban heat islands". In comparison with surrounding rural areas, urban climate is characterized by higher mean temperature, especially during heat waves and during nights. This results in a higher energy requirement for air conditioning in buildings and in a greater bioclimatic discomfort for urban populations. The reasons of this phenomena are ascribable principally to the increase of solar radiation storage and to the decrease of dissipation of water by evapotranspiration in urban environment respect to rural ones. The aim of this paper is to give a quantification of the air temperature increase due to an urbanization process. This quantification is conducted by comparing surface energy balance (incoming and outcoming radiation and turbulent fluxes) in urbanized area versus rural areas. This quantitative approach will be validated using a fluidodynamic model (Envi-Met) in a case study area representative of one among the various regional models of urban area growth. In particular, the model of expansion of small towns around big cities (2003–2008 land use changes) of a plain near-urban area in the Po Valley region (Italy) was used.

**Keywords:** urban heat island; urbanization; urban surface energy balance; fluidodynamic modeling; Envi-Met

#### **1. Introduction**

The urbanization effect refers to a general increase in population and in the amount of industrialization of a settlement, due principally to the increase, in number and extent, of cities and the movement of people from rural to urban areas. The "urban sprawl" is used to define the increase in spatial scale or in the peripheral area of the cities. At the present, more than half of the global population lives in cities and cities themselves are growing to unprecedented sizes [1].

The high density of population and the consequent use of primary resources by urban residents, especially in the North Hemisphere, make cities and their inhabitants key drivers of global environmental changes [2,3].

Land-use and land-cover changes are recognized as causes of local, regional and global warming: The urban areas are the major sources of anthropogenic carbon dioxide emissions from the burning of fossil fuels for heating, from industrial processes, from transportation of people, etc. [4,5]. The main effects of land use changes (from rural to urban) can be found on the surface energy and water balances changes. The partitioning of sensible and latent heat fluxes is a function of varying soil water content and vegetation cover [6,7]. Water balance is strongly influenced by the soil sealing: Urbanization makes the surface permanently covered by impermeable artificial material (e.g. asphalt and concrete), for example through buildings and roads.

An increase of settlement areas over time is defined as land take, also referred to as land consumption. This process includes the development of scattered settlements in rural areas, the expansion of urban areas around an urban nucleus (including urban sprawl), and the conversion of land within an urban area (densification) [8,9]. Depending on local circumstances, a greater or smaller part of the land take will result in actual soil sealing. Urban sprawl can be defined as the unplanned incremental urban development, characterised by a low-density mix of land uses on the urban fringe. It is important to underline that even planned urban development may result in land take and soil sealing.

Soil sealing and realization of buildings on it alter the surface energy balance with two different main mechanisms:


The combination of these two different mechanisms results in a different thermal trend that is observed in the built environment compared to the surrounding rural areas.

One of the most prominent feature of the urban climate is the urban heat island (UHI) effect, which is strongly tied to the geometry and dimensions of building, land use patterns, vegetation cover and the intensity of the anthropogenic heat release [11,12]. The UHI makes the city warmer than its rural surroundings and it is stronger at night than during the day. This effect also decreases with increasing wind speed and cloud cover and it's less pronounced in summer and winter [11].

On average, urban temperatures may be 1–3 ◦C warmer than surrounding urban environment [6].

The aim of this work is to quantify the change in temperature range that can be expected on the basis of different energy balances resulting from the transformation of rural areas in urban residential or productive areas. The changes in the soil use (from natural surface to no permeable surface) returns an environment where the storage of heat during the day by building materials is released during the night, increasing the urban heat island. This effect is the main subject of the present work.

#### **2. Materials and Methods**

#### *2.1. Energy Balance Method*

The Earth surface radiation balance is described by the following equation:

$$\mathbf{R\_n} = \left(\mathbf{S}\mathbf{w\_{in}} - \mathbf{S}\mathbf{w\_{out}}\right) + \left(\mathbf{L}\mathbf{w\_{in}} - \mathbf{L}\mathbf{w\_{out}}\right) \tag{1}$$

where Rn is the net radiation, Swin is the incoming shortwave (visible) radiation, Swout is the outcoming shortwave radiation, Lwin is the incoming longwave (infrared) radiation and Lwout is the outcoming longwave radiation.

The term (Swin − Swout) is the net shortwave radiation and it can be also described as (1 − α)Swin, where α is the surface albedo used to quantify the solar radiation reflected by a surface. Albedo is a characteristic of the specific surface and depends on the optical characteristics of the reflecting surface.

Lwin and Lwout depend on the atmosphere and surface temperature following the black body equation:

$$\mathbf{L}\mathbf{w} = \boldsymbol{\varepsilon} \,\,\mathbf{\sigma} \,\,\mathbf{T}^4 \tag{2}$$

where ε is the body infrared emissivity, σ is the Stefan–Boltzmann constant (5.67 10−<sup>8</sup> W m−<sup>2</sup> K<sup>−</sup>4) and T is the body temperature expressed in K.

The net radiation given by (1) is then utilized to the partitioning processes at the surface, following the surface energy balance equation [12]:

$$\mathbf{R\_n} = \mathbf{H} + \mathbf{L}\mathbf{E} + \mathbf{G} \tag{3}$$

where H is the sensible heat flux, LE is the latent heat flux and G is the ground heat flux.

Usually the greatest part of the net radiation is used for sensible and evapotranspiration processes while the G term represent only the 10% of the total available energy (Rn).

To estimate the lower dissipation of latent heat flux due to the transformation of cultivated land in urbanized area, it is necessary to quantify the evapotranspiration (ET) before and after the transformation. The ET of a cultivated area depends on climatic characteristics of the site and on eventual water contributions from irrigation, on vegetable cover and on soil type. As vegetable cover, for this study, it was assumed the cultivation of winter wheat, not irrigated and growing in optimal pedological and climatic conditions to ensure a good water supply during vegetation season (November to June). Wheat is one of the more widespread dry crops in the Emilia-Romagna Region plain.

From a sealed surface, the evaporation of rain water can be considered negligible. On the other hand, from a winter wheat crop, in the case of optimal soil and climatic conditions, about 700 lm–2 of water per year changes from liquid into vapor phase.

A very simplified, but amply adopted method to calculate potential evapotranspiration (ETc) of a crop consists on multiplying reference evapotranspiration (ET0) per crop cultural coefficient (Kc):

$$\text{ETc} = \text{Kc } \text{ET}\_0 \tag{4}$$

In this study ET0 was calculated according to Hargreaves method [13] using the data of a meteorological station relatively close to the case study area, San Pietro Capofiume (44.648993 ◦N, 11.650055 ◦E, 11 m a.s.l.). The values of the cultural coefficient Kc for wheat were adopted following the [14] suggestions and reasonable Kc's for spontaneous grasses, growing after crop harvesting, were chosen.

If a land use change occurs, such as a significant urbanization, the terms of this budget equation change significantly and other terms, depending on different processes, must be taken into account. The Equation (3) of the surface energy balance can therefore be rewritten in this way:

$$\mathbf{R\_n} = \mathbf{H} + \mathbf{L}\mathbf{E} + \mathbf{G} + \mathbf{Q}\mathbf{s} \tag{5}$$

where Qs is the storage heat flux, which is strongly dependent on the ratio between green and sealed areas and on the building geometrical, optical and thermal characteristics.

Oke [6] proposed these equations for Qs:

$$\text{t day Qs} = \left(0.20\lambda\text{v} + 0.33\lambda\text{p}\right)\text{R}\_{\text{n}} + 3\lambda\text{v} + 24\lambda\text{p} \tag{6}$$

$$\text{migh}\,\text{Qs} = (0.54\lambda\text{v} + 0.90\lambda\text{p})\,\text{R}\_\text{n} \tag{7}$$

where λv is the green area fraction and λp is the building area fraction.

The reconstruction of the thermal change given by the lower amount of evapotranspiration of a sealed surface was performed following the Fick law [12]:

$$
\Delta \mathbf{T} = \mathbf{Q} \,\mathrm{D} \mathbf{z} \,\mathrm{k}^{-1} \tag{8}
$$

where

Q = Latent Heat Flux + Storage Heat Flux (W m<sup>−</sup>2)

Dz = reference height (2 m)

k = air thermal conductivity (0.026 W m−<sup>1</sup> K−<sup>1</sup> )

#### *2.2. Fluidodynamic Simulation (Envi-Met): The Case Study*

ENVI-met [15] is a three-dimensional non-hydrostatic microclimate model designed to simulate the surface–plant–air interactions within daily cycle in urban environment with a typical resolution of 0.5 to 10 m in space and 10 s in time.

The model has been widely used in many previous studies to simulated flow around and between buildings, exchange processes of heat and vapor at the ground surface and at the walls, turbulence exchange of vegetation and vegetation parameters, bioclimatology, and particle dispersion [16]. The model was validated as reported in [17].

In order to run the model, the detailed data on soil characteristics, buildings, vegetation, and initial atmospheric conditions for the area of interest were inserted.

After this input phase, the model is ready to run and the desired variables have been selected and saved into the output files [16].

ENVI-met can be used for several studies to test various urban canyon aspects as well as ratios and orientation effects on outdoor thermal comfort, the role of vegetation in the mitigation of the urban heat island effect, and other factors.

#### **3. Results and Discussion**

Hourly measurements of the surface radiation balance components (Swin, Swout, Lwin, Lwout) are available in stations relatively close to the case study area. San Pietro Capofiume was used as the representative station for rural open space. In this site CNR-IBE performed radiation components measurements for the period 1 January 2002–31 December 2003.

Bologna urban ARPAe (HydroMeteorological Service of the Emilia-Romagna Regional Agency for Environmental Protection) station (44.500754 ◦N, 11.328789 ◦E, 78 m a.s.l.) has been considered as the representative station for sealed conditions (measurements period 1 January 2006–31 December 2009). Although this data represents only few years, furthermore non-coincident, from Table 1 it can be inferred that the values of the two stations are substantially super imposable with regard to the Swin and the Swout values, while the differences between the values of Lwin and Lwout are attributable to the different surface types (rural and urban), principal aim of this paper.



Throughout the hourly measurements, in Table 2 are reported the diurnal and nocturnal Rn values obtained for the two stations and for each month of the year.


**Table 2.** Annual mean of net radiation values in nocturnal and diurnal hours for rural (San Pietro Capofiume) and urban (Bologna urban) sites.

Table 3 shows the values of mean monthly latent heat flux which would be missed in case of transformation of a wheat field into a sealed surface.

**Table 3.** Monthly mean values of reference evapotranspiration (ET0), crop coefficient (Kc), crop potential evapotranspiration (ETc) and corresponding latent heat flux.


The area chosen as a case study is located in San Giovanni Persiceto (44.6283 ◦N, 11.1992 ◦E, 21 m a.s.l.), a small town close to Bologna city. The urbanization occurred in this area from 2003 to 2008 was estimated comparing aerial images (orthophoto by AGEA for 2008 and QuickBird satellite orthoimages by Digital Globe™–Telespazio for 2003) (Figure 1). The built-up surface varied from 5.5% in the 2003 to 30% in 2008 and the natural cover varied from 94.5% in 2003 to 70% in 2008. Inserting these values as λv and λp in (6) and (7), the difference of the storage heat flux Qs between 2003 and 2008 was obtained, separately for night and daytime hours. The annual difference in the cumulated storage heat flux is quite similar between night and day (+80 W m–2 for the diurnal hours and <sup>−</sup>84 W m–2 for nocturnal ones) which means that all the heat stored during the day by urban materials is re-emitted into the atmosphere during the night. On a monthly basis, this effect is visible in the winter months and during months with high energy availability (May, June, and July) the day energy accumulation is greater than the night emission (Table 4).

**Figure 1.** QuickBird satellite orthoimages (®Digital Globe™–Telespazio, 2003) (left) and aerial orthophotographs (AGEA) for 2008 (right).

**Table 4.** Monthly mean values of storage heat flux during diurnal and nocturnal hours for 2003 and 2008 and relative differences for each year.


Since the evapotranspiration processes occur only during diurnal hours, the difference between the two reference years in terms of latent heat flux is assigned entirely to daylight hours. The greater energy availability for 2008 for heat transfer processes was obtained summing the diurnal latent heat flux differences (proportioned according to the fractions of cultivated and urban areas in 2003 and 2008) with the differences in diurnal storage heat flux. Nocturnal thermal variation depends only on variations in storage heat flux reported in Table 4.

Table 5 illustrates the total monthly differences (2008 vs 2003) in energy budget for diurnal hours and the corresponding diurnal, nocturnal and mean temperature variations calculated according to (8).

The largest increases in temperature due to an urbanization process occur during the summer months, when the energy available for exchange processes is greater.

This approach is an estimate of changes in energy balance and this leads at uncertainties, but in general the results show an annual mean increment in air temperature of 0.36 ◦C due to this urbanization process. Again, the single months show the differences due to energy availability, giving greater values of air temperature increment during summer months.



The same area was studied with the Envi-Met fluidodynamic model in order to have a modeling feedback of estimation method based on the surface energy balance.

The two input areas (2003 and 2008), shown in Figure 2, were inserted through the program ENVI-met Eddie, taking into account the geographical location, building dimensions, land use patterns. The vegetation was assumed to be winter wheat. The domain model consisted of a 140 × 140 × 20 grid with a spatial resolution of 5 m <sup>×</sup> 5 m <sup>×</sup> 2 m, resulting in a horizontal area of 700 <sup>×</sup> 700 m<sup>2</sup> with a 40 m vertical extent.

**Figure 2.** Input maps inserted into Envi-met model for 2003 (**a**) and 2008 (**b**).

The case study reported here focused on the simulation of a typical hot day whose values of meteorological variables were obtained considering the data of the nearest ARPAe weather station (San Pietro Capofiume).

The simulation started at 9:00 a.m. and lasted for 24 hours. The configuration file contained the atmospheric values:


As first results the potential temperature (the temperature that a sample of air attains if reduced to a pressure of 1000 millibars without receiving or losing heat) to the environment difference between the two simulation years (2003 and 2008) were plotted (Figure 3) at 14:00 p.m. The urbanization process carried out from 2003 to 2008 led to a warming of some areas especially those in which it has had a greater overbuilding. In these areas the 2008 potential temperature at 1.6 m height is about 0.4–0.6 ◦C higher than 2003.

**Figure 3.** Potential temperature difference (2008–2003) at 14:00 pm at 1.60 m height for the considered case study.

During the night (Figure 4 at 2:00 a.m.) the potential temperature difference between the two years reaches 1 ◦C. This is typical of the urban heat islands, where the greatest effect is recorded during the night, so substantial increase occurs in daily minimum temperatures [18]. In addition, the total area tends to be warmer in 2008 than in 2003, and this effect is more evident during the night hours than daytime.

As a matter of fact, if the air temperature trends is plotted for a strategic point of the considered area (a point which was completely rural in 2003 and that became urban in 2008, marked as P in Figure 2) it can be seen that during the day air temperature tends to be the same for the two years and during the night the change turns out to be even 1 ◦C (Figure 5). The daily mean of air temperature differences is 0.35 ◦C that, looking at Table 5 for June month, is comparable with the night values, but not with the daytime ones. Probably the high available energy during this month means that further exchange processes come into play, that the balance method, used in this work, does not

account for. A better diurnal values estimate is surely necessary, even if during August and September (months similar to June as far as concerns air temperature values) the obtained delta temperature (0.34 ◦C and 0.36 ◦C) is very close to the model result.

**Figure 4.** Potential temperature difference (2008–2003) at 2:00 a.m. at 1.60 m height for the considered case study.

**Figure 5.** Air temperature for 2003 and 2008 in the P point marked in Figure 2.

June is the month in which the budget method assigned greater evapotranspiration values, but the fluidodynamic model does not consider specifically crop type and evapotranspiration needs. For this reason, the strongest disagreement results especially in this month.

On the other hand, the model results are in strong agreement with the annual average increase obtained with the method of the energy budget (0.36 ◦C).

This suggests a profitable use of this methodology to estimate the increase in air temperature due to urbanization processes on a regional scale.

#### **4. Conclusions**

The effect of an urbanization in a small town shows the consequences of the well-defined and studied "urban heat island": The impact of such kind of urbanization leads an increment in the air temperature, and therefore a strengthening of the urban heat island, close to 1 ◦C.

The energy balance method is verified and supported by a fluidodynamic model to have the right data interpretation. The idea is to develop a simple methodology to consider the urbanization effects in terms of temperature increase, economic cost that goes with it and biometeorological uncomfortable for the people.

This methodology could be applied at regional scale to improve and develop spatial planning but taking into considerations the mitigating effects strategies (i.e. urban green, special building materials, study of the orientation of buildings, and shadows).

The radiation and energy budget method utilized to calculate the air temperature differences due to an urbanization process showed some limits and some uncertainties.

Surely, the different approaches to compute evapotranspiration used in the budget method and in the fluidodynamic model and some hypothesis adopted (absence of water deficit in budget method), strongly influenced the results.

Anyway, the fluidodynamic model confirms an air temperature increment of the same magnitude order of the energy budget method, suggesting this last methodology can be considered a sufficiently reliable method to estimate and predict the variation in the thermal field due to changes in land-use. Its application to a regional scale should take into account the different climatic zones and the different models of urban growth (existing residential or production areas densification, expansion of urban area with new buildings, etc.). The storage heat flux computation can be improved considering further variables: volumetric ratio between sealed and unsealed surfaces and optical and thermal properties of building materials.

Land use change, and in particular urbanization process, results in air temperature increase. The methodology developed in this work can be improved to furnish grater information during a regional territory planning. Air temperature increment give a surplus of energy needs of air conditioning systems: throughout thermotechnical computations it is possible to obtain the economic cost to reinstate indoor air temperature after urbanization process.

Moreover, the bioclimatic discomfort for the population caused by urban land use increment can be computed and successively taken into account during the initial design of an urban area.

**Author Contributions:** The authors contribute at all in this study thanks to different thanks to the different skills: fluid dynamics and agronomics. The manuscript was written together.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors would like to thank IdroMeteorological Service (ARPAe) of Emilia-Romagna Region for furnishing meteorological data through DEXTER service, AGEA for orthophotographs and Digital Globe™–Telespazio for QuickBird orthoimages, acquired through internal Map Service of Emilia-Romagna Region.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Air Pollution Flow Patterns in the Mexico City Region**

#### **Alejandro Salcido \*, Susana Carreón-Sierra and Ana-Teresa Celada-Murillo**

Instituto Nacional de Electricidad y Energías Limpias, Reforma 113, Palmira, Cuernavaca 62490, Morelos, Mexico; susana.carreon@ineel.mx (S.C.-S.); atcelada@ineel.mx (A.-T.C.-M.)

**\*** Correspondence: salcido@ineel.mx; Tel.: +52-777-362-3811 (ext. 7087)

Received: 18 September 2019; Accepted: 31 October 2019; Published: 5 November 2019

**Abstract:** According to the Mexico City Emissions Inventory, mobile sources are responsible for approximately 86% of nitrogen oxide emissions in this region, and correspond to a NOx emission of 51 and 58 kilotons per year in Mexico City and the State of Mexico, respectively. Ozone levels in this region are often high and persist as one of the main problems of air pollution. Identifying the main scenarios for the transport and dispersion of air pollutants requires the knowledge of their flow patterns. This work examines the surface flow patterns of air pollutants (NO2, O3, SO2, and PM10) in the area of Mexico City (a region with strong orographic influences) over the period 2001–2010. The flow condition of a pollutant depends on the spatial distribution of its concentration and the mode of wind circulation in the region. We achieved the identification and characterization of the pollutant flow patterns through the exploitation of the 1-hour average values of the pollutant concentrations and wind data provided by the atmospheric monitoring network of Mexico City and the application of the k-means method of cluster analysis. The data objects for the cluster analysis were obtained by modeling Mexico City as a 4-cell spatial domain and describing, for each pollutant, the flow state in a cell by the spatial averages of the horizontal pollutant flow vector and its gradients (the divergence and curl of the flow vector). We identified seven patterns for wind circulation and nine patterns for each of NO2, O3, PM10, and SO2 pollutant flows. Their seasonal and annual average intensities and probabilities of occurrence were estimated.

**Keywords:** pollution flow patterns; wind circulation patterns; emission inventory; criteria pollutants; Mexico City

#### **1. Introduction**

In cities and large urban settlements, tropospheric ozone and particulate matter (PM10 and PM2.5) are the most dangerous air pollutants for human health [1]. Air pollution is a major environmental issue in urban areas, which can affect the well-being and quality of life of citizens. More and more frequently, there are detected chronic diseases of great importance that are associated with continuous exposures to high concentrations of air pollutants. Epidemiological studies have reported that exposure to air pollutants such as particulate matter, nitrogen oxides (NOx), sulfur dioxide (SO2), and surface ozone (O3) associates with an increase in mortality and hospital admissions predominantly related to respiratory and cardiovascular diseases [1,2]. This critical issue concerns the 20 million people (including 9 million children) living in the Mexico City Metropolitan Area (MCMA). Despite the reductions in the emissions of common air pollutants in MCMA since the early 1990s, millions of people remain exposed to concentrations above the critical levels associated with increased risks for cardiovascular and respiratory diseases. The anthropogenic sources still produce and release to the atmosphere every year large amounts of carbon monoxide (CO), nitrogen oxides (NOx), sulfur dioxide (SO2), particulate matter (PM10 and PM2.5), and volatile organic compounds (VOC) [3]. However, because of the frequency of occurrence of high levels, persistence, and spatial distribution, the most critical air pollutants in MCMA are by far ozone and PM10 [4–6].

At the MCMA, the complexity of the air pollution problem is also strongly related to other essential factors such as the geographical setting, meteorology, and topography. MCMA lies inside a subtropical basin with latitudes between 19.2 and 19.6 ◦N and longitudes between 98.9 and 99.4 ◦W, and an average altitude of 2240 m, surrounded by high mountains (Figure 1). In the north direction, the basin extends into the Mexican plateau and the arid interior of the country, with the Sierra de Guadalupe creating a small 800 m barrier above the basin floor. Its climate classification comprises two seasons: the rainy season from May to October and the dry season from November to April. This classification stems from the two main meteorological patterns on the synoptic-scale: dry westerly winds with anticyclone conditions from November until April, and moist flows from the east due to the weaker trade winds along the other six months [7]. The MCMA meteorology, however, is by far more complicated than this simple classification expresses it. The basin interacts with the Mexican plateau and the lower coastal areas. Moreover, due to the MCMA location, large-scale pressure gradients are generally weak, and intense solar radiation is registered here throughout the year [7–11]. These conditions and the presence of high mountains in the surroundings are ideal for the development of thermally driven winds.

**Figure 1.** Complex topography in the Mexico City region.

The knowledge of wind circulation and air pollutant spatial distributions constitutes a relevant element to understand the flow dynamics of urban air pollution. It allows knowing how the air pollutant emissions produced in an urban settlement may disperse in the atmosphere, how these air pollutants may be distributed spatially in the city, and how the winds may transport them towards neighboring sites [12,13]. From the standpoint of fluid dynamics, the pollutant-flow vector, defined as the product of the fields of the pollutant concentration and wind velocity,

$$\mathbf{J}(\mathbf{x},t) = \mu(\mathbf{x},t)\mathbf{v}(\mathbf{x},t),\tag{1}$$

describes the transport of a given pollutant in the atmosphere. Then, the surface flow condition of a pollutant depends on both factors: the surface distribution of its concentration and the local wind circulation.

Local wind circulation in Mexico City and its relation to driving forces and air pollution have been studied extensively for almost three decades. However, most of the studies have been performed on an episode-by-episode basis and using small data sets obtained from short-term experimental campaigns with different approaches [14–21]. Some of the main exceptions are the first long-term micrometeorological campaign carried out by Salcido et al. [11] in this region during 2001, and the studies of Klaus et al. [22], de Foy et al. [23], Salcido et al. [24], and Carreón-Sierra et al. [25], which were carried out using data provided by the atmospheric monitoring network of Mexico City (SIMAT: Sistema de Monitoreo Atmosférico de la Ciudad de México [26]). Klaus et al. [22] reported a principal component analysis of air quality and meteorological data for the period February–December 1995. The study of de Foy et al. [23] reported an examination of wind transport during the experimental campaign of the Megacity Initiative: Local and Global Research Observations (MILAGRO) project; they used cluster analysis for making a comparison to climatology between the period of March 2006 and the period 1998–2006 using hourly wind data of the warm, dry season. In their work, Salcido et al. [24] used a lattice wind model at a meso-β scale to carry out a description and characterization of Mexico City local wind events of the period 2001–2006. Furthermore, Carreón-Sierra et al. [25] proposed an extension of the local wind state using a non-local description based on the wind velocity and its gradients and applied hierarchical cluster analysis to recognize and characterize the Mexico City wind circulation patterns in the period 2001–2006. As far as we know, the last paper reported the first work where a non-local approach is used to characterize the wind condition for identifying patterns with more detail than the simple results produced by the wind rose method.

In this paper, extending our previous work [25], we examine the main surface patterns of wind circulation and pollutant-flow of NO2, O3, SO2, and PM10 in the Mexico City region over the period 2001–2010. We achieved the identification and characterization of the pollutant-flow patterns through the exploitation of the 1-hour average values of the pollutant concentrations and wind data provided by the atmospheric monitoring network of Mexico City (SIMAT [26]) and the application of cluster analysis as a data mining procedure. First, we considered Mexico City as a 2D lattice domain and defined the flow condition (or flow state) of a given pollutant in each lattice cell in terms of a flow vector field and its gradients (precisely, the divergence and curl of this flow vector). Second, we used the 1-h average values of the pollutant concentrations and wind data provided by SIMAT to estimate, using Kriging interpolation, discrete representations of the pollutant-flow vector field and its gradients over the lattice domain. Finally, we used the k-means method of cluster analysis [27–29] to identify and characterize the pollution-flow patterns from the clustering modes of the flow states. We identified seven patterns for wind circulation and nine patterns for each of NO2, O3, PM10, and SO2 pollutant-flows. Their seasonal and annual average intensities and probabilities of occurrence were also estimated.

#### **2. Materials and Methods**

In this section, we detail the study domain, data, and procedure that we used to identify and characterize the wind and pollution-flow patterns.

#### *2.1. Study Domain*

We considered the part of Mexico City located at 19.3◦–19.6◦N and 99.0◦–99.3◦W as the study domain. In Figure 2, we showed this region as enclosed by a solid line square. This region contains 75% of the stations of the wind-monitoring network (REDMET) of SIMAT, 95% of the NO2 stations, 83% of the O3 stations, 93% of the PM10 stations, and 86% of the SO2 stations of the air quality monitoring network (RAMA) of SIMAT (Figure 3). Because of the topographic complexity of the mountains surrounding Mexico City, we examined the pollutant flows in this region, assuming it divided into quadrants (NE, NW, SW, and SE). We defined these quadrants using a reference frame with origin at (19.43◦N, 99.13◦W) and the axes along the west to east (W–E) and south to north (S–N) directions (dotted lines in Figure 2). The origin of the reference frame is the geometric centroid of the REDMET stations. The positions of the geometric centers of the city quadrants are: (19.5◦N, 99.1◦W), (19.5◦N, 99.2◦W), (19.4◦N, 99.2◦W), and (19.4◦N, 99.1◦W) for the NE, NW, SW, and SE quadrants. The dimensions of each quadrant were 15.66 km length in the west–east direction and 16.56 km length in the south–north direction (each quadrant is a square of 540 × 540 arcsec). In principle, we expect to detect the ventilation effects due to the openings located at the west and east sides of the Sierra de Guadalupe (north of the city) at NW and NE quadrants, respectively, and to observe in the SE quadrant the effect of the gap situated in the southeast of the Mexico basin. Also, in the SE and NE quadrants, we expect to recognize the effects of the winds blowing along the ventilation channel determined by

the volcanos (Sierra Nevada) on the east side of the city. Moreover, we expect to recognize the wind effects due to the Mexico City mountain-valley system in all quadrants, but mainly in the NW and SW quadrants.

#### *2.2. Data*

Figures 2 and 3 illustrate the study domain and the spatial distributions of the stations of the monitoring networks of SIMAT for wind, nitrogen dioxide, ozone, PM10, and sulfur dioxide. These monitoring networks provided the 1-hour average values of the meteorological and air quality variables systematically and made them publicly available through its web site [26]. For the present work, we collected the wind data (wind speed and wind direction) and the air quality data (NO2, O3, PM10, and SO2 surface concentrations) for the period 2001–2010, which constitute a database with 87,600 records, approximately, for each hourly averaged variable.

**Figure 2.** Study domain for Mexico City and spatial distribution of the stations (white dots) of the wind-monitoring network (REDMET) of the atmospheric monitoring network of Mexico City (SIMAT). The REDMET provides 1-hour average values of wind speed and wind direction, and other properties such as temperature and relative humidity.

**Figure 3.** Spatial distribution of the stations (white dots) of the NO2, O3, PM10, and SO2 monitoring networks (RAMA) of SIMAT. The RAMA provides 1-hour average values of the surface concentrations of NO2, O3, PM10, and SO2, among other chemical compounds.

#### *2.3. Study Database*

In a first step, we represented the study domain (*D*) as a square grid *G* with 289 nodes *Grs*, where the subscripts *r* = 1, 2, ... , 17 and *s* = 1, 2, ... , 17 vary along with the West–East (WE) and South–North (SN) directions, respectively. Using the wind speed and wind direction data supplied by the REDMET stations, and the NO2, O3, PM10, and SO2 concentration data supplied by the RAMA stations, we estimated (using Kriging interpolation) the horizontal components of the flow vector field at each grid node *Grs* and time *t*,

$$I\_{\rm WE}(r, s, t) = \mu(r, s, t) v\_{\rm WE}(r, s, t), \qquad I\_{\rm SN}(r, s, t) = \mu(r, s, t) v\_{\rm SN}(r, s, t) \tag{2}$$

For a given pollutant, at the node *Grs* and time *t*, μ(*r*,*s*,*t*) is the pollutant concentration, and *vWE*(*r*,*s*,*t*) and *vSN*(*r*,*s*,*t*) are the wind velocity components. Here, time *t* is an integer running from 1 to the number *H* of hours (87,600, approximately) in the period 2001–2010.

Then, we defined a 2D lattice *L* composed by 64 regular non-overlapping cells *Cij* covering *D*, with *i* = 1, 2, 3, ... , 8 and *j* = 1, 2, 3, ... , 8, such that each cell *Cij* is surrounded by nine grid nodes of G, as Figure 4 illustrates it.

**Figure 4.** Arrangement of the calculation grid nodes *Grs* and the lattice cells *Cij*. The green crosses represent the grid nodes, and the red crosses indicate the centers of the lattice cells. In gray, we show the cell *Cij*.

The estimations, for each pollutant, of the flow vector components at the grid nodes *Grs*, allow calculating (numerically) the spatial averages of the components of the pollutant-flow vector and its gradients at the cells *Cij* of the lattice *L* using their 2D discrete definitions in terms of centered finite differences [24,30]. For a given pollutant, we will denote the flow vector components at *Cij* by *jx*(*i*, *j*, *t*) and *jy*(*i*, *j*,*t*), and the associated divergence and curl of the flow vector will be denoted by Γ(*i*, *j*,*t*) and Ω(*i*, *j*,*t*). Here, *x* and *y* denote the WE and SN directions, respectively. At this level, we will be referring to this non-local discrete representation of the flow condition of the study domain as the 64-cell model. It involves 256 (= 4 × 64) parameters to define an *event* of the system (i.e., to describe the flow condition over the complete spatial domain). In this work, for simplicity in the handling and interpretation of the results, we considered the reduction of the 64-cell model to the cases of the 4-cell and 1-cell models, which we obtained using a spatial averaging procedure. In any case, the four parameters (*jx*, *jy*, Γ, Ω) define the flow state at a lattice cell; therefore, the 4-cell and 1-cell models respectively entail 16 and 4 parameters to describe an *event* of the system. It is convenient to note that the introduction of the additional variables (Γ, Ω) in the flow state definition, allows recovering some of the information lost during the spatial averaging of the flow vector components. This assumption gives a slight non-local character to this description [31].

For each cell *Cij*, we calculated, hour by hour, the average of the four state-variables over the 2001–2010 period, defining a one generic year database for each pollutant and cell model considered. The results were normalized and then organized in seasonal groups: January–March (winter), April–June (spring), July–September (summer), and October–December (autumn).

We performed data normalization as follows. We determined the maximum absolute values of the magnitude of the pollutant-flow vector |**j**|max, and the divergence |Γ|max, and curl |Ω|max of the flow vector, in the generic year database. Then the normalized and dimensionless state variables (all of them now ranging in the interval [−1, +1]) were calculated using the following expressions:

$$\hat{j}\_{\mathbf{x}}(i,j,t) = \frac{j\_{\mathbf{x}}(i,j,t)}{\left|\mathbf{j}\right|\_{\max}},\tag{3}$$

$$f\_y(i, j, t) = \frac{j\_y(i, j, t)}{\left|\mathbf{j}\right|\_{\text{max}}},\tag{4}$$

$$
\hat{\Gamma}(i, j, t) = \frac{\Gamma(i, j, t)}{|\Gamma|\_{\text{max}}},\tag{5}
$$

$$
\hat{\Omega}(i, j, t) = \frac{\Omega(i, j, t)}{|\Omega|\_{\text{max}}}.\tag{6}
$$

#### *2.4. Clustering Method*

We used the k-means method of cluster analysis as a procedure of data mining to identify and characterize the main patterns of wind and flow of pollutants in the Mexico City region.

Cluster analysis comprises a broad set of techniques for finding non-overlapping subgroups of data objects in a dataset. When clustering the data objects of a dataset, we seek to organize them into a given number of distinct subsets (groups, or clusters), so that they constitute a partition of the initial dataset; this way, the data objects within each subset are quite similar to each other, while data objects in different groups are quite different from each other [29]. Clustering is an unsupervised problem because we are trying to discover structure (distinct clusters) based on a dataset, without being trained by a response variable.

The k-means clustering method is the simplest and the most frequently used technique of cluster analysis for partitioning a dataset into a set of k groups. In k-means clustering, we seek to partition the set of data objects into a pre-specified number of clusters, so that the total within-cluster variation (TWCV) is as small as possible.

Suppose we have a dataset *<sup>X</sup>* = {*x*1, *<sup>x</sup>*2, ... , *xN*}, *xi* <sup>∈</sup> **<sup>R</sup>***d*, and we try to partition *<sup>X</sup>* into *<sup>M</sup>* disjoint clusters *C*1, *C*2, *CM*. Then, the k-means algorithm detects local optimal solutions concerning the TWCV defined as the sum of squared Euclidean distances between each data object *xi* and the centroid *mk* of the cluster *Ck* that contains *xi*. Analytically, the TWCV is given by

$$\varepsilon(m\_1, m\_2, \dots, m\_M) = \sum\_{i=1}^{N} \sum\_{k=1}^{M} \lambda\_{ik} \|\mathbf{x}\_i - m\_k\|^2 \tag{7}$$

where λ*ik* = 1 if *xi* ∈ *Ck*, and 0 otherwise. The TWCV measures the clustering goodness, and the optimization problem that defines the k-means clustering resides in minimizing TWCV. In practice, we used the software *DataLab* developed by Hans Lohninger [32] to perform the cluster analysis with the k-means algorithm.

For each air pollutant considered, we define its set of data objects as

$$H = \{ \mathbb{Q}(t) | t = 1, 2, \dots, T \}. \tag{8}$$

where T = 8760 is the number of hours of the generic year, and

$$Q(t) = \left\{ Q\_{\hat{i}\hat{j}t} \equiv (\hat{j}\_{\mathbf{x}}, \hat{j}\_{\mathbf{y}}, \hat{\Gamma}, \Omega)\_{\hat{i}\hat{j}t} \middle| i = 1, 2, \dots, N \mathbf{x}; j = 1, 2, \dots, Ny \right\},\tag{9}$$

with *Nx* = 2 (*Nx* = 1) and *Ny* = 2 (*Ny* = 1) in the case of the 4-cell (1-cell) model. We organized the dataset *H* in the seasonal subsets *Hwinter*, *Hspring*, *Hsummer*, and *Hautumn* comprising, respectively, the data objects of the periods January–March (winter), April–June (spring), July–September (summer), and October–December (autumn). We identified each data object of these subsets with a label of the form MMDDHH that specifies the date (relative to the generic year) and hour of occurrence.

To carry out the cluster analysis for each seasonal period, we build a data matrix structured as described in Table 1.


**Table 1.** Structure of the data matrix for the cluster analysis.

For each pollutant, using the respective data matrixes of the seasonal subsets, we can apply the k-means algorithm for clustering the associated data objects considering the parameters of the 1-cell or 4-cell model as the clustering variables. Although we can expect that the 4-cell model variables will produce a physically more detailed clustering of the wind and pollutant-flow events, in this work, for easiness in handling the results, we used only the 1-cell model parameters as clustering variables. The application of the clustering procedure to the 1-cell model data objects will organize them, assigning the cluster number to the date-time labels that identify the events. Consequently, the 4-cell model data objects also receive cluster numbers correspondingly. It means that the data objects of the 4-cell model were also organized into clusters by the same clustering procedure. Therefore, for each pollutant and seasonal period, we have clusters of events characterized by the sixteen parameters of the 4-cell model:

$$(\hat{j}\_{\mathbf{x}}, \hat{j}\_{\mathbf{y}}, \hat{\Gamma}, \hat{\Omega})\_{\mathbf{C}0'} (\hat{j}\_{\mathbf{x}}, \hat{j}\_{\mathbf{y}}, \hat{\Gamma}, \hat{\Omega})\_{\mathbf{C1'}} (\hat{j}\_{\mathbf{x}}, \hat{j}\_{\mathbf{y}}, \hat{\Gamma}, \hat{\Omega})\_{\mathbf{C2'}} (\hat{j}\_{\mathbf{x}}, \hat{j}\_{\mathbf{y}}, \hat{\Gamma}, \hat{\Omega})\_{\mathbf{C3}} \tag{10}$$

where C0, C1, C2, and C3 are the cells of the model.

#### *2.5. The Number of Clusters*

The k-means clustering is a simple and fast algorithm that can efficiently deal with large datasets. However, the k-means approach requires pre-specifying the number of clusters; and preferably, we would like to use an optimal number of clusters that we defined according to some given criterion. In this work, because of the motivation expressed in the following two paragraphs, we applied the k-means clustering method to organize each seasonal period in six clusters of data objects (events).

#### 2.5.1. Physical Motivation

In Figure 5, we presented the mean diurnal behavior of solar radiation, temperature, and wind speed observed during 2001. We obtained these results from data we measured at an urban site southeast Mexico City (Xochimilco) [11,25]. The plots evidence that the meteorological events comprised in the six time-periods 0–4, 4–8, 8–12, 12–16, 16–20, and 20–24 h, are quite different from each other. The periods 0–4 and 20–24 h show the cooling phase of the atmosphere during the night. In the period 4–8 h, we observe the sunrise occurrence, and that temperature and wind speed reach their minimum values. The period 8–12 h shows the growing phases of solar radiation, temperature, and wind speed, with solar

radiation reaching its maximum. The period 12–16 h depicts the temperature reaching its maximum, while wind speed keeps growing, and solar radiation starts to decrease. Finally, the period 16–20 h shows the sunset occurrence, the wind speed reaching its maximum, and temperature decreasing. Appreciation of these different behaviors suggests searching for six groups during the application of the k-means algorithm for clustering the data objects defined by the events of wind circulation and pollutant-flow in Mexico City.

**Figure 5.** Average diurnal behavior of solar radiation (red line), temperature (green line), and wind speed (blue line) in an urban site (Xochimilco) of Mexico City during 2001 [11,25]. We observe that the meteorological events comprised in the six intervals of 0–4, 4–8, 8–12, 12–16, 16–20, and 20–24 h are different from one period to another. This appraisal suggests searching for six clusters during the cluster analysis of the wind and pollutant-flow events.

#### 2.5.2. The Elbow Method

From the standpoint of cluster analysis, one of the most popular methods for determining the optimal clusters is the *elbow method*. It provides a simple solution. The basic idea is to compute the k-means clustering algorithm using different numbers, *k*, of clusters. The next step is drawing the total within-cluster variation (or total within-cluster sum of squares) versus the number of clusters. Commonly, the analysts consider the location of a bend in the plot as an indicator of the appropriate number of clusters. However, it is not always clear where the bend point position is. We preferred to calculate the percentage reduction (δ) of the total within-cluster sum of squares of a given *k*, relative to its value for the previous number of clusters *k* − 1, and then to plot δ as a function of *k*. As our optimal number of clusters, we selected the value of *k* from which on the reduction δ is less than 10% each unit step in increasing *k*. For example, according to this procedure, Figure 6 shows that we must select *k* = 6 when clustering the winter wind data objects of the present work.

**Figure 6.** Percentage reduction (δ) of the total within-cluster sum of squares, relative to its previous value, as a function of the number of clusters. Here, for the clustering of the winter wind data objects, we observe that from 6 clusters on, the reduction δ is less than 10% each unit step in the increase of the number of clusters: 0% ≤ δ ≤ 10%.

#### **3. Results and Discussion**

Figure 7 summarizes the application of the k-means clustering method (with k = 6) to the seasonal subsets, *Hwinter*, *Hspring*, *Hsummer*, and *Hautumn*, of the data objects (events) of wind circulation and pollutant-flows. This figure shows an arrangement with four columns (one per seasonal period) and five rows (one for wind circulation, and four for the pollutant-flows of NO2, O3, PM10, and SO2). Each graph shows six plots (one per cluster) differentiated by the line color. The algorithm enumerates the clusters according to their sizes, from the larger to the smaller, so colors mean no more. Each plot presents the *hourly population* of one cluster, that is to say, the number of data objects belonging to the cluster as a function of the hour of the day (Mexico City local time, UTC-6 h).

**Figure 7.** *Cont*.

**Figure 7.** Application of the k-means algorithm (with k = 6) to the seasonal subsets of the data objects of the wind circulation (first row) and pollutant-flows (next four rows) events. Each graph shows six plots, one per cluster, individualized by the color. Each plot expresses the *hourly population of the cluster*, i.e., the number of data objects (or events) of the cluster as a function of the hour of the day.

For wind circulation or any given pollutant-flow, a qualitative comparison of graphs in the respective row allowed recognizing some seasonal similarities and differences between the clusters. For example, for the wind circulation (first row), if we pay attention to the clusters represented by a blue-line plot in the four seasons, we can observe an evident regularity: the hourly populations have similar trends and indicate that the majority of their events occurred during the night, mainly from midnight until dawn. We interpreted this regularity as the possible existence of a wind pattern. The plots of the pollutant-flows also show a similar regularity.

The trends of the hourly populations of the clusters provide useful but insufficient information for discovering patterns in sets of events of wind circulation or pollutant-flows. In [25], the authors detected two wind clusters with similar trends in the hourly populations (both in the number of events and in the times of occurrence), but in one case the winds were blowing from the north and in another one from the south, indicating different patterns of wind circulation. However, we can use the information provided by the values of the 16 parameters of the 4-cell model to complement the pattern identification process.

Eight of the sixteen parameters of the 4-cell model are the components of the flow vectors in the quadrants (NW, NE, SW, and SE) of the MCMA. This complementary information and the hourly population plots allowed recognizing the wind circulation patterns shown in Figure 8, and the pollution-flow patterns shown in Figures 9–12. We found seven patterns for wind circulation and nine patterns for each of the pollutant-flows (NO2, O3, PM10, and SO2). In Figure 8, we show the hourly population of each wind pattern and a graphic representation of the wind velocity at the quadrants. Figures 9–12 show the hourly population of each pattern, a graphic representation of the flow vector at the quadrants, the surface distribution of the pollutant concentration, and the mean spatial pollutant concentration (expressed in ppb for NO2, O3, and SO2, and in μgm−<sup>3</sup> for PM10). In these figures, the hourly population of a pattern is an average over the seasons the pattern occurred; the vectors (the wind velocity and the pollutant-flow vectors) are averages over the elements of the pattern and the seasons; the spatial distribution of the pollutant concentrations is the average over the pattern elements and seasons; and the mean spatial concentration of a pollutant is the average over the positions of the average pollutant distribution. The lengths of the arrows that represent the wind and flow vectors were scaled to get the larger one fitted in the square that represents the quadrant, allowing a qualitative comparison of the vector magnitudes among the patterns of wind circulation or the patterns of any of the pollutant flows.

#### *3.1. Wind Circulation Patterns*

The proposed methodology for examining the database of the wind events that occurred in Mexico City from 2001 to 2010 allowed discovering there the presence of seven wind circulation patterns and the estimation of their seasonal and annual occurrence frequencies, which we summarized in Table 2. We denoted these patterns as WIND-Pn, with n = 1, 2, ... , 7, and enumerated them trying to follow

the order of their appearance during the day. In Figure 8, we presented the hourly population of each wind pattern, including a sketch out of the corresponding wind velocity in the city quadrants.


**Table 2.** Seasonal and annual frequencies of the wind circulation patterns (%).

In the following paragraphs, we provide brief descriptions of the wind circulation patterns. We used the Mexico City local time (UTC-6 h) in all figures and descriptions.

WIND-P1: early morning downslope winds. WIND-P1 was the wind pattern with the highest annual frequency (30%). It was observed systematically throughout the year during the night, showing a large peak around hour 6 (Figure 8) and a small peak around hour 15. The mean wind velocities at the quadrants suggest downslope winds from the surrounding mountains, converging towards Mexico City basin. The converging winds represented by the small peak reflect the presence of convective upwards winds due to the urban heat island (UHI) in the Mexico City region. Jauregui [7] and Klaus et al. [33] pointed out that the Mexico City downslope winds are reinforced by UHI.

WIND-P2: northeasterly and easterly winds. The annual average frequency of this wind pattern was 21%. It is observed systematically throughout the year, especially during winter (25%) and summer (24%). The mean velocities correspond to winds blowing from NE in almost all the quadrants of the city (trade winds). The wind events of this pattern occurred from the sunrise (hour 7) up to the hour 17. Its hourly population presents a high peak between hours 10 and 11.

WIND-P3: midday northerly winds. This pattern had an annual frequency of 9%. It was observed during spring (12%), summer (9%), and autumn (14%) from hour 7 to hour 21, with a peak around hour 13. The mean wind velocities in the quadrants indicate winds blowing from the north and northeast sectors.

WIND-P4: westerly and southerly winds. This wind pattern was observed only winter season with a seasonal frequency of 20%; its events occurred throughout the day, but mainly from noon to midnight, reaching its maximum population around hour 17. The mean wind velocities represent winds with a westerly component in the west quadrants and southerly winds in the east quadrants. These velocities correspond to winds blowing from the south and southeast city sectors through a gap between the mountains of the Sierra Ajusco-Chichinautzin and the Sierra Nevada, following the ventilation channel S–N located at the east side of the city (Figure 1). The local winds of this pattern could be related to the subtropical jet stream of winter [34] or with the westerlies that are permanently occurring in subtropical and middle latitudes, but coupled with the winds driven by the gap between the mountains located at the southeast. Doran and Zhong [35] described the main characteristics of a gap wind system in the southeastern corner of Mexico City that produces low-level jets occurring regularly during the late winter.

WIND-P5: afternoon northerly winds. This pattern was observed throughout the year with an annual frequency of 12%. It occurred from noon to midnight, with its maximum around hour 19. The mean wind velocities at the quadrants indicate winds blowing from the north sectors.

WIND-P6: midnight downslope winds. This pattern was observed systematically throughout the year, mainly during the nighttime hours, from the sunset (hour 18) to the sunrise (hour 6 of the next day), with its maximum around midnight. It got an annual average frequency of 18%. The population plot of this pattern has a small peak close to hour 15, which could be related to the UHI effect of the city. The wind velocities in the western quadrants of the city represent downslope flows from Sierra de las Cruces (at west) and Sierra del Ajusco-Chichinautzin (located at SW and S) toward the town. In the eastern quadrants, otherwise, late afternoon northerly gap winds are observed guided by the S–N ventilation channel of the city.

WIND-P7: UHI-driven winds. This wind pattern was observed only during the spring (9%) and autumn (10%) seasons, with a low annual occurrence frequency of 5%. The wind events of this pattern occurred from the hour 14 to the hour 22, and its population reached a maximum value close to the hour 17. The mean wind velocities of the city quadrants indicate cyclonic winds converging towards Mexico City downtown, which seems to be related to an upwards flow driven by the UHI.

Summarizing, the wind patterns with the most significant frequencies were WIND-P1 and WIND-P2, with annual values of 30% and 21%, respectively. These two patterns, however, correspond to small velocities in comparison with the other five wind patterns. An outstanding feature of WIND-P1 is that, although its events occurred mainly during the night with converging winds, it also comprises daylight events of winds with the same characteristic of convergence towards the city downtown. We cannot interpret these last events as katabatic winds produced by the mountain-valley system; instead, we must understand them as winds associated with convective upwards flows driven by the urban heat island (UHI) effect. The pattern WIND-P6 also occurred during the night, but it peaked at midnight. This wind pattern comprises the downslope winds of the first hours of the night and also reflects the effect of the late afternoon northerly gap winds. The patterns WIND-P3 and WIND-P5 reflect the well-known prevailing winds from N and NE of the Mexico City region, and WIND-P2 evidence the occurrence of the trade winds.

**Figure 8.** *Cont*.

**Figure 8.** Wind circulation patterns in the Mexico City region during the period 2001–2010. For each pattern, we presented the hourly population (left) and the mean wind velocity (relative to the pattern) at the city quadrants (right). The size of red arrows that stand in for the velocity vectors indicates the velocity magnitude in a scale relative to the edge size of the squares that represent the quadrants.

The wind patterns with the most significant frequencies were WIND-P1 and WIND-P2, with annual values of 30% and 21%, respectively. These two patterns, however, correspond to small velocities in comparison with the other five wind patterns. An outstanding feature of WIND-P1 is that, although its events occurred mainly during the night with converging winds, it also comprises daylight events of winds with the same characteristic of convergence towards the city downtown. We cannot interpret these last events as katabatic winds produced by the mountain–valley system; instead, we must understand them as winds associated with convective upwards flows driven by the urban heat island (UHI) effect. The pattern WIND-P6 also occurred during the night, but it peaked at midnight. This wind pattern comprises the downslope winds of the first hours of the night and also reflects the effect of the late afternoon northerly gap winds. The patterns WIND-P3 and WIND-P5 reflect the well-known prevailing winds from N and NE of the Mexico City region, and WIND-P2 evidence the occurrence of the trade winds.

Our wind patterns agree quite well with the Mexico City wind patterns obtained in [25] for the period 2001–2006 with a slightly different methodology. In Table 3, we described the correspondence between these two sets of patterns. The first column contains the wind patterns WP1-WP7 reported in [25]; the second column presents the wind patterns obtained in the present work; and, the last column provides a brief description of the main characteristics of these patterns.


**Table 3.** Comparison against the wind patterns reported by Carreón-Sierra et al. [25].

The main differences between both sets of patterns were:

• The pattern WIND-P4 includes the wind patterns WP4 (southerly winds) and WP5 (westerly winds) reported in [25]. However, we detected the pattern WIND-P4 only during the winter season, while the pattern WP4 occurred throughout the year (although with tiny frequencies in spring, summer, and autumn), and the pattern WP5 was only observed during the first semester of the year (winter and spring seasons).

• We recognized a pattern (WIND-P7) associated with cyclonic winds strongly converging towards Mexico City downtown during the daylight hours. The winds in this pattern indicate an upwards flow driven by the UHI. No report exists in [25] about a similar pattern.

#### *3.2. Pollutant Flow Patterns*

In Tables 4–7, we enumerated the pollutant-flow patterns (NO2, O3, PM10, and SO2) that we identified for the Mexico City region, their seasonal and annual frequencies and flow intensities (the magnitudes of the flow vectors), and the associated wind circulation patterns. In Figures 9–12, for each pollutant-flow pattern, we presented the average population plot, the mean pollutant-flow vectors at the quadrants NW, NE, SW, and SE of the city, the pollutant surface distribution averaged over the pattern, and its mean spatial value. We denoted the pollution-flow patterns as POL-Pn, where POL indicates the pollutant (POL = NO2, O3, PM10, and SO2), and n is a consecutive integer number from 1 to 9. The cells with the value 0.00 in the tables indicate the seasonal periods where some flow patterns were not detected.

In general, although the pollutant concentrations modulate the magnitudes of the flow vectors, the pollutant-flow patterns reflected the directional characteristics inherited from the wind patterns.

#### 3.2.1. The NO2 Flow Patterns

We identified nine NO2 flow patterns for the Mexico City region. Table 4 shows the seasonal and annual average frequencies and average flow intensities of the patterns, expressed in percent (%) and μgm<sup>−</sup>2s−1, respectively. As shown in Figure 9, these flow patterns closely resemble the wind circulation patterns, and Table 4 summarized the relationship between the pollutant-flow and wind circulation patterns.


**Table 4.** Frequencies (%) and flow intensities (μg/m2s) of Mexico City NO2 flow patterns (2001–2010).

The NO2 flow patterns with the most significant annual frequencies were NO2-P1 (30%), NO2-P4 (20%), and NO2-P7 (12%), but the patterns with the most substantial annual flow intensities were NO2-P5, NO2-P7, and NO2-P9 with 25.4, 25.3, and 24.5 μgm<sup>−</sup>2s−1, respectively.

The flow patterns NO2-P4, NO2-P7, NO2-P8, and NO2- P9 carry nitrogen dioxide from the north to the south quadrants of the city, particularly the pattern NO2-P9, although it occurred only during the second semester of the year. The events of these patterns take the ozone precursor to the SW and SE quadrants, contributing actively to the ozone formation in this area throughout the year.

The patterns NO2-P2 and NO2-P5 carry the pollutant from the eastern to the western quadrants of the city following the trade winds. The events of the NO2-P3 and NO2-P6 patterns, differently, take nitrogen dioxide from west to east on the west side of the city, but from south to north on the eastern side. However, while the pattern NO2-P3 reflects a coupling of the westerly winds and the afternoon southerly gap winds, the pattern NO2-P6 reflects cyclonic transport driven by the UHI.

During the night, the events of the pattern NO2-P1, follow the downslope winds from the surrounding mountains and carry NO2 to the city downtown. In comparison with the patterns of other pollutants (O3, PM10, and SO2), the NO2-P1 flow pattern is the only one that revealed considerable nocturnal transport due to the downslope winds. The main reason is that NO2 accumulates during the night since there are no photochemical reactions that consume it, and some of the mobile sources, which are responsible for approximately 86% of nitrogen oxide emissions in the Mexico City region [3], remain active during the night.

The NO2 surface distributions, which are averages over the events of the NO2 flow patterns, reveal a North–South part of the city as the zone with the higher NO2 concentrations. This zone extends around the N–S axis that separates the west from the east quadrants (although slightly shifted to the eastern quadrants). The highest NO2 levels occurred, in general, in the south sector, close to the Huizachtepetl (Cerro de la Estrella).

The surface NO2 distributions with the highest concentrations were in connection with the NO2-P5 and NO2-P6 flow patterns, with mean spatial levels of 46 and 41 ppb. These high levels of NO2 were detected mainly during the winter and autumn. We note that, in these two cases, the high levels of NO2 extended spatially, covering almost completely the city. In the first case, it was a consequence of blocked transport by the mountains of the Sierra las Cruces, while in the second case, the high levels were a consequence of the cyclonic wind convergence driven by the UHI.

**Figure 9.** *Cont*.

**Figure 9.** *Cont*.

**Figure 9.** Flow patterns of nitrogen dioxide in the Mexico City region during the period 2001–2010. For each pattern, we included the hourly population (second column), the mean flow vectors (relative to the pattern) at the city quadrants (third column), the surface concentration distributions (fourth column) averaged over the events of the flow pattern, and the mean spatial concentration (last column). The size of red arrows that stand in for the NO2 flow vectors indicates the flow intensity in a scale relative to the edge size of the squares that represent the quadrants.

#### 3.2.2. The O3 Flow Patterns

In Table 5 and Figure 10, we summarized the main characteristics of the O3 flow patterns that we identified. In Table 5, we included the average seasonal and annual occurrence frequencies (%) and flow intensities (μg m−2s−1) of the flow patterns, and the wind patterns from which they inherited their circulation characteristics.

The flow vectors of the diurnal flow patterns O3-P6, O3-P8, and O3-P9 have considerable northerly components and convey ozone from the NW and NE to the SW and SE city quadrants, providing substantial contributions to the high ozone concentrations frequently observed in the south sectors (especially in SW) of the city throughout the year.


**Table 5.** Frequencies (%) and flow intensities (μg/m2s) of the Mexico City O3 flow patterns (2001–2010).


**Figure 10.** *Cont*.

**Figure 10.** Flow patterns of ozone in the Mexico City region during the period 2001–2010. For each flow pattern, the figures show the hourly population (second column), the mean flow vectors (relative to the pattern) at the city quadrants (third column), the surface concentration distributions (fourth column) averaged over the events of the flow pattern, and the mean spatial concentration (last column). The size of red arrows that stand in for the O3 flow vectors indicates the flow intensity in a scale relative to the edge size of the squares that represent the quadrants.

The flow patterns O3-P6, O3-P8, and O3-P9 revealed the highest flow intensities (43, 82, and 74 μgm<sup>−</sup>2s−1, on annual average). Seasonally, we observed the highest O3 flow intensities during the spring, summer, and autumn.

The flow patterns O3-P3 and O3-P5 bring ozone from East to West in the city following the trade winds (pattern WIND-P2). These patterns reveal flow vectors from NE in the west quadrants of the city, and because of the mountains barrier (Sierra Las Cruces and Sierra Ajusco Chichinautzin), they contribute to the ozone accumulation at SW sector of the Mexico City region. The O3 surface distributions, averaged over the events of these flow patterns, reveal, in fact, an evident and significant ozone accumulation (ranging from 70 to 100 ppb, approximately) at the SW quadrant.

The flow patterns O3-P1 and O3-P7 comprised nighttime flow events and had the first and second-largest occurrence frequencies. However, their flow intensities are tiny because of the low ozone concentrations (12 ppb and 22 ppb on average, respectively) in the nocturnal atmosphere (ozone production is of photochemical nature). It is interesting to note that the patterns O3-P1 and

O3-P7 and the patterns NO2-P1 and NO2-P4 reveal opposite behaviors (during the night, one observes NO2 accumulation while O3 decreases to its lowest levels) due to its relationship through the atmospheric photochemistry.

The O3-P2 is a flow pattern detected during the evening, driven by the UHI winds. It brings ozone from the surrounding parts of the city towards downtown, where an upwards convective flow takes ozone out there again, keeping the ozone levels relatively small in the city.

During winter, the flow events of the O3-P4 pattern bring ozone from west to east in the west side quadrants with small flow intensities, but the same flow pattern carries a considerable amount of ozone from south to north in the east side quadrants of the city, following the ventilation channel located at the west side of Sierra Nevada.

In general, the O3 surface distributions show relative small concentrations where the NO2 surface distributions show large emissions and vice versa. It is an evident behavior because of the NOx-O3 photochemical interactions.

#### 3.2.3. The PM10 Flow Patterns

We recognized nine PM10 flow patterns, which we enumerated in Table 6 and sketched out in Figure 11. In Table 6, we presented the seasonal and annual averages of the occurrence frequencies and the flow intensities of the PM10 flow patterns, expressed in % and μg/m2s, respectively. The directional characteristics of wind inherited by the PM10 patterns are also briefly indicated in this table.

**Table 6.** Frequencies (%) and flows (μg/m2s) of the Mexico City PM10 flow patterns (2001–2010).


**Figure 11.** *Cont*.


**Figure 11.** *Cont*.

**Figure 11.** Flow patterns of PM10 in the Mexico City region during the period 2001–2010. For each flow pattern, the figures show the hourly population (second column), the mean flow vectors (relative to the pattern) at the city quadrants (third column), the surface concentration distributions (fourth column) averaged over the events of the flow pattern, and the mean spatial concentration (last column). The size of red arrows that stand in for the PM10 flow vectors indicates the flow intensity in a scale relative to the edge size of the squares that represent the quadrants.

The PM10 flow patterns with the highest annual flow intensities were PM10-P7 (43 μgm−2s−1), PM10-P6 (39 μgm−2s−1), PM10-P5 (35 μgm−2s−1), and PM10-P8 (26 μgm−2s−1); but their annual frequencies were of the smaller (3%, 4%, 2%, and 5%, respectively). All these four patterns reveal a remarkable flow intensity at the NE city quadrant throughout the year, although with different flow directions. The patterns PM10-P6, PM10-P7, and PM10-P8 exhibited also flow vectors with intense northerly components in the NW quadrant, which convey particle PM10 from the north to the south in the west side of the city. The events of these patterns also carried particulate matter from the northeast to the south sectors of the city, particularly in the NE quadrant. The flow intensities of the flow pattern PM10-P6 were 71, 56, and 27 μgm<sup>−</sup>2s−<sup>1</sup> during the winter, spring, and summer seasons; for PM10-P7 were 110 and 63 μgm<sup>−</sup>2s−<sup>1</sup> during the spring and summer, and for PM10-P8 were 43 and 55 μgm<sup>−</sup>2s−<sup>1</sup> during the summer and autumn, respectively.

The events of the PM10-P4 and PM10-P5 flow patterns carried particles from south to northwest in the quadrants of the east side of the city: PM10-P4 with flow intensities of 29 and 27 μgm−2s−<sup>1</sup> during the winter and spring, and PM10-P5 with flow intensities of 67 and 71 μgm−2s−<sup>1</sup> during the winter and autumn, respectively.

The events of the flow patterns PM10-P1 and PM10-P2 occurred throughout the year with the highest frequencies of occurrence (44% and 22%, on an annual average). The flow intensity of the PM10-P1 (a nocturnal pattern) was too small all year, but the flow intensity of the PM10-P2 (a diurnal pattern) was similar to those of the other flow patterns. The patterns PM10-P1 and PM10-P3 reflect the influence of the downslope winds driven by the mountain-valley system. These patterns, furthermore, reflect the effect of the urban heat island, which is revealed in the hourly population by the small peak around the hour 15 (daylight winds converging to the downtown).

For all the PM10 flow patterns, the average surface distributions of the pollutant reveal the NE quadrant as the zone of the city with the highest PM10 concentrations, and the SW quadrant as the zone with the lowest levels of PM10. It suggests that, at the NE of the town, in the surrounding area of Xalostoc on the east side of the Sierra de Guadalupe, there exists considerable sources of particulate matter, which release PM10 to the atmosphere all year long.

It is interesting to observe that the events of the pattern PM10-P5 revealed flow vectors at the NE and SE quadrants with intense southerly components, which produced, on average, the highest mean spatial concentration of PM10 (81 μgm<sup>−</sup>3) in the city. This result displays the east side quadrants of the city as the most polluted areas by PM10 during the winter and autumn.

#### 3.2.4. The SO2 Flow Patterns

Table 7 and Figure 12 enumerate and sketch out the SO2 flow patterns that we identified. Table 7 presents the average seasonal and annual occurrence frequencies and flow intensities of these patterns, including their relationship with the wind circulation patterns. Figure 12 shows the hourly population of the patterns, the corresponding flow vectors, and the surface distribution of the SO2 emissions produced by the events of these flow patterns.


**Table 7.** Frequencies (%) and flows (μg/m2s) of the Mexico City SO2 flow patterns (2001–2010).

We observed, in general, that all the SO2 flow patterns presented small flow intensities and produced small surface concentrations (we underline that the Mexican air quality standard for SO2 is 110 ppb on a 24 h average, and 25 ppb on an annual average). Nevertheless, it is particularly interesting to observe that the flow patterns SO2-P3 (winter, spring, and autumn), SO2-P4 (winter, summer, and autumn), SO2-P5 (winter and autumn), and SO2-P9 (spring and summer), which were detected mainly during the night hours, exhibit the flow vector at the NW quadrant with a flow intensity relatively larger than in the other quadrants. However, the corresponding wind patterns show winds of similar magnitudes in all the quadrants.

**Figure 12.** *Cont*.

**Figure 12.** *Cont*.

**Figure 12.** Flow patterns of sulfur dioxide in the Mexico City region during the period 2001–2010. For each flow pattern, the figures show the hourly population (second column), the mean flow vectors (relative to the pattern) at the city quadrants (third column), the surface concentration distributions (fourth column) averaged over the events of the flow pattern, and the mean spatial concentration (last column). The size of red arrows that stand in for the SO2 flow vectors indicates the flow intensity in a scale relative to the edge size of the squares that represent the quadrants.

The events of the flow patterns SO2-P7 (spring, summer, and autumn) and SO2-P8 (spring and summer) occurred during daylight hours and show the NW-quadrant flow vector larger than in the other quadrants. Consistently with these observations, the mean spatial SO2 surface distributions reveal the NW quadrant of the city as the zone with the highest levels of concentration, particularly in the case of the patterns SO2-P4 and SO2-P5. This situation seems to indicate that in the NW quadrant of the city, there are intense activities that release sulfur dioxide to the atmosphere during the night.

The SO2-P1 and SO2-P2 flow patterns, as was the case with the other pollutants, were the most frequent patterns (33% and 19%, respectively, on annual average) and the only ones detected all over the year; however, their flow intensities were tiny.

#### **4. Conclusions**

We used a simple model to define the pollutant-flow conditions in terms of the horizontal components of the flow vector and its gradients (divergence and curl). We applied this model to the hourly data of wind and pollutant concentration measured during the period 2001–2010 by the atmospheric monitoring network of Mexico City. We obtained the main flow patterns of NO2, O3, PM10, and SO2 using a 4-cell model of the city and the k-means algorithm of cluster analysis. Estimations of the seasonal and annual occurrence frequencies and flow intensities of the flow patterns were obtained.

The pollution flow patterns reflect some of the wind circulation conditions of the Mexico City wind patterns; however, since the pollutant-flow conditions also depend on the pollutant concentration, the flow patterns revealed situations of pollutant transport that cannot be inferred simply from the wind circulation modes. Different pollutants under the same wind patterns produced different pollutant-flow patterns. For example, the SO2 flow patterns (SO2-P3, SO2-P4, and SO2-P5) display large flow intensities at the NW quadrant bringing this pollutant towards the Mexico City region, mainly during winter and autumn seasons; however, the corresponding wind patterns exhibit wind velocities with similar magnitudes in almost all the quadrants. The events of the flow patterns SO2-P8 occurred with similar conditions, but during the spring and summer periods. Similarly, the PM10-P4 and PM10-P5, and PM10-P6 and PM10-P7 flow patterns revealed large flow intensities in the NE quadrant, although winds have similar magnitudes in other quadrants.

The pollutant-flow patterns of NO2, O3, PM10, and SO2 represent the main scenarios of the pollutant transport in the Mexico City region. They provide information that can guide, enrich, and enlighten the results of Mexico City air quality studies based on modeling techniques.

**Author Contributions:** Conceptualization, A.S.; Formal analysis, A.S.; Investigation, S.C.-S. and A.-T.C.-M.; Methodology, A.S., S.C.-S. and A.-T.C.-M.; Software, A.S. and S.C.-S.; Supervision, A.S.; Writing—original draft, A.S., S.C.-S. and A.-T.C.-M.; Writing—review & editing, A.S. and A.-T.C.-M.

**Funding:** This research received no external funding.

**Acknowledgments:** We thank Ana Laura Colín Aguilar (Instituto Nacional de Electricidad y Energías Limpias) for beneficial comments to make the paper more comprehensible.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
