Past and Future Spatial Growth Dynamics of Chihuahua City, Mexico: Pressures for Land Use

Abstract

In this study, the transitions of land use that occurred in the urban and peripheral areas of Chihuahua City, Mexico, were determined for the period 1989–2014. Landsat TM and OLI scenes, as well as the method of Markov Chains (MC) were used. Grasslands and Shrublands were the land uses that experienced the highest pressures for land use. Grasslands occupied 23.5% of the area in 1989, decreasing to 16.01% in 2014. Likewise, Shrublands were reduced from 54.53% to 48.06%. The areas occupied by Croplands, Oak forest, Water bodies and Riparian vegetation stayed in general constant. In contrast, the urban area increased from 13.6% to 28.6% of the total area studied. In addition, projections of land use for 2019 and 2024 were generated through the method of MC and Cellular Automata (CA). According to the projections, validated with an agreement of 0.90, the Human settlements would continue to expand, occupying 38.57% by 2019 and almost half of the studied territory (47.33%) by 2024. The ecosystems with the highest pressure for land use change will continue to be the Grasslands and Shrublands. By 2024, the former would lose 15.8% while the latter would lose 16.7% of the area. These methods are valuable for urban planning and the results could support growth plans for Chihuahua City, Mexico, with a sustainable approach.

1. Introduction

The increase in population and urbanization is one of the most complex processes because it involves changes in land use and vegetation at local, regional and global scales [1,2]. Although urban areas cover only 2% of the planet’s surface, they have significantly altered the natural landscape [3,4,5,6]. During the last decade, urban sprawl has become a topic of particular interest due to the accelerated growth of human settlements on the planet and the great impact involved in the phenomenon [7,8,9,10].

Cities are responsible for the production of 78% of the greenhouse gases, contributing significantly to global climate change [11]. Other effects of urbanization include alteration of the biogeochemical cycles [12] and the reduction of areas dedicated to agricultural crops, grasslands, forests and in general of the ecosystems located nearby. This has resulted in land fragmentation and degradation [13]. Therefore, the understanding of the growth dynamics of urban areas is of great importance to elaborate better and more environmentally friendly urban growth plans, and to take actions for the preservation of the natural resources [14].

To analyze the structure and growth dynamics of urban systems it is necessary to link the spatial patterns with the landscape to quantify the causes and consequences of their evolution [15]. Several methods for detecting changes in the urban area are based on remote sensing [16,17,18]. Such methods either employ multi-temporal analyses of satellite images using algebra of maps [19] or apply imaging spatial regression techniques [20]. The latter are the ones most recently employed to estimate land use through the variation of a regression model [21]. However, they have limitations for the quantification of changes on a temporal basis [22].

Markov Chains (MC) and Cellular Automata (CA) are stochastic models that incorporate the interaction of spatial and temporal dynamics [22,23,24,25,26]. These methods can serve to analyze the dynamic behavior of land use in a time-space pattern and provide forecasts of future changes that can help in decision-making [23,27]. Some studies have shown the strong capabilities of traditional Markov models to describe trends in land use change [28,29,30]. Even though the Markov analysis itself cannot simulate and predict changes in land use, MC together with CA have the capability of determining the spatiotemporal dynamics and project future scenarios when fed with appropriate susceptibility and limitations criteria [31,32,33]. Therefore, the integration of MC and CA give complementary results [34]. The method of MC quantifies the transition changes based on the past while CA uses this parameter to estimate changes in the future and their location [35].

Chihuahua City, Mexico, has experienced rapid growth in past decades. From 8489 ha occupied in 1980, Chihuahua City grew to 19,024 ha by 2005 [36]. This urban growth has caused a process of fragmentation and loss of biodiversity, resulting in significant losses of area for the natural ecosystems that were once located in the peripheral areas of the city. Such ecosystems included mainly Grasslands and Shrublands. These Grasslands are immersed in the Chihuahuan desert and they belong to the North America Grasslands Priority Conservation Areas [37]. Besides that, Grasslands are one of the most threatened ecosystems on the planet [13], and, specifically in this region, they possess a great biodiversity and a high degree of endemism [38].

Chihuahua City requires high inputs of water for domestic and industrial operations. It has been reported that a total of 150.2 × 10⁶ m³ of water is spent by the city on an annual basis. This concentrates great pressure over the aquifers due to the amount of water extracted from them. Besides that, some of the city’s growth has occurred over their recharged zones [38]. In addition, the growth of the city is expected to continue at high rates in the coming years. The lack of local policy on this topic in Chihuahua is threatening the sustainability of the water governance system on a long-term basis, with serious externalities on other areas such as agriculture [39]. If this is regulated, urban growth would occur by taking into account both the space demand and the impact over the natural resources. However, the magnitude and direction of such growth is not known precisely, limiting urban managers for making effective growth plans to mitigate environmental impacts.

The objective of this study was to analyze the growth dynamics and the pressures for land use change in the urban and peripheral areas of Chihuahua City, Mexico. The analysis was based on the methodologies of MC and CA. The land use transitions for the period 1989–2014 were determined through MC. In addition, projections of land use for the years 2019 and 2024 were generated by CA. Analysis and discussions on the effects of the future growth of the city over the nearby ecosystems are presented.

2. Materials and Methods

2.1. Study Area

The urban and peripheral areas of Chihuahua City, Mexico, were studied. The city lies at the geographic coordinates of 28°40′N and 106°05′W (Figure 1). The topography of the area has elevations ranging from 1306 to 2665 m above the sea level. The land uses of the nearby city areas are Grasslands, Shrublands, Oak forest, Water bodies, Croplands and Riparian vegetation. In 2010, Chihuahua had a total population of 819,543 [40].

2.2. Collection and Pre-Processing of the Data

Four scenes covering the study area and taken by the Landsat sensor (Path 32, Row 40) were used. The spatial resolution of the scenes was 30 m × 30 m. The four scenes corresponded to the years 1989, 1999, 2009 and 2014 and they were obtained from the United States Geological Survey [41]. The characteristics of each of the scenes are presented in

Table 1. Characteristics of the scenes corresponding to the urban and peripheral areas of Chihuahua.

Satellite	Capture Data	Characteristics	Path/Row
Landsat (TM)	1989	7 spectral bands, 30 m resolution	32/40
Landsat (TM)	1999	7 spectral bands, 30 m resolution	32/40
Landsat (TM)	2009	7 spectral bands, 30 m resolution	32/40
Landsat (OLI)	2014	8 spectral bands, 30 m resolution; panchromatic band, 15 m resolution	32/40

.

The scenes were radiometrically corrected. The conversion from digital numbers (DN´s) to reflectance values was performed with the Top of the Atmosphere (TOA) process, which allows making comparisons among images from different dates. The radiometric conversion for the Landsat TM sensor was performed by following Equations (1) and (2), where the spectral radiance ($L_{\lambda }$) and the TOA reflectance (${\rho }_{\lambda }$) were obtained:

$$L_{\lambda }=\left(\left(Lma{x}_{\lambda }-Lmi{n}_{\lambda }\right)/\left(QCALmax-QCALmin\right)\right)\ast \left(QCAL-QCALmin\right)+Lmi{n}_{\lambda }$$

(1)

$${\rho }_{\lambda }=\frac{\pi \ast L_{\lambda }\ast d^2}{ESU{N}_{\lambda }\ast cos{\theta }_s}$$

(2)

where $QCAL$ is DN; $QCALmin$ and $QCALmax$ are the minimum and maximum quantized calibrated pixel value, respectively; $Lmi{n}_{\lambda }$ is the spectral radiance scales to $QCALmin$; $Lma{x}_{\lambda }$ is the spectral radiance scales to $QCALmax; d$ is the distance from the earth to the sun; $ESU{N}_{\lambda }$ is the mean solar exoatmospheric irradiance; and ${\theta }_s$ is the solar zenith angle.

In the case of the data from the Landsat OLI, the radiometric conversion was performed applying Equation (3).

$${\rho }_{\lambda }^{\ast }=\frac{{\rho }^{\lambda ´}}{sin{\theta }_{SE}}$$

(3)

where ${\rho }^{\lambda ´}$ is the TOA planetary reflectance, with a correction for the solar angle, and ${\theta }_{SE}$ is the local sun elevation angle.

For the reflectance normalization of the images from 1989, 1999 and 2009, the image from Landsat OLI was used. This process allowed an improvement on the histograms by modifying the brightness values in the images from 1989, 1999 and 2009, taking as a reference the image from 2014. With this, spectral variations of the land use covers were minimized [42].

As a final step for data processing, the scenes were edited and the study area was defined on the images by using the software ArcMap 10.2^®. The edges were made similar to those of the watersheds of the area. Such edges were taken from the digital elevation model of the state of Chihuahua. All scenes were ensured to cover the same area after the edition process.

2.3. Land Use Classification

Image layer stacking was performed with the software ERDAS^®. This procedure allowed generating false/true colored images, which were required for the analysis of land use classification. For the case of the Landsat TM sensor, band 5 corresponding to the infrared channel (1.55–1.75 µm), band 4 in the near infrared range (0.76–0.90 µm), and the band 3 in the red range (0.63–0.69 µm) were used. These combinations were made based on the recommendations by Lillesand and Kiefer [43] and applied to the images from the years 1989, 1999 and 2009. Likewise, band 6 corresponding to the medium infrared channel (1.57–1.65 µm), band 5 in the near infrared range (0.85–0.88 µm), and band 4 in the red range (0.64–0.67 µm) were used for the Landsat OLI. These combinations were applied to the image from the year 2014.

A classification based on the method of maximum likelihood was applied to obtain the information of land use (Equation (4)). This method employed Gaussian probability. As a result, thematic maps of land use were obtained for each year indicated in the second column of Table 1. The classification areas included the following land use types: (1) Croplands; (2) Human settlements; (3) Shrublands; (4) Grasslands; (5) Oak forest; (6) Water bodies; and (7) Riparian vegetation. The image carries properties allusive to each type of land use as shown in

Table 2. Definition of land uses for the study area and properties of corresponding training areas.

Land Use	Properties of the Training Areas from Each Land Use Type
Croplands	Irregular shape composed of pixels colored “Peridot Green” and “Ultramarine”
Human settlements	Irregular shape composed of pixels colored “Glacier Blue”
Shrubland	Irregular shape composed of pixels colored “Lime Dust”
Grasslands	Irregular shape composed of pixels colored “Cantaloupe”
Oak forest	Irregular shape composed of pixels colored “Green Leaf”
Water bodies	Irregular shape composed of pixels colored “Dark Navy”
Riparian vegetation	Irregular shape composed of pixels colored “Quetzal Green”

.

$$gi \left(x\right)=In p \left({\omega }_i\right)-\frac{1}{2} In \left|\sum i\right|-\frac{1}{2} {\left(x-m_i\right)}^T{\sum i}^{-1} \left(x-m_i\right)$$

(4)

where $g$i is the class, x represents the n-dimensional data (where n is the number of bands), p (${\omega }_i$) is the probability that the class ${\omega }_i$ appears in the image and that is assumed for all classes, |∑i| is the determinant of the co-variance matrix with data from class ${\omega }_i$, T is the transposed matrix, ∑i⁻¹ is the inverse matrix and m_i is the vector.

2.4. Land Use Classification Accuracy

The cartography of land use land cover from the government of the state of Chihuahua (scale 1:50,000) was used for accuracy validation [44]. The cartography corresponding to land use land cover from the Mexican Institute of Statistics, Geography and Informatics (INEGI) (scale 1:250,000) was also employed [40]. In addition, data from field sampling and photointerpretation were considered during the validation process.

The statistical index K_APPA [45] was used to determine the accuracy of the land use classification on the maps. K_APPA is a discrete multivariate technique for comparing classes through a matrix [46]. It can be used to establish the degree of similarity between mapped and actual or real values of land use [47]. A K_APPA value equal to one indicates a 100% similarity between mapped and real values. Conversely, a value equal to zero suggests a similarity of 0%. K_APPA is represented by Equation (5).

$$K_{APPA}=\frac{N{\sum }^K{X}_{ii} {\sum }^k\left(X_{i+}\ast X_{+i}\right)}{N^2-{\sum }^k\left(X_{i+}\ast X_{+i}\right)}$$

(5)

where K_APPA is the Kappa index; k is the total number of matrix rows; X_ii is the observation number on row i and column i (along the diagonal); X_i+ and X_+i are total marginal for row i and column i, respectively; and N is the total number of observations.

To estimate the K_APPA index, a group of sample points in the peripheral area was employed. The resulting accuracy is shown in

Table 3. Accuracy of the land use classification.

Land Use	Precision
	1989	1999	2009	2014
Croplands	0.82	1.00	1.00	1.00
Human settlements	1.00	1.00	0.83	1.00
Shrublands	0.61	0.85	0.57	0.44
Grasslands	0.83	0.80	1.00	0.51
Oak forest	0.87	1.00	1.00	1.00
Water bodies	1.00	1.00	1.00	1.00
Riparian vegetation	1.00	1.00	1.00	0.81
General Precision KAPPA	0.87	0.95	0.91	0.82

.

2.5. Models of Markov Chains and Cellular Automata

To elaborate projections of land use change for future years, the geosimulation techniques of MC and CA were employed [48]. They account for the changes in land use between two dates by extrapolating them assuming constant changes [49]. The CA technique includes a simulation model where space and time are discrete variables while the assigned interactions are local variables [30]. It is fed with the results from the MC methodology to simulate land use in a future time. In this study, the classifications of 1989 and 1999 were used to estimate the land uses of 2009. Likewise, the land uses of 2014 were estimated based on the classifications of 1999 and 2009.

The MC methodology was implemented in the MARKOV module of the software Idrisi Selva^®. This technique is based on a stochastic model that describes the probability of change from one state to another through a transition probability matrix [28]. The results of the MC methodology include a probability matrix, a matrix of transition areas and transition probability maps. The probability matrix includes the probability of one land use to change from one category to another. This matrix is the result of the crossing between the images after setting a proportional error. The matrix of transition areas states the number of pixels that are expected to change from each land use to the others during the period of time analyzed. The transition probability maps are generated based on the projections of possible changes during the period analyzed. In this study, changes in land uses of 1989, 1999, 2009 and 2014 were used to develop the transition probability matrix that helped to develop the land use projections for 2019 and 2024. The mathematical expression of the transition probability is:

$${\displaystyle \sum _{I=1}^m}Pij=1 i=1,2\dots \dots m$$

(6)

$$P=\left(Pij\right)=\begin{array}{ccc}P11& P12\dots & P1m\\ P21& P12& P2m\\ Pm1& Pn2& Pmm\end{array}$$

(7)

where Pij is the the probability of transition from one land use to another, and m is the total number of land use types of the study area. Pij values stay within the range 0–1.

According to the non-after effect of the Markov methodology and the condition equations of Bayes, the Markov prediction model is obtained (Equation 8):

$$P\left(n\right)=P_{\left(n-1\right)}{P}_{ij}$$

(8)

where $P\left(n\right)$ is the probability $P$ in time $n$, $\left(n-1\right)$ is the probability of the previous time $n$.

The combination of MC and CA was implemented through the module of CA_Markov available in the software IDRISI Selva^®, which allows simulating the dynamics of growth based on the increase of the number of pixels. Each pixel can take a value from a finite set of states [34]. All pixels are affected by a transition function that takes as arguments the measured values of the pixels and the values of the neighboring pixels as a function of time. Cellular Automata and MC were implemented to simulate the land uses of 2009 and 2014. Since this period comprehends five years, periods of the same number of years were used to estimate land uses of the future. That is, land uses were simulated for the years 2019 and 2024.

For the simulations, it was assumed that the probabilities of change were low and constant during the periods analyzed. Thus, the transition probability matrix created from the changes observed between 2009 and 2014 was used to simulate the land uses of 2019. Likewise, the transition probability matrix of 2014 and 2019 was used to estimate the land uses of 2024.

In an iterative process, the module of CA_Markov uses the transition probability of each land use to establish the susceptibility of each pixel, based on their properties, to be occupied by each of the other types of land use. While that is performed, a spatial filter restricts the susceptibility of pixels located away from the class being processed, which is done by assigning a value of greater preference to neighboring areas. In this study, spatial filters of 5 × 5 were applied (Figure 2).

2.6. Suitability Parameters and Limitations for Urban Growth

The parameters to define the suitability of a given pixel to change from one type of land use to another were defined. Such parameters serve to represent the susceptibility of the land to be occupied by each of the other land uses. These parameters were assigned to the variables of elevation, slope, distance to rivers and distance to roads. It was assumed that these variables remained unchanged over the 25 years represented by the dates of the oldest and most recent satellite images analyzed in this study. The suitability parameters for the variables of elevation and slope were defined as zero for not suitable and one for suitable (Figure 3). The suitability parameters for distance to rivers and distance to roads were assigned on a scale of 0 to 255 to represent a minimum and maximum convenience, respectively (Figure 4).

2.7. Validation of the Land Use Change Projections

The transition probability matrix, matrix of transition areas and transition maps for 1989 and 1999 were created with the Markov methodology. Then, we used the CA_Markov module to simulate the land use of 2009. Likewise, the transition probability matrix, the matrix of transition area and the transition maps of 1999 and 2009 were generated and then used to simulate the change in land use for 2014.

The validation of the model to simulate land use change was conducted by comparing the results of the estimated land use changes with the land uses verified through supervised classifications for 2009 and 2014. For this, a randomly stratified sampling design was used and the K_APPA index was employed. Once the simulated land uses were validated, estimations of land use for 2019 and 2024 were made. Finally, the procedure used in this study is summarized in Figure 5.

3. Results

3.1. Detection of Land Use Changes

Results from the analysis of land use show a remarkable gain for the surface area corresponding to Human settlements (

Table 4. Occupation percentages of seven land use types in the urban and peripheral areas of Chihuahua City in 1989, 1999, 2009 and 2014.

Land Use	Occupation (%)
	1989	1999	2009	2014
Croplands	4.12	3.76	3.66	3.48
Human settlements	13.57	17.01	24.63	28.58
Shrublands	54.53	53.64	49.35	48.06
Grasslands	23.50	21.67	18.18	16.01
Oak forest	2.89	2.86	2.86	2.86
Water bodies	0.07	0.08	0.11	0.11
Riparian vegetation	1.31	0.97	1.21	0.90

). Chihuahua City increased more than twice its occupied area in the past 25 years. The expansion occurred mainly to the north and southeast directions (Figure 6). In these directions, the lowest elevations exist and this condition makes the terrain desirable for residential development. Of the categories analyzed, the only one that showed continuous growth is the Human settlements, with 13.57%, 17.01%, 24.63% and 28.50% for 1989, 1999, 2009 and 2014, respectively. In contrast, the land uses of Shrublands and Grasslands showed a reduction in their area during the same period (Figure 7). Regarding the classes of Croplands, Oak forest, Water bodies and Riparian vegetation, their occupied surface areas stayed in general constant. Each of these classes presented a change lower than 1% for the period 1989–2014.

The dynamics of land use changes are presented in

Table 5. Change dynamics of seven types of land use in the urban and peripheral areas of Chihuahua City. Positive and negative numbers indicate gains and losses of surface, respectively.

Category	Surface (ha)			Total
Year	1989–1999	1999–2009	2009–2014
Croplands	−290	−82	−146	−518
Human settlements	2772	6143	3182	12,097
Shrublands	−26	0	0	−26
Grasslands	12	21	0	33
Oak forest	−720	−3458	−1041	−5219
Water bodies	−1471	−2820	−1746	−6037
Riparian vegetation	−273	195	−249	−327

. The Grasslands presented the greatest loss of surface with 1471 ha during the period 1989–1999. This land use increased its loss in surface to 2820 ha during the period 1999–2009 and further reduced its area to 1746 ha during the period 2009–2014. The area occupied by Shrublands was, after the Grasslands, the category that lost most of the surface area. By contrast, the Human settlements experienced the largest growth with a total increase of 12,097 ha for the period of 1989–2014. The location of housing projects has markedly contributed to the increase in the area occupied by Human settlements.

The transition probabilities of the land use corresponding to the periods 1989–1999, 1999–2009, and 2009–2014 are shown in

Table 6. Transition probability matrix for the periods 1989–1999, 1999–2009, and 2009–2014.

Land Use	Year	CL	HS	S	G	OF	WB	R
CL	1989–1999	0.81	0.05	0.13	0.00	0.00	0.00	0.00
	1999–2009	0.83	0.16	0.00	0.00	0.00	0.00	0.00
	2009–2014	0.85	0.14	0.00	0.00	0.00	0.00	0.00
HS	1989–1999	0.00	0.89	0.10	0.00	0.00	0.00	0.00
	1999–2009	0.00	0.89	0.10	0.00	0.00	0.00	0.00
	2009–2014	0.01	0.90	0.01	0.01	0.01	0.01	0.01
S	1989–1999	0.00	0.11	0.87	0.00	0.00	0.00	0.00
	1999–2009	0.00	0.15	0.82	0.00	0.00	0.00	0.01
	2009–2014	0.00	0.12	0.87	0.00	0.00	0.00	0.00
G	1989–1999	0.00	0.16	0.00	0.82	0.00	0.00	0.00
	1999–2009	0.00	0.23	0.01	0.75	0.00	0.00	0.00
	2009–2014	0.00	0.20	0.00	0.79	0.00	0.00	0.00
OF	1989–1999	0.00	0.00	0.11	0.00	0.88	0.00	0.00
	1999–2009	0.01	0.01	0.01	0.01	0.90	0.01	0.01
	2009–2014	0.01	0.01	0.01	0.01	0.90	0.01	0.01
WB	1989–1999	0.01	0.01	0.01	0.01	0.01	0.90	0.01
	1999–2009	0.01	0.01	0.01	0.01	0.01	0.90	0.01
	2009–2014	0.01	0.01	0.01	0.01	0.01	0.90	0.01
R	1989–1999	0.02	0.00	0.29	0.00	0.00	0.01	0.66
	1999–2009	0.00	0.16	0.00	0.00	0.00	0.00	0.83
	2009–2014	0.01	0.01	0.01	0.01	0.01	0.01	0.90

CL, Croplands; HS, Human settlements; OF, Oak forest; WB, Water bodies; S, Shrublands; G, Grasslands; R, Riparian vegetation.

. Bold numbers along the diagonal show the transition probabilities among the study periods. The probability of change from the land use of agriculture to the one of Human settlements was 16% for the period 1999–2009 and decreased to 14% during 2009–2014. Shrublands had a probability of change of 11% during 1989–1999, 15% during 1999–2009 and 12% during 2009–2014. Moreover, Grasslands presented the highest increase in the probability of change during the same period: 16% for the period 1989–1999, 23% for the period 1999–2009 and 20% for the period 2009–2014.

This transition probability matrix (Table 6) indicates that the classes of Croplands, Shrublands, Oak forest, Water bodies and Human settlements had been stable with a tendency to stay in the same land use during the periods 1989–1999 to 2009–2014, as indicated by the probabilities close to 1.0 in the transition matrix. Regarding Grasslands, there was a decrease in the transition probability of 0.82 from the period of 1989–1999 to 0.75 during 1999–2009 and to 0.79 for the period 2009–2014. The decrease in the transition probability indicates an increased likelihood for change of the Grasslands to another type of land use.

The algorithm of CA_Markov determines the exact location of the changes. Thus, the probability of change of any pixel depends on the previously assigned filter and on the restriction and suitability parameters applied to neighboring pixels. The Human settlements showed an increase in the transition probability, indicating a stabilization of this class due to population growth [40], as shown in Figure 8.

The previous results indicate three main findings from the analysis of the growth dynamics of Chihuahua City: (1) the class that had the greatest loss of land use and the one that had the greatest transition change was Grasslands; (2) it is expected that Croplands continue to change to urban areas; and (3) Shrublands are another type of land use that will increase its probability of change to urban area (Table 6). This shows a dominance of the urban land use over the Croplands, Grasslands and Shrublands. In addition, this indicates a continuous change and a pressure on the neighboring natural resources due to the expansion of the urban area through time.

3.2. Validation of the Land Use Change Projections

For the validation of land use change projections, the land uses of 2009 and 2014 were simulated with CA_Markov. The simulated land uses were then compared with the actual land uses resulting from the supervised classifications of the same years. The comparison of simulated and actual values of land use are shown in Figure 9, where deviations can be visually assessed.

It can be observed in Figure 9 that for both, 2009 and 2014, the differences between the actual and simulated land uses were small. The overall accuracy was 90% for 2009 and 91% for 2014. The greatest accuracy was obtained for the class of Oak forest, while the class that showed the lowest accuracy was Human settlements with 0.79 and 0.70 for 2009 and 2014, respectively (

Table 7. Agreement between the real and simulated values of seven land uses of 2009 and 2014 in the urban and peripheral areas of Chihuahua City.

Land Use	Agreement between Real and Simulated
	2009	2014
Croplands	0.88	0.86
Human settlements	0.79	0.70
Shrublands	0.90	0.95
Grasslands	0.85	0.85
Oak forest	0.99	0.99
Water bodies	0.99	0.99
Riparian vegetation	0.92	0.99
General Precision KAPPA	0.90	0.91

).

The 100% agreement between the simulated and real values for the class of Oak forest were obtained because the probability of change of this class was unaffected by other classes. In addition, this class is not geographically close to the urban area. In the case of Human settlements, the resulting lower agreement could be due to the growth patterns influenced by different causes such as socio-economic, political, and other larger contextual factors. Still, all agreement values estimated during the validation process could be considered as very good. The values of K_APPA close to 1.0 indicate a high similarity between the simulated and real values, as well as a spatial distribution of the simulated land uses close to reality. Based on the above, CA_Markov model can be used to predict future changes in land uses in the study area.

3.3. Simulations of Land Use

The model of CA_Markov was used to simulate the areas occupied by seven land uses for the years 2019 and 2024, as an effect of the urban development. For the simulation of the land uses of 2019, the transition probability matrix of the period 2009–2014 was used, setting 2009 as the base year and considering the following five years of change. Afterwards, the year of 2019 was established as the starting point to generate the land uses of 2024. For that, the transition probability matrix of the period 2014–2019 was employed.

The transition probabilities for the years 2019 and 2024 are shown in

Table 8. Transition probability matrix for 2019 and 2024.

	Year	CL	HS	S	G	OF	WB	R
CL	2014–2019	0.74	0.25	0.00	0.00	0.00	0.00	0.00
	2019–2024	0.68	0.31	0.00	0.00	0.00	0.00	0.00
HS	2014–2019	0.00	0.89	0.10	0.00	0.00	0.00	0.00
	2019–2024	0.00	0.89	0.10	0.00	0.00	0.00	0.00
S	2014–2019	0.01	0.18	0.79	0.00	0.00	0.00	0.00
	2019–2024	0.01	0.18	0.79	0.00	0.00	0.00	0.00
G	2014–2019	0.00	0.32	0.02	0.64	0.00	0.00	0.00
	2019–2024	0.00	0.38	0.03	0.58	0.00	0.00	0.00
OF	2014–2019	0.01	0.01	0.01	0.016	0.90	0.01	0.01
	2019–2024	0.01	0.01	0.01	0.01	0.90	0.01	0.01
WB	2014–2019	0.01	0.01	0.01	0.01	0.01	0.90	0.01
	2019–2024	0.01	0.01	0.01	0.01	0.01	0.90	0.01
R	2014–2019	0.01	0.01	0.01	0.01	0.01	0.01	0.90
	2019–2024	0.01	0.01	0.01	0.01	0.01	0.01	0.90

CL, Croplands; HS, Human settlements; OF, Oak forest; WB, Water bodies; S, Shrublands; G, Grasslands; R, Riparian vegetation.

. The class of Grasslands, which had the lowest probability represented by a 0.64 in 2014–2019 and 0.58 in 2019–2024, is the class which would be subjected to the greatest probability of change with a value calculated of 0.79 for 2014 (Table 6) and an estimated and reduced value of 0.58 for the next 10 years (Table 8).

It is estimated that Human settlements is the land use that will show the greatest percentage of change in coverage beginning with 28.57% in 2014 (Table 4) and reaching 38.57% in 2019 (

Table 9. Simulated coverage areas and percentage of occupation of seven land uses of the urban and peripheral areas of Chihuahua City by 2019 and 2024.

Year	2019		2024
Land Use	Surface (ha)	Occupation (%)	Surface (ha)	Occupation (%)
Croplands	2565.23	3.22	2328.43	2.92
Human settlements	30,666.91	38.57	37,636.71	47.33
Shrublands	33,667.13	42.34	30,075.10	37.82
Grasslands	9244.07	11.62	6105.60	7.67
Oak forest	2305.23	2.89	2305.26	2.89
Water bodies	87.02	0.10	86.90	0.10
Riparian vegetation	975.01	1.22	972.73	1.22
Total	79,511	100	79,511	100

). This type of land use will end up covering almost half of the study area (47.33%) in 2024. The classes of Water bodies, Oak forest and Riparian vegetation will be nearly unaffected by the urban growth. However, according to the results from the projections for 2019 and 2024, Human settlements will be extended over adjacent and sensitive areas such as Grasslands and Shrublands (Figure 10).

The estimated dynamics of land use gains and losses for the periods 2014–2019 and 2019–2024 are shown in

Table 10. Dynamics of land use changes for seven land use types of the urban and peripheral areas of Chihuahua City by 2019 and 2024.

Year	2014–2019	2019–2024
Land Use	Surface (ha)	Surface (ha)
Croplands	−241.76	−236.79
Human settlements	7629.91	6969.79
Shrublands	−5070.86	−3592.03
Grasslands	−3662.92	−3138.47
Oak forest	−0.76	0.02
Water bodies	0.02	−0.12
Riparian vegetation	248.01	−1.27

. During these two periods, it is estimated that Shrublands and Grasslands will suffer the greatest losses of areas. Human settlements, represented by gains of 7,629.92 and 6,969.79 ha for the periods of 2014–2019 and 2019–2024, respectively, is the class that would experience the greatest expansion.

4. Discussion

The model of CA_Markov has been widely used for simulations of land use changes and its impact on the landscape by projecting possible trends. In this study, the model was used to project the land uses of the urban and peripheral areas of Chihuahua City by 2019 and 2024. To first understand the urban growth dynamics, this research determined the land uses of the past years as a basis for simulating future changes in the study area. This tool is an alternative means of support for urban planners.

Remote sensing produces valuable data with quick acquisition, which can be used for analyses of land use; for example, the data from the Landsat sensor, which provides images taken from 1972 to date. This satellite has a worldwide coverage with a medium spatial resolution. The Markov prediction method employs the historical data from the Landsat sensor to analyze the dynamic behavior of land use in a time-space pattern. Based on that, forecasts of future changes are estimated.

The methodology employed in this study showed a good level of precision, with values of the K_APPA index above 0.82. This precision is comparable with the ones estimated in other studies employing similar methodologies [2]. The high precision in this case was due to a clear spatial distribution of the land uses in the study area. These land uses are strongly related to the topography where the City is located. The plain areas are clearly occupied by ecosystems of Grasslands while surfaces conformed by terrains with slight slopes are dominated by communities of Shrublands. The results of this study show the feasibility and validity of the CA_Markov based model for simulating urban land use change.

From 1989 to 2009, the city of Chihuahua grew mainly to the north and southeast directions. One of the main reasons is the increasing manufacturing industry present in those parts of the city. In these directions, the lowest elevations exist and these conditions make the terrain desirable for industrial development. In its growth stage, this industry has been settled on plain lands with access to the main roads. The most important road in Chihuahua is the one running in the north-south direction, which connects the city with the rest of the country and with the United States of America. The latter represents the main market of the manufacturing items produced in the city.

Another reason for this growth could be attributed to the increase on the number of small houses, which were constructed for people with low incomes. Many of these people work on the manufacturing industry. Thus, the location of housing projects has markedly contributed to the increase in the area occupied by Human settlements. Together, these two factors have influenced the city growth dynamics, the urban structure, its geographical expansion, and the location of the jobs generated in the city.

Before the 1970s, the number of jobs generated by the manufacturing industry was small and their location was scattered around the city. In those days, jobs were related to mining or logging activities. This scenario changed with the installation of the industrial parks called “Complejo Industrial Chihuahua”, “Las Americas” and “Saucito”. Thereafter, Human settlements had a remarkable growth, with areas where jobs related to the industry sector are concentrated [38]. The location of industrial parks near the main roads and the houses of social interest are factors that have influenced the growth of the city, as it can be visually verified on Figure 6. Other factors include the requirement of labor force, the access to highways, and the proximity of the edge of the urban area to the rural areas.

With the model of CA_Markov the growth of the urban area and the land use change on the suburban areas were assessed and quantified. Due to the absence of a conservation policy for suburban land uses, it is expected that a number of both economic and social factors cause alterations on the land uses of Grasslands, Shrublands and Riparian vegetation areas. The establishment of buffer zones could improve conditions for the use of land surrounding the city.

Even though the MC and CA methodologies have been criticized for its inability to incorporate social factors such as human decision [28], this study simulated land use changes for the years of 2009 and 2014 with a high degree of accuracy. One of the reasons for that could be the period between the dates of the images used, which was in general consistent (10-year period), compared to other studies that employed only three dates [50] or dates with varied periods among the dates of the images [51]. This gave confidence about the results from the simulations of land use for 2019 and 2024. It is possible that the estimated changes, in the absence of policy intervention, become a reality and mainly affect Grasslands and Shrublands, as indicated in the results of this study. Information on land use changes generated in this study could be useful for decision-making and for the creation of public policies focused on urban planning.

The probability matrices revealed that Grasslands were the least stable land use. This suggests that urban development will mainly occur on the plains and small slopes. Grasslands are one of the most threatened ecosystems worldwide [13]. This ecosystem possess a high degree of endemism in the region [38] and provide us with ecosystem services such as water harvest, carbon sequestration, soil retention, and contributions to weather stability, just to name a few [52]. This class lost the biggest surface area. All this area has been converted to the urban use. Urban planners in Chihuahua should take these findings into account and promote a more equilibrated growth.

Meanwhile, Shrublands have also been affected by the expansion of Human settlements due to the construction of commercial and residential buildings. This land use is distributed in lands with slopes generally greater than those of the lands where Grasslands are located. The lack of urban planning has led to a non-organized growth of Chihuahua City, with a relatively large urban area with a small population density. The increase in population, the demand for residential buildings, and the introduction of industrial parks are additional factors causing the change of the landscape.

The population of Chihuahua City has increased in the past 20 years from 530,783 to 867,910 inhabitants in the year 2014, representing an increment of 337,127 inhabitants [40]. The population growth rate for the periods 1989–1999, 1999–2009 and 2009–2014 was of 27%, 21% and 7%, respectively. These percentages mean an economic growth of the city that promotes population migration from small towns, especially the ones located nearby. People from these small towns move to the city looking for job opportunities. This produces an economic diversification demanding more labor force and space. Given the amount of territory reserves declared by the “Instituto de la Vivienda del Estado de Chihuahua” (State Institute for Housing), the city might continue growing towards the North direction unless new industrial developments occur in other directions. Lands with good characteristics for industrial development are also located south of the city.

5. Conclusions

Markov Chains and Cellular Automata, applied to remote sensing data, showed their potential as a tool for urban planning. This study established the change dynamics of seven land uses of the urban and peripheral areas of Chihuahua City. Chihuahua is experiencing a rapid urban growth regardless of the land use types of the surroundings and the urban area is becoming the main land use. In contrast, the land uses of Shrublands and Grasslands were the ones experiencing the greatest pressures from land use change.

The methodology of CA_Markov allowed describing the future behavior of the areas occupied by seven land uses in the study area. The urban growth of Chihuahua City will be mainly directed towards the North and East. Housing projects and the establishment of manufacturing industries are trigger factors for urban growth. This condition is expected to persist for over the next 10 years. The growth of the urban area indicated from this study, will cover 50% of the surface area by 2024, mainly affecting the ecosystems of Grasslands and Shrublands located nearby.

Urban planning through public policies, accompanied by projections of urban growth, could contribute to mitigate the impact over the ecosystems located nearby the City. The methods employed in this study, which identified land use transitions, represent an alternative tool for urban and territory planning. Furthermore, these results could support the elaboration of urban growth plans for Chihuahua City, Mexico, with a sustainable approach.

The model of CA_Markov has some limitations for this application. The model does not integrate socio-economic data, such as population growth, social demand, political decisions, the willingness of landowners to sell their property, or the policy changes regarding land use during the study period. It is considered that these factors notoriously influence the urban expansion. Therefore, the inclusion of these variables can improve the accuracy of the simulations; however, such variables have to be first generated in a spatiotemporal basis for the study area.