1. Introduction

Extreme climate events (ECE) must be periodically monitored and analyzed in detail, due to their role as high impact agents on society, environment, and ecosystems. An extreme climatic or meteorological event refers to the occurrence of a climatic or meteorological variable which value is above or below a threshold that is close to the upper (or lower) limits of the range of observed values of the variable (^{Seneviratne et al., 2012}). ECE are important because of their impacts, but they are difficult to quantify statistically, as they are infrequent and occur at multiple scales (^{Palmer and Räisänen, 2002}). The definitions of “rare” vary, but an ECE would normally be as rare as, or rarer than, the 10th or 90th percentile of the observed probability density function.

It is likely that different ECE affect specific regions and increases in their frequency and intensity have been detected in several regions of the world (^{Brown et al., 2008}; ^{Almazroui et al., 2014}; ^{Chen et al., 2015}; ^{Wypych et al., 2017}; ^{Caloiero, 2017}). It is expected that these ECE will intensify in the future in response to global climate changes caused by the emission of greenhouse gases (^{Beniston et al., 2007}; ^{Gao et al., 2012}; ^{Lau and Nath, 2012}; ^{Easterling et al., 2016}; ^{Grotjahn et al., 2016}; ^{Schoof and Robeson, 2016}).

There are basically two fundamental approaches to study ECE: global circulation models (GCMs), and statistical models using the extreme value theory (EVT). Contemporary GCMs, such as those used for the 5th Coupled Model Intercomparison Project (^{Taylor et al., 2012}) are a key component of regional climate change projections, but their limited spatial resolution reduces their utility in estimating local or regional extremes without substantial post-processing (^{Schoof and Robeson, 2016}). On the other hand, the central issue of EVT is the modeling of extreme events, and the main purpose of this theory is to provide asymptotic models for the distribution tails (^{Furió and Meneu, 2011}). Therefore, EVT aims at deriving a probability distribution of events at the far end of the upper or lower ranges of the probability distributions (^{Coles, 2001}); its main advantage is that it allows for estimating and analyzing the probability of occurrence of events that are outside of the observed data range (^{Raggad, 2018}).

For these reasons, EVT is the approach that has been chosen in this research due to its wide applicability in different fields that are related to extreme weather and climate events and their impact: ecology (^{Moritz, 1997}; ^{Meehl et al., 2000}; ^{Dixon et al., 2005}; ^{Katz et al., 2005}; ^{Jentsch et al., 2007}; ^{Burgman et al., 2012}); extreme temperatures and heat waves (^{Meehl and Tebaldi, 2004}; ^{Della-Marta et al., 2007}; ^{Parey et al., 2007}; ^{García-Cueto et al., 2010}; ^{Waylen et al., 2012}; ^{Tanarhte et al., 2015}; ^{Liu et al., 2015}; ^{Shen et al., 2016}); extreme rainfall (^{Katz et al., 2002}; ^{Koutsoyiannis, 2004}; ^{Friederichs, 2010}; ^{Papalexiou and Koutsoyiannis, 2013}; ^{Kim et al., 2015}; ^{Boucefiane and Meddi, 2019}); and damages to the communities that affect agroecosystems through changes in soil moisture and evapotranspiration rates (^{Miralles et al., 2014}; ^{Whan et al., 2015}; ^{Guan et al., 2015}; ^{Hatfield and Prueger, 2015}).

Changes in extremes of temperature and precipitation have been evaluated in different regions of the world. However, until the Fourth Assessment Report of the Intergovernmental Panel on Climate Change (^{Trenberth et al., 2007}), the cities had been treated as “noise-generating” entities in globally studied climatic signals. In the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (^{IPCC, 2013}), a special chapter was dedicated to the subject of cities, and their role in climate change. This is not surprising, as although cities are very important contributors to social and economic well-being, they require an uninterrupted source of energy for all their activities. Cities consume approximately 75% of global primary energy and emit between 50-60% of the greenhouse gases (GHG) on the planet (^{Rosenzweig et al., 2011}). This figure can be raised to 80% when indirect emissions generated by inhabitants of the cities are included (^{UN-Habitat, 2011}; ^{Kraussmann et al., 2017}). Thus, cities promote global warming, and contribute to an increase in average surface temperature at the planetary level.

In addition to the global effects of GHG, the worldwide upward trend of urban growth must be considered, both in population and in its areal extension. This growth has generated environmental problems, not only air pollution and solid waste, but also those concerning a different environmental product, such as the genesis of an urban climate. In particular, the main connotation of urban climate is the formation of an urban heat island, which in turn requires additional water and energy to maintain thermal comfort through air-conditioned spaces (^{Coutts et al., 2012}; ^{Wang et al., 2016}; ^{de Munck et al., 2018}, ^{Skelhorn et al., 2018}). Thus, cities that are already particularly vulnerable to ECE caused by global climate change must now also consider the effects caused by local climate change. Therefore, it can be inferred that because a city is the geographical space that brings together a majority of the population and provision of services, it will be where the greatest vulnerabilities associated with the impacts of climate change manifest themselves (^{UN-Habitat, 2011}).

In Mexico approximately three out of four people (72.3%) live in cities according to Fundación Centro de Investigación y Documentación de la Casa and Sociedad Hipotecaria Federal (^{CIDOC-SHF, 2011}). This percentage is expected to increase in the medium term, as according to projections of the National Population Council the number of people in 384 localities of the National Urban System will increase by 16.6 million (from 82.6 million in 2010 to 99.3 million in 2030) as a result of an annual average growth rate of 0.92% (^{Hernández et al., 2014}). The urban proportion of the national population will increase to 77.9% (18.1 million new urban inhabitants). This trend of the geographic dynamics of cities is inequitable with low levels of quality of life and urban sustainability, and not all cities have the same development potential. Thus, the challenges in facing climatic risks such as heat waves and floods will be massive if quantitative scientific studies are not carried out.

Very little research has been conducted in Mexico on extreme climate values in urban environments (^{Magaña et al. 2003}, ²⁰¹²; ^{Cavazos and Rivas, 2004}; ^{Ríos-Alejandro, 2011}; ^{García-Cueto et al., 2013}, ²⁰¹⁴, ²⁰¹⁸; ^{Martínez-Austria and Bandala, 2017}). This can be explained by a limited access to measured climate data, limitations in the geographical coverage of the networks of stations, and interruptions in climate series due to missing data. In view of this and given the quantitative uncertainty of the climate extremes mentioned at the urban level and their great importance for the assessment of risks and adaptation proposals, this study selects some cities in Mexico with important increases in population and in areal extension, which have recently been affectated by ECE.

Thus, this study has three main objectives: (a) to detect thermal behavior in some growing cities of Mexico, (b) to model extreme temperatures through the EVT, and (c) to make projections of return levels for extreme temperatures in a future changing climate.

This paper is organized as follows: section 2 presents the climatology of extreme values in Mexico; section 3 presents the study area and the climate data, and section 4 provides a theoretical outline of the methodology. The results and discussion are presented in section 5. In section 6, conclusions are drawn.

2. Climatology of extreme temperature values of Mexico

Figures 1 and 2 show the 10th and 90th percentiles of minimum and maximum temperatures, respectively (^{Cavazos et al., 2013}), for the period 1961-2000; the methodology used in these figures is described in ^{Colorado-Ruiz et al. (2018)}. According to extreme climate indices derived from reliability ensemble averaging (REA) (Fig. 1), the coldest winters and their 10th percentile (P10) during the winter months of December, January, and February (DJF) for Mexico are characterized by minimum temperatures below 0 ºC in the highlands of the Sierra Madre Occidental and the Mexican High Plateau, and between 0 and 5 ºC in much of northern Mexico. Figure 2 shows that the extreme values of the 90th percentile (P90) of the maximum temperatures during the summer months of June, July, and August (JJA) expand the area covered by the isotherms of 25 and 40 ºC to a narrower strip in the border region. This region of semiarid climate is the most extreme in Mexico; therefore, it is the most susceptible to negative impacts caused by increases in temperature.

Fig. 1 Thresholds of P10 of minimum temperature for winter, obtained with the ensemble of the Reliability Ensemble Averaging (REA) for 1961-2000 (^{Cavazos et al., 2013}). 

Fig. 2 Thresholds of P90 of maximum temperature for summer, obtained with the ensemble of the Reliability Ensemble Averaging (REA) for 1961-2000 (^{Cavazos et al., 2013}). 

The projected climate changes at a national level (^{Cavazos et al., 2013}) agree with those obtained globally (^{Alexander et al., 2006}; ^{Caesar et al., 2011}; ^{Song et al., 2014}). Regarding those related to the increase in the frequency of heat days, heatwaves and variability in precipitation will exacerbate problems related to natural and human systems (^{Ummenhofer and Meehl, 2017}). If we analyze the 10th and 90th percentiles of extreme temperatures at the national level, there is a notable lack of reliable and timely information at the regional and urban levels for risk assessment, which justifies the study carried out herein.

3. Data and study area

A digital historical collection of daily temperature data for several cities in Mexico was possible using information from climatological stations operated by the Servicio Meteorológico Nacional (SMN, National Weather Service). Unfortunately, the selection of urban areas of interest faced some limitations, as not all climate stations have the same record length, and none had a strict data quality control. Based on a previous analysis of the growth of some cities in Mexico and the occurrence of climatic events that have affected them in important ways, 21 cities were selected in the first instance. However, an analysis of the quality of climate information limited the study to only 12 cities. These cities are heterogeneously distributed in the country because, as mentioned above, they were selected for their population growth and urban development. Data quality control was carried out for the daily information for each city, consisting of the following synthetic process.

The SMN database from the Climate Computing Project (CLICOM) was extracted in the .csv format, using Matlab and RClimdex software. The climatic information was explored, and quality control was performed for the original database. The daily meteorological values of maximum temperature (TXX) and minimum temperature (TNN) were selected with a computation routine; at the same time, a continuity of the time series was sought by adding missing dates and values, if possible, and adding a -99.99 label for missing data. Subsequently, the database was imported into the RClimdex program (^{Zhang and Yang, 2004}), ensuring the internal and temporal consistency of the daily climatological information. This quality control validated the following: (1) internal coherence, by verifying that the maximum temperature was always greater than the minimum temperature; (2) identification of atypical values and changes in the seasonal cycle or variability of the data through visual inspection of the time series of TXX and TNN, and (3) identification of values located more than four standard deviations (σ) from the mean as outliers and possible errors. Finally, outliers were verified individually to determine if they had been caused by an atypical event or if the measurement was incorrect and had to be discarded.

Temporal homogeneity of the data was then evaluated using the RHtest V3 software (^{Wang and Feng, 2010}) to identify abrupt jumps or change-points. This homogeneity test is based on a two-phase regression model with a linear trend for the entire series, applied to selected series for each of the cities.

The final selection included the following cities and/or intra-urban regions for the detailed analysis of extreme temperatures: Aguascalientes, Mexicali, Tijuana, Tuxtla Gutiérrez, Mexico City (with analysis of the climatological stations of Úrsula Coapa, Gran Canal, and Milpa Alta), León, Guadalajara, Monterrey, Puebla, Tlaxcala, Veracruz and the metropolitan area of the Valley of Mexico, including the urban areas of Toluca, Aculco and Chapingo (Table I). Their locations are presented in Figure 3 and, as can be seen, the analysis includes cities distributed in central Mexico (18-25º N latitude), a city in southeastern Mexico (Tuxtla Gutiérrez), another in northeastern Mexico (Monterrey), and two more in northwestern Mexico (Tijuana and Mexicali).

Fig. 3 Locations of intra-urban cities or regions according to Table I. 

Our study focused on ECE using a parametric non-stationary generalized extreme value (GEV) distribution model, which tacitly assumes that extreme events are changing over time as climate changes. Instead of the climate change indices used in ^{García-Cueto et al. (2018)}, the current study directly applies extreme temperature values. In addition, the non-stationary GEV model is applied to investigate how the return levels of extreme temperatures might change in the future.

It is worth mentioning that for the city of Mexicali, due to the importance of extreme temperatures for human comfort and energy consumption for use in air conditioning, other studies have been carried out (^{García-Cueto et al., 2010}, ²⁰¹³). Unlike current research, in ^{García-Cueto et al. (2010)} warm days were modeled with the GEV and the maximum temperature was included as a covariate, without performing a temporal trend analysis of extreme temperatures. Their difference with respect to ^{García-Cueto et al. (2013)} is that the trends of extreme temperatures and their return periods are updated with a methodological process that requires a re-parameterization of the GEV model, which gives greater reliability in its estimation.

4. Methodology

In this section we present the techniques used for trend analysis and a brief review of the GEV distribution, to provide a basis for the modeling of extreme temperature events. As described in section 1, the GEV distribution is used to model extremes in atmospheric science and in many other scientific fields. Using the trends of annual temperature series, we describe the non-stationary models, the criteria for deciding which GEV model to use, the estimation of the parameters of the GEV distribution, and the return levels.

5. Results and discussion

In this section we present and analyze the results of the parametric approach based on the GEV distribution for modeling the maximum annual value of the TXX and TNN for each of the selected urban areas in Mexico. The data series are analyzed in each climatological station, the use of stationary and non-stationary models is evaluated, and a statistical evaluation of the changes in TXX and TNN is presented.

We modeled the data series through the GEV distribution of three parameters, using stationary and non-stationary models for periods that varied from 83 years (Veracruz) to 31 years (Tuxtla Gutiérrez); the other cities had intermediate periods (Table I). The inclusion of non-stationarity is plausible for our modeling approach, as can be visualized for many locations. As an example, time series of TXX (Figs. 4, 6, 8, and 10) and TNN (Figs. 5, 7, 9, and 11) are shown for four of the 12 urban areas (Tijuana, Guadalajara, Toluca, and Puebla), with a trend line using Sen’s slope estimator. The graphs for several cities show a clear trend in the maximum and minimum annual data. The graphical analysis and results of the Mann-Kendall trend test justify the use of non-stationary GEV models over time as a covariate. The trends shown are statistically significant with a p-value < 0.05, except for the TXX of Toluca, whose trend is significant with a p-value < 0.1, and the TXX of Tijuana, whose trend is not significant.

Fig. 4 Trend of extreme maximum temperature (TXX) in Tijuana from 1950 to 2012. 

Fig. 5 Trend of extreme minimum temperature (TNN) in Tijuana from 1950 to 2012. 

Fig. 6 Trend of extreme maximum temperature (TXX) in Guadalajara from 1957 to 2013. 

Fig. 7 Trend of extreme minimum temperature (TNN) in Guadalajara from 1957 to 2013. 

Fig. 8 Trend of extreme maximum temperature (TXX) in Toluca from 1974 to 2009. 

Fig. 9 Trend of extreme minimum temperature (TNN) in Toluca from 1974 to 2009. 

Fig. 10 Trend of extreme maximum temperature (TXX) in Puebla from 1955 to 2013. 

Fig. 11 Trend of extreme minimum temperature (TNN) in Puebla from 1955 to 2013. 

The time series of TXX and TNN were used to estimate the parameters in the distribution of GEV with and without trends, as well as the periods of return. The models were adjusted to the maximum annual temperature and minimum annual temperature of each of the 16 climate stations using the maximum likelihood (ML) method.

An attempt was made to improve the modeling approach by allowing the location parameter (μ) to depend on time. As mentioned, three model selection criteria (AIC, BIC, and LRT) were used to select the best model from a collection of nested models. The best adjustment models for TXX and TNN through the GEV distributions, as selected by the criteria separately listed for each urban area, are summarized in Table II.

In total, 58 models were generated. Table II shows the resulting estimators and the criteria for selecting the best model to be used in the estimation of return periods. According to this table, for the temporal trend of TXX, significant at a 95% confidence level, it can be seen that: (a) 44% of locations (seven urban areas) show a significant trend, (b) 50% of locations (eight urban areas) have no significant trend, and (c) 6% of locations (one urban area) show no trend. For the temporal trend of TNN, also significant at the 95% confidence level, it was found that: (a) 82% of locations (13 urban areas) show a significant trend, (b) 12% of locations (two urban areas) show no significant trend, and (c) 6% of locations (one urban area) show no trend.

According to the shape parameter (ξ), and by performing a detailed analysis, it was found that for the case of TXX, 63% of the climatic stations were adjusted to the Weibull distribution (ξ < 0), 25% were adjusted to the Gumbel distribution, 6% were adjusted to the Fréchet distribution, and 6% did not fit any of the three distributions. In the case of TNN, 94% of the analyzed climatic stations were adjusted to the Weibull distribution, and only 6% did not adjust to any of the distributions.

In the case of TXX, which exhibits a significant temporal trend, it was found that 72% was adjusted to the Weibull distribution, 14% to the Gumbel distribution, and 14% to the Fréchet distribution. For TNN, it was found that those that show a significant temporal trend (82%) and those that are stationary (12%) conformed to the Weibull distribution.

Overall, for both extreme temperatures TXX and TNN, it was found that the 57% following the Weibull distribution are statistically significant to the temporal trend at the 95% confidence level. One of the properties of this distribution, which may even be controversial because of the trend found, is that both extreme temperatures have an upper limit that cannot be exceeded.

Figures 12 and 13 describe the temporal trend pattern of decadal trends for the location parameter (μ) of the maximum and minimum temperatures according to non-stationary GEV models (Table II). About the TXX (Fig. 12), the most significant positive trend occurred in the center and east regions (Toluca, Chapingo, Puebla and Veracruz), whereas the lowest trends occurred in the northwest regions (Tijuana and Mexicali).

Fig. 12 Spatial patterns of trends of the location parameter (μ) for the maximum temperature. Red (blue) triangles mean positive (negative) values. Full triangles mean significant trends at the 5% level. The symbol ⊗ indicates without change. 

Fig. 13 Spatial patterns of trends of the location parameter (μ) for the minimum temperature. Red (blue) triangles mean positive (negative) values. Full triangles mean significant trends at the 5% level. 

Respect to the TNN (Fig. 13), the most significant positive trends were most relevant in different geographical areas of the country, i.e., higher in the west center (Guadalajara and León), southeast (Tuxtla Gutiérrez), some urban areas located in the central part (Úrsula Coapa, Toluca, Chapingo, Tlaxcala and Puebla), and in the northeast and northwest (Monterrey and Mexicali), whereas the lowest positive trend occurred in the urban area of Tijuana (with a value of 0.30). The most significant negative trend occurred in an urban area.

Once the best models were selected, the return levels of the extreme, maximum, and minimum temperatures were estimated. Tables III and IV present the estimators and confidence intervals for return levels to 10, 20, 50, and 100 years, for both stationary and non-stationary models. The stationary return levels for TXX in Table III, indicate that the values increase for increasingly larger periods (10, 20, 50, and 100 years), except for the urban area of Tuxtla Gutiérrez. In addition, the confidence intervals were increasingly wider as the periods of return increased. In Aguascalientes the TXX can be expected to progressively exceed 37.1 ºC on average every 10 years, 37.9 ºC on average every 20 years, 38.9 ºC on average every 50 years, and 39.7 ºC on average every 100 years. The 95% confidence intervals for these return periods, in ºC, were 36.4-37.8, 37.0-38.9, 37.5-40.4, and 37.8-41.5, respectively.

Among the locations considered in that block, and in the context of the stationary return levels, Monterrey (in the northeast of the country) was associated with the highest return levels, and Milpa Alta (in the metropolitan area of the Valley of Mexico in the center of the country) had the lowest return levels. Based on the 95% confidence interval and according to Table III, we can expect that the largest TXX event recorded for Milpa Alta could reappear for the next 10 years. Moreover, there is a high probability that the annual TXX will exceed the maximum historical value in the next 50 years, except in the sites Gran Canal and Monterrey.

The estimated return levels assume stationarity, meaning that the level of return for a return period is the same for the successive years. This implies that the statistical properties of the parameters μ, σ, and ξ are constant.

In a non-stationary case, the parameters of the GEV distribution vary over time, and the return levels of extreme temperatures will also follow that temporal trend. Under a changing climate, the return value can be interpreted as an extreme quantile of a temperature distribution that varies over time (for example, a return value of 20 years can be interpreted as a value that has a 5% probability to be exceeded in a particular year).

It is now possible to estimate return levels for any year, which are also presented in Table III. For example, note that the average return levels (in ºC) for Mexicali for 10, 20, 50, and 100 years, are 48.5, 49.4, 50.9 and 52.7, respectively. The confidence intervals at 95% (in º C) are 48.0-49.1, 48.8-50.1, 50.1-51.7, and 51.7-53.7, respectively. The differences in return levels between stations is remarkable, both in the stationary GEV models and in the non-stationary models. Mexicali, in the northwest of Mexico, is associated with the highest return levels, whereas Milpa Alta, a metropolitan area in Mexico City, has the lowest return levels.

With respect to TNN, Table IV shows that only stations Milpa Alta and Veracruz have stationary return levels, and that the values increase slightly for increasingly large return periods (10, 20, 50, and 100 years). In addition, the confidence intervals remain nearly constant as the periods of return increase. For example, Milpa Alta could expect the TNN to exceed 1.9 ºC on average every 10 years, 2.3 ºC on average every 20 years, 2.6 ºC on average every 50 years, and 2.8 ºC on average every 100 years. The confidence intervals, at 95% (in ºC) are 1.4-2.3, 1.8-2.7, 2.2-3.1, and 2.3-3.3, respectively.

As already mentioned, the estimated return levels assume stationarity, meaning that the level of return for a return period is the same for successive years. This also implies that the statistical properties, as mentioned for TXX, keep the parameters μ, σ, y, and ξ constant.

In a non-stationary case (as with TXX), the parameters vary in time and the return levels of extreme temperatures will follow a similar temporal trend. It is possible to estimate return levels for any year of interest in a time period. The return levels for the non-stationary TNN series are presented in Table IV. Note that the average return levels of 10, 20, 50, and 100 years (in ºC) for Aguascalientes are -1.3, -0.3, 2.1, and 5.9, respectively, whereas their confidence intervals, at 95% and in ºC, progress from -1.9 to -0.7, -1.0 to 0.4, 1.3 to 3.0, and 4.9 to 7.0, respectively.

In the case of a non-stationary series, for both TXX and TNN, a positive trend of the location parameter (μ ₁ ) will be reflected in the positive trend of the return levels (^{Chen and Chu, 2014}). This means that with greater magnitudes of the trend return levels will increase considerably, and threshold values of return levels will contain greater variation over time. As examples of estimators of return levels, the time series of return levels for TXX for the city of Puebla (Fig. 14) and for TNN for the city of Guadalajara (Fig. 15) are plotted, according to the GEV adjustment non-stationary model. The difference in the return level for 2 years for the city of Puebla from the beginning of the period (1955) to the end (2013) is approximately 5 ºC, whereas the difference in the return level of 2 years for the city of Guadalajara from the beginning of the period (1957) to the end (2013) is approximately 3 ºC.

Fig. 14 Time series of extreme maximum temperature (TXX) return levels in Puebla in accordance with the non-stationary generalized extreme value (GEV) model. The red, green, and blue lines represent the return levels of 2, 20, and 100 years, respectively. The black line represents observed values. 

Fig. 15 Time series of extreme minimum temperature (TNN) return levels in Guadalajara in accordance with the non-stationary GEV model. The red, green, and blue lines represent the return levels of 2, 20, and 100 years, respectively. The black line represents observed values. 

The application of the stationary and non-stationary GEV distributions has allowed this study to model the behavior of extreme temperatures in some cities of Mexico. Due to the lack of reliable data, only 16 climate stations were used, focusing on 12 urban areas. The latitudinal dispersion of these climatic stations ranged from 16º 45’ (Tuxtla Gutiérrez) to 32º 39’ (Mexicali), and with altitudes above sea level varying from 4 m (Mexicali) to 2726 m (Toluca). This resulted in several types of climates according to the geographical location (^{García-Cueto et al. 2018}), making the analysis by chosen city even more interesting. This is particularly true for those cities that show significant trends and periods of return that could put the adaptation of people and the urban ecosystem at risk, particularly the fauna and flora. We found that the TXX is increasing in Mexicali, Tijuana, Toluca, Chapingo, Puebla and Veracruz, and no significant changes were detected in the other 10 urban zones. A larger number of climatological stations indicate increasing TNN compared to TXX, namely Aguascalientes, Mexicali, Tijuana, Tuxtla Gutiérrez, Úrsula Coapa, León, Guadalajara, Aculco, Toluca, Chapingo, Monterrey, Puebla and Tlaxcala. Non significant increases were detected in only two stations, i.e., Milpa Alta and Veracruz.

A comparison of the obtained results with those from studies using the ETCDDI indices shows similarities and some differences. For example, ^{García-Cueto et al. (2018)} used the same database herein and found a statistically significant increasing trend in the TX90 (warm days) and TN90 (warm nights) indices for Aguascalientes, whereas for Milpa Alta, both indices (TX90 and TN90) decreased significantly. The difference is between the TXX (stationary in this study) and TX90 (increasing trend) for Aguascalientes; in the case of Milpa Alta, there are coincidences in both indices (TXX with TX90 and TNN with TN90) that show decreasing trends.

At the country level, according to ^{Gosling et al. (2011)}, Mexico has experienced a generalized warming since 1960. The frequency of cold days has decreased, and the frequency of warm nights has increased. There has also been a general increase in average winter temperatures in the country as a result of human influences on climate. Thus, during the winter season, the occurrences of warm temperatures are more frequent, and those of cold temperatures are less frequent. For the A1B emissions scenario, the projected temperature increases of the Coupled Model Intercomparison Project Phase 3 (CMIP3) over Mexico are approximately 4 ºC near the US border, with increases in the rest of the country between 2.5 and 3.5 ºC. This implies that the influence of GHGs will have an additive effect on the potential warming caused by urbanization. This can become a major problem affecting city dwellers both positively and negatively, depending on their latitudinal location.

Significant trends in extreme maximum and minimum temperatures found in many parts of the world (^{Heim Jr., 2015}; ^{Easterling et al., 2016}; ^{Houngninou et al., 2017}) reinforce the results obtained here. In the US, extremely hot maximum and minimum temperatures have shown increasing trends between the 20th and 21st centuries and during the last four decades (1971-2013), whereas the presence of extremely low and minimum temperatures has decreased during these periods (^{Hartmann et al., 2013}).

The application of a theory of extreme events to the TXX and TNN in 16 climatological stations located in 12 urban areas of Mexico, in time series with different periods, showed general trends in urban warming, especially at night. However, values varied from city to city. A less clear detection signal was the diurnal heating. Although the causes for this behavior of heterogeneous trends for TXX and TNN may be of purely local origin and caused by the dynamics of urbanization (land use change, anthropogenic activities), it cannot be ruled that they could be caused by some other climatic phenomena of a greater scale. For example, Mexico could be exposed to blocking patterns, or modes of climate variability, such as the El Niño Southern Oscillation, Pacific Decadal Oscillation, North Atlantic Oscillation, or the Madden-Julian Oscillation, among others. These could be causal factors that contribute to thermal trends since in other regions of the world they have significantly influenced the behavior of extreme temperatures (^{Arblaster and Alexander, 2012}; ^{Parker et al., 2014}; ^{Burgess and Klingaman, 2015}; ^{Matsueda and Takaya, 2015}). All these factors are considered important for possible attribution but are not the object of this study.

It is well known that cities will remain the most vulnerable geographic spaces, especially in developing countries such as Mexico, as they often concentrate large populations without adequate infrastructure. Thus, knowledge of regions where the climate is unstable (changing), and the knowledge that extreme climate events are not controllable, but have showed signals of an increasing trend, will serve to generate some potential urban adaptation strategies for facing current and future extreme temperatures. For example, such temperatures can be addressed by redesigning the urban space (making it safe, ecological and acceptable) to reduce vulnerability and strengthen urban water security. It should be noted that each city is a case of analysis, and that an adequate scheme of sustainable urban planning can only be accomplished through the participation of governmental authorities, productive sectors, academics, and local actors.

6. Conclusions

This study modeled the TXX and TNN recorded in 16 climatological stations corresponding to 12 cities in Mexico, for periods that varied according to the series of available data. The longest record corresponded to Veracruz with 85 years (1930-2014), and the shortest to Tuxtla Gutiérrez with 31 years (1980-2010). Through the application of the non-parametric Mann-Kendall test and the Sen method, a trend towards urban warming was detected, but no homogeneous behavior in all cities analyzed was found. An important observation was the significant prevalence of non-stationary series with the TXX in half of the cities analyzed; only Guadalajara, located in the center-west of the country, presented a significant negative trend, perhaps due the effect of “diming” by increases in aerosol pollution (^{Fonseca-Hernández et al., 2018}). Therefore, the annual trends in daytime warming, represented by TXX, are not necessarily identical in all selected urban areas of Mexico. In the case of the TNN, which is related to night warming, the behavior was more uniform. Specifically, 90% of the cities are non-stationary with a significant positive trend, and only two areas presented a stationary series: one urban area located to the east of the metropolitan area of the Valley of Mexico (Milpa Alta), and another on the coast of the Gulf of Mexico (Veracruz). With the adjustment of the non-stationary GEV distribution to the data set and by incorporating the trend of the location parameter, the trends of the return levels of 10, 20, 50, and 100 years were estimated, while keeping the shape and scale parameters constant.

The results are very similar to those from non-parametric trend detection methods, both for TXX and TNN. This confirms the non-stationary behavior of half of the stations for TXX, and 90% of the stations for TNN. Return periods of the thermal extremes estimated in a changing climate for many of the cities analyzed in this study vary significantly, and the regularity with which extreme temperatures are occurring is becoming more frequent. Thus, statistical modeling should consider this behavior, due to its importance for risk assessments for human health, flora and fauna, and urban infrastructure. The proposed non-stationary GEV model has provided new and important information concerning changes in extreme temperatures, and by estimating return periods, it has provided their probability of recurrence. Some main modes of climate variability as causal mechanisms of attribution of extreme temperatures have not been studied, making it a pending task.