DOI: 10.20937/ATM.53103

Received: July 15, 2021; Accepted: April 25, 2022

Application of geostatistical models for aridity scenarios in northern Mexico

Javier de Jesus Correa-Islas1, Juan Manuel Romero-Padilla1*, Paulino Pérez-Rodríguez1 and Antonio Vázquez-Alarcón2

1 PSEI-Estadística, Colegio de Postgraduados Campus Montecillo, Km. 36.5, México-Texcoco, 56230 Montecillo, Estado de México.

2 Departamento de Suelos, Universidad Autónoma Chapingo, km 38.5 Carretera México-Texcoco, Chapingo, 56230, Estado de México.

* Corresponding author:


Se calculó un mapa de temperatura media anual mediante el método de interpolación de Kriging para la zona centro-norte de México con la finalidad de obtener la condición actual de aridez, así como posibles escenarios en el futuro cercano y lejano. El gradiente altitudinal se estimó mediante regresión lineal y se usó en la estimación de la temperatura media. Las Áreas de Influencia Climática (CIA) se obtuvieron superponiendo la capa de precipitación oficial y la capa de temperatura media anual con la ayuda de herramientas de Sistemas de Información Geográfica. Se generaron bases de datos mensuales de variables climáticas para cada CIA y se estimó la evapotranspiración potencial utilizando la metodología de Thorthwaite. El Índice de Aridez (IA) se calculó y mapeo para un escenario base (1970-2000). Posteriormente, se proyectó y mapeó el comportamiento de aridez para algunos escenarios, utilizando los modelos de clima global HADGEM 2.0, GFDLCM 3.0, MIP_ESM y CRNMCM5. Se pronosticaron algunos escenarios, en el mejor escenario la aridez debilitará los ecosistemas húmedos y en el peor escenario aparecerán climas hiper áridos en el área de estudio.


An annual mean temperature map was calculated using the Kriging interpolation method for the north-central zone of Mexico to obtain the current aridity, as well as possible scenarios for the near and distant future. The altitudinal gradient was estimated by linear regression, and it was used to estimate the mean temperature. Climate Influence Areas (CIA) were obtained by superimposing the official precipitation layer and the annual mean temperature layer using Geographic Information Systems tools. Monthly databases of climatic variables were generated for each CIA and potential evapotranspiration was estimated using the Thorthwaite methodology. The Aridity Index (AI) was calculated and mapped for a base scenario (1970-2000). Subsequently, the aridity behavior of some scenarios was projected and mapped using the global climate models HADGEM 2.0, GFDLCM 3.0, MIP_ESM, and CRNMCM5. Under the best scenario projected, aridity will weaken the humid ecosystems and in the worst scenario, hyper-arid climates will appear in the study region.

Keywords: Kriging, Climate Change, Potential Evapotranspiration, Linear Regression.

1. Introduction

The United Nations Framework Convention on Climate Change (UNFCCC) defines climate change as “a change attributed directly or indirectly to human activity that alters the composition of the global atmosphere and affects the natural variability of the observed climate over comparable periods”. In Mexico, aridity is promoted by climate change, which directly impacts population development since it affects the production of goods and services. The advance in aridity can contribute enormously to the loss of land productivity and become a severe problem of desertification. According to data from the United Nations Department of Economic and Social Affairs (ONU, 2018), desertification affects approximately one-sixth of the world’s population and 70% of all drylands. Furthermore, desertification decreases food production, increases infertility and soil salinization, increases flooding in lower areas, lead to water scarcity, increases health problems due to wind-borne dust (e.g., eye infections, respiratory diseases, and allergies), changes biological cycles, and reduces livelihoods, which can contribute to stimulate migration (Hori et al. 2011).

Arid, semi-arid, and dry sub-humid lands are found in Mexico mainly in the central regions influenced by the Western and Eastern Sierra Madre mountain ranges. Drylands occupy approximately 101.5 million hectares in Mexico, just over half of the territory; arid zones represent 15.7%, semi-arid regions 58%, and dry sub-humid regions 26.3% (UACh, 2011), so, it is important to study their behavior considering that climate is being modified by humans.

In the present work, aridity is addressed through variables such as temperature, precipitation, and evapotranspiration; statistical tools, geostatistical concepts, model theoretical foundations, theory of spatial autocorrelation, statistical assumptions, and interpolation are described. Different projected scenarios are compared and estimates of current and future weather conditions are presented. The world climate is a whole, and what happens in a particular area may be replicated in other places with similar characteristics; the present study contributes to validate results from other sites worldwide.

2. Study region

Eight economic regions of Mexico are described in reports by the National Commission for the Exploitation and Use of Biodiversity (Bassols, 1985; CONABIO, 2010). Each region can include more than one state that share characteristics of the biophysical environment and with common economic activities in the primary, secondary, and tertiary sectors. In the present work, region II is considered given its current characteristics as dry climates and limiting surface of water distribution, conditions that make it especially vulnerable to climate change from the point of view of drought (CONAGUA, 2021). Region II (Figure 1) covers the states of Chihuahua, Coahuila, Durango, San Luis Potosí and Zacatecas and it is located between 110 º 0 ‘0’ and 95 º 0 ‘0’ west longitude and 35 º 0 ‘0’ and 30 º 0 ‘0’ north latitude.

Figure 1

Fig. 1. Economic regions of Mexico (Bassols, 1985 and CONABIO, 2010).

3. Materials and methods

Precipitation and temperature data were obtained from stations grouped in the Rapid Climate Information Extractor (ERIC) of the National Institute of Water Technology (IMTA) version III (IMTA, 2010). Data quality of many climatological stations is not good, with large gaps of missing data due to little maintenance and even permanent deactivation. Only 496 stations with 80% or more viable information were considered. The analyzed variables were maximum temperature (Tmax) and minimum temperature (Tmin) for the period 1970-2000. Precipitation values were obtained from the Informatics Unit for Atmospheric and Environmental Sciences (UNIATMOS) of the National Autonomous University of Mexico from 1950 to 2000 (UNAM, 2017).

Additional information was obtained from cartography of climate change models, scenarios of climatology (Hijmans, et al., 2005), climate global models, and emission scenarios of the Intergovernmental Panel on Climate Change (IPCC) that were processed and interpolated for Mexico and Central America, based on the Model for the Assessment of Greenhouse-gas Inducted Climate Change – MAGICC (UCAR, 2005).

3.1 Geostatistical analysis

The concept of spatial autocorrelation is supported by Tobler’s principle, which states that in a geographic space everything is related and the closest areas are more related (Tobler, 1970). Spatial autocorrelation helps to understand the degree of similitude of one object with another (ESRI, 2019). In general, there is spatial autocorrelation (SA) whenever there is a pattern in the behavior of the variable according to the geographical location of the data.

Geostatistical techniques require second-order stationarity data, that is, at least the variance must be the same in different sub-regions. The lack of stationarity is commonly linked to anomalies in space, or to the existence of a trend or spatial gradient whose dimension is greater than the study region (Cressie, 1993). Stationarity can be a problem when interpolating points in space, but it does not justify abandoning geostatistics in favor of other interpolation techniques such as inverse distance weighted (IDW) since it is equally sensitive to lack of stationarity as described by Isaaks and Srivastava (1989).

3.2 Moran’s Global Index and G general of Getis-Ord

The global Moran Index (Moran, 1950) measures spatial autocorrelation based on locations and values of the entities simultaneously. Given a set of entities and an associated attribute, it evaluates whether the expressed pattern is clustered, sparse, or random (ESRI, 2019). The Moran index (I) for spatial autocorrelation is given by equation (1):

Eq1 (1)

where, zi = (xi) is the deviation of an attribute or characteristic from its mean; xi is an attribute value for characteristic i; wi,j is a spatial weight between feature i and j, n is the total number of features, and S0 is the aggregate of all spatial weights which is computed as follows, Eq. To test if the spatial configuration occurs at random (H0) vs the alternate hypothesis (H1) that the spatial configuration does not occur at random, we use the standardized Moran index, eq.

The Getis-Ord general statistic (G) (Getis and Ord, 1992) measures the concentration of high or low values in a given study region and is given by equation (2):

Eq2 (2)

where xi and xj are attribute values for characteristics i and j, wi,j is the spatial weight between characteristic i and j, n is the number of characteristics in the data set. The null hypothesis for the High / Low Clustering statistic (General G) states, H0 : There is no spatial grouping of entity values vs the alternate hypothesis H1: There is a spatial grouping of entity values. The test statistic is based on the standardized Getis-Ord index given by Eq with Eq. The zG statistic is important for the interpretation when the null hypothesis is rejected. If zG is positive, high values are clustered, otherwise low values are clustered. The software ArcMap was used to calculate Moran’s Index (I) and Getis-Ord index (G), the z-scores, and the p-values to assess the significance (ESRI, 2019).

3.3. Kriging method

Kriging is an advanced geostatistical procedure that generates an estimated surface from a set of scattered points with z-values. The Kriging tool performs an interactive investigation of the spatial behavior of the phenomenon represented by the z-values before selecting the best estimation method to generate an output surface. The kriging technique assumes that data are correlated in space. The kriging method is considered to be BLUE (Best Linear Unbiased Estimator) and is similar to Inverse Distance Weighted. The general equation is given by:

Eq3 (3)

where: Z(Si) is the value measured at location i; λi is an unknown weight for the value measured at location i; s0 is a location and n is the number of measured values.

In ordinary kriging, λi, depends on a fitted model, the distance to the location, and the spatial relationships between the measured values around the location. The kriging interpolation method describes dependency rules with variograms and covariance functions (spatial autocorrelation). The kriging methodology assumes that data follow a Gaussian distribution, present stationarity, are spatially continuous, and the autocorrelation is known, only then, kriging is considered an optimal predictor (Gallardo, 2006). In this work, ordinary kriging was chosen due to the nature of the data, in addition, it offers the possibility of making different adjustments through pre-established models.

For realizations of a variable Z at points xi, i = 1,2,… n, say, Z (x1), Z (x2),… Z (xn), we want to estimate Z (x0), and there was no measurement at x0. In this circumstance, the ordinary kriging method estimates the variable as a linear combination of the n random variables as follows:

Eq4 (4)

where the λi represent weights of the original values. (x0) is the best linear estimator because the weights are obtained in such a way that they minimize the variance of the estimation error; other interpolation methods such as inverse distances or polygonal do not guarantee minimum estimation variance (Samper and Carrera, 1990). Weights are estimated by minimizing V [(x0) – Z(x0)] subject to eq, using the method of Lagrange multipliers. The weights λ, also, can be estimated through the semi-variance function, for which it is required to know the relationship between the co-variogram and semi-variance functions. The modeling of semi-variograms includes two fundamental stages (Xie and Myers, 1995): estimation of the empirical semi-variogram, and fitting of a theoretical model, to determine the descriptive parameters of the semi-variogram that will later be used in the estimation (Journel and Huijbregts, 1976).

3.4 Aridity index

Temperature varies as a function of altitude; the empirical method developed by Fries et al. (2012) estimates temperature in mountainous regions, calculates the ranges of variation of the annual temperature at different elevations by a simple linear regression model:

Eq5 (5)

where yi is the temperature (in ºC), xi is the elevation (in meters above the sea level), i = 1,…, n the number of data points, α is the intercept or temperature at sea level, while β is the temperature variation for each meter that increases the elevation, as in nature there is an inverse relationship between elevation and temperature, the value of β is negative, ei is the random error, we assume ei~ NIID (0, σe2), where “NIID” stands for Normal Independent and Identically distributed and σe2 is the variance associated to the errors. The elevation was estimated with models of geographical elevation provided by the National Institute of Geography and Statistics (INEGI, 2020).

The mean annual temperature layer was obtained with ordinary kriging interpolation; ranges of 1ºC were generated. The Climate Influence Areas (Gómez et al., 2008) were obtained by cross correlating the cartography of mean annual temperature with the estimated annual precipitation. Potential Evapotranspiration (E) was calculated with Thornthwaite method as follows. First, the monthly heat index was calculated from the monthly mean temperature (Ti, º C), as Eq; subsequently, the annual heat index (I) was calculated with equation 6:

Eq6 (6)

The monthly potential evapotranspiration (mm/month) for months of 30 days and 12 hours of sun (theoretical), E “uncorrected” was calculated also based on temperature data, as:

Eq7 (7)

where a = 675(10–9) (I3) – 771(10–7) (I3) + 1792 (10–5)(I) + 0.49239 is a theoretical adjustment coefficient. Finally, the calculation of the corrected E for month i is obtained through equation 8,

Eq8 (8)

where, Ni is the maximum number of hours of sunshine in month i and di is the number of days in month i. The Aridity Index (AI) is a qualitative characteristic of the climate, which makes it possible to measure the degree of sufficient or insufficient precipitation. The Aridity Index (AI) was established in in Mexico through the Official Gazette of the Federation, on 1 June 1995 to define an aridity classification based on the Precipitation (P) and Evapotranspiration (E) as AI=P/E (DOF, 1995). The AI allows the classification into eight categories: hyper-arid (< 0.051), arid (0.051-0.201), semi-arid (0.201-0.501), dry subhumid (0.501-0.651), humid subhumid (0.651-0.751), humid (0.751-1.251), very humid (1.251-2.5) and per humid (> 2.5).

3.5 Setting scenarios in the near and distant future

Climate change scenarios are representations of future climate, based on an internally coherent set of climatologic relationships, which are constructed to investigate the potential consequences of climate change. Scenarios were constructed by superimposing layers of annual mean temperature and precipitation; and then databases were generated for the potential evapotranspiration and the AI. The centroids were obtained with automatic tracing of isotherms and isohyets; finally, global climate models were used to project data for near future (2015-2039) and far future (2075-2099), considering two values (4.5 and 8.5) of the Representative Concentration Pathway (RCP).

4. Results and discussion

4.1 Global Moran Index and Getis Ord G

The z-score of the Moran Global Index and G Index of Getis Ord, constructed with a weighting matrix under the K-nearest neighbor’s method, were 3.7797 and 5.5738 respectively; in both cases there is less than 1% probability that the pattern is the result of a coincidence, so there is indeed a spatial autocorrelation grouped in clusters for the case of the calculated mean annual temperature.

4.2 Predicted temperature

Fries, et al. (2012) indicate that it is necessary to consider the altitudinal gradient to obtain a correct temperature layer once the interpolation has been carried out, so the fitted linear regression model was:

Eq9 (9)

The residuals were distributed randomly and with a presumably normal histogram, so we can conclude that the model complies with the assumptions of normality and homoscedasticity. Once the linear regression model was fitted, the equation proposed by Fries, et al. (2012) was used to predict the temperature:

Eq10 (10)

where, Tmonth corresponds to the mean annual temperature in º C; Γ corresponds to the altitudinal gradient calculated in the regression model (β̂ = Γ = –0.0042); TDet corresponds to a reference altitude (in this case 2000 meters above the sea level, used in the kriging methodology); and Zest corresponds to the elevation of each weather station in meters above the sea level.

4.3 Data normality

The Kriging procedure uses transformations in case data do not comply with the normality assumption, but it is difficult to use transformed values since it will later be necessary to work with the precipitation variable. A Shapiro-Wilk test was performed to verify data normality, with the hypothesis H0: The population is normally distributed vs. H1: The population is not normally distributed. A test value W = 0.97707 was obtained, with a p-value < 1 × 10–5, so setting the significance of the test at α = 0.05, we conclude that the population does not come from a normal distribution and a type of transformation is required. To correct the problem of normality, a subdivision of the study region was proposed, considering the physiographic provinces proposed by CONABIO, (CONABIO 2001) as shown in Figure 2.

Figure 2

Fig 2. Physiographic sub regions and distribution of meteorological stations.

The physiographic provinces reduce extreme values and have homogeneous distributions. The Shapiro-Wilk goodness of fit test is shown in Table I, indicating that none of the sub-regions has normality problems. Once the assumptions were corroborated, the variograms of each physiographic province were obtained, as well as the adjustment of different models depending on the dispersion of Gamma values, as shown in Table II.

Table I. Results of the Shapiro Wilk goodness of fit test.

Shapiro Wilk’s
W value
P-value for
normality test
R1 0.9893 0.3298
R2 0.9832 0.1325
R3 0.9836 0.2689
R4 0.9824 0.1659

Table II. Summary of adjustment of variograms, models marked in bold have the best fit and were selected.

Model SSR Variogram and Fitted Model Variogram attributes
Psill: 0.9
Psill: 4.0
Psill: 2.0
Psill: 3.2

4.4. Model validation

Model validation through different error indicators was carried out by two methodologies: Cross-validation and validation. Table III presents the summary of cross-validation, obtained through a random sample made up of 80% of the data, the variable temperature was estimated and validated with the stations not included in the sample. In general, values for mean, root mean square, standardized mean, root mean square standardized, and average standard error are small and some of them close to zero, therefore it can be concluded that the adjusted model seems to be correct.

Table III. Summary of the cross-validation indicators, where R1, R2, R3 and R4 represent the previously defined sub regions.

Topic R1 R2 R3 R4
Sample 117 142 173 98
Mean –0.0008 –0.0157 0.0366 0.0245
Root Mean Square (RMS) 0.905 1.5362 0.0258 1.2813
Standardized mean –0.0010 –0.0082 0.0209 0.0174
RMS Standardized 1.0858 1.0823 0.9856 0.8987
Average Standard Error 0.8319 1.3954 1.1772 1.4239
Regression function –0.006*x + 17.17 –0.52 * x + 8.70 0.30 * x + 11.51 0.54 * x + 7.43

For the validation methodology, the data were divided into 2 subsets: a training set and a test set, the subsets are large enough to be representative and to generate statistically significant results. Again the estimations for average error, root mean square, and average standard error are close to zero, which indicates a good fit (Table IV). Once the validation of the models is completed, the temperatures are interpolated and the altitude gradient is adjusted.

Table IV. Summary of the values of the errors obtained in the validation.

Topic R1 R2 R3 R4
Average error 0.075 –0.333 –0.173 0.614
Root mean square 0.829 1.420 1.193 1.420
Average standard error 0.089 –0.240 –0.194 0.429

4.5 Isotherms

Sub-region R1, mainly presents temperatures between 14 ºC and 19 ºC, with some surfaces between 8 ºC and 10 ºC in the elevated regions, as well as temperatures of 22 ºC to 24 ºC in the northwest. For sub-region R2, there are low average temperatures throughout the summit of the western Sierra Madre with temperatures ranging between 8 ºC and 10 ºC, with the highest temperatures also occurring towards the east of the mountain, ranging between 24 and 26 º C. In sub-region R3, the highest temperatures are obtained in the southeast, coinciding with the warmest region of the state of San Luis Potosí, with average temperatures of 25 ºC to 26 ºC, as well as the coldest in the center with isotherms of 9 ºC to 10 º C. In sub-region R4, temperatures are generally observed ranging between 16 ºC and 18 ºC with the warmest surfaces to the east with temperatures ranging between 20 ºC and 23 º C (Fig. 3).

Figure 3

Fig. 3. Temperature map of region II.

4.6 Precipitation

The total annual precipitation was obtained using map algebra

Eq11 (11)

where i corresponds to the i-month evaluated. The result, reclassified to every 100 mm of rainfall, is presented in Figure 4. Twenty-one classes of precipitation were obtained, the most intense colors of blue being the regions with the highest rainfall.

Figure 4

Fig. 4. Precipitation map of region II.

The precipitation is distributed as expected, the most humid regions are those corresponding to the windward part of the western Sierra Madre and the southern part of the state of San Luis Potosí. On the other hand, the regions with less precipitation correspond to the leeward part of both the western and eastern Sierra Madre, coinciding with the theory of Gómez et al. (2008).

4.7 Aridity index of base scenario

Monthly databases of climatic variables were generated for each CIA. Centroids of each polygon were obtained using geometric methods. The aridity index was calculated with estimations of annual mean temperature (Tmed), potential evapotranspiration (E), and precipitation (P). The baseline scenario was built with climate data from the 1970-2000 period. The data for the period 2000-2020 were not used because it shows a considerable decrease in quality in most of the stations considered. The calculated aridity baseline scenario is presented in Figure 5. The percentages for each classification were as follows: hyper-arid (0%), arid (1%), semi-arid (54%), dry sub-humid (15%), humid sub-humid (9%), humid (15%), very humid (6%) and per humid (0%). The baseline scenario for region 2 presents a predominance of a semi-arid climate with more than half of the surface region (in km2), the percentages of Per-humid and Hyper-arid climates are close to 0%. The highest concentration of humid regions is concentrated in the western highlands, more specifically in the windward region, and the driest regions are in the north-central region.

Figure 5

Fig. 5. Aridity Base Scenario of region II (1970-2000).

4.8 Global climate models

Table V and Figure 6 show the estimates of four global climate models considering two Representative Concentration Pathway (RCP) and two horizons (H1: 2015-2039 and H2: 2075-2099), in general, there is a transition from left to right, from humid classes to arid classes, the most optimistic scenario is CRNMCM5 where the “hyperarid” class is not reached and in the worst scenario, the closest horizon only occupies 6245.40 km2 of the “arid” class. The most pessimistic models are MPI-ESP and GFDL-CH3, since in both cases the “hyperarid” class is reached with 5639.26 km2 (0.9%) and 5621.52 km2 (0.9%) respectively. The class “perhumid” represented by just 1 km2 only exists in the base scenario and is extinguished in the rest of the scenarios. A model that is neither very pessimistic nor very extreme is HADGEM 2.0, which predicts an “arid” surface of 349325.99 km2 (53.8%) without reaching the “hyperarid” class.

Table V. Results of estimates in percentage for different horizons in eight climate classes. Four models are presented.

Model RCP H Per
Wed Sub
humid Wed
humid Dry
Arid Hyper
base scenario - 0 0.0% 6.3% 14.6% 9.5% 15.1% 53.5% 1.0% 0.0%
HADGEM 2.0 4.5 1 0.0% 4.8% 14.1% 2.0% 13.5% 64.7% 1.0% 0.0%
2 0.0% 2.8% 15.6% 2.3% 5.1% 73.1% 1.0% 0.0%
8.5 1 0.0% 7.6% 12.9% 3.4% 13.7% 61.5% 1.0% 0.0%
2 0.0% 1.9% 7.0% 1.0% 2.1% 34.3% 53.8% 0.0%
CRNMCM5 4.5 1 0.0% 9.9% 12.6% 6.4% 9.1% 61.1% 0.9% 0.0%
2 0.0% 3.6% 16.4% 1.1% 11.3% 66.7% 1.0% 0.0%
8.5 1 0.0% 9.6% 11.4% 4.6% 12.6% 60.9% 0.9% 0.0%
2 0.0% 2.8% 7.3% 10.5% 8.8% 69.6% 1.0% 0.0%
GFDL-CH3 4.5 1 0.0% 8.8% 15.4% 5.6% 8.3% 61.1% 0.9% 0.0%
2 0.0% 8.6% 15.3% 5.0% 9.1% 61.2% 0.9% 0.0%
8.5 1 0.0% 3.0% 15.4% 1.1% 5.8% 72.9% 1.7% 0.0%
2 0.0% 2.8% 5.6% 8.8% 6.6% 20.9% 54.4% 0.9%
MPI-ESP 4.5 1 0.0% 7.3% 12.9% 2.6% 9.7% 66.5% 1.0% 0.0%
2 0.0% 3.5% 7.0% 10.4% 9.9% 68.3% 1.0% 0.0%
8.5 1 0.0% 3.6% 16.8% 2.4% 9.1% 67.2% 0.9% 0.0%
2 0.0% 1.9% 7.9% 0.3% 10.5% 24.4% 54.2% 0.9%

H: Horizons, 0: 1970-2000; 1: 2015-2039; 2: 2075-2099.
RCP: Representative Concentration Pathway.
HADGEM2.0: Met Office Hadley Centre (MOHC) (United Kingdom model).
CRHMCM5: Centre National de Recherches Météorologiques (French model).
GFDL-CH3: Geophysical Fluid Dynamics Laboratory (USA model).
MPI-ESP: Max Planck Institute for Meteorology (Germany model).

Figure 6

Fig. 6. Mapping of aridity evolution scenarios for models: a) HADGEM2.0; b) MIP-ESP; c) CRNMCM5; d) GFDL-CM3; I: RCP=4.5, Horizon 2015-2039; II: RCP=4.5, Horizon 2075-2099; III: RCP=8.5, Horizon 2015-2039; IV: RCP=8.5, Horizon 2075-2099.

In summary, the results indicate that all the assumptions of the Kriging model were fulfilled and, therefore, the best unbiased linear estimator was developed. The validation model showed that estimates were correct. A base aridity index and projections from global climate models were obtained for near and far futures. Such projections are useful to make decisions in the present and change the projected trend. All the global climate models predict an increase in aridity for the far future, and RCP = 8.5, so we should take action to avoid large RCP levels. The information generated in this study has technical foundations.

For future research, the study can be extended to develop projections for all regions of Mexico. Having reliable information helps us make the best decisions and, in terms of climate, it allows us to reduce climate change effects.

5. Final remarks

The availability of climate information in Mexico is severely affected by different factors, such as meteorological stations with no maintenance and less than 70% of the data. The aridity in regions where land has severe degradation problems can lead to a loss of productive capacity practically irreversible. The use of geostatistical interpolation methodologies for unsampled points is very useful to model the behavior of different climatic variables. Temperature and precipitation detailed maps were obtained for economic region II of Mexico. Different scenarios were predicted with global climate models. Results show adverse scenarios: in the best prediction, aridity will strongly weaken the humid ecosystems, and in the worst scenario hyper-arid climates will appear, practically inhospitable for many of the species with limited environmental resilience that today inhabit these regions. It is urgent to carry out specialized studies of territorial order, natural resources management, regeneration, and adaptation work for the regions where ecosystems are most vulnerable. Some proposals that can be highlighted are: well-managed reforestation plans, constant cooperation with multilateral organizations, transition towards clean and renewable energies, reducing emissions into the atmosphere, and reducing the potential of the highest Representative Concentration Pathways.

Although Mexico does not have an atmospheric circulation model, some efforts have been made to represent the country’s environmental conditions, and this study, supported by geostatistics, is one of them.


Bassols BA. 1985. Recursos naturales de México: teoría, conocimiento y uso. Editorial Nuestro Tiempo. Segunda Ed.

CONABIO. 2001. Provincias Fisiográficas de México, México DF.

CONABIO. 2010. Catálogo de metadatos geográficos. Comisión Nacional para el Conocimiento y Uso de la Biodiversidad. Available at: (accessed 2021 January 3)

CONAGUA. 2021. Monitor de Sequía en México. Available at: (access 2021 May 19)

Cressie N. 1993. Statistics for Spatial Data. Revised Edition. John Wiley and Sons. New York.

DOF. 1995. Diario Oficial de la Federación, (DOF: 07/06/1995). Available at: (accessed 2021 May 25)

ESRI. 2019. Cómo funciona Autocorrelación espacial. Available at: (accessed 2021 January 3)

Fries A, Rollenbeck R, Naub T. 2012. Near surface air humidity in a megadiverse Andean mountain ecosystem of southern Ecuador and its regionalization. Agricultural and Forest Meteorology, 152: 17-30.

Gallardo A. 2006. Geostadística. Ecosistemas, 15(3): 48-58.

Getis A, Ord JK. 1992. The Analysis of Spatial Association by Use of Distance Statistics. Geographical Analysis , 24: 189-206.

Gómez JD, Monterroso-Rivas A, Gay C. 2008. Spatial estimation of mean temperature and precipitation in regions of scarce meteorological information. Atmósfera, 21: 35-56.

Hijmans RJ, Cameron SE, Parra JL, Jones PG. 2005. Very high resolution interpolated climate surfaces for global land areas. International Journal of climatology, 25: 1965-1978.

Hori Y, Stuhlberger C, Simonett O. 2011. Desertification: a visual synthesis. United Nations Convention to Combat Desertification (UNCCD), First edition, GRAPHI 4 Press, Bresson France.

IMTA. 2010. Extractor Rápido de Información Climatológica, México: Instituto Mexicano de Tecnología del Agua.

INEGI. 2020. Continuo de Elevaciones Mexicano (CEM). Available at: (accessed 2021 January 3).

Isaaks E, Srivastava R. 1989. An introduction to applied geostatistics. Oxford University Press.

Journel AG, Huijbregts CJ. 1976. Mining Geostatistics. United Kingdom.

Moran PAP. 1950. Notes on Continuous Stochastic Phenomena. Biometrika, 37: 17–23.

ONU. 2018. Ordenación de los ecosistemas frágiles: Lucha contra la desertificación y la sequía. Departamento de Asuntos Económicos y Sociales. Available at: (accessed 2021 January 3)

Samper CFJ, Carrera RJ. 1990. Geoestadística. Barcelona: Centro Internacional de Métodos Numéricos en Ingeniería.

Tobler WR. 1970. A computer movie simulating urban growth in the Detroit region. Economic Geography, 46: 234-240.

UACh. 2011. Actualización y delimitación de zonas áridas, semiáridas y subhúmedas de México. Reporte final de proyecto de investigación. Departamento de Suelos. Universidad Autónoma Chapingo.

UCAR. 2005. Model for the Assessment of Greenhouse-gas Inducted Climate Change. A Regional Climate Scenario Generator. Available at: (accessed 2021 January 3)

UNAM. 2017. Portal de datos de escenarios de cambio climático. Responsorio Institucional Instituto de Ciencias de la Atmósfera y Cambio Climático, UNAM. Available at: (accessed 2021 January 3)

Xie T, Myers DE. 1995. Fitting Matrix-Valued Variogram Models by Simultaneous Diagonalization (Part I: Theory). Mathematical Geology, 27(7): 867-875.