Monitoring and Modeling of Spatiotemporal Urban Expansion and Land-Use / LandCover Change Using Integrated Markov Chain Cellular Automata Model

Spatial–temporal analysis of land-use/land-cover (LULC) change as well as the monitoring and modeling of urban expansion are essential for the planning and management of urban environments. Such environments reflect the economic conditions and quality of life of the individual country. Urbanization is generally influenced by national laws, plans and policies and by power, politics and poor governance in many less-developed countries. Remote sensing tools play a vital role in monitoring LULC change and measuring the rate of urbanization at both the local and global levels. The current study evaluated the LULC changes and urban expansion of Jhapa district of Nepal. The spatial–temporal dynamics of LULC were identified using six time-series atmospherically-corrected surface reflectance Landsat images from 1989 to 2016. A hybrid cellular automata Markov chain (CA–Markov) model was used to simulate future urbanization by 2026 and 2036. The analysis shows that the urban area has increased markedly and is expected to continue to grow rapidly in the future, whereas the area for agriculture has decreased. Meanwhile, forest and shrub areas have remained almost constant. Seasonal rainfall and flooding routinely cause predictable transformation of sand, water bodies and cultivated land from one type to another. The results suggest that the use of Landsat time-series archive images and the CA–Markov model are the best options for long-term spatiotemporal analysis and achieving an acceptable level of prediction accuracy. Furthermore, understanding the relationship between the spatiotemporal dynamics of urbanization and LULC change and simulating future landscape change is essential, as they are closely interlinked. These scientific findings of past, present and future land-cover scenarios of the study area will assist planners/decision-makers to formulate sustainable urban development and environmental protection plans and will remain a scientific asset for future generations.


Introduction
Policy makers in developing countries face unprecedented challenges with regard to governing, urban planning and land-use/land-cover (LULC) management because of the prevailing high dynamic growth.Therefore, knowledge concerning past, current and future growth plays an important role in the decision-making process [1].The conversion of natural habitats into agricultural land for the maintenance of human livelihoods has been identified as the greatest driver of global environment change [2].The spatial distribution of agricultural land for urban development [3][4][5] and heterogeneous settlement patterns over agricultural land is a universal trend [6].Changes in LULC and rapid urban growth are subjects of great concern worldwide.Rural-urban migration, natural population growth and administrative reclassification of rural areas to urban areas are the major components [7,8] of global patterns of urbanization [9,10].World population crossed the first billion mark in 1830, reached 2 billion in 1930, 3 billion in 1960, 5 billion in 1987, 6 billion in 1999 [11], 7.6 billion in 2017 and is projected to reach 8.6 billion in 2030, 9.8 billion in 2050 and 11.2 billion in 2100 [12].The world urban population has increased much faster than with rural population and it was 14% in 1900 to 29.1% in 1950, 47% in 2005 [11] and 2014 accounted for 54% and projected about 66% by 2050 [13].As a result of various anthropogenic and natural factors, developing countries experience relatively more rapid urban growth and LULC changes than developed countries [14] and this growth triggers several environmental and ecological problems at multiple spatial scales [1,[14][15][16][17].
In recent decades, Southeast Asia and other developing countries under transitional economies (e.g., Cambodia, Laos, Myanmar and Vietnam) [18] as well as their cities have recorded a higher rate of urbanization [14].There was an increase in urban population from 9%, 12%, 24% and 19% in 1980 to 20%, 32%, 33% and 35% in 2012 regarding Cambodia, Laos, Myanmar and Vietnam, respectively [18].Despite its short urban history, Nepal has recorded the highest urban growth rate (6.6%) in South Asia exceeding Sri Lanka (2.2%), India (2.9%), Pakistan (4.4%) and Bangladesh (5.3%) [19].In the last 50 years, Nepal's population has increased nearly threefold, but the area of cultivated land has only increased twofold during the same period [20].The urban population of Nepal was 2.9% in 1952/54 and increased to 40.49% in 2015 [21].Such trend of urban expansion is a serious issue for sustainable urban development in developing countries [22]; since it is directly associated with sustainability issues such as food security and a secured urban future.
Remote sensing and a geographical information system (GIS) have been applied extensively to help understand the spatiotemporal patterns of urbanization [11] and urban planning [23].They have proven to be significant, powerful [16,17], cost effective and accurate [14,18,24,25].Different simulation and empirical models have been applied to predict changes in LULC [26,27]; however, few have used system dynamics models to determine the main drivers [28].Spatial simulation models are efficient tools for quantitative simulation [29] and integrating the Markov-chain model and cellular automata (CA) is considered by some to be one of the best options for the analysis of LULC on different spatial scales [30].One of the basic assumptions of the Markov chain LULC change is a stochastic process which describes the probability of one state i changing to another state j [25,29,31].The CA model predicts spatial change over a specific time period using a probability matrix [29].
Economic growth is indicated by the occurrence of urban sprawl and population growth.However, human-induced deforestation and development of the natural landscape to establish building, water supply, sewage, transport networks and large scale of construction work create multiple unused impacts on land and soil, biodiversity, air, water qualities as well as other environmental activities.Similarly, major construction works transform cities into impervious surfaces with the formation of heat islands which contribute to climate changes including increased occurrences of multiple natural disasters [32][33][34][35].
Thus, it is very important to access, monitor and model past, present and future LULC scenarios using RS/GIS tools and simulation models.This is because obtaining overtime spatial-temporal characteristics is a fundamental prerequisite for the formulation of effective urban policies, as well as economic, demographic and environmental plans to ensure sustainable development and maintain environmental equilibrium [1,13,30,36,37].
The study area of this work is the Jhapa District of Nepal, located in the southern plain of the country.This area has been experiencing accelerated urban growth since 1996.The urbanization process is influenced by migration and population growth has recorded a remarkable impact on other LULC, mostly prime agricultural land.Settlements are scattered and unmanaged and are likely to influence the natural/social environment, thereby resulting in an unbalanced ecosystem.Urban planning and management has been affected due to the lack of reliable urban databases, spatiotemporal information of LULC change and environmental concerns.There exist poor planning, weak institutional capacity and insufficient investment in the urban sector which hinder effective remedies which address urban poverty, disaster risk as well as social and natural environment issues [3,5,6,14,22].The sustainable development and urban future of the study area is skeptical.
The major objective of this study was to explore the spatiotemporal analysis of LULC, monitor past/present urban expansion trends and future change simulation by 2026 and 2036, using the CA-Markov Model.More specific objectives are to (a) analyze the spatiotemporal dynamics of LULC from 1989 to 2016, (b) evaluate the urban expansion ratio and spatiotemporal extent, and (c) predict future LULC change by 2026 and 2036.

Study Area
The study area, Jhapa district is located between 26.36     in the southern part of Nepal (Figure 1).The district shares its eastern and southern boundaries with India and its northern and western boundaries with Ilam and Morang districts, respectively.Administratively, it is divided into eight municipalities and seven rural municipalites.There was an increase in the district population from 593,737 in 1991 to 688,109 in 2001 and 812,650 in 2011 [20].
ISPRS Int.J. Geo-Inf.2017, 6, 288 3 of 21 spatiotemporal information of LULC change and environmental concerns.There exist poor planning, weak institutional capacity and insufficient investment in the urban sector which hinder effective remedies which address urban poverty, disaster risk as well as social and natural environment issues [3,5,6,14,22].The sustainable development and urban future of the study area is skeptical.
The major objective of this study was to explore the spatiotemporal analysis of LULC, monitor past/present urban expansion trends and future change simulation by 2026 and 2036, using the CA-Markov Model.More specific objectives are to (a) analyze the spatiotemporal dynamics of LULC from 1989 to 2016, (b) evaluate the urban expansion ratio and spatiotemporal extent, and (c) predict future LULC change by 2026 and 2036.

Study Area
The study area, Jhapa district is located between 26.36°, 26.80°, 87.63° and 88.20° in the southern part of Nepal (Figure 1).The district shares its eastern and southern boundaries with India and its northern and western boundaries with Ilam and Morang districts, respectively.Administratively, it is divided into eight municipalities and seven rural municipalites.There was an increase in the district population from 593,737 in 1991 to 688,109 in 2001 and 812,650 in 2011 [20].Jhapa district has a subtropical monsoon climate with summer temperatures ranging from 32 to 35 °C, winter temperatures ranging from 8 to 15 °C and an average annual rainfall of 270 mm.It is low lying, at 58 to 792 m in elevation and has highly fertile soil where rice, wheat, maize, oilseed, pulses, and cash crops such as tea, jute and nuts are grown.The area is rich in biodiversity, containing subtropical evergreen hardwood forests to decideous forests depending on the altitude Jhapa district has a subtropical monsoon climate with summer temperatures ranging from 32 to 35 • C, winter temperatures ranging from 8 to 15 • C and an average annual rainfall of 270 mm.It is low lying, at 58 to 792 m in elevation and has highly fertile soil where rice, wheat, maize, oilseed, pulses, and cash crops such as tea, jute and nuts are grown.The area is rich in biodiversity, containing subtropical evergreen hardwood forests to decideous forests depending on the altitude [38].Geologically, the area has sediment of the Indo-Gangetic plain, rocks from the Siwalik range as well as small parts of the Lesser Himalaya.The Himalayan frontal thrust, Intra Chure thrust and main boundary thrust are located in the southern to northern parts.The Mechi, Deuniya, Biring, Kankai, Ratuwa and Mawa Khola are major rivers experiencing bankfull discharges during the monsoon season, resulting in river bank erosion and loss of farmland.The district is well supplied with water and has a good climate that makes it a major agricultural center.Several proximate and underlying causes are contributing to the migration of people to the area and have increased the population density in the district and urbanization, often at the expense of prime agricultural land.

Satellite Data and Pre-Processing
For spatial-temporal analysis of LULC change and monitoring of urban expansion of the study area, six time-series of atmospherically-corrected surface reflectance Landsat images from 1989 to 2016 were obtained from the United States Geological Survey (USGS) website [39].All scenes were verified for geometric accuracy and all data were projected on WGS 1984, UTM zone 45N.Surface reflectance data are ready for application and these data sets have been post-production processed and included atmospheric correction, geometric correction as well as other adjustments (for more details visit: http://landsat.usgs.gov/CDR_LSR.php).The maximum cloud-free images were included in the study.All scenes were verified for geometric accuracy and all data was projected on to UTM WGS 1984, zone 45N.Overall, the OLI bands are spectrally narrower than the ETM+ band, particularly in the near infrared (NIR) region [40,41].The ultra-blue band (0.43-0.451 µm) of Landsat 8 was removed for consistency and renamed for other images (e.g., 1, 2, 3, 4, 5 and 7).Band 1 of Landsat 8 was used for the coastal areas and shallow water and aerosol, dust and smoke detection studies [40].High spatial resolution images from Google Earth [42] and digital topographical maps (scale 1:25,000) published by the Survey Department of the Government of Nepal in 1995 (sheet numbers: 2687 (04D, 07B, 07D, 08A-08D, 11B, 12A, 12B, 12D) and 2688 (01C, 01D, 05A, 05D, 09A, 09C)) [43] were used to validate the results.Region of interest (ROI) boundaries representing district, municipal and village level study areas were delineated for analysis.The Jhapa district profile [38] was also collected for refinement of the study.

Extraction of LULC Maps and Analysis
Images were stacked, subset and analyzed in the ENVI [44] environment and classified using the maximum likelihood algorithm and change analysis for 1989, 1996, 2001, 2006, 2011 and 2016.Supervised approaches using a maximum likelihood classifier algorithm [23,30,[45][46][47][48] were applied for the extraction of LULC.A modified land-cover classification system was used for remote sensing data as recommended by Anderson et al. [49].For the extraction of LULC, seven classes were identified: urban (built-up), cultivated land, forest, shrub, sand, water bodies and tea farming area (Table 1).The study primarily concentrated on mapping of urban growth and sprawl.Because Jhapa is one of the few important districts for tea farming, this land use has been classified separately.The land change modeler (LCM) of TerrSet software [50] was used to assess changes and transitions in LULC in different years.Assessment of classification accuracy is necessary when LULC maps are derived from remote sensing techniques [24].Although the accuracy is limited by satellite image resolution [6], accuracy assessment of LULC maps was achieved using a random sampling method.After land-cover classification, 210 sample points were randomly selected for the evaluation of classification accuracy for seven land-cover classes.Accuracy assessment was based on the calculation of the overall accuracy, user's accuracy, and producer's accuracy.Due to the unavailability of updated land-cover data, topographical maps (scale 1:25,000) of 1995, field experience and Google Earth images were used as supporting details for accuracy assessment.

Quantification of LULC Based Transition Analysis
Land-cover change was calculated for six time periods (1989 to 1996, 1996 to 2001, 2001 to 2006, 2006 to 2011 and 2011 to 2016).For each time period, the transition matrix consists of rows and columns of landscape categories at times T 1 and T 2 .

Measuring Urban Expansion Rate and Analysis of Urban Extent
The urban expansion growth rate [51] of the study area was measured by calculating the total new urban area.The urban expansion rate indicates the average annual urban area growth in succeeding years.
where MUER refers to the urban expansion rate (km 2 /year) and U 1 , U 2 represent the urban area (km 2 ) between times T 1 and T 2 in years.
We quantified the centrality of the land-use spatial pattern using the ring-based analysis approach.Ring-based buffer analysis [10,[52][53][54][55] was used to explore the urban expansion characteristics based on landscape orientation.The boundary of the buffer zone outward from the city center (Birtamod) was drawn in Arc GIS 10.1 at equal intervals of 2 km.The buffered area from 1989 to 2016 was mapped.

Simulation of LULC Change
Modeling systems are intended to simplify the complexity of urban systems and make them easier to understand [56].Markov chain models have been widely applied for ecological modeling [57] as they show the descriptive power and simple trend projection [30].The CA model addresses spatial allocation and location of change [1] and as it has many advantages for modeling urban phenomena [36] including (a) cell (b) neighborhood, (c) rules, (d) time, and (e) state [29,58,59].Mainly the modeling system consisted of the following procedures: (a) land-cover maps of 1996, 2006 and 2016 were assessed, (b) the transition area matrix was calculated using the Markov process, (c) transition suitability images were prepared using multi-criteria evaluation (MCE), analytic hierarchy process (AHP) models and fuzzy membership functions in TerrSet and applied in CA-Markov (d) the actual and predicted maps of the year 2016 were evaluated, and finally (e) the LULC maps for the year 2026 and 2036 were simulated using the CA-Markov Model.

CA-Markov Modeling
The Markov model itself does not provide the spatial location of future LULC, so the hybrid CA-Markov was used to achieve the objectives.First, the transition matrix file was created and applied to the model for the specified time of 20 years.Next, CA filters were determined as a standard 5 × 5 contiguity filter in the modeling process.The filter used analysis is based on: The LULC map of 2016 was used as the base map to simulate LULC maps for the years 2026 and 2036 by calculating the transition area matrix of 2006-2016 and 1996-2016, respectively.

Generating Transition Potential Maps
In this study, GIS algorithms, MCE and fuzzy membership function were applied to extract the transition potential maps of land-cover types.This step determined the status of change.Transition potential maps represent the ability of a pixel to change to a new category or remain unchanged in each transition based on driving factors [60].In this study, slope, distance to main roads, distance to water bodies and distance to built-up areas were set as driving factors.They were chosen based on their use in similar previous studies [1,29,60].Fuzzy membership functions (e.g., J-shaped monotonic decrease functions) were used to rescale driver maps into the 0-255 range.AHP was then run to determine the weight of the driving factors using pairwise evaluation.The weighting parameters and control points for each driving factor were chosen based on expert knowledge, interviews and recent similar studies e.g., [1,29,61].The individual weights and control points determined are listed in Table 2. Evaluation of the model was done by comparing the predicted map for 2016 with the real map of 2016 based on the kappa variation.An accuracy higher than 80% infers some confidence in the simulation [29,61].A study framework was presented to achieve the research objectives presented in Figure 2.

LULC Dynamics
The overall classification accuracy for the extracted LULC maps were: 1989 (80.66%), 1996 (85.2%), 2001 (82.33%), 2006 (85.88%), 2011 (84.8%) and 2016 (86.2%).These results indicate the suitability of the classified remote sensing images for effective and reliable LULC change analysis and modeling.Also, analysis of user's and producer's accuracies shows that the classes with the highest producer's accuracy were those of tea (100%, 2001), water (96.43%,2011), and built-up areas (96.31%, 2001) followed by forest (95.59%, 2006).The lowest producer's accuracy was obtained for sand (76.67%, 2001).User's accuracy was higher for forest (93.44%, 2016), cultivated lands (93.33%, 1996), and built-up areas (92.53%, 2016), followed by water bodies (91.00%, 1996).The lowest user's accuracy was found for the sand (73.33%, 1989) (Table 3).Between 1989 and 2016, a wide range of change in land-use was observed, including a linear increase in urban area, the continuous decline of cultivated land, minimal change in forest cover and shrubs and non-linear transformations of water bodies and sand.The smallest fluctuations were observed for the area of tea farming.Although the water bodies decreased in some time periods,

LULC Dynamics
The overall classification accuracy for the extracted LULC maps were: 1989 (80.66%), 1996 (85.2%), 2001 (82.33%), 2006 (85.88%), 2011 (84.8%) and 2016 (86.2%).These results indicate the suitability of the classified remote sensing images for effective and reliable LULC change analysis and modeling.Also, analysis of user's and producer's accuracies shows that the classes with the highest producer's accuracy were those of tea (100%, 2001), water (96.43%,2011), and built-up areas (96.31%, 2001) followed by forest (95.59%, 2006).The lowest producer's accuracy was obtained for sand (76.67%, 2001).User's accuracy was higher for forest (93.44%, 2016), cultivated lands (93.33%, 1996), and built-up areas (92.53%, 2016), followed by water bodies (91.00%, 1996).The lowest user's accuracy was found for the sand (73.33%, 1989) (Table 3).Between 1989 and 2016, a wide range of change in land-use was observed, including a linear increase in urban area, the continuous decline of cultivated land, minimal change in forest cover and shrubs and non-linear transformations of water bodies and sand.The smallest fluctuations were observed for the area of tea farming.Although the water bodies decreased in some time periods, there was no overall decline of permanent water resources.Instead, the area covered by water was seasonal and determined by the changes in rainfall.The Biring, Kankai, Mechi, Deuniya and Ratuwa rivers as well as some streams which cross the district all originate in the mid-hills.The silty Churiya hill is located in the northern belt and during the monsoon season, the swollen rivers spill over causing bank erosion, erosion of agricultural land and accumulation of sand and silt.All these natural activities result in a strongly unpredictable arrangement of sand, water bodies and cultivated land in areas affected by floods.
The study period showed an overall seven-fold increase in urban area for Jhapa district with an area of 12.35 km 2 in 1989 increasing to 13 [25] states that Jhapa is a Tarai district that experienced over 25% inter-district lifetime migration, a high rate of immigration (24.5%) and a low rate of emigration (9.2%).The population density of the district was 299 people per km 2 in 1981, which reached 506 per km 2 in 2011, for which the national level population density was 180 in 2011.All these factors influenced the population growth of the district.The enormous mass migration, mainly from the adjacent mid-hill districts of Taplejung, Ilam and Panchthar along with the Terhathum and Dhankuta districts was due to political upheaval between 1996 and 2006.After 2006, as the political insurgency settled down, the pace of infrastructure development such as construction of road networks, hospitals, bus parks and academic institutions accelerated.This condition resembles to the urban scenario of Sri Lanka, particularly the Colombo metropolitan area where the development activities were accelerated after the settlement of three decades of civil war in 2009 [62].
The area of cultivated land increased from 1259.34 km 2 in 1989 to 1290.34 km 2 in 1996 at the cost of shrubs, forest, water, sand and tea (Tables 4 and 5; Figures 3 and 4).Significant decreases in forest and shrub cover were the other main changes from 1989 to 1996.The shrub area increased from 32.31 to 42.88 km 2  from 1989 to 2016.The region consists of the plain arable land which is strongly affected by floods, especially during monsoon season.Most of the sandy area along the riverbanks in the district has been transformed to cultivated land, which has led the statistical increase in cultivated land; however, fertile agricultural land was converted into urban land.There is an intricate relationship between anthropogenic factors, LULC change and sustainable future environment.Uncontrolled urbanization is mainly associated with prime farmland loss and is closely related with sustainable food security.Furthermore, unplanned road networks are also associated with the secure future and decreased environmental equilibrium.The region is seriously affected by seasonal floods and inundation.Each year, several catastrophes such as flood, inundation and riverbank erosion due to seasonal rainfall, result in the loss of lives and property.The recent flood and inundation killed 12 people, affected 2404 families, and destroyed huge physical properties [63].
Although almost part of the district is plain with fertile cultivated land, some parts of the northern territory belong to Churiya (Siwalik) region which is composited of gravel and infertile land.This region is a center for the extraction of sand and gravel required for construction activities in the urban areas which could become a threat in the future.There is an intricate relationship between anthropogenic factors, LULC change and sustainable future environment.Uncontrolled urbanization is mainly associated with prime farmland loss and is closely related with sustainable food security.Furthermore, unplanned road networks are also associated with the secure future and decreased environmental equilibrium.The region is seriously affected by seasonal floods and inundation.Each year, several catastrophes such as flood, inundation and riverbank erosion due to seasonal rainfall, result in the loss of lives and property.The recent flood and inundation killed 12 people, affected 2404 families, and destroyed huge physical properties [63].
Although almost part of the district is plain with fertile cultivated land, some parts of the northern territory belong to Churiya (Siwalik) region which is composited of gravel and infertile land.This region is a center for the extraction of sand and gravel required for construction activities in the urban areas which could become a threat in the future.

Spatial Transitions
Various natural and anthropogenic factors have resulted in LULC change and the transition ratios differ for each time period.Figures 5a-e shows the gains and losses of LULC during the period under investigation.LULC change for the district can be grouped into four stages.The first stage represents the period from 1989 to 1996 and saw a rise in cultivated area from 1259.34 to 1290.34 km 2 mainly from conversion of 13.22 km 2 of shrub, 12.45 km 2 of forest, 9.51 km 2 of water, 6.6 km 2 of sand and 4.13 km 2 of tea farming area.Two other major transformations were an increase in sand area from 48 to 73.76 km 2 and a remarkable decline in water bodies from 71.72 to 34.40 km 2 mainly due to the conversion of 29.34 km 2 of water bodies to sand (Figure 5a).In the second stage from 1996 to 2001, the urban area doubled from 13.70 km 2 in 1996 to 27.12 km 2 in 2001 at the cost of 13.55 km 2 of cultivated land.Further decreases in cultivated land to 1219.26 km 2 were due to loss of forest (9.44 km 2 ), shrubs (12.81 km 2 ), sand (35.77 km 2 ), water (5.34 km 2 ) and tea (3.21 km 2 ).Of this, 4.66 km 2 of forest changed to cultivated land and 3.42 km 2 to shrubs (Figure 5b).

Spatial Transitions
Various natural and anthropogenic factors have resulted in LULC change and the transition ratios differ for each time period.Figure 5a-e shows the gains and losses of LULC during the period under investigation.LULC change for the district can be grouped into four stages.The first stage represents the period from 1989 to 1996 and saw a rise in cultivated area from 1259.34 to 1290.34 km 2 mainly from conversion of 13.22 km 2 of shrub, 12.45 km 2 of forest, 9.51 km 2 of water, 6.6 km 2 of sand and 4.13 km 2 of tea farming area.Two other major transformations were an increase in sand area from 48 to 73.76 km 2 and a remarkable decline in water bodies from 71.72 to 34.40 km 2 mainly due to the conversion of 29.34 km 2 of water bodies to sand (Figure 5a).In the second stage from 1996 to 2001, the urban area doubled from 13.70 km 2 in 1996 to 27.12 km 2 in 2001 at the cost of 13.55 km 2 of cultivated land.Further decreases in cultivated land to 1219.26 km 2 were due to loss of forest (9.44 km 2 ), shrubs (12.81 km 2 ), sand (35.77 km 2 ), water (5.34 km 2 ) and tea (3.21 km 2 ).Of this, 4.66 km 2 of forest changed to cultivated land and 3.42 km 2 to shrubs (Figure 5b).
In the third stage (2001 to 2006), 10.34 km 2 of cultivated land was urbanized, 8.32 km 2 of sandy area changed to cultivated land and 11.49 km 2 to water, but 6.42 km 2 of water changed to sand, which resulted in an overall decrease in sand from 121.26 to 107.87 km 2 and increase in water from 28.61 to 33.32 km 2 between 2001 and 2006 (Figure 5c).In the fourth period (from 2006 to 2011), urbanization was intense with the conversion of 20.73 km 2 of cultivated land.Although the transformation of 4.48 km 2 of shrubs to cultivated land occurred, the change of 6.77, 6.34 and 3.83 km 2 of cultivated land to sand, water and tea, respectively, along with urban change caused an overall decrease in cultivated land.Increased population was accommodated in scattered settlements and/or the centers developed at the cost of agricultural land (Figure 5d).The final period of the study (from 2011 to 2016) was characterized by an increase in urban areas from 58.65 to 70.90 km 2 with a loss of 11.45 km 2 of cultivated land to urban land.The conversion of 6.26 km 2 of tea farms to cultivated land demonstrated a pronounced decrease in tea farming area from 28.95 to 22.52 km 2 (Figure 5e). Figure 5f shows the generalized pattern of change between 1989 and 2016 with urbanization moving outwards from the core.The darker red areas represent the more intensively changed areas which gradually decrease towards the outskirts of the city.This cubic analysis was presented using the LULC map of 1989 and 2016.The darker red areas represent the more intensively changed areas which gradually decrease towards the outskirts of the city.The eastern and western parts have more urban areas than the southern and central areas.In the third stage (2001 to 2006), 10.34 km 2 of cultivated land was urbanized, 8.32 km 2 of sandy area changed to cultivated land and 11.49 km 2 to water, but 6.42 km 2 of water changed to sand, which resulted in an overall decrease in sand from 121.26 to 107.87 km 2 and increase in water from 28.61 to 33.32 km 2 between 2001 and 2006 (Figure 5c).In the fourth period (from 2006 to 2011), urbanization was intense with the conversion of 20.73 km 2 of cultivated land.Although the transformation of 4.48 km 2 of shrubs to cultivated land occurred, the change of 6.77, 6.34 and 3.83 km 2 of cultivated land to sand, water and tea, respectively, along with urban change caused an overall decrease in cultivated land.Increased population was accommodated in scattered settlements and/or the centers developed at the cost of agricultural land (Figure 5d).The final period of the study (from 2011 to 2016) was characterized by an increase in urban areas from 58.65 to 70.90 km 2 with a loss of 11.45 km 2 of cultivated land to urban land.The conversion of 6.26 km 2 of tea farms to cultivated land demonstrated a pronounced decrease in tea farming area from 28.95 to 22.52 km 2 (Figure 5e). Figure 5f shows the generalized pattern of change between 1989 and 2016 with urbanization moving outwards from the core.The darker red areas represent the more intensively changed areas which gradually decrease towards the outskirts of the city.This cubic analysis was presented using the LULC map of 1989 and 2016.The darker red areas represent the more intensively changed areas which gradually decrease towards the outskirts of the city.The eastern and western parts have more urban areas than the southern and central areas.

Urban Expansion
As mentioned in Section 2.4, the total extent of urban area of the region in 1989 was 12.35 km 2 with an average annual growth rate of 1.56% per year in the succeeding seven years.In 1996, the urban coverage reached 13.38 km 2 , but sharply increased to 27.12 km 2 in 2001 with the noticeable annual urban growth rate of 19.57%.It exponentially increased to 37.09 km 2 , but relatively decreased to 7.35% by 2006.The dramatic expansion of 58.65 km 2 with an average annual growth rate of 11.63% created a historical transformation in the urban history of the district and the continuous urban expansion to 70.90 km 2 with an average growth rate of 4.18% can be observed between 2011 and 2016.
Driving factors, which are defined as the underlying elements [14,35] trigger urban expansion are responsible for the urban expansion which occurs at specific temporal and spatial levels.Contextual reference is required for a better understanding of the drivers of urbanization, since distinct factors play different roles in specific time and space.Physical, socioeconomic and neighborhood factors are the main factors which directed the urbanization of Beijing; of which socioeconomic factors contributed the highest, except in the first 12 years of the study period-1972 to 2010.Also, the magnitude of physical and neighborhood factors declined as the socioeconomic factors increased [64].Factors such as physical setting, population growth and economic development influenced urbanization in Greater Dhaka [14].
Proximate and underlying factors of public service accessibility, economic opportunity, population growth, globalization, political conditions, government plans and policies as well as social, cultural and rural urban linkages largely contributed to the migration of the population in the Nepalese context [6,[65][66][67][68][69].According to the CBS report 2014 [20], migration for agriculture (21%) was predominant over all other economic migration in rural areas, whereas service oriented migration was predominant in urban areas (17%).
The accelerated urbanization of Jhapa district is mainly driven by economic opportunity, political condition and development activities in specific periods which are further accompanied by access to public services, population growth, topography, government plans/policies and globalization factors.Such conditions are also seen in other Nepalese cities.The process of urbanization in Pokhara is mainly governed by greater economic opportunities, public service accessibility and globalization [66,68,69].Economic opportunities in core, population growth in the fringe and political situation in the rural areas highly contributed to the urbanization process of Kathmandu valley [19].Similarly, urban sprawl is additionally driven by physical conditions, public service accessibility and globalization as well as plans and policies [19,21,66].The major drivers of urbanization in the rapidly urbanizing city Biratnagar are explored in Rimal, 2011 study [68] which properly reflects and represents the urbanization scenario of the Tarai region, particularly the adjoined Jhapa District.

Urban Extent from City Center Outwards
In our analysis, the "center" generally refers to the urban center [52,55] and the assumed urban area gradually decreased outwards.Exploring the spatiotemporal change dynamics, orientation, cause and consequences of urban sprawl is highly essential.Nevertheless, very limited studies have outlined the contexts in terms of developing world.
The city center is usually regarded as the socioeconomic center [64].The direction of urban expansion in the study area has been highly influenced by its socioeconomic process than its physical setting.Kakarbhitta, Birtamod, Damak and Bhadrapur are the four cities and Birtamod was assumed as the central point in our study, since it represented the largest part of urban area from the very beginning.
Ring-based analysis identified the greatest variation in the urban expansion characteristics from 1989 to 2016 (Figures 6 and 7).The areas near the East-West highway and Mechi highway showed higher rates of urbanization compared to the south and south-west peripheries.There was a gradual increase in urban area within 2-4 km of the center of Birtamod and dense settlements at Charpane and Buttabari along with the central point itself are located within these rings.The four emerging cities of Chandragadhi, Surunga, Dhulabari and Budhabare are located within ring six and Kakarbhitta in ring nine.Thus, the area from 12 to 18 km from Birtamod has experienced considerable urbanization.Consequently, the area from 28 to 32 km hold important settlements at Charpane and Buttabari along with the central point itself are located within these rings.The four emerging cities of Chandragadhi, Surunga, Dhulabari and Budhabare are located within ring six and Kakarbhitta in ring nine.Thus, the area from 12 to 18 km from Birtamod has experienced considerable urbanization.Consequently, the area from 28 to 32 km hold important urban coverage because the densely-settled Damak and Gauradaha are located within the specified area.There was a gradual increase in urban area within 2-4 km of the center of Birtamod and dense settlements at Charpane and Buttabari along with the central point itself are located within these rings.The four emerging cities of Chandragadhi, Surunga, Dhulabari and Budhabare are located within ring six and Kakarbhitta in ring nine.Thus, the area from 12 to 18 km from Birtamod has experienced considerable urbanization.Consequently, the area from 28 to 32 km hold important urban coverage because the densely-settled Damak and Gauradaha are located within the specified area.

Analysis of Transition Matrix
The transition potential matrix was computed based on LULC conditions during the periods of 1996-2006, 2006-2016 and 1996-2016 to show how each land type was projected to change.Data elements lying on the diameter of the matrix indicate that a phenomenon will probably remain constant over time, while off-diagonal data indicate that various phenomena are likely to be converted to one of the other phenomena [16].For example, the transition probability matrix shows that the probability of future loss of cultivated land to built-up area from 1996 to 2006 was 4.3%.This probability of change increased reasonably to 9% in 2016.Table 6 shows that, for both periods, cultivated land possessed the highest likelihood of transforming into built-up areas.Moreover, urban, cultivated, forest and water bodies were more likely to remain stable in the second period (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016) compared to the first one (1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006), while other phenomena exhibited higher probabilities of changing.Tea land experienced the biggest decrease in both periods and water bodies showed the greatest increase between the first and the second projections.The likelihood of stability of land under tea cultivation decreased from 84% to 70% and the probability for most of these classes to be used for other agricultural purposes was 26% (Table 6).Comparison of the ground-truth map from 2016,the map simulated by the CA-Markov model using kappa variations and the Kstandard, Kno, and Klocation indicators resulted in measures of 87.5%, 90% and 94.5%, respectively.Visual comparison also shows great similarity between the two maps of 2016 (Figure 8); therefore, based on the kappa values, the CA-Markov model can be used to simulate future land-cover conditions.
The present research showed that "distance from built-up areas" played a key role in urban development.Urban development principles, access to facilities (including hospitals, schools, stores, etc.), and security priorities cause this factor to take on special importance."Distance from main roads," as an index of access, and "slope," as a limiting factor, are the next effective factors in urban development, respectively (Table 2).These key drivers are usually considered by managers when solving key problems related to land-use development.About 4.43% of the total land area in the region was urban in 2016 and is predicted to reach 8.11% and 10.5% by 2026 and 2036, respectively (Table 7).In contrast, cultivated land will decline from 1181.27 km 2 (73.7%) to 993.61 km 2 (62%), and forestland will increase from 143.01 km 2 (8.9%) to 146.2 km 2 (9.1%) by 2036.

Analysis of Simulation Results
About 4.43% of the total land area in the region was urban in 2016 and is predicted to reach 8.11% and 10.5% by 2026 and 2036, respectively (Table 7).In contrast, cultivated land will decline from 1181.27 km 2 (73.7%) to 993.61 km 2 (62%), and forestland will increase from 143.01 km 2 (8.9%) to 146.2 km 2 (9.1%) by 2036.Despite the urban development and population increase recorded in recent years, the expansion of agricultural lands has stopped completely in the region since about 1996 (Table 4) and expansion of urban lands has taken place instead.As shown in Figure 9, urban land will expand in the city outskirts, because agricultural land is very suitable for urban development with respect to access (proximity to built-up area) and topography [29].Modeling further indicated an increase in natural land-cover (Natural land = Total land area − (Farmland area + Built-up area)) [70].Natural areas are expected to cover 25.9% of the study area by 2036, which is an increase of 5.5% in comparison with the current distribution.These results are in contrast to those of Keshtkar and Voigt [60], who have shown a decrease in natural land-cover and increase in landscape fragmentation with the expansion of urban land.The increase of natural lands is mainly due to the increase of sand and shrubs.This increase is principally a result of the massive rural-urban migration and the abandonment of agricultural land.This migration can be caused by socio-economic, ecological and mismanagement drivers [71].First, the abandoned agricultural sites converted to the bare lands and show spectral behavior like sand.Finally, the bare lands were either planted with shrubs and trees or gradually replaced by spontaneous growth of grass, shrubs, and trees via secondary succession [72,73].Based on the authors' local and field visit knowledge, it is explored that newer settlements have expanded onto private farmland with the direct involvement of real estate agents who purchase the land from the local owners, divide the land into smaller plots and resell them at higher rates.The rapid development of the real estate business has intensified agricultural land loss on the urban fringe because investment in the land is considered to be a secure sector.It has resulted in land fragmentation and farmland loss and fosters incompatible and disorganized land use.
The government of Nepal has introduced different land-use acts, policies, laws and regulations regarding LULC management.The National Land Use Policy of 2015 [77] emphasizes land consolidation and compact settlement to reduce land fragmentation and scattered settlements.The Settlement Development, Urban Planning and Building Construction Guidelines in 2015 provide guidance for sustainable urban planning.However, the Ministry of Urban Development (MoUD) [78] has reported that their implementation has been hindered by the lack of long-term strategies and political commitment and poor understanding of the concept of sustainability.The lack of a proper land-use plan has resulted in urban growth and loss of productive cultivable land for development [21].Random urbanization and scattered settlements are the common features of the study area due the lack of effective urban planning.Another harmful effect of urban development is the excessive absorption of solar radiation that worsens the urban heat island effect.Satellite studies showed that all cities in the world (especially the metropolises), are faced with this problem as a result of the materials used, especially dark-colored building materials [79].Heat islands can directly or indirectly influence the health and social well-being of citizens.In the United States, 1000 people die annually due to heat intensity [80].Higher temperatures, especially in summer, increase energy demand for cooling purposes, and this increase generally increases the emission of particles that pollute the air, such as greenhouse gases from power plants  The increase of natural lands confirms that management decisions were taken with the goal of protecting the natural ecosystems in the region.Although population increases and the consequent urban growth in the study region will not reduce the natural areas, they will increase landscape fragmentation.Landscape fragmentation and land holding are high in Nepal with an average of more than three parcels per holding [74].Previous studies have shown that landscape fragmentation plays a fundamental role in reducing the connection between ecosystems [75], decreasing the possibility of migration for organisms especially for plant species [60], and increasing competition between species and their extinction [76].
Based on the authors' local and field visit knowledge, it is explored that newer settlements have expanded onto private farmland with the direct involvement of real estate agents who purchase the land from the local owners, divide the land into smaller plots and resell them at higher rates.The rapid development of the real estate business has intensified agricultural land loss on the urban fringe because investment in the land is considered to be a secure sector.It has resulted in land fragmentation and farmland loss and fosters incompatible and disorganized land use.
The government of Nepal has introduced different land-use acts, policies, laws and regulations regarding LULC management.The National Land Use Policy of 2015 [77] emphasizes land consolidation and compact settlement to reduce land fragmentation and scattered settlements.The Settlement Development, Urban Planning and Building Construction Guidelines in 2015 provide guidance for sustainable urban planning.However, the Ministry of Urban Development (MoUD) [78] has reported that their implementation has been hindered by the lack of long-term strategies and political commitment and poor understanding of the concept of sustainability.The lack of a proper land-use plan has resulted in urban growth and loss of productive cultivable land for development [21].Random urbanization and scattered settlements are the common features of the study area due the lack of effective urban planning.
Another harmful effect of urban development is the excessive absorption of solar radiation that worsens the urban heat island effect.Satellite studies showed that all cities in the world (especially the metropolises), are faced with this problem as a result of the materials used, especially dark-colored building materials [79].Heat islands can directly or indirectly influence the health and social well-being of citizens.In the United States, 1000 people die annually due to heat intensity [80].Higher temperatures, especially in summer, increase energy demand for cooling purposes, and this increase generally increases the emission of particles that pollute the air, such as greenhouse gases from power plants [81].Moreover, population increase and urban development naturally result in an increase in the number of vehicles.The addition of each vehicle to the urban environment means that a new source of pollution has been added to the city [82].The aforementioned problems, along with others such as global warming and scarcity of water resources, emphasize the importance of predicting the trend of urban development and the possible damage to the natural environments, so that managers can make preventive decisions if necessary.

Conclusions
Understanding the change in the spatial pattern of land-cover and urban growth dynamics of any area over time is important for effective land management and sustainable urban planning.Here we have described the spatiotemporal pattern of LULC and the urban expansion scenario of Jhapa district of Nepal from 1989 to 2016 with the help of multi-date images and predicted future change by 2026 and 2036 using an integrated CA-Markov model.First, the analysis explored linear urban expansion, continuous decline of cultivable land, minimum change in forest cover and shrubs, non-linear transformations of water bodies and sand and slight fluctuations in tea farming areas.Urban areas, which accounted for a total of 12.35 km 2 in 1989, increased to 13.7 km 2 in 1996, 27.12 km 2 in 2001, 37.69 km 2 in 2006, 58.65 km 2 in 2011 and 70.91 km 2 in 2016.This rapid urban expansion and decline of fertile agriculture land occurred mainly after 2001.Second, we explored the urban extent and orientation using ring-based analysis to show the spatiotemporal pattern of urban expansion from the urban center outwards.The rate of urban expansion has not been uniform throughout the district.Higher levels of urbanization have occurred mostly on former agricultural land near road networks and in the areas which provide quick access from the mid-hills.Urbanization was highly concentrated along the East-West and Mechi highway.Finally, using an integrated CA-Markov model, we showed that the urban area will increase to 130.12 km 2 in 2026 and 167.57km 2 in 2036 while cultivated land will decline to 1085.99 and 993.61 km 2 in the respective years.All other land-use classes will experience some increase.
Because massive urban expansion is ongoing at the cost of prime cultivated land and is likely to continue in future decades, food insecurity and environmental disequilibrium are probable.Developing and implementing proper urban plans for the protection of high-productive agricultural land is urgently required.Careful urban planning ensuring the preservation of cropland, green outskirts and open spaces is essential to create a resilient urban environment and sustainable development.Adopting preventive measures to dissuade high-volume immigration through various migration controls and redirection schemes should be the foremost priority of the government.The CA-Markov model is recommended as an appropriate tool for further research of the complex nature of urban and LULC.Overall, our research should serve as an important benchmark for urban planners, policy makers and other researchers.

Figure 1 .
Figure 1.Location and elevation model of study area.

Figure 1 .
Figure 1.Location and elevation model of study area.

Figure 2 .
Figure 2. Detailed framework of the study.

Figure 2 .
Figure 2. Detailed framework of the study.

Figure 4 .
Figure 4. Trend of LULC change in study area.

Figure 4 .
Figure 4. Trend of LULC change in study area.

Figure 6 .
Figure 6.Buffer partitioning of study area.

Figure 7 .
Figure 7. Urban land measured from urban centers outwards.
3.5.Markov Model 3.5.1.Analysis of Transition Matrix The transition potential matrix was computed based on LULC conditions during the periods of 1996-2006, 2006-2016 and 1996-2016 to show how each land type was projected to change.Data elements lying on the diameter of the matrix indicate that a phenomenon will probably remain constant over time, while off-diagonal data indicate that various phenomena are likely to be

Figure 7 .
Figure 7. Urban land measured from urban centers outwards.
[81].Moreover, population increase and urban development naturally result in an increase in the number of vehicles.The addition of each vehicle to the urban environment means that a new source of pollution has been added to the city [82].The aforementioned problems, along with others such as global warming and scarcity of water resources, emphasize the importance of predicting the trend of urban development and the possible damage to the natural environments, so

Table 2 .
Extracted weights based on analytic hierarchy process (AHP) and fuzzy standardization for urban areas.
.70 km 2 in 1996, 27.12 km 2 in 2001, 37.09 km 2 in 2006, 58.65 km 2 in 2011 and 70.90 km 2 in 2016.The Central Bureau of Statistics (CBS) 2014 report and forest cover from 138.57 to 141.23 km 2 from 1996 to 2001.Cultivated land experienced a decrease from 1219.26 km 2 in 2001 to 1215.43 km 2 in 2006 to 1183.34 km 2 in 2011 to 1181.26 km 2 in 2016.Overall, cultivated land decreased by 78.07 km 2 while urban area increased by 58.56 km 2

Table 5 .
Magnitude of land-use/land-cover change (km 2 ).The tea farming area of 22.25 km 2 area in 1989 had decreased to 20.38 km 2 by 1996; however, it expanded to 23.10 km 2 in 2001, 26.58 km 2 in 2006 and 28.95 km 2 in 2011 before experiencing a sharp decline to 22.52 km 2 in 2016, likely due in part because of a lack of interest in tea farming and in part to escalation of the real estate business.The tea farming area of 22.25 km 2 area in 1989 had decreased to 20.38 km 2 by 1996; however, it expanded to 23.10 km 2 in 2001, 26.58 km 2 in 2006 and 28.95 km 2 in 2011 before experiencing a sharp decline to 22.52 km 2 in 2016, likely due in part because of a lack of interest in tea farming and in part to escalation of the real estate business.