Land Use / Land Cover Dynamics and Modeling of Urban Land Expansion by the Integration of Cellular Automata and Markov Chain

This study explored the past and present land-use/land-cover (LULC) changes and urban expansion pattern for the cities of the Kathmandu valley and their surroundings using Landsat satellite images from 1988 to 2016. For a better analysis, LULC change information was grouped into seven time-periods (1988–1992, 1992–1996, 1996–2000, 2000–2004, 2004–2008, 2008–2013, and 2013–2016). The classification was conducted using the support vector machines (SVM) technique. A hybrid simulation model that combined the Markov-Chain and Cellular Automata (MC-CA) was used to predict the future urban sprawl existing by 2024 and 2032. Research analysis explored the significant expansion in urban cover which was manifested at the cost of cultivated land. The urban area totaled 40.53 km2 in 1988, which increased to 144.35 km2 in 2016 with an average annual growth rate of 9.15%, an overall increase of 346.85%. Cultivated land was the most affected land-use from this expansion. A total of 91% to 98% of the expanded urban area was sourced from cultivated land alone. Future urban sprawl is likely to continue, which will be outweighed by the loss of cultivated land as in the previous decades. The urban area will be expanded to 200 km2 and 238 km2 and cultivated land will decline to 587 km2 and 555 km2 by 2024 and 2032. Currently, urban expansion is occurring towards the west and south directions; however, future urban growth is expected to rise in the southern and eastern part of the study area, dismantling the equilibrium of environmental and anthropogenic avenues. Since the study area is a cultural landscape and UNESCO heritage site, balance must be found not only in developing a city, but also in preserving the natural environment and maintaining cultural artifacts.


Introduction
Land-use/land-cover (LULC) classification, generally based on remote sensing images, is one of the major research themes in terms of global environmental change [1].LULC change is the prompt influence and response of human activities upon nature followed by significant consequences [2][3][4][5].Anthropogenic factors play a vital role in land-use changes of any area [6,7] and LULC is commonly based on urban growth [8].Urbanization is defined as a process which incorporates increased population, modernization, and the imperative socioeconomic phenomenon [9,10].It is a continuous process that is and has become a global trend under the combined influence of anthropogenic activities and biological factors [11].Accelerated urban expansion not only influences socioeconomic change, but also influences farm land loss [7,12], impacts ecology and the environment [13], and often threatens sustainable urban development [14,15].Urban expansion is rife in the cities of developing countries and discrete and inconsistent with the local plans and policies.Scientific analysis in these areas is highly essential to better understand the growth patterns and processes.Generally, cities in developing countries surpass the coping capacities, resulting in squatter settlements and shanty townscapes.Therefore, the precise projection of future urban growth and its management based on reliable statistics and a complete understanding of the patterns and urbanization trends is essential [16].
Accurate, consistent, and updated information on the urbanization trend is crucial for the need analysis and policy formulation to ensure a sustainable urban future.Over time, land-use change maps provide essential information for land-use planning [17,18] that can help to understand the drivers and dynamics of land-cover transformation and predict the future economic and environmental influences [3].GIS and remote sensing are the appropriate tools for land-cover monitoring, urban/regional planning [19,20], and exploring spatiotemporal changes of LULC from a local to a global scale [8].These are now capable of real time evaluation on mobile devices [20,21].
Spatiotemporal analysis of urban expansion is of paramount significance as detailed information on spatial location, characteristics, and consequences of urban growth are the fundaments for the (i) formulation of urban development plans; (ii) development and modification of the theories of urban morphology; and (iii) defining the boundary between the urban area and environment for some environment models [21,22].Similarly, amongst the various geospatial and statistical models, urban development models tend to predict the events and actions and their effects in a specific time and space [23][24][25].Urban growth is intrinsically associated with the change in LULC and leaves its abrupt consequences upon the overall ecology.Hence, understanding urban growth is a fundamental prerequisite for a future change simulation [26].
LULC models consist of multiple methodological approaches, can be static and spatial, and explore the change vs. rates of change [27].Land-cover change models generally comprise a change demand sub-model, transition potential sub-model, and change allocation sub-model, which determine the amount and spatiotemporal location of LULC change and the changing land-cover type from one to another [3].Estimating urban growth, analyzing spatiotemporal pattern, and assessing the characteristics and consequences of urbanization have been more possible due to the utilization of urban growth models [17,28].An urban growth simulation model (UGSM) was developed post 1950s [29] under the assumption that the dual process towards the development of agglomerations and dispersions was well-balanced.These gained a wide range of attention after the 1990s, after the acquisition of computing ability.LULC models are intended to amply support urban development planning and sustainable growth management [30] and are effectively used due to the capability to make effective predictions about actions and their effects in a particular time and space [31].
The Markov chain is a popular model which calculates future change based on the past, and Cellular Automata (CA) detects the spatial location of change [22,32,33].Integrated CA-Markov is recognized as a more powerful and effective modeling technique for the change simulation [34][35][36] of any area.It is highly applicable to predict the status of an rapidly urbanizing area [35] and identify the conversion from one state to another at each time period [1].It has been widely utilized for the study of the urbanization of cities in Nepal [20,33] and abroad, (i.e., Wuhan, China [36]; Hua Hin, Thailand [34]; and Saga, Japan [1]; Setúbal and Sesimbra, Portugal [37]; the central part of Germany [32]; London, UK [11]; Ahmedabad, India [38]; the Tehran metropolitan area, Iran [22]; the Santiago metropolitan area, Chile [39]; and the Foshan metropolitan city, China [40]).Although various additional models have been applied to simulate the urbanization process and LULC change including SLEUTH [41], DINAMICA [42], CLUE [43,44], SERGoM [45], and LUCAS [46], the present study selected the integrated CA-MC model to simulate the spatio-temporal dynamics of the study area.
Nepal is a rapidly urbanizing country in Southeast Asia [47], and is expected to become significantly more urbanized by 2050 [48].Urban population accounted for 2.9% in 1952/54 that increased to 3.6% in 1961, 4% in 1971, and gradually increased to 6.4%, 9.2%, 13.9%, and 17.1% in 1981, 1991, 2001, and 2011, respectively, while urban centers increased from 10 to 58 over time [49].The Government declared new urban centers in 2015 and urban centers totaled 217; whereas the urban population accounted for 40.49% [50].By 2017, the country acquired 292 urban centers [51], with more than 50% of the total population inhabiting the urban centers [49].Despite the pattern being rapid and discrete, the consequence of urbanization favors economic development, modernization, and social advancement of the country, contributing about 33.1% in national domestic product (GDP) [52].
The Kathmandu valley, the capital of the nation, is an economic, sociocultural, political, and educational hub with extant fertile land for agriculture [53].It is one of the fastest growing urban regions of South Asia [47] and the urbanization of the Kathmandu valley has gained complete momentum since the late 1950s [54].The accelerated pace of urbanization has explicitly contributed to the substantial change in LULC of the region [25,50,55,56], mainly the exponential increase in urban/built up areas and the sharp decline in cultivated land.Therefore, urban dynamics and LULC trajectories of the valley have been the subject of various scientific researches.Haack and Rafter [53] explored the LULC dimension between 1978 and 2000; Haack and Rafter [57] compared the maps of 1955 and 2000 to detect the change; Thapa and Murayama explored the overall LULC scenario between 1967 and 2000 [25] and simulated the urban scenario of Kathmandu Valley by 2050 [20].Rimal et al. [50] observed the urban expansion from 1976 to 2015 in relation to the earthquake hazard.Ishtiaque et al. [55] monitored the LULC dynamics of the valley from 1989 to 2016.Pradhan and Parera [58] examined the LULC dynamics and resolved the unsystematic urban sprawl at the expense of prime farm land.All these previous studies have explored the rapid, unplanned, and unmanaged urban sprawl, high population concentration, LULC change and complexities in the valley and the trend is likely to be more intense in the future.Furthermore, previous studies [25,50,55] employed maximum likelihood algorithm for LULC classification including the three districts of Kathmandu valley.The neighboring district, Kabhrepalanchok has similar socioeconomic activities and features of urbanization.It is not only the new destination of possible urban sprawl, but also the core region in terms of cultivated land and forest resources.Hence, the region still contains great opportunities to develop planned and managed urban settlements unlike the Kathmandu valley.However, previous studies have paid less attention to addressing future urban expansion simulation of adjoining areas as well as the expansion direction towards the contiguous district.
As seen in earlier studies from Kathmandu and its fringes, urban planning and management becomes urgent right before the situation irreversibly deteriorates.Hence, there appears to be a strong need to investigate the spatiotemporal dynamics of LULC to support sustainable urban planning.This study aims to identify the LULC change from 1988 to 2016 by using the support vector machines (SVM) procedure; analyze spatiotemporal dynamics of urban expansion and orientation; and simulate the LULC and urban area expansion by the years 2024 and 2032, using Markov Chain and Cellular Automata (MC-CA) models.Specifically, the following two questions are addressed in this study: (1) What was the trend of LULC change within the study area in the past , and which LULC classes were mostly affected by urbanization and where?(2) What growth and change patterns can be expected in the future?

Study Area
The study area is located in the central east part of Nepal.Administratively, it incorporates 24 cities of the Kathmandu valley, (Kathmandu: 11 cities, Lalitpur: three cities, and Bhaktapur: four cities) and the Kabhrepalanchok (6 cities) district (CKVAKD) from four districts covering a total area of 1215.23 km 2 .It is geographically enclosed between the northern latitudes of 27 • 31 40" and 27 • 49 10" and the eastern longitudes of 85 • 11 18" and 85 • 43 44" (Figure 1).The elevation of the area ranges from 546 to 2825 m above sea level (masl).Population and economic growth are the two major underlying elements that play vital roles in urban land expansion [17,19,33,38,[59][60][61][62][63][64] of the region.Economic opportunities in the city's core areas, population growth in the fringes, and the political situation in the rural areas have mainly contributed to the urbanization of the Kathmandu valley [20,25,65].Furthermore, rural urban migration, economic centrality, sociopolitical factors, and real estate boom are the drivers of the valley area [55].According to Central Bureau of Statistics (CBS), 2011, the total population of these four districts was 1,431,699 in 1991, which increased to 2,032,764 in 2001, and 2,900,971 in 2011 [49].

Study Area
The study area is located in the central east part of Nepal.Administratively, it incorporates 24 cities of the Kathmandu valley, (Kathmandu: 11 cities, Lalitpur: three cities, and Bhaktapur: four cities) and the Kabhrepalanchok (6 cities) district (CKVAKD) from four districts covering a total area of 1215.23 km 2 .It is geographically enclosed between the northern latitudes of 27°31′40″ and 27°49′10″ and the eastern longitudes of 85°11′18″ and 85°43′44″ (Figure 1).The elevation of the area ranges from 546 to 2825 m above sea level (masl).Population and economic growth are the two major underlying elements that play vital roles in urban land expansion [17,19,33,38,[59][60][61][62][63][64] of the region.Economic opportunities in the city's core areas, population growth in the fringes, and the political situation in the rural areas have mainly contributed to the urbanization of the Kathmandu valley [20,25,65].Furthermore, rural urban migration, economic centrality, sociopolitical factors, and real estate boom are the drivers of the valley area [55].According to Central Bureau of Statistics (CBS), 2011, the total population of these four districts was 1,431,699 in 1991, which increased to 2,032,764 in 2001, and 2,900,971 in 2011 [49].

Data
In this study, freely available maximum cloud-free time-series Landsat images for the years 1988,1992,1996,2000,2004,2008,2013, and 2016 (TM, ETM, and OLI) were obtained from the United State Geological Survey (USGS) website https://earthexplorer.usgs.gov(Table 1) [66].All the acquired images were post-monsoonal from 23 October to 3 April.In addition, all the images and Shuttle Rader Topographical Mission (SRTM) digital elevation model (DEM) with 30 m resolution were projected in UTM Zone 45N.Furthermore, high spatial resolution images from Google Earth (http://earth.google.com) of multiple dates and topographical maps published by Survey Department, Government of Nepal, 1995 (scale 1:25,000) [67] were collected.The administrative boundary vector data at the district and municipal levels were acquired from the Survey Department of Nepal, and district level population information of 1991, 2001, and 2011 from the Central Bureau of Statistics (CBS) [49].Major land-cover change information and road networks were gathered using Global Positioning System (GPS) in course of the field visit and field verification campaigns during June 2015 and July 2016.TM, ETM+, and OLI Landsat images were verified for geometric accuracy and all images were processed in an ENVI environment [33].These multispectral images were atmospherically and geometrically corrected in course of pre-processing [68,69].In this study, L1T (terrain corrected) images were converted from digital number (DN) to radiance, and image processing functions were conducted applying the radiometric calibration model.ENVI included radiance calibration, geometric correction, and atmospheric correction.The FLAASH atmospheric correction model was employed using radiance image by applying the appropriate model on the location of the study area.The positional root mean square (RMS) error of geometric rectification was not more than 0.5 pixels.A supervised approach support vector machines (SVM) classifier [18,70] was employed for the classification of LULC and change analysis.Several previous authors have evaluated the performance of SVM algorithms and concluded that it is one of the flexible supervised classifier options with high accuracy [68,[70][71][72][73][74].SVM is a supervised, non-linear, non-parametric classification technique that is widely used in the remote sensing field due to its specific capability to draw conclusions even with limited training samples [75].Previously, Lee [76] evaluated the SVM algorithms for landslide vulnerability mapping in the Gangwon province of Korea and the results suggested that SVM was quite suitable for a wide range of classification problems, even if the problems were of high dimension and were non-linear separable.Moreover, the approach was applied by Mahmoud et al. [71] to monitor urbanization in Abuja, Nigeria.Schneider [70] evaluated the maximum likelihood classifier (ML), decision tree (DT), and SVM algorithms for monitoring land-cover change in urban and peri-urban areas using Landsat satellite data.This research illustrated that the overall accuracy result, the DT and SVM performed better than the ML classifier.Waske and Benediktsson [77] compared the SVM classification outputs with other algorithms like DT and ML.Based on the results, SVMs achieved the highest overall accuracies of all approaches when classifying the multispectral data.Nevertheless, SVMs formulations are not completely free from drawbacks.The major problem concerns the selection of kernels [76,78].Despite this limitation, SVMs are still more popular and produce more accurate classification results than the traditional methods [76,78,79].SVMs can be generally grouped into four kernel functions-namely, linear, polynomial, radial basic function, and sigmoid [78].The mathematical representations of each kernel are mentioned in the help section of ENVI version 5.3; whereas, the same mathematical representations are presented by Lee et al. [76] using ENVI version 4.4.
In this current study, a communally used radial basis function (RBF) kernel was recognized during image classification [18].The RBF kernel is frequently used for its advantage in the classification of remotely sensed images to achieve better results than other kernels (polynomial, linear, and sigmoid) [68,80].The penalty parameter maximum value 100 was assigned and the land-cover classification scheme was adopted as recommended by Anderson et al. [81] and Thapa and Murayama [25].After the extensive field survey, six main land-use classes were identified: urban/built-up, cultivated land, vegetation cover (dense forest, scattered forest, shrub, and grass), sand area, water body, and open field (Table 2).

Accuracy Assessment
In this study, overall accuracy (OA), user's accuracy (UA), and producer's accuracies (PA) were observed based on field reference data; topographical maps published by Survey Department, Government of Nepal, 1995 (scale 1:25,000) [67]; and high spatial resolution images of Google Earth (http://earth.google.com)from multiple dates.A total of 1200 (i.e., 200 for each class) stratified random sample points were selected for the validation of classified images for each year.

Measuring Urban/Built up Area Expansion Rate
To calculate the urban expansion growth rate of the study area, the total transformation of the urban area was taken into consideration [50].The urban expansion rate refers to the average annual urban area growth in the following years.
where BER measures the urban expansion rate in (km 2 /year); B 1 , B 2 represents the urban area (km 2 ) in the year; and T 1 , T 2 indicates the time.

Quantification of LULC Based Transition Analysis
The cross function has been applied to create the transition matrix of LULC using TerrSet software developed by Clark Lab (clarklabs.org).The transition matrix consisted of rows and columns of landscape categories at time T 1 and T 2 [50].In this study, the Land Change Modeler (LCM) system of TerrSet software was applied to create the transition of a LULC map for seven time periods : 1988 and 1992; 1992 and 1996; 1996 and 2000; 2000 and 2004; 2004 and 2008; 2008 and 2013; and 2013 and 2016.

Urban Expansion and Orientation
Despite its wide range of significance, few studies [82,83] have quantified the landscape structure using landscape metrics and explored the orientation of urban expansion [50,84,85].Cushman et al. [83] applied the FRAGSTATS software and principal component analysis (PCA) for the evaluation of landscape structure, and cluster analysis was conducted.The same software package FRAGSTATS was employed to detect the gradient of landscape pattern along two directions: W-E and SW-NE [82].To analyze the spatial changes of urban expansion, Singadurbar of Kathmandu, Lagankhel of Lalitpur, Kamalbinayak of Bhaktapur, and Dhulikhel of Kavrepalanchok districts were assumed (Figure 1) as centers.Identification of the assumed city centers of each district were based on the field visit experience of the authors and classified satellite images from 1988.The outcomes should answer the questions of where and in which direction the changes mostly occurred from the city core [50].District level urban orientation was observed through the linings created by utilizing ArcGIS 10.1.The waves [84] consisted of eight subdivisions: North, North-East, East, South-East, South, South-West, West, North-West, North (N-NE, NE-E, E-SE, SE-S, S-SW, SW-W, W-NW, and WN-N), each representing 45 • from the central point (Figure 1).

Simulation of LULC Change
CA-Markov is an appropriate technique for modeling LULC change in places where understanding and describing landscape relations are difficult.This model allows the future status of an ecosystem based on its pre-existing status to be predicted.The CA-Markov model has been widely used to effectively understand and measure urban expansion [33] and landscape dynamics [86].Typically, the implementation process of the CA-Markov model consists of: ( 1 the evaluation of landscape structure, and cluster analysis was conducted.The same software package FRAGSTATS was employed to detect the gradient of landscape pattern along two directions: W-E and SW-NE [82].To analyze the spatial changes of urban expansion, Singadurbar of Kathmandu, Lagankhel of Lalitpur, Kamalbinayak of Bhaktapur, and Dhulikhel of Kavrepalanchok districts were assumed (Figure 1) as centers.Identification of the assumed city centers of each district were based on the field visit experience of the authors and classified satellite images from 1988.The outcomes should answer the questions of where and in which direction the changes mostly occurred from the city core [50].District level urban orientation was observed through the linings created by utilizing ArcGIS 10.1.The waves [84] consisted of eight subdivisions: North, North-East, East, South-East, South, South-West, West, North-West, North (N-NE, NE-E, E-SE, SE-S, S-SW, SW-W, W-NW, and WN-N), each representing 45° from the central point (Figure 1).

Simulation of LULC Change
CA-Markov is an appropriate technique for modeling LULC change in places where understanding and describing landscape relations are difficult.This model allows the future status of an ecosystem based on its pre-existing status to be predicted.The CA-Markov model has been widely used to effectively understand and measure urban expansion [33] and landscape dynamics [86].Typically, the implementation process of the CA-Markov model consists of: ( 1  To produce transition potential maps, multi-criteria evaluation (MCE), analytic hierarchy process (AHP), and fuzzy membership function were applied.The AHP method was applied to determine the weights of driving factors with the use of pairwise evaluations [87].The AHP method allows weighting of land-cover transition potential on the basis of a set of potential maps (e.g., magnitude of slope), and incorporates growth constraints.This GIS-based AHP is a strong tool because of its high ability to incorporate different types of heterogeneous variables and its simplicity To produce transition potential maps, multi-criteria evaluation (MCE), analytic hierarchy process (AHP), and fuzzy membership function were applied.The AHP method was applied to determine the weights of driving factors with the use of pairwise evaluations [87].The AHP method allows weighting of land-cover transition potential on the basis of a set of potential maps (e.g., magnitude of slope), and incorporates growth constraints.This GIS-based AHP is a strong tool because of its high ability to incorporate different types of heterogeneous variables and its simplicity to gain the weights of suitable variables [88,89].The transition potential maps indicate the ability of a pixel to turn from one state to another or remain unchanged.According to the conditions of each region, different driving factors are used to produce transition potential maps.Here, based on previous studies [18,43,81], the following driving factors were chosen: (a) distance to water bodies; (b) distance to main roads; (c) distance to built-up areas; and (d) slope.Fuzzy membership functions (e.g., linear monotonic increase function) were used to standardize the driver maps into a 0-1 range, where 0 represents unsuitable locations and 1 represents ideal locations.The AHP technique was then run for weighting and specifying priorities regarding the relative importance of the driving factors.The details of the weights and control points are listed in Table 3.Then, Kappa variations were employed to evaluate the model by comparing the LULC map of 2016 with the simulated map of 2016.The standard accuracy of the models, which was not less than 80%, determined that they were potent predictive tools [32,37].In the present study, Kno, Kstandard, and Klocation were used to validate the model.Descriptions of these Kappa indices are summarized in Keshtkar and Voigt [32].Finally, the LULC map of 2016 was set as the base map, and the transition area matrices of 2000-2008 and 2008-2016 were used to forecast the LULC maps of 2024 and 2032, respectively.

LULC Change
In this study, the overall classification accuracy (OA) of the LULC maps for 1988 (88.75%), 1992 (90.83%), 1996 (90.33%), 2000 (89.67%), 2004 (91.92%), 2008 (88.92%), 2013 (90.92%), and 2016 (92.25%) were derived.The user's accuracy (UA) and producer's accuracy (PA) for the LULC maps of each year are presented in Figure 3a,b to gain the weights of suitable variables [88,89].The transition potential maps indicate the ability of a pixel to turn from one state to another or remain unchanged.According to the conditions of each region, different driving factors are used to produce transition potential maps.Here, based on previous studies [18,43,81], the following driving factors were chosen: (a) distance to water bodies; (b) distance to main roads; (c) distance to built-up areas; and (d) slope.Fuzzy membership functions (e.g., linear monotonic increase function) were used to standardize the driver maps into a 0-1 range, where 0 represents unsuitable locations and 1 represents ideal locations.The AHP technique was then run for weighting and specifying priorities regarding the relative importance of the driving factors.The details of the weights and control points are listed in Table 3.
Then, Kappa variations were employed to evaluate the model by comparing the LULC map of 2016 with the simulated map of 2016.The standard accuracy of the models, which was not less than 80%, determined that they were potent predictive tools [32,37].In the present study, Kno, Kstandard, and Klocation were used to validate the model.Descriptions of these Kappa indices are summarized in Keshtkar and Voigt [32].Finally, the LULC map of 2016 was set as the base map, and the transition area matrices of 2000-2008 and 2008-2016 were used to forecast the LULC maps of 2024 and 2032, respectively.
Major changes observed were the significant increase of urban/built-up areas and a sharp decline of cultivated land area.The built-up area increased by 103.82 km 2 from 40.53 km 2 (3.33%) in 1988 to 144.35 km 2 (11.88%) in 2016 with an average annual growth rate of 9.15%.In contrast, cultivated land decreased by 122.91 km 2 from 764.87 km 2 (62.94%) in 1988 to 641.96 km 2 (52.83%) in 2016 (Table 5).
Major changes observed were the significant increase of urban/built-up areas and a sharp decline of cultivated land area.The built-up area increased by 103.82 km 2 from 40.53 km 2 (3.33%) in 1988 to 144.35 km 2 (11.88%) in 2016 with an average annual growth rate of 9.15%.In contrast, cultivated land decreased by 122.91 km 2 from 764.87 km 2 (62.94%) in 1988 to 641.96 km 2 (52.83%) in 2016 (Table 5).
Major changes observed were the significant increase of urban/built-up areas and a sharp decline of cultivated land area.The built-up area increased by 103.82 km 2 from 40.53 km 2 (3.33%) in 1988 to 144.35 km 2 (11.88%) in 2016 with an average annual growth rate of 9.15%.In contrast, cultivated land decreased by 122.91 km 2 from 764.87 km 2 (62.94%) in 1988 to 641.96 km 2 (52.83%) in 2016 (Table 5).

Spatiotemporal Transition of LULC
Despite the overall expansion in urban/built-up areas, the expansion pattern varied across individual time periods as briefly outlined below: Period one (1988-1992): The speed of urban expansion was slower before 1992 and accelerated remarkably in the preceding periods due to socioeconomic and political factors.Increase in urban/built-up areas and cultivated land and the decline in vegetation were the major transformations observed.Urban area increased by 2.92 km 2 from 40.53 km 2 to 43.45 km 2 with an annual average growth rate of 1.8%.This was due to the transformation of remaining land-use classes, mainly 2.73 km 2 of cultivated land area into urban/built-up areas (Figure 6a).Vegetation cover declined by 5.69 km 2 , from 396.32 km 2 to 390.63 km 2 due to the conversion of 11.97 km 2 of vegetation cover into cultivated land.With this, cultivated land increased by 1.85 km 2 , from 764.87 km 2 to 766.72 km 2 .
resulting in the land price increasing by more than 300% since 2003 [55].The flourishing land market influenced aggressive urban area expansion, and cultivated land declined by 5.3 km 2 , from 718.08 km 2 to 712.78 km 2 (Figure 6d).
Period five (2004-2008): Remarkable increase in the urban/built-up areas and vegetation cover, as well as a sharp decline of cultivated land were witnessed during this period.Urban area increased by 16.05 km 2 , from 84.30 km 2 to 100.35 km 2 with an annual average growth rate of 4.75%.The continuous shift of people towards the Kathmandu valley as well as city centers from rural areas was also driven by access to public services and economic opportunities [20,25,65].Despite the alteration of 9.8 km 2 vegetation cover into cultivated land, it dropped by 24.21 km 2 , from 712.78 km 2 to 688.57km 2 .This increase was mainly driven by the conversion of 15.7 km 2 of cultivated land into urban areas and 20.36 km 2 into vegetation cover.The Shivpuri Nagarjun wild life reserve area was upgraded to a national park status in 2002.Additionally, strongly managed community forests have played key roles in the protection of the vegetation cover [55] (Figure 6e).
Period six (2008-2013): Continuous expansion of urban/built-up areas and a decline of cultivated land and vegetation cover were investigated during this period (Figure 6f).The urban area increased by 21.24 km 2 , from 100.35 km 2 to 121.59 km 2 , with an average annual growth rate of 4.23%.Cultivated land declined by 14.96 km 2 , from 688.57km 2 to 673.61 km 2 because 20.63 km 2 of cultivated land changed into urban areas and 6.47 km 2 into vegetation cover.The Kathmandu valley is a gateway and central touristic hub with several heritages which receives 90% of the tourists in the country.Hence, the infrastructure to cope with the huge influx of tourists was developed at the expense of cultivated land in the valley [92].
Period seven (2013-2016): Increases in urban/built-up areas and vegetation cover, and a decrease in cultivated land were the major changes observed during this period.Urban areas increased by 22.76 km 2 , from 121.59 km 2 to 144.35 km 2 with an average annual growth rate of 6.23%.In contrast, cultivated land declined by 31.65 km 2 , from 673.61 km 2 to 641.96 km 2 .Vegetation cover increased by 8.75 km 2 , from 408.42 km 2 to 417.17 km 2 .All these changes were associated with the conversion of 22.45 km 2 of cultivated land into urban areas and 14.86 km 2 into vegetation cover (Figure 6g).Period two (1992)(1993)(1994)(1995)(1996): Increase in urban/built-up areas and vegetation cover and decrease of cultivated land were the three important changes witnessed during this period.Urban area increased by 11.14 km 2 from 43.45 km 2 to 54.59 km 2 with an average annual growth rate of 6.41%.Urban parks, forest reserves, and forest resorts maintained the vegetation cover within the Kathmandu valley [55].The expansion of urban/built-up areas and the sharp decline of cultivated land was mainly because 10.98 km 2 of cultivated land was converted into urban/built-up areas (Figure 6b).The increase in forest area was closely associated with the community-based forest management biodiversity practice from 1975, the program launched across in Nepal after 1995 intended to conserve and rehabilitate degraded forest resources [90], which has been successful in achieving a significant increase in forest cover and density not only in the study area, but in the hill region of Nepal [91].
Period three (1996)(1997)(1998)(1999)(2000): Gradual increase in urban/built-up and continuous decline of cultivated land are the two major transformations observed during this period (Figure 6c).Urban areas increased by 23.11 km 2 , from 54.59 km 2 to 77.7 km 2 with an average annual growth rate of 10.58%.The country confronted a decade-long political upheaval from 1996 to 2006, which enforced internal mass migration to the cities [47].Urban in-migration contributed around 40% of the population growth of the Kathmandu valley in the 1990s.Due to high migration and population growth, the real estate market culminated after 1999 and vigorously continued in the preceding years [55].This led to intense pressure upon the cultivated land resulting in a decrease of 23.23 km 2 from 741.31 km 2 to 718.08 km 2 over the period.
Period four (2000)(2001)(2002)(2003)(2004): The trend of urban increase and decline of cultivated land also continued during this period.Urban areas increased by 6.6 km 2 , from 77.7 km 2 to 84.3 km 2 with an average annual growth rate of 2.12%.Of the expanded urban area, 6.05 km 2 was sourced from cultivated land.Since the political turmoil peaked during this period, the trend of in-migration continued.More than 500,000 refugees were displaced and accommodated inside the valley [92], resulting in the land price increasing by more than 300% since 2003 [55].The flourishing land market influenced aggressive urban area expansion, and cultivated land declined by 5.3 km 2 , from 718.08 km 2 to 712.78 km 2 (Figure 6d).
Period five (2004)(2005)(2006)(2007)(2008): Remarkable increase in the urban/built-up areas and vegetation cover, as well as a sharp decline of cultivated land were witnessed during this period.Urban area increased by 16.05 km 2 , from 84.30 km 2 to 100.35 km 2 with an annual average growth rate of 4.75%.The continuous shift of people towards the Kathmandu valley as well as city centers from rural areas was also driven by access to public services and economic opportunities [20,25,65].Despite the alteration of 9.8 km 2 vegetation cover into cultivated land, it dropped by 24.21 km 2 , from 712.78 km 2 to 688.57km 2 .This increase was mainly driven by the conversion of 15.7 km 2 of cultivated land into urban areas and 20.36 km 2 into vegetation cover.The Shivpuri Nagarjun wild life reserve area was upgraded to a national park status in 2002.Additionally, strongly managed community forests have played key roles in the protection of the vegetation cover [55] (Figure 6e).
Period six (2008-2013): Continuous expansion of urban/built-up areas and a decline of cultivated land and vegetation cover were investigated during this period (Figure 6f).The urban area increased by 21.24 km 2 , from 100.35 km 2 to 121.59 km 2 , with an average annual growth rate of 4.23%.Cultivated land declined by 14.96 km 2 , from 688.57km 2 to 673.61 km 2 because 20.63 km 2 of cultivated land changed into urban areas and 6.47 km 2 into vegetation cover.The Kathmandu valley is a gateway and central touristic hub with several heritages which receives 90% of the tourists in the country.Hence, the infrastructure to cope with the huge influx of tourists was developed at the expense of cultivated land in the valley [92].
Period seven (2013-2016): Increases in urban/built-up areas and vegetation cover, and a decrease in cultivated land were the major changes observed during this period.Urban areas increased by 22.76 km 2 , from 121.59 km 2 to 144.35 km 2 with an average annual growth rate of 6.23%.In contrast, cultivated land declined by 31.65 km 2 , from 673.61 km 2 to 641.96 km 2 .Vegetation cover increased by 8.75 km 2 , from 408.42 km 2 to 417.17 km 2 .All these changes were associated with the conversion of 22.45 km 2 of cultivated land into urban areas and 14.86 km 2 into vegetation cover (Figure 6g).

Urban Area Expansion and Orientation
Exploring where urbanization has occurred is important when determining how urbanization may be better managed in the future, and the city center is usually active in socioeconomic activities [61].In our analysis, the "center" generally referred to the urban center [33,85] and the assumed urban area gradually decreased outwards.Although, the essentiality to explore the orientation and extent of urban expansion exists, few studies have attempted it [50,84].For the orientation and extent analysis, we assumed Dhulikhel in Kabhrepalanchok, Kamalbinayak in Bhaktapur, Lagankhel in Lalitpur, and Singhadurbar in Kathmandu as the urban centers since they were the socioeconomic hubs where urbanization took place since the initial phase of our study, which was 1988.
According to the study, every direction from the city core has gained a remarkable urban increase, however, expansion and density has gradually decreased outside the assumed centers.The urban area of Kathmandu was limited to 25.43 km 2 in 1988, which gradually increased to 27.25 km 2 in 1992, 34.77 km 2 in 1996, and aggressively increased to 50.6 km 2 in 2000.The expansion trend continued and reached 55.06 km 2 in 2004, 66.53 km 2 in 2008, 78.44 km 2 in 2013, and 90.97 km 2 in 2016.The N-NE (Tokha, Budhanilkantha areas), NE-E (Shankarapur and Kageshwori Manohara areas), SW-W (Chandragiri and Kirtipur areas), NW-N (Tarakeshwor area), and E-SE (Baneshwar and Airport area) directions were highly urbanized while the W-NW (Nagarjun area) direction was relatively less urbanized.Regarding the SE-S direction, almost the entire area was urban, and the boundary was adjoined with the densely settled Lagankhel of Lalitpur district.The S-SW (Dakshinkali area) direction was comparatively less urbanized (Figure 7a).boundary was adjoined with the densely settled Lagankhel of Lalitpur district.The S-SW (Dakshinkali area) direction was comparatively less urbanized (Figure 7a).
Peripheral areas of the assumed city center, Lagankhel in Lalitpur district have been densely settled for several years.The total urban coverage of the district during 1988 was 7.73 km 2 , which increased to 8.33 km 2 in 1992 to 28.12 km 2 in 2016.The SW-W, W-NW, NW-N, N-NE, and NE-E directions were highly urbanized, but were relatively transected as smaller urban clusters during the course of the study.East sprawl was relatively complacent (Figure 7b).The areas at the proximity to the road networks were more urbanized.Higher urban enlargement was reported in the Lubhu, Lamatar, and Chapagaun areas of Lalitpur [50].
Bhaktapur is one of the most important gateways to the Kathmandu valley and consists of a large plain area where urban sprawl is more prevalent.The total urban coverage accounted for 4.8 km 2 in 1988 that gradually increased to 5.1 km 2 in 1992, and to 18.27 km 2 in 2016.The majority of the urban area expanded towards the SW-W (Balkot, Katunje, Sirutar, and Dadhikot areas) and W-NW (Madhyapur Thimi, Duwakot) directions; whereas the N-NE (Chhaling area) and E-SE (Tathali and Chipatol area) directions gained comparatively lower urban expansion (Figure 7c).
Peripheral areas of the assumed center, Dhulikhel in Kabhrepalanchok district have also been densely settled for several years.The total urban coverage of district during 1988 was 2.55 km 2 , which increased to 6.93 km 2 in 2016.The highest urban expansion occurred in the SW-W and W-NW (Panauti, Dhulikhel, and Banepa areas) directions.Similarly, medium expansion was witnessed towards the E-SE, NE-E, and N-NE (Panchkhal, Mandandeupur, and Namobudda areas), and relatively lower expansion has taken place towards the SE-S, S-SW, and NW-N directions (Dapcha and Balthali areas) (Figure 7d).

CA-Markov Model
In this study, the transition probability matrices produced by the Markov model prepared the details on the probability of transition between LULC types during the periods 2000-2008 and 2008-2016 (Table 6).
Models with accuracies in excess of 80% are typically considered very strong predictive tools [32,37].In evaluating the overall accuracy, using the value of Kno is better than the Kstandard [93].The Kno value was 0.91, which verifies the accuracy of this model.The model performed very well in its overall ability to predict LULC map of 2016 (Kstandard = 0.87), and the Klocation value of 0.88 indicates that the model provides a reasonable representation of location.Thus, based on the results, the CA-MC model is a reliable and powerful tool for predicting LULC changes.
These results showed that the probability of future changes of vegetation to cultivated land from 2000 to 2008 was 13.71%.This likelihood of transition decreased reasonably to 8.5% in 2016.Of all the phenomena, built-up areas experienced the largest increase between the first and second projections.After eight years (during 2024), there was an 89.5% chance a current built-up pixel would remain part of the status quo, however, after sixteen years (during 2032), this value would increase to 91.1%.Cultivated areas possessed the highest likelihood of changing into built-up areas for both periods.Peripheral areas of the assumed city center, Lagankhel in Lalitpur district have been densely settled for several years.The total urban coverage of the district during 1988 was 7.73 km 2 , which increased to 8.33 km 2 in 1992 to 28.12 km 2 in 2016.The SW-W, W-NW, NW-N, N-NE, and NE-E directions were highly urbanized, but were relatively transected as smaller urban clusters during the course of the study.East sprawl was relatively complacent (Figure 7b).The areas at the proximity to the road networks were more urbanized.Higher urban enlargement was reported in the Lubhu, Lamatar, and Chapagaun areas of Lalitpur [50].
Bhaktapur is one of the most important gateways to the Kathmandu valley and consists of a large plain area where urban sprawl is more prevalent.The total urban coverage accounted for 4.8 km 2 in 1988 that gradually increased to 5.1 km 2 in 1992, and to 18.27 km 2 in 2016.The majority of the urban area expanded towards the SW-W (Balkot, Katunje, Sirutar, and Dadhikot areas) and W-NW (Madhyapur Thimi, Duwakot) directions; whereas the N-NE (Chhaling area) and E-SE (Tathali and Chipatol area) directions gained comparatively lower urban expansion (Figure 7c).
Peripheral areas of the assumed center, Dhulikhel in Kabhrepalanchok district have also been densely settled for several years.The total urban coverage of district during 1988 was 2.55 km 2 , which increased to 6.93 km 2 in 2016.The highest urban expansion occurred in the SW-W and W-NW (Panauti, Dhulikhel, and Banepa areas) directions.Similarly, medium expansion was witnessed towards the E-SE, NE-E, and N-NE (Panchkhal, Mandandeupur, and Namobudda areas), and relatively lower expansion has taken place towards the SE-S, S-SW, and NW-N directions (Dapcha and Balthali areas) (Figure 7d).

CA-Markov Model
In this study, the transition probability matrices produced by the Markov model prepared the details on the probability of transition between LULC types during the periods 2000-2008 and 2008-2016 (Table 6).Models with accuracies in excess of 80% are typically considered very strong predictive tools [32,37].In evaluating the overall accuracy, using the value of Kno is better than the Kstandard [93].The Kno value was 0.91, which verifies the accuracy of this model.The model performed very well in its overall ability to predict LULC map of 2016 (Kstandard = 0.87), and the Klocation value of 0.88 indicates that the model provides a reasonable representation of location.Thus, based on the results, the CA-MC model is a reliable and powerful tool for predicting LULC changes.
These results showed that the probability of future changes of vegetation to cultivated land from 2000 to 2008 was 13.71%.This likelihood of transition decreased reasonably to 8.5% in 2016.Of all the phenomena, built-up areas experienced the largest increase between the first and second projections.After eight years (during 2024), there was an 89.5% chance a current built-up pixel would remain part of the status quo, however, after sixteen years (during 2032), this value would increase to 91.1%.Cultivated areas possessed the highest likelihood of changing into built-up areas for both periods.
According to the analysis, 11.88% of the entire study area was occupied by urban land in 2016 (Table 7), which will increase to 19.6% by 2032.The results also showed that the rapid development of urban land would be accompanied by a sharp decline in agricultural and vegetation lands(Table 7) at the expense of other land use types, but mostly at the cost of cultivated land (loss of 15.6%).However, minor loss has been predicted for vegetation cover, from 417.12 km 2 (34.3%) to 405.97 km 2 (33.4%) over the respective years (Figure 8).According to the analysis, 11.88% of the entire study area was occupied by urban land in 2016 (Table 7), which will increase to 19.6% by 2032.The results also showed that the rapid development of urban land would be accompanied by a sharp decline in agricultural and vegetation lands(Table 7) at the expense of other land use types, but mostly at the cost of cultivated land (loss of 15.6%).However, minor loss has been predicted for vegetation cover, from 417.12 km 2 (34.3%) to 405.97 km 2 (33.4%) over the respective years (Figure 8).
The results showed that the rapid development of urban land would be accompanied by a sharp decline in agricultural and vegetation land (Table 7).Consequently, in addition to the gradual decline of cultivated land area, this is closely associated with endangering food security and health, and the integrity of the ecosystem will be heavily affected [94].Urban development is also associated with an increase in population and personal and public transportation [33,94].Therefore, we can expect an increase in environmental pollutants in the future; and due to the destruction and degradation of natural purifiers (plants and soil), environmental hazards may appear as strong threats to the residents of this region.Gradual urban expansion in the four individual districts of the study area was predicted from the simulation analysis.Newer settlements will emerge, and existing settlements will be in-filled.The urban area of Kathmandu occupied 90.97 km 2 in 2016, which is expected to expand to 124.56 km 2 by 2024 and 148.47 km 2 by 2032.Furthermore, the urban area is expected to increase significantly in The results showed that the rapid development of urban land would be accompanied by a sharp decline in agricultural and vegetation land (Table 7).Consequently, in addition to the gradual decline of cultivated land area, this is closely associated with endangering food security and health, and the integrity of the ecosystem will be heavily affected [94].Urban development is also associated with an increase in population and personal and public transportation [33,94].Therefore, we can expect an increase in environmental pollutants in the future; and due to the destruction and degradation of natural purifiers (plants and soil), environmental hazards may appear as strong threats to the residents of this region.
Gradual urban expansion in the four individual districts of the study area was predicted from the simulation analysis.Newer settlements will emerge, and existing settlements will be in-filled.
The urban area of Kathmandu occupied 90.97 km 2 in 2016, which is expected to expand to 124.56 km 2 by 2024 and 148.47 km 2 by 2032.Furthermore, the urban area is expected to increase significantly in the N-NE, SW-W, W-NW, and NW-N directions (Figure 9a).The outskirt boundaries of Tarakeshwor, Nagarjun, Chandragiri/Kirtipur areas, and the eastern boundary of Shankarapur and Kageshwori Manohara areas are adjoined to hills where less urban expansion can be assumed in the future.Meanwhile, the urban area of Lalitpur district covered 28.12 km 2 during 2016, which is expected to grow continuously and reach 38.67 km 2 by 2024 and 41.84 km 2 by 2032.Higher urban expansion is predicted to occur in the E-SE, SE-S, and S-SW directions (Figure 9b).The trend of urbanization is expected to be outflanked in the city outskirts and hinterland adjacent to road networks.The urban coverage of Bhaktapur totaled 18.27 km 2 , which is predicted to reach 27.77 km 2 and 37.18 km 2 by 2024 and 2032, respectively.The majority of urban expansion will take place in the SW-W and W-NW directions (Figure 9c).During 2016, the urban coverage of Kabhrepalanchok district totaled 6.93 km 2 , which is expected to occupy 9.75 km 2 and 12.08 km 2 by 2024 and 2032, respectively.The N-NE, NE-E, S-SW, and SW-W directions will experience relatively higher expansion (Figure 9d). the N-NE, SW-W, W-NW, and NW-N directions (Figure 9a).The outskirt boundaries of Tarakeshwor, Nagarjun, Chandragiri/Kirtipur areas, and the eastern boundary of Shankarapur and Kageshwori Manohara areas are adjoined to hills where less urban expansion can be assumed in the future.Meanwhile, the urban area of Lalitpur district covered 28.12 km 2 during 2016, which is expected to grow continuously and reach 38.67 km 2 by 2024 and 41.84 km 2 by 2032.Higher urban expansion is predicted to occur in the E-SE, SE-S, and S-SW directions (Figure 9b).The trend of urbanization is expected to be outflanked in the city outskirts and hinterland adjacent to road networks.The urban coverage of Bhaktapur totaled 18.27 km 2 , which is predicted to reach 27.77 km 2 and 37.18 km 2 by 2024 and 2032, respectively.The majority of urban expansion will take place in the SW-W and W-NW directions (Figure 9c).During 2016, the urban coverage of Kabhrepalanchok district totaled 6.93 km 2 , which is expected to occupy 9.75 km 2 and 12.08 km 2 by 2024 and 2032, respectively.The N-NE, NE-E, S-SW, and SW-W directions will experience relatively higher expansion (Figure 9d).

Conclusions
The study investigated the spatial-temporal pattern of LULC change, in particular, the urban expansion dynamics of CKVAKD using maximum cloud-free time-series Landsat images between 1988 and 2016.As SVM tends to result in relatively higher accuracy than the ML classifier [70,77], the current study has utilized the SVM algorithm and is able to achieve an overall accuracy of LULC maps between 85% and 93%.The LULC analysis explored the remarkable increase in urban/built-up areas and forest cover and a significant decline of the cultivated land area.Increase in the urban area was assumed to be closely associated with the rise in infrastructure to accommodate the rising population, which in turn has caused losses of prime farm land.Urban expansion has sparked towards the western part of the study area when compared to the eastern section.Urban sprawl was multifaceted and largely attributed to the proximity and accessibility of road networks.Simulation analysis conducted for the years of 2024 and 2032 using a CA-Markov model predicted that urban area will gradually increase in future years while cultivated and vegetation land will continuously remain under pressure.Unmanaged settlement patterns and the continuous urban development process, leading to land fragmentation and pressure in land resources, continue the urban development process [20].Such unplanned urban sprawl and its detrimental consequences are likely to result in overall environmental disequilibrium in the valley [53,55].

Conclusions
The study investigated the spatial-temporal pattern of LULC change, in particular, the urban expansion dynamics of CKVAKD using maximum cloud-free time-series Landsat images between 1988 and 2016.As SVM tends to result in relatively higher accuracy than the ML classifier [70,77], the current study has utilized the SVM algorithm and is able to achieve an overall accuracy of LULC maps between 85% and 93%.The LULC analysis explored the remarkable increase in urban/built-up areas and forest cover and a significant decline of the cultivated land area.Increase in the urban area was assumed to be closely associated with the rise in infrastructure to accommodate the rising population, which in turn has caused losses of prime farm land.Urban expansion has sparked towards the western part of the study area when compared to the eastern section.Urban sprawl was multifaceted and largely attributed to the proximity and accessibility of road networks.Simulation analysis conducted for the years of 2024 and 2032 using a CA-Markov model predicted that urban area will gradually increase in future years while cultivated and vegetation land will continuously remain under pressure.Unmanaged settlement patterns and the continuous urban development process, leading to land fragmentation and pressure in land resources, continue the urban development process [20].Such unplanned urban sprawl and its detrimental consequences are likely to result in overall environmental disequilibrium in the valley [53,55].
The study assumed the eastern part of the study area as the potential zone for future urban expansion where the opportunities to develop planned and managed urban settlements still exist.Discouraging the haphazard use of urban land, controlling the encroachment of public land, and preparation of better guidelines in terms of physical infrastructure development is highly essential to secure a sustainable urban future [56].In addition, the government authorities need to enforce strict urban growth policies in coming days [20].The urbanization process of the study area was driven by population growth, migration, socioeconomic, and political factors.The continuous urban expansion trend in the preceding decades is likely to result in severe environmental disequilibrium.Thus, careful urban planning is always emphasized.

Figure 1 .
Figure 1.Location map of the study area.Figure 1. Location map of the study area.

Figure 1 .
Figure 1.Location map of the study area.Figure 1. Location map of the study area.
) the preparation of LULC maps with the same time interval (here, 2000, 2008, and 2016); (2) the calculation of transition area matrices based on LULC maps; (3) the generation of transition potential maps using driving factors; (4) evaluating the model's ability to simulate future changes based on kappa indices; and (5) simulating the LULC maps for the coming years (here, 2024 and 2032).The transition matrix file was calculated using the Markov model in TerrSet software.A transition area matrix records the number of cells that are expected to change from one land-cover class to another in a given period in the future.In this part, transitional area matrices were generated based on the historical LULC information from 2000-2008, 2008-2016, and 2000-2016 to explore how each land-cover was projected to change.Detailed study framework has been presented in Figure 2.
) the preparation of LULC maps with the same time interval (here, 2000, 2008, and 2016); (2) the calculation of transition area matrices based on LULC maps; (3) the generation of transition potential maps using driving factors; (4) evaluating the model's ability to simulate future changes based on kappa indices; and (5) simulating the LULC maps for the coming years (here, 2024 and 2032).The transition matrix file was calculated using the Markov model in TerrSet software.A transition area matrix records the number of cells that are expected to change from one land-cover class to another in a given period in the future.In this part, transitional area matrices were generated based on the historical LULC information from 2000-2008, 2008-2016, and 2000-2016 to explore how each land-cover was projected to change.Detailed study framework has been presented in Figure 2.

Figure 2 .
Figure 2. Framework of the study.

Figure 2 .
Figure 2. Framework of the study.

Figure 5 .
Figure 5. Trend of LULC change (a) Urban built-up, cultivated land, and vegetation cover; (b) sand and others, water body, and open field.

Figure 5 .
Figure 5. Trend of LULC change (a) Urban built-up, cultivated land, and vegetation cover; (b) sand and others, water body, and open field.

Figure 5 .
Figure 5. Trend of LULC change (a) Urban built-up, cultivated land, and vegetation cover; (b) sand and others, water body, and open field.

Table 3 .
Extracted weights based on analytic hierarchy process (AHP) and fuzzy standardization for urban areas.

Table 3 .
Extracted weights based on analytic hierarchy process (AHP) and fuzzy standardization for urban areas.

Table 7 .
Distribution of LULC categories (in km 2 ) and percentage of changes for 2016-2032.

Table 7 .
Distribution of LULC categories (in km 2 ) and percentage of changes for 2016-2032.