Monitoring and Predicting Spatio-Temporal Land Use/Land Cover Changes in Zaria City, Nigeria, through an Integrated Cellular Automata and Markov Chain Model (CA-Markov)

Monitoring land use/land cover (LULC) change dynamics plays a crucial role in formulating strategies and policies for the effective planning and sustainable development of rapidly growing cities. Therefore, this study sought to integrate the cellular automata and Markov chain model using remotely sensed data and geographical information system (GIS) techniques to monitor, map, and detect the spatio-temporal LULC change in Zaria city, Nigeria. Multi-temporal satellite images of 1990, 2005, and 2020 were pre-processed, geo-referenced, and mapped using the supervised maximum likelihood classification to examine the city’s historical land cover (1990–2020). Subsequently, an integrated cellular automata (CA)–Markov model was utilized to model, validate, and simulate the future LULC scenario using the land change modeler (LCM) of IDRISI-TerrSet software. The change detection results revealed an expansion in built-up areas and vegetation of 65.88% and 28.95%, respectively, resulting in barren land losing 63.06% over the last three decades. The predicted LULC maps of 2035 and 2050 indicate that these patterns of barren land changing into built-up areas and vegetation will continue over the next 30 years due to urban growth, reforestation, and development of agricultural activities. These results establish past and future LULC trends and provide crucial data useful for planning and sustainable land use management.


Introduction
Land use/land cover (LULC) changes are presently becoming issues of great concern regionally and globally due to numerous environmental challenges. The rapid increase in population and the expansion of urban areas due to urbanization has contributed immensely to the change in land use/land cover patterns and the conversion of fertile land and vegetation to land use of different functions [1][2][3]. The continuous transformation of land uses and widespread deforestation has continuously led to multiple scales of emerging environmental issues with severe effects on the ecosystem. These changes have directly and indirectly affected most urban centers [3][4][5][6]. It alters climatic conditions, biodiversity, and physical systems, serving as a significant threat to natural resource management [7]. Various studies

Remotely Sensed Satellite Data
In this study, we utilized three remotely sensed satellite images at an interval of 15 years to analyze LULC change dynamics. These images comprised Landsat TM (4) for 1990, Landsat ETM+ (7) for 2005, and Landsat OLI (8) for 2020. The Landsat image scenes were acquired from the freely accessible data portal (http://earthexplorer.usgs.gov/) of the United States Geological Survey (USGS). The downloaded images were already geo-referenced and projected to the Universal Transverse Mercator (UTM) map zone 32 N, with a datum and ellipsoid of WGS84. Therefore, Landsat TM and Landsat ETM+ satellite images were geometrically corrected using ground control points (GCPs) as used by Ogunjobi et al. [59]. The various pre-processing operations were conducted in ENVI 5.3 and ArcGIS 10.7.1 image processing software. Table 1 present the different Landsat datasets used in the study.

Remotely Sensed Satellite Data
In this study, we utilized three remotely sensed satellite images at an interval of 15 years to analyze LULC change dynamics. These images comprised Landsat TM (4) for 1990, Landsat ETM+ (7) for 2005, and Landsat OLI (8) for 2020. The Landsat image scenes were acquired from the freely accessible data portal (http://earthexplorer.usgs.gov/) of the United States Geological Survey (USGS). The downloaded images were already geo-referenced and projected to the Universal Transverse Mercator (UTM) map zone 32 N, with a datum and ellipsoid of WGS84. Therefore, Landsat TM and Landsat ETM+ satellite images were geometrically corrected using ground control points (GCPs) as used by Ogunjobi et al. [59]. The various pre-processing operations were conducted in ENVI 5.3 and ArcGIS 10.7.1 image processing software. Table 1 present the different Landsat datasets used in the study. We conducted a random field survey of the different parts of Zaria city in Nigeria to identify the city's LULC classes. The identified land cover classes were matched with similar types observed in the satellite images to interpret the different spectral signatures of the LULC class on each image. The prominent land cover types identified in both the satellite images and a field survey conducted included built-up areas that are spares and dense, un-vegetated lands with bare soils, vegetation cover of different types, and water bodies. These observations were utilized to classify and map the 4 broad land cover types: built-up areas, vegetation, water bodies, and barren land.

Method
The methodological approach employed in this study was sub-divided into 5 stages. These stages are as follows: (i) pre-processing of satellite data, (ii) image classification, (iii) accuracy assessment, (iv) land use/land cover change detection, and (v) modelling and prediction of LULC using the CA-Markov model. Figure 2 illustrates the flowchart of the methodology.  We conducted a random field survey of the different parts of Zaria city in Nigeria to identify the city's LULC classes. The identified land cover classes were matched with similar types observed in the satellite images to interpret the different spectral signatures of the LULC class on each image. The prominent land cover types identified in both the satellite images and a field survey conducted included built-up areas that are spares and dense, un-vegetated lands with bare soils, vegetation cover of different types, and water bodies. These observations were utilized to classify and map the 4 broad land cover types: built-up areas, vegetation, water bodies, and barren land.

Method
The methodological approach employed in this study was sub-divided into 5 stages. These stages are as follows: (i) pre-processing of satellite data, (ii) image classification, (iii) accuracy assessment, (iv) land use/land cover change detection, and (v) modelling and prediction of LULC using the CA-Markov model. Figure 2 illustrates the flowchart of the methodology.

Pre-Processing of Satellite Images
The pre-processing of satellite images is vital before image classification to monitor and predict land use/cover changes. Satellite images are subjected to either company or user pre-processing operations to establish the relationship between acquired data and the various biophysical phenomena [60]. Therefore, we used ENVI 5.3 software to pre-process satellite images for various operations such as geometric correction, radiometric calibration, atmospheric correction, and image enhancement. The geometric correction process involved geo-referencing images to remove geometric distortions in the acquired data and position the satellite images into a geographical coordinate system. It was performed with ground control points (GCPs) aided using aerial photographs, topographic maps, and Google Earth images of the study area. The satellite images were rectified to remove atmospheric distortions Sustainability 2020, 12, 10452 6 of 21 due to differences in orientation sensor parameters and noise from the acquiring platform. This helped considerably in reducing the misinterpretation of data during image classification. The images were further enhanced to improve the satellite images' visual quality and clarity for further analysis.

Image Classification
Image classification is one of the most effective methods of processing data from satellite imageries [61]. It helps detect, identify, and classify an image's different features on the basis of the actual classes it represents on the ground or earth surface [62]. Image classification is divided into supervised and unsupervised classification methods in remote sensing [63]. We adopted the supervised classification method using the maximum likelihood classifier (MLC) in classifying the acquired satellite images. Although several methods are devised in implementing the supervised classification, the maximum likelihood classifier, also known as the Bayesian decision rule, is one of the most widely used supervised classification algorithms due to its availability and simple training process [64]. In this method, satellite image pixels are classified on the basis of their probability of belonging to a specific LULC class. The method concept assumes equal probabilities for all classes and normal distribution for all input bands. The input bands used in the study to produce false-color composite maps consisted of bands 4, 3, 2, for Landsat 4 TM and 7 ETM+ and bands 5, 4, 3, for Landsat 8 OLI. The spectral signature of each image pixel was matched with the training samples of the study area and the satellite images classified into 4 broad land use/land cover (LULC) types. These LULC types were composed of built up/urban area, vegetation, barren land, and water bodies, as described in Table 2. The image classification procedure involved the generation of 150 training samples for each satellite image using the region of interest (ROI) tool in ENVI 5.3 image processing software. This was done to classify land use/land cover spectral classes successfully. In addition, we employed existing land use maps obtained from the city's urban planning agency and Google Earth maps to aid the creation of 50 validation samples. These samples were used to evaluate the accuracy of the image classification. According to [65][66][67], the algorithm for calculating the maximum likelihood (L i ) of unknown measurement vector (x), belonging to one of the known classes (M c ), is determined on the basis of the Bayesian equation shown in Equation (1).
where the discriminant function in the maximum likelihood algorithm is L i (x). The class is a c , i = 1, 2, 3, 4, M. M is the total class number. X is the n-dimensional pixel of a vector. n represents the number of bands. p(a c ) is the probability of the exact class at position X in a c for a pixel. The determinant of the covariance matrix of the data in class a c is |Cov c |. The inverse of the covariance matrix is Cov c and the mean vector is M c . Table 2. Land use/land cover (LULC) classification scheme.

1.
Built-up/urban areas Areas that include residential, industrial, and commercial areas, mixed-use buildings, roads, and other transport facilities.

Vegetation
It comprises agricultural and vegetable areas, crop and fallow lands, forest areas, scrubs, conifers, and other plantation of different varieties.

Barren lands
Includes areas with exposed soils, un-vegetated lands, landfill sites, and active excavation lands.

4.
Water bodies Such areas cover the city's permanent open water, rivers, streams, lakes, ponds, and various reservoirs.

Accuracy Assessment
An important process in determining image classification methods' reliability and accuracy is evaluating results [68,69]. The two broad methods of assessing a classified image's accuracy are the confusion (error) matrix technique and the receiver operating characteristic (ROC) curve. We utilized the error matrix method and considered an accuracy level above 85% to be an excellent and reliable LULC classification as recommended by previous studies [70][71][72][73]. The confusion (error) matrix technique consists of an array of numbers displayed in rows and columns showing the number of pixels or polygons attributed to a given land cover class relative to the actual class on the ground. [74,75]. It is used for obtaining descriptive and analytical statistics in classification accuracy assessment [72,74,76,77]. Producer and user accuracy were computed as suggested by [78] to obtain the average accuracies. The producer accuracy is the probability of classifying a reference pixel correctly, while user accuracy is the probability of classified pixel on each map representing the actual class on the ground or real-world location [74,76,79]. The confusion matrix comprises the overall accuracy and kappa coefficient as its two key evaluation indexes [80,81]. These indexes are computed by Kohavi and Provost [75], as presented in Equations (2) and (3). The overall accuracy (OA) is the ratio of the sum of correctly classified pixels (i.e., the total amount of major diagonal entries) to the sum of pixels in the confusion matrix. The kappa coefficient (K) is the percentage of correctly classified pixels extracted from the actual percentage expected by chance. It uses a discrete multivariate method in assessing the level of agreement or accuracy [70]. The kappa coefficient has a value that varies between −1 and 1 but usually falls between 0 and 1. A coefficient of kappa between (i) 0.00 and 0.20 signifies slight correlation, (ii) 0.21 and 0.40 signifies fair agreement, (iii) 0.41 and 0.60 signifies moderate correlation, (iv) 0.41 and 0.60 signifies substantial agreement, and (v) 0.81 to 1 indicates an almost perfect correlation [82].
Overall Accuracy (OA) = (Pc) + (Nc) where OA is the overall accuracy and percentage of correctly classified cases. Pc is the number of positive cases that are correctly identified. Nc is the number of negative cases classified correctly. Fp is the number of negative cases classified as positive incorrectly. Fn is the number of positive cases incorrectly classified as negative. P(e) is defined as the anticipated ratio by chance (i.e., the proportion of the sum of marginal probabilities multiplication per class to the total class entries).

Land Use/Land Cover Change Detection
We employed the post-classification change (PCC) detection method, previously utilized by various urban planning and environmental studies, in detecting the nature and extent of LULC change [83]. An overlay procedure was adopted using the GIS technique to acquire the land use/land cover spatial distribution and change dynamics of the study area. The overlay procedure result was a two-way cross-tabulated matrix, which showed the various changes in LULC classes. The change matrix produced in ENVI 5.3 extracted the overall quantitative land use/land cover changes and the gains and losses in each land cover type from 1990 to 2020. The cross-tabulate matrix analysis using pixel comparison facilitated the determination of the amount of change from one LULC type to others and their corresponding land cover area over the period evaluated. The change analysis was conducted to detect the LULC changes between the different periods. We monitored and analyzed the various land cover changes between 1990 and 2005, 2005 and 2020, and 1990 and 2020 defined as periods 1, 2, and 3, respectively. The analysis identified the various areas that changed from one LULC type to another in a spatial and quantitative approach during the specified period. The study presented the statistical and graphical analysis of LULC class losses and gains; the LULC transitions for periods 1, 2, and 3; and the overall contribution of other land cover classes to the built-up areas.

LULC Prediction Using the CA-Markov Model
For predicting LULC change dynamics, we utilized cellular automata and the Markov model, popularly known as the CA-Markov. The model has been efficiently used in recent studies to simulate land use changes due to its combination of the Markov model's contiguity and cellular automata's advantages. It recognizes the potential spatial distribution of transitions [84,85]. The Markov model is a stochastic model that forecasts change probability from one particular class to another, taking the LULC changes at different periods into consideration [38,39]. It has random variation whose probabilities depend on a period's historical data using time-series techniques [86]. The outcome of the Markov model relies on transition probability [87]. The transition probability (P ij ) between state (j) and (i) is the probability in which land cover (i in pixels) in time (x) changes to land cover class (j) in time (x + 1). In this model, the change dynamics for any study area depends on the previous or current land cover condition [88] and mathematically computed as adopted by [36] using Equations (4) and (5).
The transition probabilities are obtained from the transition samples occurring during a specific time interval [87] and presented in the transition matrix (P).
where L (x+1) and L (x) are the land use/land cover conditions at time (x + 1) and (x), respectively. 0 ≤ P ij < 1 and m j=1 P ij = 1, (i, j = 1, 2, 3, . . . , m) is the transition probability matrix. We performed a Markov chain analysis to generate the transition matrix of the LULC change and the probabilities of change from 1990 to 2005, 2005 to 2020, and 1990 to 2020. The transition matrix serves as the basis for projecting future LULC change dynamics. However, the Markov chain model is not explicit spatially. It lacks the scientific explanation of the processes of change and neglects the spatial distribution of LULC, which is exceptionally significant in simulating land cover patterns [89]. The cellular automata (CA) model is also widely used for LULC prediction due to its spatial capacity to alter and control processes of complex distributed systems. The CA model comprises the cell, cell space, neighbor, time, and rule. The model describes the new pattern of LULC, considering the state of previous neighborhood cells [39,86]. The distance between the neighbor and the cell defines the weight factor of changing to a particular land cover. The weight factor of the change to a specific land cover is determined by the distance between the cell and the neighbor. The weight factor was then combined with the transition probabilities to estimate neighborhood cells' condition so that change prediction is not only based on a random decision. The cellular automata (CA) model considers Markov's previous or current state of LULC and utilizes the neighborhood cells' condition for its transition rules [87]. It provides a suitable environment for dynamic modelling in GIS and remote sensing due to its analytical engine [85]. However, the CA model has limitations in defining transition rules and modelling structures. Therefore, the combination of different empirical and dynamic models, such as the CA-Markov model, is vital in achieving a dynamic LULC spatial modelling [90].
The study used the land change modeler (LCM) in IDRISI-TerrSet geospatial monitoring and modelling software to monitor and simulate the potential LULC change dynamics from 2020 to 2050. The procedures used to run the CA-Markov in LCM comprised the following procedures. The first stage involved running three different models using the land covers maps of 1990-2005, 2005-2020, and 1990-2020. This procedure was used to generate the transition probability matrix (TPM) and transition area matrix (TAM), and transition suitability maps (TSM). The TPM is obtained by cross-tabulating two multi-temporal images. It contains the probability of each land cover class changing to another over a pre-determined time. The TAM contains the estimated number of pixels that might change over a specific number of time units from a particular land cover class to another, while the conditional probability map expresses the probability of each land cover class pixel having a place with the assigned class after a specific time [91]. The transition suitability maps (TSM) were generated for the 4 land cover classes used in the study (i.e., built up, vegetation, water bodies, and barren soil).
The model of different historical LULC scenarios of 1990-2005 and 2005-2020 was used to produce the transition probability matrix of periods 1 and 2. The second stage involved using a standard contiguity filter of 5 × 5 to define each cell's neighborhoods and generate the spatially explicit weighing factors. After calibrating the model, we used the scenario-bound approach to simulate the potential LULC pattern in the final stage. The process involved using the classified land cover maps of 1990 and 2005 to calibrate and optimize the Markov chain algorithm. The earliest year (i.e., 1990) was used as time 1, while the later year (i.e., 2005) was used as time 2. The transition probabilities between time 1 and 2 were used to simulate the LULC pattern in 2020. The study validated the CA-Markov model to determine the accuracy of the 2020 prediction. The process requires a statistical tool to differentiate between location errors and quantities' errors in analyzing the similarities between 2 images [92]. We employed the kappa statistics index to determine the level of agreement between the projected and the actual land cover map of 2020. The result indicated an acceptable kappa index, signifying a reliable LULC modelling and prediction. Therefore, the 2020 classified land cover map was used as a base map to forecast the potential LULC in 2035 and 2050, using the transition probabilities of 2005-2020.

Land Use/Land Cover Change
The LULC change analysis results illustrate the historical land use/land cover changes and the change dynamics in each land cover class. The LULC results would significantly help urban planners and decision-makers determine the nature, position, and rate of change in the study area. The classified land use/land cover maps of 1990, 2005, and 2020 were generated using the supervised maximum likelihood classification. The values of the overall accuracies and kappa coefficients for all the classified LULC maps of 1990, 2005, and 2020 were found to be above 80% and 0.75, as shown in Table 3. This indicates a reliable and accurate classification of images for analyzing land use/land cover change. It was observed that the land use/land cover pattern of Zaria city has changed dramatically between 1990 and 2020, with the built-up areas experiencing continuous growth over the last three decades. The city's vegetation cover has also increased, while barren lands and water bodies decreased during the period. Figure 3 displays the spatial distribution of each LULC in 1990, 2005, and 2020. The three maps depict the various stages each land cover class has undergone, with each stage having a different rate and quantity of change. The quantitative result is presented in Table 4, which shows the statistics of the four LULC classes and the change dynamics during the three-time nodes.
decades. The city's vegetation cover has also increased, while barren lands and water bodies decreased during the period. Figure 3 displays the spatial distribution of each LULC in 1990, 2005, and 2020. The three maps depict the various stages each land cover class has undergone, with each stage having a different rate and quantity of change. The quantitative result is presented in Table 4, which shows the statistics of the four LULC classes and the change dynamics during the three-time nodes.    The result shows steady growth in the city's built-up area and significant changes in other LULC classes between 1990 and 2010. It indicates an increase in built-up area from 1857.51 ha in 1990 to 2784.24 ha in 2005. This is an increase from 4.72% to 7.08% of the total area during the first period (1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005). The built area further increased to 5444.55 ha (13.84% of the total area) in 2020, indicating a growth of 2.36% and 6.76% during the first and second periods, respectively. The result reveals an increase from 4.72% of the total area in 1990 to 13.84% in 2020, which indicates the expansion of the built-up area with about 9.12% between 1990 and 2020. The city's vegetation cover and barren land witnessed an inconsistent change during the study period. The vegetation decreased from 12,937.77 ha in 1990 to 11,606.94 ha in 2005 and increased rapidly to 18,209.25 ha in 2020. This indicates an increase from 32.89% of the total area in 1990 to 46.29% in 2020. The growth in Zaria city's built-up area and vegetation cover can be attributed to the city's rapid increase in population due to rural-urban migration and the increased demand for agricultural produce that contributes to the rise in cash crop farming. The barren land increased from 22,298.13 ha in 1990 to 23,152.77 ha in 2005, and decreased drastically to 13,674.87 ha in 2020. The result indicates a decrease from 56.68% of the total area in 1990 to 34.76% in 2020. The water bodies in Zaria city did not witness much change during the entire study period; it decreased from 2247.03 ha (5.71% of the total area) in 1990 to 2011.77 ha (5.11% of the total area) in 2020. The historical land use/land cover results in 1990, 2005, and 2020 are illustrated graphically in Figure 4. Change analysis was further conducted using the LULC maps of 1990, 2005, and 2020 to demonstrate the losses and gains in each land cover class during the three-time nodes.
The change analysis result clearly shows both losses and gains in all LULC classes during the different study periods.
The result shows steady growth in the city's built-up area and significant changes in other LULC classes between 1990 and 2010. It indicates an increase in built-up area from 1857.51 ha in 1990 to 2784.24 ha in 2005. This is an increase from 4.72% to 7.08% of the total area during the first period (1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005). The built area further increased to 5444.55 ha (13.84% of the total area) in 2020, indicating a growth of 2.36% and 6.76% during the first and second periods, respectively. The result reveals an increase from 4.72% of the total area in 1990 to 13.84% in 2020, which indicates the expansion of the built-up area with about 9.12% between 1990 and 2020. The city's vegetation cover and barren land witnessed an inconsistent change during the study period. The vegetation decreased from 12,937.77 ha in 1990 to 11,606.94 ha in 2005 and increased rapidly to 18,209.25 ha in 2020. This indicates an increase from 32.89% of the total area in 1990 to 46.29% in 2020. The growth in Zaria city's built-up area and vegetation cover can be attributed to the city's rapid increase in population due to ruralurban migration and the increased demand for agricultural produce that contributes to the rise in cash crop farming. The barren land increased from 22,298.13 ha in 1990 to 23,152.77 ha in 2005, and decreased drastically to 13,674.87 ha in 2020. The result indicates a decrease from 56.68% of the total area in 1990 to 34.76% in 2020. The water bodies in Zaria city did not witness much change during the entire study period; it decreased from 2247.03 ha (5.71% of the total area) in 1990 to 2011.77 ha (5.11% of the total area) in 2020. The historical land use/land cover results in 1990, 2005, and 2020 are illustrated graphically in Figure 4. Change analysis was further conducted using the LULC maps of 1990, 2005, and 2020 to demonstrate the losses and gains in each land cover class during the threetime nodes. The change analysis result clearly shows both losses and gains in all LULC classes during the different study periods.  Vegetation cover experienced the most significant loss during this period. This decline can be mainly linked to the city's urban expansion and development due to population growth. Other factors contributing to the reduction in vegetation cover between 1990 and 2005 include inadequate legislation, institutional weakness, unaccountability, and uncontrolled human activities. The city's inhabitants often engaged in forest encroachment and destruction by removing vegetation cover for agricultural purposes and cutting down forest trees for firewood used in cooking. This served as a major threat to the city's natural habitats and contributed immensely to the vegetation cover loss. The result of this period aligns with the study on urban sprawl pattern and its implications for urban management, which revealed Zaria city's inhabitants to often engage in forest encroachment and destruction by removing vegetation for agricultural purposes and cutting down forest trees for firewood used in cooking [54]. This contributes immensely to the loss of vegetation cover and serves as a threat to the city's environmental sustainability.
During This result is in contrast with recent studies that indicate a continuous decrease in vegetation cover due to urban growth and development [3,18,19]. In this case, reforestation and massive investment in Fadama agricultural projects by federal, state, and various non-governmental organizations in Zaria city has contributed to the significant increase in agricultural/vegetable areas and the consequent decrease in barren land.
During the third phase between 1990 and 2020 (i.e.,  Table 5 and Figure 5. The spatial changes in the four LULC classes and the unchanged areas during the study period are shown in Figures 6 and 7.

Markov Transition Matrix Analysis
We produced the change dynamics of the four LULC classes using a Markov chain analysis as presented in the transition probability matrix (TPM), computed for the periods between 1990 and 2005, and 2005 and 2020. The TPM presented in Table 6 has values that range between 0 and 1. The higher the value, the more likelihood of changing from the horizontal land cover class to the vertical and vice versa. The diagonal values in each period's transition matrix represent the likelihood of each land use/land cover class persisting. In contrast, the off-diagonal values represent the likelihood of a land cover class changing to another over the period. The transition probability matrix presents the potential LULC change trajectories for periods 1 and 2. It can be observed from the matrix that the future change likelihood of barren land to built-up areas increased from 3.1% in 1990-2005 to 12.1% in 2005-2020. The likelihood of water bodies changing to built-up areas increased reasonably from 8.0% in 1990-2005 to 11.6% in 2005-2020. This increase in the possibility of barren land and water bodies changing to built-up areas can be ascribed to the city's historical setting and rapid population growth, mainly due to agricultural activities. The rapid increase in the city's population also led to an increased demand for urban infrastructures and services, which contributes to the expansion of built-up areas, as similarly observed in the growth and development of prominent cities such as Lagos, Uyo, and Abuja in Nigeria [17,51,52].
The Markov transition matrix result indicates the change probability of barren land to vegetation cover as 15.8% from 1990 to 2005, while the probability of changing from water bodies to vegetation was 36.2%. The change probability from barren land to vegetation increased to 34.4% for the period between 2005 and 2020, while water bodies to vegetation increased to 64.0% within the same period. This remarkable increase in the likelihood of barren land and water bodies changing to vegetation may have been unconnected to agricultural activities' rapid growth and various irrigational schemes. The probability of vegetation cover changing to the built-up area also increased over the last three decades, from 3.4% to 4.1%. The change probabilities between the two periods suggest rapid urban growth and a gradual decrease in agricultural lands due to urbanization. Therefore, the modelling of future scenarios is required to monitor and analyze LULC trends in 2035 and 2050.

Modelling and Predicting Land Use/Land Cover Changes
We employed the integrated CA-Markov model and satellite-derived land cover maps to model the future spatial and quantitative LULC scenarios in Zaria city, Nigeria. However, the model's validation is necessary to evaluate the simulation process's consistency in terms of quantity and location [93]. The transition area files for the years 1990-2005 were used to model the 2020 LULC map. The simulated 2020 land cover map was then compared with the satellite-derived LULC map, which served as a reference map representing the city's actual ground conditions. This procedure was used to evaluate the validity of the CA-Markov model. The validation module comprises four statistical kappa parameters: the overall kappa, kappa location, kappa histogram, and accuracy percentage [94]. The overall kappa was utilized in evaluating the model's overall success, while kappa location was utilized in evaluating the model's ability to identify location accurately. These values are essential in evaluating the overall accuracy of the simulation model. A kappa statistic of 0% indicates no agreement between land cover maps and their probability, while 100% indicates a perfect agreement. The validation statistics present a kappa no (K no ) value of 0.7698, kappa location (K location ) value of 0.8336, kappa locationStrata (K locationStrat ) value of 0.8336, and kappa standard (K standard ) value of 0.6997. The results indicate good kappa index values, verifying the model's efficiency for future land cover simulations [95]. Therefore, we used the classified 2020 LULC map alongside the transition probability matrix and the transition area matrix between 2005 and 2020 to forecast the city's land cover pattern in 2035 and 2050, as shown in Figure 8. The corresponding areas of the modelled land cover classes are presented in Table 7.  [95]. Therefore, we used the classified 2020 LULC map alongside the transition probability matrix and the transition area matrix between 2005 and 2020 to forecast the city's land cover pattern in 2035 and 2050, as shown in Figure 8. The corresponding areas of the modelled land cover classes are presented in Table 7.     13.05% between 2020 and 2050. Water bodies will increase slightly from 2011.77 ha in 2020 to approximately 2559.87 ha in 2050. This result suggests that the increase in built-up areas and the continuous reforestation programs will contribute to the expected decrease in non-vegetated areas, exposed soils, landfills, and excavations sites. The city's expansion of the built-up areas can be associated with various educational, socio-economic, and agricultural activities [53,54]. These activities have contributed to different change dynamics observed in the city's land cover.
This forecast aligns with previous studies in India, which revealed socio-economic and biophysical factors to have significantly influenced the projected growth and development of built-up areas in Warangal and Allahabad Urban district [19,35]. The availability of infrastructural, educational, and medical facilities in these districts will significantly increase the urban population and alter the city's land uses, affecting the natural environment. Similarly, the city of Ravansar in Iran is further expected to experience an unsustainable natural environment changing into built-up areas in the future due to excessive human activities [3]. Studies in China highlight rapid urbanization as a driving factor for its extensive construction activities [18]. This development will lead to a future reduction in the country's agricultural land and counter the state's ability to achieve sustainable urban expansion unless the government achieves a reduced urban construction.
Therefore, in order to balance the expansion of urban areas and conserve the natural environment in Zaria city, certain policies that necessitate rapid urban growth and expansion without providing a harmonious relationship between human activities and the ecological environment need to be amended. This will help urban planners and other local authorities effectively manage land uses and maintain their complex change dynamics.

Conclusions
The study monitored and predicted the spatio-temporal land use/land cover change in Zaria city through an integrated CA-Markov that combined the advantages of cellular automata and the Markov chain model. Multi-temporal satellite data and GIS techniques were employed to monitor the city's historical land use/land cover pattern using the land cover maps of 1990, 2005, and 2020. In the study, the land cover maps were classified into four major LULC classes (built-up area, vegetation, barren land, and water bodies). This was to comprehensively understand the city's LULC changes and predict the future scenario. First, we simulated the LULC map of 2020 and evaluated its accuracy using kappa statistical parameters. Its validation showed a strong correlation between the simulated land cover map and satellite-derived map, which proved the simulation model's reliability. We utilized the calibrated CA-Markov model to predict the future LULC change dynamics in 2035 and 2050. The results presented in this study showed noticeable changes in the spatial and quantitative distribution of land use/land cover. It revealed the city's built-up area and vegetation to have grown continuously between 1990 and 2020, while barren land decreased significantly during the study period. The outcome of the LULC change prediction suggests a continuous increase in the extent of built-up and vegetation at 2.60% and 5.81%, respectively, from 2020-2035, while barren land will decrease by 9.35% and 13.05% from 2020 to 2035 and 2020 to 2050, respectively. The findings of the study revealed the historical land use/land cover changes and highlight the future change scenarios in the next 30 years, relying on the continuity of the historical LULC trends. The study successfully demonstrated the CA-Markov model's efficiency, using remotely sensed data and GIS techniques to monitor spatio-temporal LULC changes and predict future land cover patterns. Such data are vital for informed decision-making in urban planning, providing the potential information required to monitor growth and improve environmental sustainability. Although the present study effectively predicted the expected LULC of 2035 and 2050, it lacks up-to-date census data on the city's population growth. Therefore, further research needs to examine more closely the relationship between land cover changes and population growth. This is to utilize demographical data in having a holistic understanding of the long-term dynamics of LULC change and its implications for achieving sustainable development.

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