Quantifying the Spatiotemporal Pattern of Urban Expansion and Hazard and Risk Area Identification in the Kaski District of Nepal

The present study utilized time-series Landsat images to explore the spatiotemporal dynamics of urbanization and land use/land-cover (LULC) change in the Kaski District of Nepal from 1988 to 2016. For the specific overtime analysis of change, the LULC transition was clustered into six time periods: 1988–1996, 1996–2000, 2000–2004, 2004–2008, 2008–2013, and 2013–2016. The classification was carried out using a support vector machine (SVM) algorithm and 11 LULC categories were identified. The classified images were further used to predict LULC change scenarios for 2025 and 2035 using the hybrid cellular automata Markov chain (CA-Markov) model. Major hazard risk areas were identified using available databases, satellite images, literature surveys, and field observations. Extensive field visits were carried out for ground truth data acquisition to verify the LULC maps and identify multihazard risk areas. The overall classification accuracy of the LULC map for each year was observed to be from 85% to 93%. We explored the remarkable increase in urban/built-up areas from 24.06 km2 in 1988 to 60.74 km2 by 2016. A majority of urban/built-up areas were sourced from cultivated land. For the six time periods, totals of 91.04%, 78.68%, 75.90%, 90.44%, 92.35%, and 99.46% of the newly expanded urban land were sourced from cultivated land. Various settlements within and away from the city of Pokhara and cultivated land at the river banks were found at risk. A fragile geological setting, unstable slopes, high precipitation, dense settlement, rampant urbanization, and discrete LULC change are primarily accountable for the increased susceptibility to hazards. The predicted results showed that the urban area is likely to continue to grow by 2025 and 2035. Despite the significant transformation of LULC and the prevalence of multiple hazards, no previous studies have undertaken a long-term time-series and simulation of the LULC scenario. Updated district-level databases of urbanization and hazards related to the Kaski District were lacking. Hence, the research results will assist future researchers and planners in developing sustainable expansion policies that may ensure disaster-resilient sustainable urban development of the study area.


Introduction
Land-cover refers to the physical and biological materials over the Earth's surface [1].Land-cover change occurs as a natural process, whereas land use change occurs due to human activities [2,3].Urbanization is the proportion of a population living in urban areas [4].Rapid urbanization affects land use worldwide and creates major challenges for planners and policy makers in making cities for more sustainable [5].Hazards are single, consecutive, or collective physical events that cause the loss of lives or property, socioeconomic disruption, or environmental deterioration.Vulnerability is the condition generated by physical, social, economic, and environmental factors that influence the susceptibility of a community to the impacts of hazards [6].These phenomenon are either introduced via natural causes and/or by humans [7,8] and are globally recognized due to their widespread occurrences [9][10][11][12].Numbers of natural hazards and the loss and damages due to such events have increased globally in the past few decades [12].A total of 330 disasters were triggered by the beginning of 2014, and this increased to 376 in 2015 [13].In 2015, various disasters killed a total of 22,765 people worldwide, victimized 110.3 million, and caused loss and damage worth $70.3 billion.The period between 2010 and 2015 was characterized by the severe earthquakes, particularly in Chile, Haiti, New Zealand, Japan, and Nepal [14].The continent of Asia was the site of 62.7% of the deaths and 49.15% of the total damage caused by disasters.Seven Asian countries are in the top ten list of countries with the highest total death count [13].These are mainly developing countries with low-income and lower-middle economies [8,15]; they experience greater loss and damage due to hazards.Low income and lower-middle economies account 69.9% of the total victims.Nepal, DPR Korea, Myanmar, and Somalia have acquired the highest number of victims since 2005 [13].Nepal is listed as the 20th most disaster-prone country in the world and the 30th most flood-vulnerable country [16].Landslide, mass movements, and floods are common occurrences in hilly and mountainous regions of Nepal [17].In 2014, there were 588 disaster-related mortalities of which 294 were landslide and flood related [18].More than 4300 landslides were recorded in Nepal after the Gorkha earthquake (of magnitude 7.8) that took place on 25 April 2015 and its subsequent aftershocks.Flashfloods and landslides that took place on July 2016 killed 120 people, 18 remained missing, 151 were injured, and 2100 homes were damaged [19].
Remote sensing (RS) and geographic information system (GIS) technologies are globally applied key tools for the accurate characterization of LULC change [20][21][22][23][24][25] and monitoring and mapping multiple hazard areas [11,15,[26][27][28][29][30][31].RS and GIS-based risk-area mapping is recognized as an imperative process to obtain effective hazard risk reduction and mitigation output [9].Monitoring over-time changes requires consistency among maps, both ontologically and statistically [22], and time-series remote sensing data has made it possible to monitor changes in urban structure [32] and LULC [33,34].Application, GIS, and remote sensing along with LULC modeling provides an wider understanding of human and environmental systems and improved knowledge of future changes cities are to confront [35,36].Land change-related research works are cooperative for land managers and policy makers to predict potential future landscape changes; simulation modeling tends to function as the experimental tool that helps governments to characterize hazard-exposure consequences and their incorporation into hazard mitigation plans [37,38].LULC models are capable of simplifying the complexity of LULC, analyzing past changes and simulating future changes effectively [39,40].Urban growth modeling and simulation history dates back to the 1950s [41].It is gradually gaining a wide range of attention [35,42,43].A cellular-automata Markov (CA-Markov) model is appropriate to simulate the LULC change in multiple time scales [35,40] and determine the actual amount of change in land use classes [44][45][46][47].Furthermore, few studies has integrated LULC scenarios with hazards exposure and population [37].Consequently, exploring multiple hazard risks is a fundamental prerequisite to establishing sustainable disaster risk management plans.
The study area, the Kaski District, contains huge natural caves.Deep gorges exist from the northernmost to the southernmost part of the area where water often disappears underneath and only its sound can be heard (e.g., the Mahendra bridge area within Pokhara City) [48].The walls of the Sabche Gang consists of sedimentary rocks (mostly limestone) belonging to the Tibetan series of the Higher Himalaya [49], which primarily contain metamorphic with some sedimentary and igneous rocks [50].Similarly, the Higher Himalaya is prone to Glacial lake outburst floods GLOFs and avalanches, while the Fore Himalaya receives high rainfall and is vulnerable to geodisasters.Rock fall Land 2018, 7, 37 3 of 22 is more frequent because of the reduced soil deposits on the steep slopes.Due to the orographic effects, the south Pokhara area gets higher rainfall than other parts of the country [51].The mean annual rainfall at the Lumle and Pokhara stations is 5700 and 3967 mm, respectively [52].Various geohazards in the district are associated with the geophysical setting of the Kaski District.A total of 5000 houses in various villages of the district are at risk of hazards [53].Land collapse, land subsidence, and sinkhole hazards are most prominent within Pokhara City [54].Thulibesi phant in Armala, where 117 sinkhole spots were observed, is exposed to severe sinkhole hazards [55].Sarangakot and Kaskikot within the Phewa watershed [56] and the Rupakot area of the Rupa watershed [57] are prone to landslide hazards.The Phewa Lake catchment area in Pokhara witnessed a shallow seated landslide due to extreme rainfall in 2007 [51].The area has an earthquake hazard probable intensity of VIII MMI, [58] while several settlements and river banks are at high risk of liquefaction [8,58].Consequently, this area experiences bare rock cliffs and landslides [50].Studies [49,59] have explored the historical occurrences of mass wasting [60] and flash floods from 400 to 1100 CE.The latest hyper-concentrated flash floods occurred on 5 May 2012 in the Seti River within Pokhara City.It claimed [60,61] 40 lives, while 32 went missing; economic loss was worth NRS 85 million [62].Hazards of this kind may have low recurrence rates but leave important consequences and are a substantial threat for the growing population [60].Unplanned cities [6,63] with rampant urban sprawl within the unregulated municipal system [64], like the unplanned city centers and fringes of the Kaski District, are at high risk of hazards.Settlements within the core city have weak structures and also are extended in the hazard risk areas.This trend is predicted to continue in the future [8].
Additionally, the region has been undergoing high population growth and an unprecedented rate of urbanization in the recent decades.Rural-urban migration, high population growth increase, and the reclassification of rural areas into urban areas are significantly influencing an urban population increase [8].According to the Central Bureau of Statistics (CBS) in 2011 [52], the total number of annual immigrants in the district was 8000.The population of the district was 221,272 in 1981, and this increased to 292,945 in 1991, 380,527 in 2001, and 492,098 in 2011.The average annual growth rate as per the periodical analysis was 3.23%, 2.98%, and 2.93% during 1981-1991, 1991-2001, and 2001-2011, respectively.Based on new and local-level reconstruction, the total urban population of the district accounted for 84.16% (in Pokhara-Lekhnath) of the total district population of 2011.The new constitution of Nepal-promulgated in 2015 through the Constitution Assembly-legally and institutionally mandated the reconstruction of the administrative units within a federal state.During the reconstruction, former Village Development Committees were integrated into what was the Pokhara Sub metropolitan; this is now known as the Pokhara-Lekhnath Metropolitan City.The urban land area of the Seti watershed was 24.03 km 2 in 1990; it has increased by 60% and reached 54.20 km 2 in 2013.A total of 29.19 km 2 prime cultivated land was transformed into urban area [8].
The land resource mapping project (LRMP 1986) and the topographical data prepared by the Cadastral Survey Branch for the Survey Department in 1998 [65] are the two official databases for the Kaski District.These databases cover a limited time period and, hence, do not show the long-term LULC change of the study area.A detailed geographical study of Pokhara, Kaski was carried out by Gurung in 1965 [66].Yet, there are few studies have assessed the hazard risk scenario [54][55][56][57][58] and urbanization trends [67,68] of Pokhara City and the Seti watershed [8].However, long-term time-series studies of urbanization and the hazard risk of the entire Kaski District is lacking.Disaster-resilient sustainable development and urban planning is hindered due to the lack of reliable district-level databases.Disasters capture widespread attention only in relation to the extreme loss of lives and properties.With this in mind, this research aims to analyze the spatiotemporal dynamics of urban expansion on an annual basis from 1988 to 2016 using Landsat time-series images.It also aims to predict the urban expansion scenario by 2025 and 2035 using the CA-Markov model and identify the major hazard risk areas of the district.

Study Area
The study area (Figure 1), Kaski District in the Gandaki Zone of Nepal's Western Development Region, is geographically situated in the southern lap of the Annapurna Mountain range, the Mahabharat Great Himalayan range, and the Midlands.The district is located between 28 • 6 to 28 • 36 N and 83 • 40 to 84 • 12 E.It covers an area of 2088.25 km 2 and ranges from 450 to 8091 MASL (i.e., Annupurna, the tenth highest mountain peak of the world).
The land surface of the district can be divided into the High Himalayan Massit, the Middle Mountain, and hills.The climatic condition varies according to the elevation of the area.It is subtropical in altitudes of up to 1500 m in the southern part, temperate in altitudes of 1500-2000 m, temperate-cold in altitudes of 2000-3000 m, alpine in altitudes of 3000-4500 m, and tundra in altitudes of 4500 m and above.The region above 4500 m gets occupied by a thick snow cover year-round and has 11 integrated mountain peaks above 7000 m [69].The Seti, Mardi, Madi, and Modi rivers (and their tributaries) as well as various lakes (Phewa, Rupa, Begnas, Dipang, Maidi, etc.) are the major water sources of the district.
The southern part of the Annapurna Himal presents one of the sharpest contrasts in landscape; and the highlands in the north have the highest relief value than southern hills [66].The mountain ranges in the north, the deep gorges, Seti River, Davi's Falls, several natural caves, and the large lakes are attractions for the large numbers of domestic and global tourists; this fosters the socioeconomic activities that influence urbanization.Additional drivers of urbanization processes of the area include its plentiful water resources, suitable climatic condition, promising livelihood opportunities, and access to various public services [70].Administratively, the Kaski District is divided into one metropolitan (Pokhara-Lekhnath) and four rural municipalities (Rupa, Madi, Machhapuchhre, and Annapurna).

Study Area
The study area (Figure 1), Kaski District in the Gandaki Zone of Nepal's Western Development Region, is geographically situated in the southern lap of the Annapurna Mountain range, the Mahabharat Great Himalayan range, and the Midlands.The district is located between 28°6′ to 28°36′ N and 83°40′ to 84°12′ E. It covers an area of 2088.25 km 2 and ranges from 450 to 8091 MASL (i.e., Annupurna, the tenth highest mountain peak of the world).
The land surface of the district can be divided into the High Himalayan Massit, the Middle Mountain, and hills.The climatic condition varies according to the elevation of the area.It is subtropical in altitudes of up to 1500 m in the southern part, temperate in altitudes of 1500-2000 m, temperate-cold in altitudes of 2000-3000 m, alpine in altitudes of 3000-4500 m, and tundra in altitudes of 4500 m and above.The region above 4500 m gets occupied by a thick snow cover year-round and has 11 integrated mountain peaks above 7000 m [69].The Seti, Mardi, Madi, and Modi rivers (and their tributaries) as well as various lakes (Phewa, Rupa, Begnas, Dipang, Maidi, etc.) are the major water sources of the district.
The southern part of the Annapurna Himal presents one of the sharpest contrasts in landscape; and the highlands in the north have the highest relief value than southern hills [66].The mountain ranges in the north, the deep gorges, Seti River, Davi's Falls, several natural caves, and the large lakes are attractions for the large numbers of domestic and global tourists; this fosters the socioeconomic activities that influence urbanization.Additional drivers of urbanization processes of the area include its plentiful water resources, suitable climatic condition, promising livelihood opportunities, and access to various public services [70].Administratively, the Kaski District is divided into one metropolitan (Pokhara-Lekhnath) and four rural municipalities (Rupa, Madi, Machhapuchhre, and Annapurna).

Data
Geometrically, radiometrically, and topographically corrected surface reflectance (SR) Landsat thematic mapper (TM), enhanced thematic mapper (ETM+), and operational land image (OLI) (with Path/Row 142/40) images from 1988 to 2016 were collected from the United States Geological Survey (USGS) website (https://earthexplorer.usgs.gov/)[71] (Table 1).Landsat 5 TM and Landsat 7 ETM with SR data product was generated from the Landsat Ecosystem Disturbance Adoptive Processing System (LEDAPS) developed by NASA [72].The Landsat 8 OLI-SR data was generated from the Landsat Surface Reflectance Code (LaSRC) [73], which is commonly used in land-cover monitoring and mapping [20].Most images were clear and nearly free of clouds.However, some were covered by seasonal snowfall in northern part of the district.The Scan Line Corrector (SLC)-off data was excluded.All images were projected in UTM zone 45 of WGS 1984.The administrative boundary data at district, municipal, and village levels were acquired from the Survey Department of Nepal [74].Shuttle Rader Topographical Mission (SRTM) digital elevation model (DEM) with 30-m resolution (https:/lta.cr.usgs.gov/SRTM) was used for the verification of slope of the study area.In addition, high-resolution images from Google Earth (http://earth.google.com) of topographical maps from multiple dates published by the Survey Department of the Government of Nepal in 1998 (scale 1:25,000) [74] were used for references.Socioeconomic data was collected from the CBS [52,75].Additional datasets such as hazard risk area database for the Seti watershed was collected from Department of Mine and Geology (DMG), Government of Nepal [76], disaster reports from Ministry of Home Affairs (MoHA) for 2015 [16] and 2017 [77] and District profile (district information report) 2016 was collected from District Development Committee (DDC) office website (e.g., www.ddckaski.gov.np)[69].Furthermore, major land-cover change, road networks, and hazard area information were gathered using GPS during the data collection and field verification campaign in 2015 and 2016.

Image Processing and Analysis
Image processing was carried out in the ENVI environment.A supervised approach using a support vector machine (SVM) classifier [21] was employed for the abstraction of LULC and the change analysis.There are multiple approaches for LULC classification [21,22].In this study, we have used a nonparametric machine learning SVM algorithm, which has been applied successfully with sufficient accuracy in several studies [78][79][80].The outputs of SVMs varies depending on the choice of the Kernel function and its parameters [81].In the current study, the default work situation (DWS), a radial basis function (RBF) Kernel has been accepted during image classification.In the current version of ENVI (version 5.3), SVM classifier provides four types of Kernels: e.g., linear, polynomial, radial basis function (RBF) and sigmoid.The RBF Kernel, which works well in the most cases [82].The penalty parameter maximum value 100 was allocated during the classification.

LULC Extraction, and Analysis
A modified version of the land-cover classification scheme was used for remotely-sensed data, as recommended by Anderson et al. [83].After extensive field observation and using topographical map developed by Government of Nepal in 1998 (scale 1:25,000) [74], 11 land-use classes were identified namely urban/built-up, cultivated, forest, shrub, grass, barren land, sand, water body, swamp, ice and snow cover, and open (Table 2).The northern part of the study area gets covered by permanent and seasonal snow.Meanwhile, cloud-free images of the region during summer and autumn seasons were unavailable.Due to the snowfall in the winter, snow-cover area fluctuates according to season each year.This also determines the area of barren land, forest, and grass.Thus, to achieve more reliable outputs, LULC analysis was conducted on an annual basis.

LULC Based Transition Analysis
After the final production of an LULC maps, land change modeler (LCM) of the IDRISI/TerrSet software was taken into consideration to explore the LULC transition.For the close analysis of cultivated-to-urban land conversion, LULC statistics in six time periods : 1988-1996, 1996-2000, 2000-2004, 2004-2008, 2008-2013, and 2013-2016 were analyzed using the LCM.The computed transition matrix consisted of rows and columns of landscape categories at times t 1 and t 2 [35,84].

Identification of Multihazard and Risk Analysis
Historical events of landslides, floods, and sinkholes and their loss and damage information were collected based scientific research works and other literatures including disaster reports of Ministry of Home Affairs (MoHA) for 2015 [16] and 2017 [77].Hazard risk area database for the Seti watershed area with the scale of 1:50,000 was collected from Department of Mine and Geology (DMG), Government of Nepal [76] and District profile (district information report) 2016 [69].The field visit data, hazards risk areas database and LULC information developed by Rimal et al. [8] were also utilized for the analysis.Additional landscape or risk area observations and primary data acquisition were conducted during the field data collection in 2015 and field verification campaign in 2016.Highly risky areas in terms of flood, landslide, sinkhole, and edge fall were captured with the help of GPS.High-resolution Google Earth images were used for the further verification of hazard risk zones.Furthermore, local experts' knowledge was also taken into consideration.The collected information of risk areas was plotted, and final hazard risk layer was developed using Geographic Information System (GIS) to prepare multiple hazard risk map.The map was overlaid with the land-use map of 2016 and the predicted map of 2035 to develop the risk sensitivity land-use map of 2016 and hazard risk area of 2035.

The CA-Markov Model
In this section, the trend of land-use transformation was observed to predict transformations that would occur by 2025 and 2035.For this, a hybrid model (i.e., CA-Markov) has been widely used to effectively understand and measure urban expansion [35,85] and landscape dynamics [42].The model mainly comprises the following steps: (a) calculation of the transition area matrix using a Markovian process; (b) generation of transition potential maps; (c) model evaluation based on the Kappa index; and (d) simulation of future land-use using the CA-Markov model.Description of the CA-Markov model and its assumptions are summarized in Keshtkar and Voigt [40].In this study, the Markov model has computed transitional area matrices based on the historical land-use information from 1995-2005, 2005-2015, and 1995-2015.To produce transition potential maps of urban areas, four general agents (namely distance to main roads, water body, urban areas, and slope) were set as driving factors.The ancillary data was chosen based on similar previous studies [35,40].Fuzzy membership functions were applied to rescale drive maps into the range of 0 to 1, where 0 represents unsuitable locations and 1 represents ideal locations.In addition, an analytic hierarchy process (AHP) model was run to regulate the weight of driving forces with the use of pair-wise evaluations.The distinct weights and control points are summarized in Table 3.Then, the performance of the CA-Markov model was evaluated using actual (i.e., derived from satellite image) and predicted (i.e., simulated using the transition area matrix of 1995-2005) maps from the year 2015 based on the Kappa index.Since the CA-Markov model tends to predict LULC changes better within the same time span for which the model is calibrated [40]

Evaluation of Land-Cover Modeling
Kappa variations were applied to evaluate the model by comparing the land-use map of 2015 with the simulated map of 2015.The accuracy of the models, which was more than 80%, determined them to be potent predictive tools [40,45].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 [40].The detailed study design is presented in Figure 2.

Evaluation of Land-Cover Modeling
Kappa variations were applied to evaluate the model by comparing the land-use map of 2015 with the simulated map of 2015.The accuracy of the models, which was more than 80%, determined them to be potent predictive tools [40,45].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 [40].The detailed study design is presented in Figure 2.

Accuracy Accessment
Assessment of classification accuracy becomes highly essential for the LULC maps obtained using remote sensing technology [22,84,86].Due to the unavailability of updated land-cover data, high spatial resolution images of Google Earth from multiple dates, topographical maps published by the Survey Department of the Government of Nepal in 1998 (scale of 1:25,000) [74], land-cover maps from 1990 and 2013 developed by Rimal et al. [8] and GPS points were used as references.In this study, a stratified random sampling was implemented to select a total of 660 ground truth data for the validation of classified images.At least 60 sample points were represented for each class.Finally, accuracy assessment was based on the calculation of the overall accuracy, user's accuracy, and Land 2018, 7, 37 9 of 22 2.9.Accuracy Accessment Assessment of classification accuracy becomes highly essential for the LULC maps obtained using remote sensing technology [22,84,86].Due to the unavailability of updated land-cover data, high spatial resolution images of Google Earth from multiple dates, topographical maps published by the Survey Department of the Government of Nepal in 1998 (scale of 1:25,000) [74], land-cover maps from 1990 and 2013 developed by Rimal et al. [8] and GPS points were used as references.In this study, a stratified random sampling was implemented to select a total of 660 ground truth data for the validation of classified images.At least 60 sample points were represented for each class.Finally, accuracy assessment was based on the calculation of the overall accuracy, user's accuracy, and producer's accuracy.

LULC Change Analysis of 1988-2016
In the study, overall accuracy, user's accuracy, and producer's accuracies of the LULC map for each year were observed in course of accuracy assessment task (Supplementary Information S1).The overall classification accuracy ranged between 85% and 93% (1988 to 2016) was observed, which refers to the suitability and reliability of remote sensing based LULC classification and modeling.The highest producer's accuracy (100%) was obtained in the swamp cover during 2013 and the lowest (74.24%) in shrub land during 1988.Consequently, the highest user's accuracy (96.67%) was acquired in urban/built up, 2013 and the lowest (76.67%) in barren land, 1989.
Exploring the spatiotemporal pattern of LULC dynamics is a strong foundation for sustainable urban planning and effective land management [35,84] As per the analysis, it was observed that urban/built-up area experienced the largest transformation during 1988-2016 (Figure 3a).Remarkable changes occurred mainly after 2000 (Supplementary Information S2).The majority of new urban areas over the period were sourced from cultivated land.During 1988-1996, the annual rate of urban growth was 1.79% (increased by 3.46 km 2 ).Of this, 3.15 km 2 (91.04%) was converted from cultivated land.Consequently, between 1996 and 2000, the annual rate of urban growth was 4% and there was a total increase of 4.41 km 2 in new urban area, of which 3.47 km 2 (78.68%) was converted from cultivated land.Between 2000 and 2004, 6.06 km 2 of new urban area was developed, and 4.6 km 2 (75.9%) of this was from cultivated land (Figure 3a).The annual rate of urban growth reached to 4.74% during this period.In the period between 2004 and 2008, with the annual growth rate of 5.16%, there was 7.85 km 2 of new urban area, and 7.1 km 2 (90.44%) of this was from cultivated land.Between 2008 and 2013, an additional 9.29 km 2 of urban area was developed, and 8.58 km 2 (92.35%) of this was converted from cultivated land alone.The growth rate dropped to just 4.05% during this period.The urban expansion rate during 2013-2016 was 3.39%.A total of 5.61 km 2 of new urban area was created, and 5.58 km 2 (99.46%) of this was transformed from cultivated land.Figure 3b presents the spatiotemporal pattern of urban/built-up areas for six periods from 1988 to 2016.and 7.1 km 2 (90.44%) of this was from cultivated land.Between 2008 and 2013, an additional 9.29 km 2 of urban area was developed, and 8.58 km 2 (92.35%) of this was converted from cultivated land alone.The growth rate dropped to just 4.05% during this period.The urban expansion rate during 2013-2016 was 3.39%.A total of 5.61 km 2 of new urban area was created, and 5.58 km 2 (99.46%) of this was transformed from cultivated land.Figure 3b presents the spatiotemporal pattern of urban/built-up areas for six periods from 1988 to 2016.Urban/built-up areas, which totaled 24.06 km 2 in 1988, increased to 60.74 km 2 in 2016 (Figure 4a).Contrarily, cultivated land gradually declined from 445.02 km 2 in 1988 to 402.55 km 2 in 2016 (Figure 4b).Despite fluctuation during different years, forest cover has gained an overall increase in the 28 years of the investigation period.Forest had occupied 841.22 km 2 in 1988, and this expanded to 887.25 km 2 by 2016 (Figure 4c).Expansion of forest cover is closely associated with privately owned as well as community based forest management programs [8].The reason behind this is the overall preservation of forest resources because of community forest management programs, and the expansion of private forests in the abandoned cultivated land.The remaining land-cover classes (shrub and water) remained almost constant over the study periods (Figure 4d).
Almost the entire northern portion of the district usually gets covered by ice and snow.However, some variation in ice/snow-cover area was found due to the variation in date and time of the captured images.This has determined the fluctuation of barren land and grass.The higher the ice/snow cover, the lower the barren land and grass and vice versa (Figure 4e).The annual distribution of sand, open field and swamp area is shown in Figure 4f.Urban/built-up areas, which totaled 24.06 km 2 in 1988, increased to 60.74 km 2 in 2016 (Figure 4a).Contrarily, cultivated land gradually declined from 445.02 km 2 in 1988 to 402.55 km 2 in 2016 (Figure 4b).Despite fluctuation during different years, forest cover has gained an overall increase in the 28 years of the investigation period.Forest had occupied 841.22 km 2 in 1988, and this expanded to 887.25 km 2 by 2016 (Figure 4c).Expansion of forest cover is closely associated with privately owned as well as community based forest management programs [8].The reason behind this is the overall preservation of forest resources because of community forest management programs, and the expansion of private forests in the abandoned cultivated land.The remaining land-cover classes (shrub and water) remained almost constant over the study periods (Figure 4d).
Almost the entire northern portion of the district usually gets covered by ice and snow.However, some variation in ice/snow-cover area was found due to the variation in date and time of the captured images.This has determined the fluctuation of barren land and grass.The higher the ice/snow cover, the lower the barren land and grass and vice versa (Figure 4e).The annual distribution of sand, open field and swamp area is shown in Figure 4f.
(shrub and water) remained almost constant over the study periods (Figure 4d).
Almost the entire northern portion of the district usually gets covered by ice and snow.However, some variation in ice/snow-cover area was found due to the variation in date and time of the captured images.This has determined the fluctuation of barren land and grass.The higher the ice/snow cover, the lower the barren land and grass and vice versa (Figure 4e).The annual distribution of sand, open field and swamp area is shown in Figure 4f.For instance, ice/snow had occupied 277.04 km 2 in 1988, while 317.36 km 2 areas was grassland.As ice/snow-covered area increased to 443.65 km 2 in 2010, grassland area dropped down to 152.53 km 2 (Supplementary Information S3).Consequently, during 2004, barren land was confined within 32.8 km 2 and 411.93 km 2 was ice/snow-covered.However, barren land increased to 80.04 km 2 in 2005 mainly because of the significant decrease of ice/snow-cover to 287.97 km 2 .LULC transition maps during 1988-1996, 1996-2000, 2000-2004, 2004-2008, 2008-2013, and 2013-2016 have been presented in Figure 5. Transition map shows that cultivated land area changed into urban/built-up as stated above.The northern part of the study area witnessed remarkable transitions between grass land and ice/snow cover.For example, during 1996-2000, ice and snow cover had occupied large part of grassland.However, during 2000-2004, remarkable areas of grass land were transformed into ice and snow cover.A similar trend was observed in the conversion of barren land and ice and snow cover area.For instance, ice/snow had occupied 277.04 km 2 in 1988, while 317.36 km 2 areas was grassland.As ice/snow-covered area increased to 443.65 km 2 in 2010, grassland area dropped down to 152.53 km 2 (Supplementary Information S3).Consequently, during 2004, barren land was confined within 32.8 km 2 and 411.93 km 2 was ice/snow-covered.However, barren land increased to 80.04 km 2 in 2005 mainly because of the significant decrease of ice/snow-cover to 287.97 km 2 .LULC transition maps during 1988-1996, 1996-2000, 2000-2004, 2004-2008, 2008-2013, and 2013-2016 have been presented in Figure 5. Transition map shows that cultivated land area changed into urban/built-up as stated above.The northern part of the study area witnessed remarkable transitions between grass land and ice/snow cover.For example, during 1996-2000, ice and snow cover had occupied large part of grassland.However, during 2000-2004, remarkable areas of grass land were transformed into ice and snow cover.A similar trend was observed in the conversion of barren land and ice and snow cover area.
mainly because of the significant decrease of ice/snow-cover to 287.97 km 2 .LULC transition maps during 1988-1996, 1996-2000, 2000-2004, 2004-2008, 2008-2013, and 2013-2016 have been presented in Figure 5. Transition map shows that cultivated land area changed into urban/built-up as stated above.The northern part of the study area witnessed remarkable transitions between grass land and ice/snow cover.For example, during 1996-2000, ice and snow cover had occupied large part of grassland.However, during 2000-2004, remarkable areas of grass land were transformed into ice and snow cover.A similar trend was observed in the conversion of barren land and ice and snow cover area.

Multihazard Risk
Study area exposed to several hazards.The main hazards and its impacted areas are summarized below (Figure 6a-d).

Multihazard Risk
Study area exposed to several hazards.The main hazards and its impacted areas are summarized below (Figure 6a-d).

Flood Hazard
The settlements and farm lands at the river banks were found to be at a particular risk of flood hazards.These farm lands and settlements local to the Seti, Bhunge and Phusre rivers, tributaries, and streams are at more high risk of monsoonal flooding.According to the Nepal Disaster Risk Reduction Report for the Kaski District [77], eight people were killed while seven remained missing,

Flood Hazard
The settlements and farm lands at the river banks were found to be at a particular risk of flood hazards.These farm lands and settlements local to the Seti, Bhunge and Phusre rivers, tributaries, and streams are at more high risk of monsoonal flooding.According to the Nepal Disaster Risk Reduction Report for the Kaski District [77], eight people were killed while seven remained missing, three were injured, and 41 households were affected by flood hazard events in the Pokhara, Machhapuchre, Sildjure, Lekhnath, Dhikurpokhari, Ghandruk, Lamachaur areas of the Kaski District from 2011 to 2017.This was along with an economic loss of NRS 7 million.The northern part of the Seti watershed in the district is highly prone to avalanches and flash floods.However, the discussion of hazard type is beyond the scope this study.The estimated annual soil erosion rate of the Phewa watershed alone is 142,959 tons/year which is 12.2 tons/ha/year.

Edge Fall and Landslides
There is high downpour every monsoon season resulting in severe landslides in several parts.Edge fall can be witnessed in the banks of various rivers and streams.Between 2011 and 2017, landslides claimed 50 lives; three were reported missing and 30 injured while 96 households were affected.Estimated loss was worth NRS 1,450,000.Pokhara, Lekhnath, Kalika, Saimarang, Dhikurpokhari, Bhaudaure Tamagi, Lumle, Hemja, Sarangkot, Salyan, and Thumki were mainly impacted by the landslide events [77].Major landslides and edge fall risk areas of Kaski District [69] are presented in Table 4.

Land Subsidence and Sinkhole
Sinkhole and land subsidence are the most common and sensitive hazard types prevalent in the district.Random land-use practices and fragile geology are considered the two major factors associated with these hazards.The study area lies in the Ghachok formation, which consists of sediments with limestone fragments.Water in the basin passes through the gravel layer to the dissolvable lime and, finally, underground ponds emerge causing sinkholes and land subsidence.The settlements in close proximity to the Seti riverbank are at high risk of these disasters [8,48].The Seti basin in Armala is a severely sinkhole-affected area where more than a hundred sinkhole spots can be witnessed.
Existing settlements are expanded to the hazard risk areas.A total of 16.63 km 2 (27.37%) of the urban/built-up area during 2016 is at risk of hazards, while 3.83 km 2 (6.30%) of the urban area of the region needs special attention (Figure 7a).According to the prediction analysis, this trend is likely to continue in the future since the new urban areas are spreading to the risk zones.Of the total urban/built-up area by 2035, 18.21 km 2 (25.72%) will expand at risk areas and 4.42 km 2 (6.24%) will be in the regions which require special attention.The areas which are at high risk to landslides, flood, sinkhole and edge fall within the district are presented in Figure 7b.

Land Subsidence and Sinkhole
Sinkhole and land subsidence are the most common and sensitive hazard types prevalent in the district.Random land-use practices and fragile geology are considered the two major factors associated with these hazards.The study area lies in the Ghachok formation, which consists of sediments with limestone fragments.Water in the basin passes through the gravel layer to the dissolvable lime and, finally, underground ponds emerge causing sinkholes and land subsidence.The settlements in close proximity to the Seti riverbank are at high risk of these disasters [8,48].The Seti basin in Armala is a severely sinkhole-affected area where more than a hundred sinkhole spots can be witnessed.These have displaced the surrounding settlements during 2014.Similarly, various locations within the Pokhara-Lekhnath Metropolitan City have experienced land subsidence.Pokhrel [54], Rijal [55], and Rimal et al. [8] have reported the highly susceptible areas to land subsidence.5).For example, this results show that the probability of future changes of shrub land to cultivated area from 1995 to 2005 is 12.9%.This likelihood of transition decreased reasonably to 2.2% in 2015.Table 5 shows that for both periods, sand and cultivated lands possessed the highest likelihood of changing into urban/built-up areas.In assessing the overall accuracy, using the value of Kno is better than Kstandard [87].The Kno and Kstandard values were 0.92 and 0.90, respectively.The Klocation value is a reasonable representation of the location by the model and had a value of 0.95.Thus, according to the results obtained from Kappa values, the CA-Markov model is a strong predictive tool for simulating future land-use changes.
According to the prediction analysis (Table 6), urban areas will occupy 64.12 km 2 by 2025 and 70.80 km 2 by 2035.Cultivated land area is estimated to decrease to 403.07 km 2 , and 399.46 km 2 , respectively (Table 6 and Figure 8).Forest-covered area will decline to 865.02 km 2 and 860.57km 2 while water body will decrease to 32.99 km 2 and 31.5 km 2 by 2025 and 2035, correspondingly.As shown in Figure 8, urban area will increase at the cost of cultivated land and also will extend in hazard risk areas.A total of 18.21 km 2 in urban/built-up areas are expected to extend in hazard risk areas within the Kaski District by 2035.

Discussion
The Kaski District possesses a long history of urbanization, and the trend has been quite rapid in the recent years.Due to its location between mountainous and Tarai regions, the district has become an important east-west stage point in the trans-Himalayan trade route [67,70] since the medieval age.It is also the center of religious gatherings and a permanent market place experiencing speedy urbanization and population growth.Nepal experienced a decade-long political, armed conflict since 1996 that peaked in 2001, caused migration within and outside country [88].The development activities undertaken after 2006 in the aftermath of political settlement further contributed to the acceleration of the urbanization process [70], and the Kaski was not exempt.Economic opportunities, public service accessibility, population growth and migration, globalization, the physical condition, tourism, plans and policies, the land market, and political factors are the major drivers resulting in accentuated urbanization in Pokhara [8,68].
The Himalayan region often experiences severe fatalities generated from various hazards [60].High precipitation in the mountain range and adjacent foothills generally results in extreme flooding in the area and the forelands [89].The hazard of flooding becomes stronger with the excessive rainfall and alteration of floodplains into impervious surfaces [90,91] and the modification of hydrological cycle due to land-use changes [91], which take place mostly in unregulated municipal systems [64] and urban poor areas [92].Landslide is one of the crucial geomorphic drivers of LULC change; it redistributes the sediments from mountains and hills to the gentler plains.Exploitation of forest cover, over steeping natural slopes, disruption in the hydrological system, high precipitation, massive earthquakes, and rampant land-use practices enhance the landslide susceptibility of any region [17].The threat is nourished by (1) rapid urbanization, haphazard land-use practice, and the accommodation of high populations in hazard prone areas and (2) a combination of structural poverty, decaying and sub-standard infrastructures, the concentration of economic assets and commercial/industrial activities, and a weak institutional capacity for the disaster risk management [91].Kaski is the district which received the highest precipitation in a year in the country.The district has been exposed to multiple hazards including avalanches, debris flow, glacial lake outburst floods (GLOFs), landslides, edge fall, floods, sinkholes, land subsidence, and fires common in the district.Geologically, the region is incorporated in the Pokhara formation which is

Discussion
The Kaski District possesses a long history of urbanization, and the trend has been quite rapid in the recent years.Due to its location between mountainous and Tarai regions, the district has become an important east-west stage point in the trans-Himalayan trade route [67,70] since the medieval age.It is also the center of religious gatherings and a permanent market place experiencing speedy urbanization and population growth.Nepal experienced a decade-long political, armed conflict since 1996 that peaked in 2001, caused migration within and outside country [88].The development activities undertaken after 2006 in the aftermath of political settlement further contributed to the acceleration of the urbanization process [70], and the Kaski was not exempt.Economic opportunities, public service accessibility, population growth and migration, globalization, the physical condition, tourism, plans and policies, the land market, and political factors are the major drivers resulting in accentuated urbanization in Pokhara [8,68].
The Himalayan region often experiences severe fatalities generated from various hazards [60].High precipitation in the mountain range and adjacent foothills generally results in extreme flooding in the area and the forelands [89].The hazard of flooding becomes stronger with the excessive rainfall and alteration of floodplains into impervious surfaces [90,91] and the modification of hydrological cycle due to land-use changes [91], which take place mostly in unregulated municipal systems [64] and urban poor areas [92].Landslide is one of the crucial geomorphic drivers of LULC change; it redistributes the sediments from mountains and hills to the gentler plains.Exploitation of forest cover, over steeping natural slopes, disruption in the hydrological system, high precipitation, massive earthquakes, and rampant land-use practices enhance the landslide susceptibility of any region [17].The threat is nourished by (1) rapid urbanization, haphazard land-use practice, and the accommodation of high populations in hazard prone areas and (2) a combination of structural poverty, decaying and sub-standard infrastructures, the concentration of economic assets and commercial/industrial activities, and a weak institutional capacity for the disaster risk management [91].Kaski is the district which received the highest precipitation in a year in the country.The district has been exposed to multiple hazards including avalanches, debris flow, glacial lake outburst floods (GLOFs), landslides, edge fall, floods, sinkholes, land subsidence, and fires common in the district.Geologically, the region is incorporated in the Pokhara formation which is related to the catastrophic debris and water flow from the upper glaciated cirque of the Sabche Pokhara basin.The urban centers, particularly downstream, are more susceptible to the adverse effects of hazards, Anthropogenic activities affect biodiversity and the natural environment [35,93,94]; they often trigger multiple hazard occurrences.Environmental problems also arise in the absence of well-defined plans and policies to dissuade the haphazard use of natural landscape.Although developing countries often lack strong policies to manage sustainable urbanization [84], Nepal has the 2015 National Land Use Policy [95], which focuses on land consolidation, compact settlement and dissuading land fragmentation and scatter settlements.However, its implementation is lacking in the absence of long-term strategies and political commitments [96].

Conclusions and Recommendations
This study examined the trend of spatiotemporal patterns of urbanization of the Kaski District between 1988 and 2016 using time-series Landsat images.Results predicted urban sprawl by 2025 and 2035 using the CA-Markov model and identified the most hazardous areas in terms of flood, landslide or edge fall, sinkhole, and land subsidence.Research explored the gradual expansion of urban land, which increased from 1.15% to 2.91% over the study period.This spatial trend of urbanization is expected to continue to grow in future years.As per the prediction analysis, urban land area will expand to 3.07% by 2024 and 3.39% by 2032.First, urbanization has occurred mainly at the expense of prime farm lands.Of the total 36.68 km 2 expansion in urban land area between 1988 and 2016, 32.48 km 2 (88.54%) was converted from cultivated land alone.Second, urban/built-up areas expanding into areas of high risk of multiple hazards.During 2016, a total of 27.37% of the urban/built-up area was determined to be at risk.Results predicted that the urban sprawl might trespass into core and fertile farmlands and hazard prone areas.These spatial trends of urbanization pose serious challenges to sustainable food security and a sustainable urban future of the region.The study confirmed the suitability and effectiveness of long term time-series data for spatiotemporal analysis of land-use change and for identifying multihazard risk areas.However, high-resolution satellite images are recommended for detailed monitoring of hazard risk for areas with complex topography.Hazard risk and management information needs to be disseminated widely and hazard prone zones should be monitored on a routine basis.Implementation of building codes in Pokhara City is praiseworthy and should be made compulsory for the entire Kaski District.
Filling the existing gap between the technical outputs and essential prevention and coping measures should help to establish sustainability and disaster resilience in the region.Thus, it is of prime importance at this time to not only gain a general overview of the area's multiple risks, but to also establish strong mechanisms for effective disaster risk management.Despite the sensitivity of the study area, urbanization has been rapidly gaining momentum.Hence, future research should focus on the influencing factors of urbanization in the hazard risk area of the Kaski District.

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.
, 2025 and 2035 LULC maps have been predicted using the opposite trend of the LULC changes during 2005-2015 and 1995-2015.Accordingly, the LULC map of 2015 was used as a base map to simulate LULC maps for 2025 and 2035 by calculating the transition area matrix of 2005-2015 and 1995-2015, respectively.
. The explored research outputs indicate a wide range of fluctuations over the study years in various land-cover classes.Overall changes include (a) a gradual increase in urban/built-up (UB), forest cover (FC) and open field (OF) land areas; (b) a linear decrease in cultivated (CL) and shrub land (SL); (c) high fluctuation in barren land (BL) grass land (GL) and ice and snow-cover (I & SC) land areas; and (d) almost constant sizes of sandy area (SA), water body (WB) and swamp cover (SC).

Figure 3 .
Figure 3. (a) LULC map for individual years from 1988 to 2016; (b) urban area from 1988 to 2016.

Figure 5 .
Figure 5. LULC transition map for six time periods.

Figure 5 .
Figure 5. LULC transition map for six time periods.

Figure 6 .
Figure 6.Multihazard risk map of the study area: (a) major flood area; (b) multihazard risk urban areas; (c) flood and edge fall and sinkhole risk areas; and (d) overlay of multiple hazards in district map.

Figure 6 .
Figure 6.Multihazard risk map of the study area: (a) major flood area; (b) multihazard risk urban areas; (c) flood and edge fall and sinkhole risk areas; and (d) overlay of multiple hazards in district map.

Table 4 .
Landslide hazard prone areas of Kaski District.Watershed LocationMadiSurrounding of Rabaidanda village of Mijure VDC and Pakhurikot, Ghartidanda and Puranokopakha of Bhachok VDC, almost areas of Saimarang VDC, eastern and western part of Yangjakot in Thumakodanda, Setikhola pakha in Narmajung VDC, Nangjung landslide area of Sildjure VDC, Madi river bank from Sodha to Sikles of Parche VDC, and the hill between Chippali to and Modi bank of Landrung village, upper belt of Kyumrong khola, Uri gaun and Kimche gaun of Ghandruk vdc, Bhurungdi khola area of Dangshigh vdc Development Committee (DDC) Kaski, 2016.

Figure 7
Figure 7 (a) Risk-sensitive land-use map of 2016; and (b) multihazard risk area in Kaski District.Figure 7. (a) Risk-sensitive land-use map of 2016; and (b) multihazard risk area in Kaski District.

Figure 7 .
Figure 7 (a) Risk-sensitive land-use map of 2016; and (b) multihazard risk area in Kaski District.Figure 7. (a) Risk-sensitive land-use map of 2016; and (b) multihazard risk area in Kaski District.

3. 3 .
Modeling and Validation of Land-Use Change from 2015 and Onwards Transitional probability matrices that are produced by Markov model prepare details on the probability of transition between LULC types during the periods 1995-2005 and 2005-2015 (Table
Land Cover Types DescriptionSwamp/Wetland Swamp land/and with a permanent mixture of water and marsh Open field Government protected area for the construction, playground, unused land Ice and snow-cover Perpetual/temporary snow-cover, perpetual ice/glacier lake Barren land Cliffs/ small landslide, bare rocks, other permanently abandoned land Sand Sand area, other unused land, river bank Water River, lake/pond, canal, reservoir Shrub Mix of trees (<5 m tall) and other natural covers Grass Mainly grass field-(dense coverage grass, moderate coverage grass and low coverage grass) Forest Evergreen broad leaf forest, deciduous forest, scattered forest, low density sparse forest, degraded forest Cultivated land Wet and dry croplands, orchard and dry croplands, orchards Urban/built-up Commercial areas, urban and rural settlements, industrial areas, construction areas, traffic, airports, public service areas (e.g., school, college, hospital)

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

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

Table 5 .
Transition probability matrix of LULC classes for the periods 1995-2005 and 2005-2015.

Table 6 .
Distribution of land-use categories (in km 2 ) and percentage of the categories for 2025-2035.

Table 6 .
Distribution of land-use categories (in km 2 ) and percentage of the categories for 2025-2035.