Urban-Induced Changes on Local Circulation in Complex Terrain: Central Mexico Basin

: Land use land cover (LULC) signiﬁcantly impacts local circulation in the Mexico Basin, particularly wind ﬁeld circulations such as gap winds, convergence lines, and thermally induced upslope/downslope wind. A case study with a high-pressure system over the Mexico Basin isolates the inﬂuence of large-scale synoptic forcing. Numerical simulations reveal a wind system composed of meridional circulation and a zonal component. Thermal pressure gradients between the Mexico basin and its colder surroundings create near-surface convergence lines as part of the meridional circulation. Experiments show that the intensity and organization of meridional circulations and downslope winds increase when LULC changes from natural and cultivated land to urban. Zonal circulation exhibits a typical circulation pattern with the upslope ﬂow and descending motion in the middle of the basin. Large values of moist static energy are near the surface where air parcels pick up energy from the surface either as ﬂuxes of enthalpy or latent heat. Surface heat ﬂuxes and stored energy in the ground are drivers of local circulation, which is more evident in zonal circulation patterns.


Introduction
The interaction of local orographic circulations with urbanization effects is the focus of many past and recent observational, theoretical, and modeling studies ( [1][2][3][4]). Urban and builtup areas continue to proliferate around the world, with many cities and industrial complexes built within complex terrain and valleys [5]. Compounding this problem, a combination of synoptic conditions, climate change, and deteriorated air quality conditions aggravate citizen's health conditions in large urban centers [6]. The conversion of extensive land use areas, from natural to urban, implies significant changes in the thermodynamic, moisture, planetary boundary layer structure (PBL), and surface energy budget fields [7]. By replacing a natural surface cover with an impervious surface, the overlying air temperature of the urban area increases relative to its natural and rural surroundings due to larger sensible heat fluxes (SHF). This phenomenon, known as Urban Heat Island (UHI), is more intense when light winds and clear skies are present. In Bornstein ([8,9]), the author was a pioneer in the study of this phenomenon in New York City. Oke [10] presented an extensive review of observational studies of the characteristics, causes, and effects of the UHI. Significant advances have been made in recent decades on the subject, particularly on the relation between synoptic conditions and the UHI development. For example in Central Europe [11], have looked at how different anticyclonic patterns influence the magnitude of the UHI while [12], and [13] find that even under disturbed anticyclonic conditions with significant cloudiness and wind the UHI can still be detected in the diurnal cycle.
Jauregui [14] and Jauregui [15] first documented the diurnal evolution of the UHI in the Mexico Basin, which is home to one of the largest megacities globally, with more than 22 million inhabitants. The consequent urban growth has promoted an extraordinary change in land cover (LULC) in the basin, especially on the piedmont and slopes of mountains generating changes in atmospheric thermodynamics whose immediate effect is the generation of very intense UHI effects. More recently, with surface data from urban and rural stations Cui and De Foy [16] finds that the most intense UHI occurs in the dry period (winter and part of spring), reaching temperature differences between urban and rural areas of about 5 • C. The UHI also gives rise to vertical motion over urban areas, which in turn produces low-level wind convergence ( [17][18][19]). This local circulation interacts with a diurnal wind system consisting of upslope and downslope winds at the Mexico Basin's mountain slopes. Solar radiation produces diurnal temperature and density changes over the mountain slopes in the morning (upslope) and radiative cooling (downslope) later at night and early morning, respectively [20]. Klaus et al. [17] discuss how Mexico Basin's UHI intensity increases during the night and early morning by downslope winds from the surrounding mountains increasing low-level convergence at or near the urban area. More recently, Ganbat et al. [2] perform two-dimensional idealized numerical experiments with a mountain and an urban area to study the interactions between the mountain upslope/downslope winds with those that the UHI generates. They find that the UHI circulation can counteract the morning upslope winds while enhancing the downslope winds later in the evening. Klaus et al. [17] find similar results from observations for the Mexico Basin during a dry period (January-February, 1995) with upper-level weak anticyclonic synoptic forcing. In addition to this orographic wind system, past studies for the Mexico Basin have uncovered and analyzed the existence of a low-level diurnal southerly jet southeast of the Mexico Basin, denoted hereafter as the Chalco jet. Usually, the Chalco jet climbs the southeastern pass from Morelos Basin (Chalco pass) and then meets a northeasterly wind amid the Mexico Basin during the day and sometimes reverses direction at night. The intensity and impact of the Chalco jet in the low-level circulation of the basin generate lines of strong wind convergence that move northeastward north of the Mexico Basin ( [19,[21][22][23]). The origin of the Chalco jet lies within the thermally induced pressure gradients from the colder air surrounding the warmer elevated Mexico Basin (see discussions by [20,23,24]). Previous studies have revealed some aspects of the diurnal three-dimensional wind system of the Mexico Basin. However, it is still an open question on how LULC impacts local circulation and near-surface wind field circulations such as gap winds, convergence lines, and thermally induced upslope/downslope wind at the foot of mountain chains. In this paper, we combine observations and WRF model to focus on the diurnal evolution of the UHI and the wind systems generated by topography under conditions of a high-pressure system over the region. Section 2 presents the materials and methods used to define optimal days where weak synoptic forcing is present with an enhanced UHI. Section 3 presents the diurnal cycle characteristics of the Mexico Basin, the control experiment validation, and the hypothesis testing experiments needed to isolate the effect of the UHI from the entire wind system. Finally, Section 4 contains a summary and conclusions.

Terrain and Observational Data
The Mexico Basin's important physiographic characteristics include its tropical location (20 • N), high elevation (2250 m amsl), and the orographic confinement on its west and south sides by an almost continuous mountain chain with tops reaching about 3000 m agl. On its southeast corner, a mountain pass extends southward along the eastern mountains (Popocatépetl and Iztaccíhuatl) and connects it with a deep southern valley (Morelos Basin) 1000 m below (Figure 1). To the north, the Mexico Basin is open to the Sierra Madre Oriental and finally to the Gulf of Mexico, where it draws some of its moisture from the easterlies during the dry (winter) and wet (summer) seasons [20]. We use hourly surface meteorological variables from a network of 24 meteorological stations belonging to the Atmospheric Monitoring Agency of Mexico City (Secretaría del Medio Ambiente, SEDEMA). Additionally, 12 automatic meteorological stations maintained by the Mexican National Weather Service (Servicio Meteorológico Nacional, SMN) complete the observational surface network (Figure 1). Concerning upper levels, data is available from a 915 Mhz Doppler radar wind profiler (RAWP) (see Lau et al. [25] for technical details) maintained by SEDEMA. Note that the RAWP is located near the northern border of Mexico City limits (green dot in Figure 1) in a downward slope between two mountain ranges where a northwesterly gap wind usually develops. The RAWP outputs horizontal wind components as function of height above the basin floor and has an approximate vertical resolution of 60 m with a temporal resolution of 6 min and a range of about 3 km above ground level. Thus the RAWP can sample air well within the atmospheric boundary layer of the Mexico Basin. In addition, a radiosonde denoted here as RAS is available twice daily at 0600 and 1800 LST is launched from the headquarters of the SMN (lat: 19.4 • long: 99.2 • W elevation: 2313 m agl) in the northwest part of Mexico City (Figure 1). The surface database is subjected to a data quality control including time series homogeneization of temperature data. Table S1 (Supplementary Material) summarizes the surface station data. The air quality network from SEDEMA provides hourly ozone data. In addition to the above data, we use ECMWF ERA 5 reanalysis to analyze prevailing synoptic conditions and to obtain initial and boundary conditions for the regional atmospheric model as described below.

Case Study and Synoptic-Scale Conditions
The research focuses on the beginning of the warm period, which corresponds to February when the UHI intensity is higher in the Mexico Basin [16]. Due to availability of wind profiler data, the year 2017 is selected. To find days in February with weak synoptic conditions, we search days with the following restrictions: (1) ozone concentration values larger than 100 ppm (parts per million) as extracted from the SEDEMA air quality network, (2) zero precipitation from SMN rain gauge stations, (3) clear skies from solar radiation measurements at the SMN stations, and (4) adiabatic lapse rate with vertical gradients of potential temperature θ such that dθ/dZ > 9.8 K km −1 from the morning radiosonde at 0600 LST (see [20]). Vertical profiles of temperature, potential temperature, and adiabatic lapse rate for 9 February to 12 February are shown in Figure S1 of the Supplementary Material. All days show strong inversions at different levels in the vertical, but only 10 February stands out from the rest with the largest dθ/dZ = 10.2 K km −1 at about 690 m agl. On 10 February, average ozone concentrations show 125 ppm, and average direct solar radiation of 785 W m −2 for most SMN stations, which is near 90% of the theoretical average maximum value for the Mexico Basin on this date. Thus, 10 February meets all the above conditions best. Next, we analyze synoptic conditions for 10 February at 700 hPa and 500 hPa levels drawn from the ERA5 reanalysis data (see Figure S2 of Supplementary Material). At around 20 • N, there is a well-defined ridge pattern precisely centered over the Mexico Basin at this time. This circulation pattern results from a mid-latitude Rossby wave train traveling westward with two deep troughs located in the Pacific and the Atlantic and a ridge in the middle over the central U.S. and Mexico. At 500 hPa, the geopotential height and wind patterns are very similar, with only minor changes in direction from its 700 hPa counterpart. While February climatology shows ridging over western North America and Mexico, the circulation pattern at hand, has its center shifted to the east relative to climatology (see https://psl.noaa.gov/cgi-bin/data/testdap/plot.comp.pl accessed on 22 May 2021) for an ERA5 rendering of geopotential height maps). One critical characteristic of this large-scale circulation pattern is the strong inversion over the Mexico Basin, which keeps air motions and humidity well confined within the atmospheric boundary layer.

WRF Modeling System
In this study, we use version 4.0 of the mesoscale atmospheric model system, Weather Research and Forecasting (WRF, [26]) to perform realistic and idealized simulations over the Mexico Basin. The one-way nested grid configuration consists of 3 domains (see Figure  S3 of Supplementary Material) using a mesh spacing of 18 km in the external domain (D01), 6 km in the second domain (D02), and 2 km in the inner most domain (D03). The latter consists of 63 cells in the west-east direction and 60 cells in the north-south direction. This domain is designed to capture the orographically forced local phenomena and to avoid problems associated with the proximity to mountain ranges. All simulations use 76 sigma vertical levels, distributed as follows: 20 levels from the valley floor to 2250 m with about 50 m resolution for the first 5 levels, 30 levels between 2250 m and 6000 m, and the remaining 26 between 6000 m and 16,000 m. The choice of large vertical resolution is guided by the need to resolve within the atmospheric boundary layer the vertical structure of the wind system at the mountain slopes and the valley floor. The initial and lower boundary conditions of the numerical experiments are obtained from the ERA5 fifth generation of ECMWF atmospheric global reanalyses data which is available with a resolution of 0.25 • × 0.25 • starting from 1979, with a time resolution of every hour. The ERA5 is obtained from the Research Data Archive, which is maintained by the Computational and Information Systems Laboratory at the National Center for Atmospheric Research ( [27]).

Model Configuration
The model is configured with the following physical parameterizations: The rapid radiative transfer model for global circulation models [28] is used to parameterize the longwave and shortwave radiation. The Kain-Fritsch scheme is selected for the representa-tion of deep and shallow convection in the outer domain (the inner two domains have it turned off). The WSM6 scheme for cloud microphysics [29] is turned on in domain D03 only in order to account for possible orographic clouds. We use the Yonsei University (YSU) Planetary Boundary Layer (PBL) model ( [30,31]), and the revised Noah land surface model LSM [32] coupled to a modified version of the Monin-Obukhov Similarity model for surface turbulent fluxes (see Section 2.6 for further discussion). The experiments start on 8 February at 0000 LST and end on 14 February at 2300 LST, with model output every 10 min. The first 24 h are not considered in the analysis since they are part of the required model spin-up. The WRF model provides two LULC alternative databases: the U.S. Geological Survey (USGS), LULC database, generated in 1993, and the Moderate Resolution Imaging Spectroradiometer (MODIS) database, generated in 2001. Both datasets are not updated and are not accurate for all regions and some land cover types, such as urban areas ( [33][34][35][36]). In atmospheric modeling, an accurate representation of LULC is key because such information impacts the correct estimate of water and energy fluxes, the atmospheric circulation, and consequently, the quality of regional climate and weather simulations. The National Institute of Statistic and Geography (Instituto Nacional de Estadística y Geografía, INEGI) is the federal agency that gathers and distributes Mexico's LULC information which can be downloaded from its web page (https://www.inegi.org.mx/temas/usosuelo accessed on 27 January 2021). Serie-VI represents the latest LULC maps produced for Mexico with a far more reliable LULC description than those currently available for the WRF model. To obtain a proper LULC for WRF, we reclassify the Serie-VI LULC into the 24 USGS defined categories while including the particular phenological and geographical characteristics of the INEGI vegetation types. Lastly, we match Serie-VI's coordinate system and spatial resolution to the original USGS's LULC (see [37] for technical details).  Table 1 shows all the LULC categories and their percent area cover.  Figure 2 shows, the geographical distribution of LULC in the hypothesis testing experiments described above, including the CTRL. The CTRL and the three hypothesis testing help define two meaningful UE: the difference (CTRL-NATX) is the UE keeping the LULC heterogenous, and the difference (URBX-VEGX) is the UE keeping the LULC homogeneous (since the DCP category is in both URBX and VEGX). The hypothesis behind these two UEs is that the difference (CTRL-NATX) is more complex than its counterpart (URBX-VEGX) since inhomogeneous LULC surface conditions can generate mesoscale eddies within the PBL that can impact the local circulation effectively. To further investigate the UE on the three dimensional local circulation patterns in the Mexico Basin we use moist static energy (MSE) defined as: MSE = C p T + gz + L v q, where L v is the latent heat of vaporization (2.5 × 10 6 J kg −1 ), q the specific humidity, and C p T is enthalpy (J kg −1 ) and gz is potential energy per unit mass. MSE is conserved following non-condensing adiabatic motion of the parcels [38]. Synoptic and local circulations for this case study restrict the conversion of water vapor into condensed water(except possibly during ascend of moist air up the mountain slopes), and the only significant sources or sinks of water vapor q and MSE lie at the lower boundary and the lateral walls of the computational volume domain of the model.

Model Parameterization Tests
In order to obtain the best configuration of physical parameterizations available from WRF we use the following statistical measures of distance between simulations and observed fields: (i) the root mean square error (RMSE), (ii) the correlation coefficient (R) and (iii) the bias (B) ( [39,40] for details). The results of the comparisons between the different parameterizations are shown using Taylor diagrams which facilitate the interpretation of statistical metrics in a compact form, with the added advantage of analyzing several variables in a single diagram. The statistical metrics for simulated and observed quantities in the Taylor diagram are: where RMSE is root mean square error, R is the Pearson pattern correlation, B is the bias, MV i and OV i are the modeling and observation variables respectively closest to the i th observational site and N is the total number of observational sites considered. The µ and σ are the mean and standard deviations of the corresponding model and observed variable, respectively.
Several PBL turbulence schemes are tested in WRF to estimate the above statistical measures (not shown). Simulations with WRF underestimate wind magnitude in the Mexico Basin and it remains a challenge and an open question that is being pursued at the moment. For this reason, the choice is to minimize the error in wind direction rather than its magnitude. Table S3 in Supplementary Material shows the Taylor metric for temperature and the horizontal wind components. The revised YSU PBL model produces, on average, a Pearson correlation of about 0.5 in wind direction for the diurnal cycle. Other choices for PBL such as Mellor-Yamada 8 [41,42]), Mellor-Yamada-Niino ( [43,44]) and Total Energy-Mass Flux [45] show lower Pearson correlation values (not shown).  Figure 1). Each area temperature is the average temperature from several SMN weather stations within each area according to elevation while model data corresponds to the closest grid-point to the station location. Stations below 2300 m amsl correspond only to the Urban area. The model captures the temperature evolution of the Urban area for the most part of the day except for the late evening hours (1900 to 2300 LST) when a cold bias of about 3 • C develops. Over the East and North areas the model temperature falls well within the spatial variability early in the morning, and shows a persistent cold bias of about 2 • C, starting at 0700 LST which grows to 5°C, in the late evening hours. A warm bias of about 4 • C for the South area is present in the early morning hours that tends to diminish later in the afternoon. Despite the fact that SMN stations have very similar elevations for each area, it is interesting to note that there are temperature deviations among them early in the morning and late at night. This behavior is particular to the South, North and East areas and not so much for the Urban area. Figure 3 lower panel shows the UHI intensity as the temperature difference between the Urban area with the East, West, and South areas. In general, the model UHI intensity falls within the spatial variability bands of the observed data with the exception of the North area which has a much stronger positive UHI than observed in the evening hours. This result is consistent with the rather strong cold bias of the North area compared to the Urban area, as shown in Figure 3 upper panel. One possible explanation for the model excessive intensification of the positive UHI between the Urban and North areas is attributed to the orographic blocking of cold low-level northerly wind by a mountain chain in the middle of the Mexico Basin (Sierra Guadalupe, see Figure 1). Another important feature to notice of the UHI is the onset and end of cold pools between the Urban area, and the South and East areas. The cold pool between Urban and South starts at 1100 LST and further along, north of the Chalco pass, the cold pool between Urban and East starts at about 1600 LST. In the late afternoon, after 1800 LST the UHI between Urban and East areas increases signaling the arrival of cold air from the south at the South area near the Chalco pass. This phenomenon has been linked in the past to the advection of cold air by low-level south-westerly jet from the surrounding basins towards the warm Mexico Basin as discussed at some length by previous authors ( [18,20,23]), and will be further discussed in the context of the Mexico Basin wind system. Figure 3 lower panel also shows the distinct positive UHI to occur from about 1800 LST of the previous day until around noon for all areas.   Figure 4 also shows the vertical profile for the observed and model specific humidity. At 0600 LST, both, SMN and model, show a shallow humid (relative to values above 2871 m agl) and cold layer of about 6 g kg −1 and 7 • C on average. At this hour, the model and the SMN specific humidity values are very close to each other near the surface, but the model shows a more gradual drop in specific humidity compared to the SMN sounding. At 1800 LST, SMN specific humidity values range from about 3.7 g kg −1 near the surface, to 2.1 g kg −1 below and near an inversion level height of about 2870 m agl. Figure 4 also shows the model to be somewhat more humid than observed by about 1.0 g kg −1 . The simulation of the horizontal wind components is well captured in amplitude and pattern by the model, with westerlies above the 400 hPa level and easterlies below that level at 0600 and 1800 LST. Notably, at 0600 LST, below the 700 hPa level, in both model and observations a northwesterly wind flows down the slope towards the city center. Neighboring surface SMN stations agree in general with this early morning low-level wind field and report a shift to northeasterlies between midday and later in the evening (not shown). In blue the corresponding CTRL simulation variables. The right-most column panels shows the vertical profile of specific humidity (g kg −1 ). Figure 5 shows selected times (1000, 1400, 1900 and 2100, LST) on 10 February 2017, for the CTRL simulation wind field at 10 m agl (black vectors) and the corresponding observed wind field at 10 m agl from 28 available weather stations (red vectors) over the Mexico Basin. Figure 6 also shows narrow bands (or lines) of surface wind convergence from model data. At 1000 LST, Figure 5a shows a low-level northerly Chalco jet in the southeastern corner of the Mexico Basin, flowing towards the Morelos Basin in its southern border. North of the basin, there are northeasterly winds towards the southern mountain chain in the model, and in the observations as well. In addition, the model shows distinct upslope winds at both flanks of the mountain chains forming convergence lines over the mountain crests. At 1400 LST the Chalco jet changes to northerly direction in the model (Figure 5b), and although there is no observational site to confirm this event directly, the reversed direction is consistent with the observed increment in the UHI intensity after 1400 LST for the East area, as discussed earlier in the previous section. The model and observed wind vector data agree in direction quite well at this time, although the model at times underestimates quite significantly the magnitude of the wind (see Taylor diagram Figure S4 in Supplementary Material). In order to emphasize wind direction over wind magnitude in the horizontal maps, the norm of all vectors (observed and modeled) is set to one unitless magnitude. The wind system shows a dominant northerly wind component and upslope winds in all mountain chains surrounding the basin at this time as well. At 1900 and 2100 LST a thermally induced pressure gradient (from a warmer elevated basin and colder surroundings) generates northeasterly winds north of the basin (Figure 5c,d). At these times, radiational cooling from the southern mountain slopes induces downslope winds which join the southerly wind from Chalco and push through the entire basin creating strong south-west, north-east oriented convergence lines in the middle of the basin at the northern city limits. Note that, previous to the arrival of convergence lines in the middle of the basin (Figure 5b), ozone level concentrations are larger than 50 ppm, with some stations reporting values of about 100 ppm (not shown), and only start to decline to 40 and 30 ppm at around 1800 LST when convergence lines finally arrive there. The observed dilution of ozone concentrations is associated not only to turbulent mixing but also to vertical transport induced by horizontal convergence, as will be discussed in the next section.  Figure 6 shows the hourly horizontal wind components as a function of height for the February 10 case study by the RAWP in the north east of the basin (See Figure 1). Figure 6 also shows the CTRL simulation for the nearest grid point to the RAWP location. Note that the RAWP reports wind field data up to different height levels during the sampling period. The CTRL simulation at the RAWP site shows a well defined transition line characterized by a sharp veering from westerlies to easterlies whose height reaches a maximum of 2200 m agl at 1700 LST. At low levels, in the early morning hours, low-level northwesterly and westerly winds are present on both, the RAWP and the CTRL simulation, however the RAWP shows a more sustained westerly wind for almost the whole sampling period. Two near weather stations (LAA and 059, see Figure 1) report wind directions which agree more with the CTRL simulation. At upper levels, both the CTRL simulations and the RAWP agree on the easterly wind direction up until midday. After 1200 LST the CTRL and RAWP differ on the wind direction, which is interesting since the SMN radiosounding does not agree with the RAWP, as well. The reason for these discrepancies in wind direction is not clearly known, at this stage of the investigation. In general, however, we deem acceptable the performance of the model CTRL simulation when compared to available observations for this case study.  Figure 1 shows three transects, two in the south-north direction: SN1 (lon: 98.87 • W), and SN2 (lon: 99.06 • W) in the east and center of the Mexico Basin, respectively. Finally, we add a third west-east transect WE1 (lat: 19.43 • N) in the middle of the basin. Based on these three transects, Figure 7 shows MSE and wind field as a function of height and latitude for the CTRL simulation and the NATX experiment at 1000, and 1400 LST. At 1000 LST both, CTRL and NATX, display a similar behavior in wind and MSE fields. In all transects, there is air subsidence that keeps relative larger values (with respect to background) of MSE below 600 m agl. Upslope winds develop in all transects except in SN1 which shows a southerly low-level jet (the Chalco jet) advecting high values (relative to background) of MSE air into the Morelos Basin, south of the Mexico Basin. At the same time, however, there is an upslope wind from the Morelos Basin underneath converging with drainage flow at the crest of the Chalco pass. At 1400 LST the PBL reaches heights of about 1100 m over the floor's basin. Upslope winds continue in all transects except on SN1 where now the Chalco jet veers to southerly wind, and changes in MSE become more discernible between the CTRL and NATX. SN1, and SN2 show, on average, air parcels with low MSE below the PBL in the NATX experiment compared to CTRL due to more intense subsidence of low MSE parcels in the NATX experiment, particularly over the Chalco pass. At this time, the WE1 transect shows a circulation pattern that reminds of idealized two-dimensional simulations [46] with sinking motions in the middle of the basin and upslope wind at the mountain chains facing the basin. The CTRL experiment shows a larger PBL height compared to NATX and less intense subsidence of parcels with low MSE values on the eastern flank. Nevertheless, there is a remarkable asymmetry in the circulation that results from larger heating on the western mountain chain compared to that in the east [4]. Furthermore, at WE1, values of MSE for the CTRL are smaller than those of NATX, particularly near the surface. This decrease in MSE is related to larger latent heat fluxes in the NATX experiment as is discussed further below. Figure 8 shows MSE and wind fields at 1800, and 2100 LST, when the UHI becomes fully developed in the CTRL simulation as discussed previously in Section 3.1.1. At 1800 LST, near-surface convergence lines develop in the middle of the basin with a southeast-northwest slant due to air masses converging from the north and the south. Transect SN2, shows a convergence line at about 19.55 • N while SN1 shows two convergence lines separated by a region of low values of MSE air at about 19.70 • N (see also Figure 5). SN1 shows the Chalco jet carrying air with low values of MSE (325 kJ kg −1 ) which lowers the PBL height to the north. From the north, air masses with higher values of MSE (335 kJ kg −1 ) have higher humidity content than the southern air masses (not shown). At this time, the WE1 transect shows a zonal circulation that has evolved from earlier hours into a full convection cell transporting air parcels downward with lower values of MSE from aloft, on the eastern side of the basin, and upslope winds on its western side. The CTRL simulation is characterized by a stronger zonal circulation with larger PBL heights than the NATX experiment. At 2100 LST, transects SN1, and SN2 show the northern air mass reaching the southern basin carrying high values of MSE particularly below the PBL which is quite shallow with PBL heights of about 200 m agl for the CTRL and about 100 m agl for NATX. Of note is the more confined nature of air parcels within the basin in the CTRL simulation compared to those in the NATX experiment. This is associated to a stronger downslope wind carrying down air parcels with low MSE values until converging at the piedmont with higher MSE air impeding further advance towards the south. On the other hand, NATX shows air parcels upslope. At this time, the WE1 transect for the CTRL simulation shows convergence in the middle of the basin over the urban built-up area with ascending motion, as is clear from the MSE and wind fields. The WE1 transect for NATX fails to develop a similar wind and MSE patterns since, in this case, air parcels are more vertically confined. At this point in time (in contrast to 1400 LST, see Figure 7), MSE values are lower in the WE1 transect for the NATX experiment compared to the CTRL simulation with larger differences closer to the surface. This change in MSE is mostly related to higher MSE values from larger enthalpy near the ground in the CTRL compared to NATX. Figures S5 and S6 in Supplementary Material shows the corresponding transects for the URBX and VEGX experiments. These are all quite similar in pattern to their CTRL, and NATX counterpars. In magnitude, however, differences in MSE between the URBX and VEGX are more discernible and less noisy, and not so much, as far as wind field is concerned.   Figure 9 shows, the UE as differences between the CTRL and NATX experiments in latent heat flux (LHF), sensible heat flux (SHF), ground flux (GRF), and net radiation at the surface (RNET). We interpret (following [47] p. 251) positive fluxes as going into the atmosphere and negative otherwise. In the morning-afternoon (1000, and 1400 LST) the RNET difference is negative implying, that RNET for NATX is larger (i.e., less negative) than RNET for the CTRL. At this time, SHF is positive for both CTRL and NATX, and very close to each other in magnitude at 1000 LST. However, at 1400 LST, the SHF of CTRL surpasses that of NATX by about 40 W m −2 as Figure 9 shows in the SHF panels. At this time, the LHF in the CTRL simulation is very small in absolute value while that of the NATX experiment is positive (not shown) leading to negative differences in LHF between the CTRL and NATX simulations. At 1400 LST this negative difference in LHF increases. This increase in LHF is in agreement with larger values of MSE near the surface as Figure 7 shows for NATX at 1400 LST. Next, we examine the night period (1800, and 2100 LST). The positive difference in SHF implies the CTRL has a larger SHF than NATX, as expected, during this part of the diurnal cycle with values larger than 45 W m −2 with maxima at 1800 LST as Figure 9 shows in the SHF column. The response in MSE can be seen in all transects at 2100 LST with larger values of MSE for the CTRL near the surface compared to NATX (see Figure 7). Finally, the UE in this very simple model of urban area is interpreted as the difference in energy storage by the urban and built-up category in the CTRL and the NATX experiment. In the morning-afternoon period, the GRF is negative, and Figure 9 shows the CTRL to be more negative than NATX as expected since the urban and built-up area is expected to store more energy than the DCP category during the day. At night, the GRF becomes positive for both CTRL and NATX experiment but this time, the difference in GRF is positive meaning that the urban and built-up area is giving off energy to the atmosphere at a greater rate than the DCP category. Figure S7 in Supplementary Material shows the homogeneous UE for the URBX and VEGX experiments.

Summary and Conclusions
The present work tests the hypothesis that changes in LULC significantly impact local circulation in the Mexico Basin, particularly near-surface wind field circulations such as gap winds, convergence lines, and thermally induced upslope/downslope wind at the foot of mountain chains. Many past modeling studies on urban heat islands support this hypothesis (e.g., [2,46]). Our results confirm this, albeit with some caveats due to the complex topography and synoptic forcing that characterizes the Mexico Basin. In order to isolate the influence of large-scale synoptic forcing, we select 10 February 2017, as a special case study with a very well-defined high-pressure system over the Mexico Basin at the center of the country. This case study selection permits a convincing simulation of the main physical processes at work. It isolates the influence that LULC has on the wind system and the thermodynamics based on hypothesis testing experiments where LULC changes from the control simulation. Thus, we interpret the urban effect regarding the differences between the control and the hypothesis testing experiments. The control simulation reveals a three dimensional wind pattern composed of: (i) meridional (south-north) circulations driven by the thermal pressure gradients between the high plateau of the Mexico Basin and surrounding coastal areas (see Figure S5 in Supplementary Material), and (ii) a zonal (west-east) upslope/downslope circulation driven by the thermal interaction with the urban and built-up area and the heating/cooling of the mountain slopes facing the basin. The sense of the meridional circulation depends on the change of pressure gradients between the Mexico Basin and its surroundings as the basin warms and cools faster than its surroundings. The replacement of the natural and cultivated land by urban and built-up areas leads to more organized meridional and zonal circulations inside the basin, particularly in the afternoon and night times when the UHI intensity is largest. Nearsurface wind convergence lines in the afternoon connect with the upper air circulation and not only help close the meridional circulation but mark the border between two distinct air masses: one with low MSE advected from the south and another with high MSE advected from the north. At this time, downslope wind from the southern mountains and Chalco jet complete the low-level wind meridional circulation branch. This type of circulation helps redistribute the MSE or passive tracers originating at the basin floor vertically and export them to the southern basin (Morelos Basin) or returning them to the basin. Results from hypothesis testing experiments with natural and cultivated areas indicate a similar meridional circulation albeit less organized and with lower PBL heights. In the afternoon, the zonal circulation exhibits a typical circulation pattern with upslope flow on both of the flanking mountains, and descending motion over the Chalco pass. Also there is a zonal asymmetry with lower MSE over the Chalco pass and larger MSE on the western side of the basin. The experiment with natural cover only, shows a more accentuated asymmetry than the control. The controls that the land use land cover exert upon the circulation result in larger values of MSE near the surface where air parcels pick up energy from the surface either in the form of fluxes of enthalpy or latent heat. Surface heat fluxes and stored energy in the ground are drivers of local circulation which is more evident in zonal circulation patterns. One surprising result from this work is the generation of convection cells at the basin scale in the afternoon. One would expect that given the low ratio of mountain height to basin horizontal length (0.03 on average), the UHI circulations would be disconnected from the upslope and downslope winds. Rendón et al. [48] in an idealized two-dimensional experiment with a 0.06 ratio find that the UHI circulations are closely linked to the upslope/downslope wind system. It would be interesting to explore in, the more complex topography context of the Mexico Basin, if there is a lower limit to this ratio where UHI circulations and upslope/downslope winds are decoupled from each other. Finally, this work opens up another important avenue of research in regards to the thermal pressure gradients as the main driver of near-surface wind divergence/convergence, and circulation within the Mexico Basin. This dynamical ingredient is required if we are to understand how large-scale variability during the dry season affects the mesoescale circulations within the basin. These type of research is of high relevance to air quality studies for the basin. The present work is a first in a series of research works dedicated to understand in detail the interactions of the Mexico Basin circulation patterns and its urban landscape.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/atmos12070904/s1, Figure S1: Vertical profiles, Figure S2: Synoptic charts, Figure S3: Computational domains, Figure S4: Day-time MSE vertical cross-sections for URBX and VEGX, Figure S5: Night-time MSE vertical cross-sections for URBX and VEGX. Figure S6: Same as Figure S5, except for 1800 and 2100 LST, Figure S7: Urban effect as the difference in energy fluxes W m −2 at the surface between the URBX simulation and the VEGX experiment for February 10 at 1000, 1400, 1800, and 2100 LST. Columns 1, 2, 3 and 4 are for latent heat fluxes (LHF), sensible heat fluxes (SHF), ground flux (GRF), and net radiation, respectively. Positive direction for fluxes implies atmosphere receives energy from the surface below. Table S1: Stations locations, Table S2: Station Taylor metrics, Table S3: Hourly Taylor metrics. Acknowledgments: We thank two anonymous reviewers, for valuable comments that significantly improved the manuscript. We thank the SEDEMA and SMN for providing the data. All figures were generated with The NCAR Command Language (Version 6.6.2) [Software]. (2019). Boulder, CO, USA: UCAR/NCAR/CISL/TDD, doi:10.5065/D6WD3XH5.

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

Abbreviations
The following abbreviations are used in this manuscript: