- Email: [email protected]
cover change in Kathmandu Valley, Nepal
Journal of Hydrology: Regional Studies 26 (2019) 100635
Contents lists available at ScienceDirect
Journal of Hydrology: Regional Studies journal homepage: www.elsevier.com/locate/ejrh
Alteration of groundwater recharge areas due to land use/cover change in Kathmandu Valley, Nepal
T
Suraj Lamichhane*, Narendra Man Shakya Department of Civil Engineering, Institute of Engineering, Pulchowk Campus, Lalitpur, Nepal
A R T IC LE I N F O
ABS TRA CT
Keywords: AHP CLUE-S model Groundwater recharge Kathmandu Valley Land use/cover change
Study region: Kathmandu Valley, Nepal. Study focus: The focus of this study is to project future LULC, delineate potential recharge areas, and evaluate encroachment in recharge areas due to future changes in LULC. New hydrological insights for this region: The consequences of urbanization in Kathmandu Valley (KV) have been observed in various forms such as change in runoff, groundwater recharge, water scarcity, and others. To sustainably utilize groundwater resources by ensuring adequate supply/ recharge to groundwater system, land use/cover (LULC) management is required. A set of models and tools such as the Conversion of Land Use and its Effects at Small regional extent (CLUE-S) model for future LULC projection; geographic information system (GIS) for spatial data management and analysis; analytical hierarchy process (AHP) to estimate appropriate weights for different layers that influence groundwater recharge; and in-situ field test and analysis for infiltration rate were used to achieve the objectives. Results showed that built-up area in the KV watershed is projected to change by +21.4%, agricultural land by -20.5%, and forest areas by -0.9%. between 2020 and 2050. In terms of recharge area, 6% of open land is projected to convert into impervious area every decade. The projected changes are expected to have implications in terms of depletion in groundwater levels and subsequent consequences in urban water environment, including base flows in rivers.
1. Introduction Subsurface water in various forms such as well, stone spout, and others are used by human being from history to the hitherto for various purposes such as domestic, agriculture, and industrial, among others (Gautam and Prajapati, 2014). As water is central to achieve the sustainable development goals (SDGs) and groundwater is a key source to supply water in many areas, it plays a vital role in achieving the SDGs. Groundwater use is governed by factors such as availability, accessibility, transportability, and cost-effectiveness. Reliable sources and affordable cost are the primary factors that prompt people to use the subsurface water (Shrestha et al., 2016). Groundwater, the renewable reserve of freshwater in many cases, is stored in aquifers. Groundwater is also the primary source of water in the river during lean seasons, particularly when rainfall is nominal (Belhassan, 2011). Groundwater impacts the socio-economic health of the urban areas. However, urbanization process affects the quality and quantity of the groundwater in the area (Wakode et al., 2018). Due to population growth, agricultural development and urban area extension means of LULC change, the volume, shape, and pattern of the groundwater recharge as well as surface runoff is affected (Ansari et al., 2016). Such a change can affect directly and/or indirectly the socio-economic conditions as well as eco-environment of
⁎
Corresponding author. E-mail addresses: [email protected] (S. Lamichhane), [email protected] (N.M. Shakya).
https://doi.org/10.1016/j.ejrh.2019.100635 Received 1 July 2019; Received in revised form 13 September 2019; Accepted 13 October 2019 Available online 01 November 2019 2214-5818/ © 2019 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/BY-NC-ND/4.0/).
Journal of Hydrology: Regional Studies 26 (2019) 100635
S. Lamichhane and N.M. Shakya
the surroundings in many ways (Han et al., 2017). This is true for the Kathmandu Valley (KV), one of the largest urban centers in Nepal, as well. Over four million people in this valley rely on groundwater as main source of water supply. Due to limited surface water resources, groundwater has been a natural choice as an alternate source for the water supply. Industries, hotels, and housing colonies in the valley are using groundwater as a safe, reliable, and cost effective source of water since the 1980s (Pandey et al., 2010). Kathmandu Upatyaka Khanepani Limited (KUKL), the major supplier of water throughout the valley, is fulfilling only 19% of the demand in the dry season and 31% in wet seasons (Thapa et al., 2018). This clearly depicts that subsurface water is being used by public to fulfil the deficit. Groundwater contribution to the river flow is estimated at the range of 30%–74% depending on the topography, soil type, climatic factor, among others (Pandey et al., 2010). However, rapid urbanization of the catchment is resulting increase in runoff and decrease in groundwater recharge (Thapa et al., 2018). The change in surface recharge will ultimately affects groundwater budget. Earlier studies have estimated groundwater recharge is the valley in the range of 4.6–14.6 Million-Cubic-Meters (MCM) per year (Pandey and Kazama, 2014) (Shrestha et al., 2017). However, groundwater level is depleting in the valley (Pandey et al., 2010) due to the combination of land use/cover (LULC) change and excessive pumping. To regulate either of the drivers, we need a scientific understanding of the extent of changes and associated impacts. Kaliraj et al. (2014) has reported that plain areas of the valley, especially, foothills, are the most suitable recharge areas, however, they are under rapid urbanization. Therefore, it is imperative to understand the extent of encroachment of recharge areas due to LULC change. It requires projection of future LULC change using a suitable model and translating the LULC change into potential groundwater recharge areas. Different types of models are available for LULC projection. They include, but not limited to, CLUE-S (Verburg et al., 2002a,b), Dinamica EGO (Soares-Filho et al., 2002), GEOMOD (Praskievicz and Chang, 2011), Landuse Sim (Pratomoatmojo, 2018) and Land Change Modeler (LCM) (Kim, 2010). The Conversion of Land Use and its Effects (CLUE) is a dynamic model which can simulate the LULC conversion with the reference of driving forces (Veldkamp and Fresco, 1996). In a small region or scale, CLUE output is not properly valid because in the local level there are many driving forces that are not built in CLUE, therefore, the Conversion of Land Use and its Effects at Small regional extent (CLUE-S) model was developed for the small extent analysis (Verburg et al., 2002a,b). CLUE-S model has been commonly used for the simulation of LULC in regional and small scale (Verburg et al., 2002a,b; Zheng et al., 2015; Zhou et al., 2016). To further translate LULC data for various scenarios into potential groundwater recharge areas, Geographic Information System (GIS) can be used as a tool, as it is capable to handle analysis of spatial data (Chowdhury et al.;, 2009; Kaliraj et al., 2014; Jhariya et al., 2016) Different thematic layers that influence groundwater recharge such as drainage density, LULC, geology, rainfall, and topography, among others, are generally used for delineating groundwater recharge areas (Chowdhury et al.;, 2009; Singh et al., 2018). The studies have integrated remote sensing and GIS techniques for the purpose. Depending upon relative importance of one layer against another, differential weights can be assigned to the thematic layers (Chaudhary et al., 2016). Analytical hierarchical process (AHP), a form of multi-criteria decision making (MCDM) techniques, can be used for deriving differential weights to the layers (Singh et al., 2018). Earlier studies have also attempted to delineate groundwater recharge areas in the valley. For example, JICA (1990) divided the KV into three zones (i.e., northern, central, and southern groundwater districts) for the purpose of groundwater development and management, and also delineated some pockets as potential recharge areas. Later on Thapa et al. (2017) also classified the five different area from the output of model. However, none of them attempted to analyze future LULC scenario and associated impact on recharge areas. This study therefore aims to evaluate the impact of LULC change in potential groundwater recharge areas in the KV by answering following research questions – i) what are projected changes in future LULC under different scenarios? ii) What are potential areas for groundwater recharge under current condition? iii) how and where the potential recharge areas are projected to encroach due to LULC change?
2. Study area Bagmati river watershed area is spread in 613 km2, and it is located between 27032'13″ and 27049'10″N and 85011'31″ and 85 31'38″E as shown in Fig. 1. It is bowl shaped valley with the elevation ranging from 1212 to 2762 m above the mean sea level. Bagmati is the major river that flows north to south within the basin originated from Shivapuri hill. Bisnumati, Hanumante, Dhobi Khola, Manohara, Balkhu, and Nakhhu are the major tributaries of the Bagmati river. Kathmandu Valley is enclosed by the Mahabharat mountain range surrounded by four hills namely Phulchowki (2762 m) in the southeast, Chandragiri/Champadevi in the southwest, Shivapuri (2722 m) in the northwest, and Nagarkot in the northeast, formerly known as the forts of the valley. Average mean monthly precipitation ranges from 4.2 mm in December to 402.1 mm in July. The total average annual rainfall of the basin is 1533 mm/yr and more than 80% of total rainfall is received during the monsoon period (June -September). The valley climate is semi-tropical and the monthly average maximum and minimum temperatures vary from 29.8 0C to 3.4 0C and average humidity of the basin is 75% (DHM, 2015). In dry season, there is minimum flow and high flow during monsoon. Within the basin, the average annual runoff of the Bagmati river at the Khokana station (Valley outlet point St. No 550.05) located at the southern tip of the valley is 456.01 million m3/year. Administratively, Kathmandu valley represents the three district Kathmandu, Lalitpur and Bhaktapur where five major cities Kathmandu, Lalitpur, Bhaktapur, Kritipur, and Thimi are situated. The outer periphery (hilly area) of KV is covered by mixed forest and the peri-urban areas are composed of a mixture of agricultural and built-up lands, and the central core of the valley is covered by built-up area (Thapa and Murayama, 2009). 0
2
Journal of Hydrology: Regional Studies 26 (2019) 100635
S. Lamichhane and N.M. Shakya
Fig. 1. Location of the study area.
3. Methodology Fig. 2 depicts the methodological framework adopted in this study. LULC map was generated for baseline case using remote sensing products. Future LULC was projected based on CLUE-S model under different projection scenarios. A suitable set of indicators/layers for delineating groundwater recharge areas were identified based on literature review, expert consultation, and data availability. Differential weights of the layers were estimated based on AHP method. Recharge areas delineated from the analysis of layers were evaluated based on infiltration rate map prepared as a part of this study. Based on these analyses, projected change in future LULC, current recharge areas, and encroachment in recharge areas were estimated and mapped. Detailed methods are elaborated in following sub-sections. 3.1. Watershed delineation The outlet of the Bagmati river basin was considered at Katuwal Daha (850 16, 41.96′' East and 270 34′ 54.6′' North) to incorporate the all urban and peri-urban areas of the KV. Watershed boundary was defined by considering the demographic area of the KV using the 30 m ASTER-GDEM (Global Digital Elevation Map). The watershed area was calculated as 613 km2.The entire basin land is divided into five categories: agricultural, built up, forest, water, and restricted area (like airport, park, governmental office area, cultural area etc.). 3.2. Data collection and processing Land sat image from 2010 to 2018 (USGS land sat 5 to 8) of the study area was downloaded for the analysis. By using GIS and land sat image, the LULC map from 2010 to 2018 were prepared. The 2010 LULC data awas downloaded from International Center for Integrated Mountain Development (ICIMOD) and the LULC was validated from ICIMOD (2010) data. Biophysical and socio-economic driving forces that influence the LULC change were selected for analysis. Topographical factors such as elevation, slope, aspect, etc. were taken from the 30 m ASTER global digital elevation map (GDEM) data. Human-geographic driving forces like river network, road network, location of settlement, etc. were taken from the GDEM and GIS database. Similarly, socio-economic data involving population was taken from CBS (2012). The geological factors which influence the LULC change like a geological formation, soil type, soil texture was taken from the geological map of Nepal. Hydro-meteorological data of the KV was taken from Department of Hydrology and Meteorology (DHM, 2015) to input in the analysis. 3
Journal of Hydrology: Regional Studies 26 (2019) 100635
S. Lamichhane and N.M. Shakya
Fig. 2. Overall analysis framework. Notes: GW:- Groundwater; AHP:- Analytical hierarchical process; MCDM:- multi-criteria decision making ; LULC:- Land use/cover; GIS:- Geographic Information System; O1:- Objective One
3.3. Past LULC change and validation of LULC map For the retrieved data as above, to classify the land category and prepare the LULC map, GIS was used as an analysis tool. Existing LULC map of the KV is shown in Fig. 3. After plotting the dataset of 2010–2018, the LULC trend of each land use type was outlined. ICIMOD also prepared the LULC map of Nepal in 2010. In this map, there are nine LULC type were identified, viz. Agricultural, Barren, Grass, River, Built Up, Needle-leaved close forest, Needle-leaved open forest, Broad-leaved close forest, Broad-leaved open forest. For the validation purpose, all the similar type of LULC was categorized in the similar group of ICIMOD data and the total LULC type was changed into five categories. The generated data of the study area was validated from the ICIMOD (2010) data by using Kappa statistical analysis (K). The Kappa values vary from 0 to 1. The maximum value (1) represents the most perfect agreement between two maps and zero means no level of agreement (Cohen, 1960). The K value can be represented by Eqn. (1) as follows:
Pr (a) − Pr (e ) ⎤ K =⎡ ⎢ 1 − pr (e ) ⎦ ⎥ ⎣
(1)
Where, Pr(a) is the observed relative agreement among all rasters and Pr(e) is the hypothetical probability of a chance of agreement. For the purpose of validation, five LULC categories were defined. Most of the forest area in the valley is situated in the peripheral 4
Journal of Hydrology: Regional Studies 26 (2019) 100635
S. Lamichhane and N.M. Shakya
Fig. 3. (a) 2010 LULC map of the study area; (b) ICIMOD, 2010 LULC map of the study area.
part and managed by community forest user groups (CFUGs). Therefore, from the past LULC data, no significant change in the forest coverage area was detected. This study therefore considered low conversion area. Areas such as Airport, Parks, Religious areas, etc were considered as a restricted area due to the governmental restriction not allowing their conversion to other LULC type. Due to the encroachment of the liner waterway of the river and water body of the KV cannot pass the peak runoff and suffer in numerous flood in the urban area (Shrestha, 2015). Currently, Bagmati Integrated River Basin project and River Road Corridor project is under implementation for the management of the waterway and river corridor. In this study, it is assumed that the LULC of the water body is not changing. LULC demand projection scenario was basically estimated from the increase in built up in map and rate of increase in population in the KV. For the sensitivity analysis of LULC demand and LULC change of the study area, five types of LULC demand scenario were analyzed based on built area and population increase rate: a) Normal land use change scenario b) Double land use change scenario c) Half of land use change scenario d) Population growth rate scenario e) Half of population growth rate scenario. 3.4. Future LULC projection The CLUE-S model was selected for future LULC projection. It simulates LULC based on an empirical analysis of location suitability combined with the dynamic simulation of competition and interactions between the spatial and temporal dynamics of LULC systems (Verburg et al., 2002a,b). The model is separated into two modules, namely, non-spatial demand module and spatially explicit allocation module. The demand module examines the demand for various LULC through past LULC trends or scenario-based LULC and then translates the demand for LULC for application by the spatial allocation module (Verburg, 2010). Ten driving forces (geology, soil type, elevation, slope, population density, distance from market, etc)were considered for the future LULC projection and the dependency of the each driving forces was analysis based on ROC value (Zhou et al., 2016). The CLUE-S uses the logistic regression model (Eq. 2) for the identification of driving factors and their probability of the LULC change. Ordinary logistic regression model, as used in other studies (e.g., Verburg et al., 2002a,b; Chen et al., 2010; Zhou et al., 2016; Shrestha et al., 2018), was used for the purpose. First, the study area was divided into grid units and the response variable (spatial distribution of land use type) was in a binary number (1 indicated the transition occur and 0 indicated that this grid did not exist and explained) (Verburg et al., 2002a,b). Then, the logistic regression model of the probability of the suitability of the LULC was analyzed based on the SPSS tools and goodness of fit of the analysis is checked by the receiver operating curve (ROC) value (Pontius and Schneider, 2001).
pi ⎞ log ⎛⎜ ⎟ = β + β x1, i -------+ β x n, i 0 1 n ⎝ 1 − pi ⎠
(2)
Where pi is the probability of specified LULC type i transition in the grid, and, βi is a coefficient to be estimated for each explanatory variable x n, i . The probability of each grid is calculated to generate the suitability map of different LULC. The model then combines the initial LULC map, suitability map and LULC transition rules to calculate total probability (Eq. 3). TPROPi,u = Pi,u + ELASu + ITERu
(3) 5
Journal of Hydrology: Regional Studies 26 (2019) 100635
S. Lamichhane and N.M. Shakya
Table 1 Indicators/Layers selected for delineating potential groundwater recharge areas. S.N.
Layer
Relation
Reference
1
Slope
2 3
River Distance Geology
Rukundo and Doğan (2019); Ansari et al. (2016) Rukundo and Doğan (2019) Rukundo and Doğan (2019)
4
Land Use/Cover (LULC)
Steeper slopes results in low contact period with the land surface and reduces infiltration and groundwater recharge Farther from the river corridor, lesser will be the infiltration and recharge Porous media (e.g., sand, gravel, etc.) has high recharge potential compared to non or less porous. Agriculture and forest areas have high recharge capacity compare the built-up areas
5
Precipitation
6
Aspect
7
Elevation
8
Population Density
9
Market Distance
10
Road Distance
Longer duration, more frequent, and less intense precipitation has higher recharge potential. East and south facing areas are drier and have less recharge potential compared to other aspects Flat/plan areas has may provide opportunity for pondage and therefore higher potential for recharge Densely populated areas are more impervious and therefore may have lesser recharge potentials Farther the distance from market centre, higher will be potential for recharge, as they are relatively less developed Farther the distance from road centers, higher will be the possibility for recharge
Rukundo and Doğan (2019); Ansari et al. (2016) Rukundo and Doğan (2019); Bashir et al. (2008) Bashir et al. (2008) Bashir et al. (2008) Bashir et al. (2008) Bashir et al. (2008) Bashir et al. (2008)
Where, TPROPi,u is the total probability of cell i is suitable for land use u. Pi,u is the spatial distribution probability obtained by logistic regression; ELASu is the transition elasticity of land use u and it varies from 0 (lesser stability of the corresponding land use) to 1 (higher stability of the corresponding land use), and ITERu is the iteration variable of land use u. From the CLUE-S model was generate the probability map and land use map. The result from the CLUE-S was found in the ASCII file and easily analysis from GIS tools. 3.5. Delineation of groundwater recharge area A set of 10 indicators (Table 1), covering three broad categories (Fig. 4), were selected based on literature review, their logical ink to recharge potential, expert consultation, data availability, and statistical test (e.g., ROC test). AHP method, as elaborated in other studies (e.g.,Chaudhary et al., 2016; Kaliraj et al., 2014; Singh and Nachtnebel, 2016), was adopted to estimate weights to the layers. The weights factor was obtained by pair-wise comparison of each layer against rest of the layers (Saaty, 2004). 3.5.1. Data standardization For using of each thematic layer of the study, each thematic layer grids have a different unique value and the standardization process converts them into the dimensionless same value of the layers. The standardization Eq. (4) and Eq. (5) are given below and obtained value of the study area was "larger the better" and "smaller the worst" of the alternative value of the grids in each thematic layer (Pei-Yue et al., 2010). Standardization used for larger the better:
yi =
xi − xi ( min) xi (max ) − xi (min)
(4)
Standardization used for larger the better:
yi =
x (max) − xi xi (max ) − xi (min)
(5)
Where, yi is the standardized value of the thematic grid and i is the index of thematic grid and xi , xi (max ) , xi (min) are the original, maximum and minimum values of the thematic grids respectively. 3.5.2. Calculation of the groundwater potential recharge area (GWPRA) index All the thematic layers were then aggregated in GIS to get a single score of potential recharge area (Eq. 6) (Malczewski, 1999). Sum of the product of each weights and its corresponding grid value gave the high potential recharge zone (Jhariya et al., 2016). All the assigned value in each layer grid and the weights factor of each layer computed the theoretical potential recharge areas. The higher value of the grid has higher potential and vice-versa. m
GWPRA =
n
∑ ∑ (wj × xi) (6)
j=1 i=1 th
Where, GWPRA is the groundwater potential recharge area of the study area; x i is the normalized weight of the i class of the thematic layer; wj is the weight derived from AHP of the jth thematic layer; and m is the total number of thematic layer and n is the total number of class in thematic layer. 6
Journal of Hydrology: Regional Studies 26 (2019) 100635
S. Lamichhane and N.M. Shakya
Fig. 4. Methodology for delineating groundwater potential.
The obtained recharge value from the analysis was categorized in the six different zone. The dimensionless output value was analysis by using the GIS application. GWPRA was classified as per Table 2 (Singh et al., 2018). 3.6. Evaluation of groundwater recharge map Theoretical recharge potential estimated as described in section 3.5 were further evaluated based on infiltration data collected from the field. In situ field test with double ring infiltration method was used to estimate information at different locations within the KV because double ring infiltrometer is likely to give better performance (Verbist et al., 2010). Infiltration rate at 86 locations as shown in Fig. 7 (c), with wide spatial coverage, was estimated from field test. By spatial interpolation, a raster map of infiltration rate was developed. Those two maps were compared using Kappa statistical analysis (K) for the evaluation. The Kappa values vary from 0 to 1, with one indicating the most perfect agreement between two maps (Cohen, 1960). After proper validation, baseline recharge area map was finalized.a Table 2 Groundwater recharge potential value. S.N.
Very Low
Low
Moderate low
Moderate high
High
Very High
Recharge potential value
2.1 - 3.0
3.1 - 4.0
4.1- 5.0
5.1 - 6.0
6.1 -7.0
7.1 - 8.0
7
Journal of Hydrology: Regional Studies 26 (2019) 100635
S. Lamichhane and N.M. Shakya
Fig. 5. (a) 2018 LULC map of the study area; (b) Comparative diagram of LULC.
3.7. Projection and analysis of encroachment of recharge area The validated pair-wise comparison of weights factor of each layer is used to generate the potential future groundwater recharge area by using the future LULC map generated by CLUE-S model. Then the future potential recharge areas were compared with baseline map to estimate encroachment in recharge areas by using GIS tools.
4. Results and discussion 4.1. Validation of satellite-driven LULC data (2010 data) LULC data was downloaded from the USGS website in Land sat image and processed by using Arc-GIS tools. The images in Figs. 3(a) and 5 (a) show LULC maps from 2010 and 2018. For the purpose of statistical data validation, the observed map of 2010 was overlapped into the ICIMOD, 2010 as shown in Fig. 3(b). Overall, the LULC map generated in this study was in good agreement with the one from ICIMOD (2010). With this validation, the Kappa coefficient (K) value 0.68 was obtained from the analysis. The result showed that the generated LULC of 2010 and the corresponding process had-considerable accuracy. The results presented in Table 3 show that the built-up area increased by 4.96% whereas the agriculture area decreased by 6.51% of the total land between 2010–2018. It can be noticed that due to the 2015 earthquake, 2016 built up area trend partially decreased but again rapidly increased during years 2017 and 2018. Due to the formation of community forest, user’s group encroachment of forest land is drastically decreased. As a result of rapid urbanization post 2015 earthquake, forest areas has been decreased partially, as reflected in the form of decrease in forest cover from 2016. Furthermore, increase in built-up area is reducing the open agriculture and barren land, thus, adversely affecting groundwater. For the future analysis, water body and the restricted area land use was assumed that there is no change in the future.
Table 3 Comparison matrix of LULC between.2010–2018. LULC
Land Use/Cover 2010 (%)
Land Use/Cover 2018 (%)
Agriculture Land Built-Up Forest Restricted Area Water Body Grand Total
Agriculture Land
Built-Up
Forest
Restricted Area
Water-Body
Grand Total
32.06 4.00 6.31 0.02 0.02 42.40
8.07 12.50 1.26 0.09 0.01 21.93
8.73 0.34 24.87 – – 33.94
– 0.11 0.00 1.31 – 1.42
0.04 0.02 0.02 – 0.21 0.30
48.91 16.97 32.46 1.42 0.25 100.00
8
Journal of Hydrology: Regional Studies 26 (2019) 100635
S. Lamichhane and N.M. Shakya
Table 4 ROC value of each LULC. LULC Type
Forest
Water Body
Built up
Agriculture
Restricted Area
ROC
0.92
0.87
0.91
0.83
0.94
4.2. Evaluation of CLUE-S model For the LULC change projection, ten different driving forces were considered as an input values which is already discussed in previous section. The ROC value validates the logistic regression model. If the ROC value is less than 0.7, the accuracy of the model is low and if the value is less or greater than 0.9 represents the accuracy is credible and high precision respectively (Pontius and Schneider, 2001). The ROC value of each LULC was greater than 0.83 and it shows that there was a strong relationship between the driving forces and the LULC type. These data indicated that the regression results are excellent and the regression coefficient from the driving forces was validated parameter for each LULC change is shown in Table 4. For the sensitivity analysis of the LULC demand of the study area was analyzed by five demand scenarios and it already describe in previous section. Normal land use demand scenario was quantified based on the LULC change from 2010 to 2018 and the population growth analysis was based on census population growth rate from 1971 to 2011 of the KV. Population growth rate demand projection is overestimated and non-realistic, and half of the normal land use change growth was underestimated. However, normal land use change is more practical or realistic compared to the other as shown in Fig. 6. According to the five different scenarios described in the previous section, the CLUE-S model was used to project the LULC change in the KV. Among these scenarios, normal land use change scenario was taken for further analysis. As water body and restricted areas are not expected to change, in all the scenarios, LULC change was considered in agriculture and built-up areas only. Simulated data from years 2010 to 2018 was used for the validation of the map and the forecasted data for year 2020 to 2050 is used for the analysis. The generated map from the CLUE-S model was validated through the Kappa (K) coefficient. The Kappa coefficient of years 2014, 2016 and 2018 was 0.75, 0.64 and 0.62 respectively. It showed that the analysis is of considerable accuracy. The model assumed that the urban development is concentrated into the existing market and dense village but the practically development is scattered due to high-cost of land and less open space in the existing market. From the simulation of CLUE-S, from 2020 to 2030, forest, water body and the restricted areas are not projected to change significantly. But the agricultural land is projected to decrease dramatically from 40.11% to 33.44% within that decade. Average annual encroachment of the fertile land is projected to be 0.67% (4.07 km2) by mostly built up area due to the construction of
Fig. 6. Future LULC change simulated by CLUE-S model (a) Normal LULC Growth rate scenario, (b) Population growth rate scenario. 9
Journal of Hydrology: Regional Studies 26 (2019) 100635
S. Lamichhane and N.M. Shakya
Table 5 Comparison matrix of LULC between 2020 and 2030, 2030 and 2040. LULC
Land Use/Cover 2030 (%) Agriculture Land
Land Use/Cover 2020 (%)
Agriculture Land Built-Up Forest Restricted Area Water Body Grand Total
33.4 – – – – 33.4 Land Use/Cover 2040 Agriculture Land 26.2 – – – – 26.2
LULC Land Use/Cover 2030 (%)
Built-Up
Agriculture Land Built-Up Forest Restricted Area Water Body Grand Total
6.2 23.4 0.7 – – 30.3 (%) Built-Up 7.2 30.3 0.5 – – 38.0
Forest
Restricted Area
Water-Body
Grand Total
0.5 – 34. – – 34.5
– – – 1.5 – 1.5
– – – – 0.3 0.3
40.1 23.4 34.7 1.5 0.3 100.0
Forest – – 34.0 – – 34.0
Restricted Area – – – 1.5 – 1.5
Water-Body – – – – 0.3 0.3
Grand Total 33.4 30.3 34.5 1.5 0.3 100.0
physical infrastructure. Meanwhile, the built area of the KV is projected to increase by 6.85% within a decade which can be clearly seen in the map and appended Table 5 and 6 as shown in Fig. 7. From the overall simulation model (2020 to 2050), agriculture and forest lands is projected to reduce by 20.4% (125.05 km2) and 0.9% (5.51 km2) of the total area respectively. The total built area is projected to increase by 21.3% (130.57 km2) of the total area. Within three decades, half of the agriculture land of the KV is projected to shift into the urban area. If the current level of urbanization continues as usual, more fertile land might convert into a constructed area. Also, the open space used as recreational area of the KV might decline. In upcoming days, the agricultural production will be decreased due to the loss of fertile land and increase the food deficiency. Meanwhile, the groundwater recharge by the encroachment of potential groundwater recharge area would lower the groundwater table. The encroachment of the open areas by the urbanization between 2010–2018 was 10.4% (63.54 km2) of the total basin area. The LULC map of the CLUE- S model (2020 to 2050) clearly shows that 21.4% (130.68 km2) of the open area is projected to change from fertile and permeable agriculture land to built-up area.
4.3. Delineation of potential recharge areas (Under current and future scenario) Ten thematic layers were used for the AHP to find the paired comparison weights factor. The grid value of each thematic layer was given by the pair-wise comparison from 1 to 9. The consistency ratio (CR) of each thematic layer was less than 8% which indicates the analysis is consistent and the generated data from each layer was suitable for the future AHP analysis. Overall ten thematic layer analysis, CR value 9.7% indicate that the pair-wise comparison was more consistent (Saaty, 2004) and the output weight age of each layer was more realistic. Comparison matrix of AHP and weights factor of each layer is as shown in Table 7. Matrix was generated from the experts’ view with the assumption that each layer influenced the rate of groundwater Table 6 Comparison matrix of LULC between 2040 and 2050, 2020 and 2050. LULC
Land Use/Cover 2040 (%)
Land Use/Cover 2050 (%)
Agriculture Land Built-Up Forest Restricted Area Water Body Grand Total
LULC
Land Use/Cover 2020 (%)
Agriculture Land
Built-Up
Forest
Restricted Area
Water-Body
Grand Total
19.7 – – – – 19.7
6.2 38.0 0.5 – – 44.7
0.3 – 33.5 – – 33.8
– – – 1.5 – 1.5
– – – – 0.3 0.3
26.2 38.0 34.0 1.5 0.3 100.0
Land Use/Cover 2050 (%)
Agriculture Land Built-Up Forest Restricted Area Water Body Grand Total
Agriculture Land
Built-Up
Forest
Restricted Area
Water-Body
Grand Total
19.7 – – – – 19.7
19.6 23.4 1.7 – – 44.7
0.8 – 33.00 – – 33.8
– – – 1.5 – 1.5
– – – – 0.3 0.3
40.1 23.4 34.7 1.5 0.3 100.0
10
Journal of Hydrology: Regional Studies 26 (2019) 100635
S. Lamichhane and N.M. Shakya
Fig. 7. (a) Land conversion map of Kathmandu Valley between 2010–2018; (b) Land conversion map of Kathmandu Valley between 2020–2050; (c) Infiltration test locations in Kathmandu Valley.
Table 7 Comparison matrix of AHP and layers weight factor. S.N.
Type
1
2
3
4
5
6
7
8
9
10
%
1 2 3 4 5 6 7 8 9 10
Slope River Distance Geology Land Use Precipitation Aspect Elevation Population Density Market Distance Road Distance
1.00 2.00 3.00 4.00 5.00 0.50 0.33 0.20 0.14 0.33
0.50 1.00 2.00 2.00 3.00 0.25 0.17 0.11 0.11 0.33
0.33 0.50 1.00 2.00 2.00 0.17 0.14 0.11 0.11 0.11
0.25 0.50 0.50 1.00 3.00 0.25 0.33 0.50 0.11 0.11
0.20 0.33 0.50 0.33 1.00 0.11 0.14 0.11 0.11 0.11
2.00 4.00 6.00 4.00 9.00 1.00 0.50 0.50 0.25 0.50
3.00 6.00 7.00 3.00 7.00 2.00 1.00 0.50 0.33 1.00
5.00 9.00 9.00 2.00 9.00 2.00 2.00 1.00 0.50 2.00
7.00 9.00 9.00 9.00 9.00 4.00 3.00 2.00 1.00 3.00
3.00 3.00 9.00 9.00 9.00 2.00 1.00 0.50 0.33 1.00
7.65 12.89 18.52 17.16 28.91 4.32 3.38 2.62 1.51 3.06
recharge. Form the analysis, precipitation (28.91%) is the main key factor for the groundwater recharge and the effect of geology (18.52%) and land use (17.16%) is relatively high. Bank of the river is mostly formed by the river deposition which has relatively high infiltration capacity and its weights factor 12.89% was also a major influencing parameter. Slope influences the contact time of runoff water and soil layer; higher slope has high velocity and less recharge and vice versa. By using GIS and AHP weights factor the potential groundwater recharge area was identified as shown in Fig. 8b. For validation of the AHP generated potential groundwater recharge area by the actual field condition, number of spouts (83 no) were identified as per the geological, hydro-meteorological, socio-economical changes location. Every variation of the geological formation, land use change in the location and change in the soil texture area is considered for the test point. The location point must be free from an anthropogenic disturbance and it must be an open space for the field test. Double ring infiltrometer (see for details: Lucke et al., 2014) test was conducted in 83 different locations of the study area. The overall location of the field test point is shown in the map (Fig. 7c). Actual field raster map was prepared by the field test infiltration rate (mm/hr) point and that was compared with the AHP generated map for the evaluation of AHP outputs. Due to lack of the information about litho-logical, moisture content, water table, soil density and so on, actual recharge rate was not accurately quantified. Due to the variable grid value unit (one infiltration value and other is unit less recharge factor), statistically, these maps could not be validated. However, from the inception of map and as shown in pervious literature (Pandey et al., 2013; Thapa et al., 2018; Adhikari et al., 2019) foothill of the mountain and the alluvial part of the river corridor were most potential area. Groundwater potential area was computed using Eq. (6) and the output of the study area was categorized into six parts as per quintile mapping in ArcGIS. Area coverage of the very low, low, moderate low, moderate high, high, and very high was 4.29%, 9.14% and 14.72%, 32.59%, 32.89%, and 7.37% respectively of the total area coverage as shown in Fig. 8. Currently, high and very high areas are still in open area in the foot hill of northern and eastern part of mountain of the valley. Central part of the valley also has the potential recharge area but now it is mostly in concrete jungle. Northern and eastern part of the valley were more groundwater recharge potential zone and the peripheral boundary of the valley are less recharge potential due to high hill and geological formation. From the comparison of the observation of field test map, AHP method map and previous literature and maps (Pandey et al., 2013; Thapa et al., 2018; Adhikari et al., 2019) potential recharge location may be a similar area. The southern and northern foothill 11
Journal of Hydrology: Regional Studies 26 (2019) 100635
S. Lamichhane and N.M. Shakya
Fig. 8. (a) Potential recharge area mapping by field test; (b) Potential recharge area mapping by AHP method.
of mountain areas are the main possible recharge location of the KV. If that potential area could be preserved and managed, and artificial recharge system will be adapted which can overcome the water security situation of the KV in the future. Form the analysis, the theoretically generated map from the different layers using MCDM techniques is realistic and useable for the location which has less groundwater recharge information like Nepal. 4.4. Change in recharge areas under projected land use/cover Rapid urbanization mostly influences the groundwater recharge. Future encroachment of recharge area was identified by using AHP weights factor and the LULC map generated from the CLUE-S model. The projected future potential recharge zones are shown in Fig. 9. From the maps shown in Fig. 9, it can be observed that every decade the open area (agriculture and forest land) of the valley is projected to encroach by the urbanization. From the data, every decade on an average of 6% of open land is projected to decrease and each year 3.66 km2 of open land will change into the urban area. The likely reason for not consistent trend across future time periods in Table 8 is due to consideration of only one factor (i.e., LULC) whereas recharge areas are affected by various other factors as well. Such decrease will lead in the reduction of groundwater recharge. Reduction in groundwater leads to the lowering of groundwater table and base flow contributions in river basin. In addition, the rate of groundwater recharge will reduce over the time primarily due to the increment of the urban population and increases the volume of the water demand. Such reduction also effects in quantity and quality of water and other urban environments. 5. Conclusion A methodological framework has been developed and implemented in this study that helps to classify the LULC type, identification of the potential recharge area, possible protection of the agriculture land for the groundwater recharge, and future possible encroachment in the groundwater recharge areas. Identification of the LULC, the future projection of LULC and encroachment of recharge area generated in this study will be useful in effective urban planning, optimal utilization of water supply demand, and water resources planning and management. The sum of the present study could be summarized in the following conclusions: a) The trend of past and present LULC analysis highlights that the agricultural land is projected to decrease by 6.51% leading to an increment of 4.96% in built-up areas. It is interesting to note that between 2010–2018, the forest area is found to be increased by 1.48% due to community forestry programs. Conversion of pervious recharge areas into impervious constructed areas would lead into the reduction of groundwater recharge as the increment of forest area does not necessarily contributes significantly into recharge. b) The predicted scenario between 2020–2050 using CLUE-S model highlights that the built-up area is projected to increase by 21.4% in the next 30 years and the agricultural and forest areas is projected to decrease by 20.5% and 0.9% respectively. The 12
Journal of Hydrology: Regional Studies 26 (2019) 100635
S. Lamichhane and N.M. Shakya
Fig. 9. (a) Potential recharge map for 2020 (b) Potential recharge map for 2050 (c) Recharge encroachment area by potential urbanization from 2020 to 2050 (d) Recharge encroachment scenario by potential urbanization from 2020 to 2050.
results outline that the fertile lands would be converted into urban areas leading to the reduction in agricultural products and productivity. Lack of irrigation water may also be a crucial case which imparts additional pressure over the groundwater. c) Precipitation, geology, land use/cover, distance from the river, slope, among others were identified to be the significant factors to contribute in delineation of potential groundwater recharge area for the AHP/MCDM technique. Notably, 7.37% of the total area was obtained to be the high recharge area, 32.89% under high, 32.59% moderately high, 14.72% moderately low, 9.14% low, and 4.29% very low recharge areas. The high and moderately high recharge areas are situated in the urban areas with greater tendency of expansion. Thus, there might be the case that due to the LULC change; more high to moderate recharge areas would be converted to urban areas leading to additional pressure on groundwater. 13
Year
2020 2030 2040 2020
S.N.
1. 2. 3. 4.
to to to to
2030 2040 2050 2050
4.14 3.15 3.42 10.71
Increase in recharge area (km2)
Table 8 Decadal encroachment area and percentage of GWPRA.
0.68 0.52 0.56 1.75
Change in percentage 40.41 43.83 39.78 124.02
Decrease in recharge area (km2)
6.62 7.81 6.51 20.31
Change in percentage
5.94 6.66 5.95 18.55
Actual decrease percentage
S. Lamichhane and N.M. Shakya
Journal of Hydrology: Regional Studies 26 (2019) 100635
14
Journal of Hydrology: Regional Studies 26 (2019) 100635
S. Lamichhane and N.M. Shakya
d) Future LULC change analysis highlighted that 6% (3.99 km2) of rechargeable land area are projected to convert to impervious land in a decade. As Kathmandu Valley is already facing water scarcity, the urban expansion would put additional pressure leading to the graver water scarcity. The application of CLUE-S, GIS and AHP/MCDM techniques are more realistic, cost-effective, simple, and reliable approaches, and scientifically sound methods for the analysis. The methodology used in the study area can also be used into other urban areas in the country as well as worldwide. The findings from the study provide valuable information for decision-makers, planners, and managers for the sustainable management of the land and water resources. The impact of LULC on groundwater and surface water would ultimately affect the river runoff pattern. To address this issue, further analyses would be carried out in the future. Declaration of Competing Interest The authors declare no conflict of interest of any kind apply to the reported work. Acknowledgments The first author acknowledges the funding provided by Institute of Engineering, Tribhuvan University for the doctoral research grant. Authors acknowledge the support provided by University of Yamanashi, Japan, NTNU, Norway, Hiroshi Ishidaira, Bhesh Raj Thapa, Dipendra Gautam, Sarita Dawadi, Saraswati Thapa, Hari Poudel, and other people who supported in various ways. References Adhikari, B.R., Shrestha, S.Das, Shakya, N.M., 2019. Future Urban Water Crisis in Mountain Regions: Example of Kathmandu Valley, Nepal. Springer, Urban Drought, pp. 169–182. Ansari, T.A., Katpatal, Y.B., Vasudeo, A.D., 2016. Spatial evaluation of impacts of increase in impervious surface area on SCS-CN and runoff in Nagpur urban watersheds. India. Arab. J. Geosci 9. https://doi.org/10.1007/s12517-016-2702-5. Bashir, A.M., Fakhry, A.A., Philip, E.L., James, W.L., Mostafa, M.S., 2008. Environmental Hydrogeology. Taylor & Francis Group, USA. Belhassan, K., 2011. Relationship between River Flow, rainfall and groundwater pumpage in Mikkes Basin (Morocco). Iran. J. Earth Sci. 3, 98–107. CBS, N., 2012. National Population and Housing Census 2011. National Report. . Chaudhary, P., Chhetri, S.K., Joshi, K.M., Shrestha, B.M., Kayastha, P., 2016. Application of an Analytic Hierarchy Process (AHP) in the GIS interface for suitable fire site selection: A case study from Kathmandu Metropolitan City. Nepal. J. Socio-Economic Plan. Sci. 53, 60–71. Chen, F., Hu, Y., Peng, X., Wang, L., 2010. Simulation of land use/cover change based on the CLUE-S model. 2010 18th International Conference on Geoinformatics. IEEE 1–5. Chowdhury, A., Jha, M.K., Chowdary, V.M., Mal, B.C., 2009. Integrated remote sensing and GIS based approach for assessing groundwater potential in West Medinipur district, West Bengal. India. Int. J. Remote Sens. 30, 231–250. Cohen, J., 1960. A coefficient of agreement for nominal scales. Educ. Psychol. Meas. 20, 37–46. DHM, 2015. Hydrological Records of Nepal, streamflow Summary, Updated Version. Department of Hydrology and Meteorology: Goverment of Nepal, Ministry of Water Resources, Kathmandu. Gautam, D., Prajapati, R.N., 2014. Drawdown and dynamics of groundwater table in Kathmandu Valley. Nepal. Open Hydrol. J. 8. Han, D., Currell, M.J., Cao, G., Hall, B., 2017. Alterations to groundwater recharge due to anthropogenic landscape change. J. Hydrol. (Amst) 554, 545–557. https:// doi.org/10.1016/j.jhydrol.2017.09.018. ICIMOD, 2010. Land Cover of Nepal 2010 [Dataset]. Available online at:. International Center for Integrated Mountain Development (ICIMOD), Kathmandu, Nepal. http://rds.icimod.org/Home/DataDetail?metadataId=9224. Jhariya, D.C., Kumar, T., Gobinath, M., Diwan, P., Kishore, N., 2016. Assessment of groundwater potential zone using remote sensing, GIS and multi criteria decision analysis techniques. J. Geol. Soc. India 88, 481–492. JICA, 1990. Groundwater Management Project in Kathmandu Valley. Nepal Water Supply Corporation, Kathmandu, Nepal. Kaliraj, S., Chandrasekar, N., Magesh, N.S., 2014. Identification of potential groundwater recharge zones in Vaigai upper basin, Tamil Nadu, using GIS-based analytical hierarchical process (AHP) technique. Arab. J. Geosci. 7, 1385–1401. Kim, O., 2010. Comparison of Two GIS based Land change modules for constructing REDD baselines in Bolivia. AAG Annu. Meet. Lucke, T., Boogaard, F., van de Ven, F., 2014. Evaluation of a new experimental test procedure to more accurately determine the surface infiltration rate of permeable pavement systems. J. Urban Plan. Transp. Res. 2, 22–35. Malczewski, J., 1999. GIS And Multicriteria Decision Analysis. John Wiley & Sons. Pandey, V.P., Chapagain, S.K., Kazama, F., 2010. Evaluation of groundwater environment of Kathmandu Valley. J. Geogr. Environ. Earth Sci. Int. 60, 1329–1342. Pandey, V.P., Kazama, F., 2014. From an open-access to a state-controlled resource: the case of groundwater in the Kathmandu Valley, Nepal. Water Int. 39, 97–112. https://doi.org/10.1080/02508060.2014.863687. Pandey, V.P., Shrestha, S., Kazama, F., 2013. A GIS-based methodology to delineate potential areas for groundwater development: a case study from Kathmandu Valley. Nepal. Appl. Water Sci. 3, 453–465. Pei-Yue, L., Hui, Q., Jian-Hua, W.U., 2010. Groundwater quality assessment based on improved water quality index in Pengyang County, Ningxia, Northwest China. J. Chem. 7, S209–S216. Pontius Jr, R.G., Schneider, L.C., 2001. Land-cover change model validation by an ROC method for the Ipswich watershed, Massachusetts, USA. Agric. Environ. Ecosyst. 85, 239–248. Praskievicz, S., Chang, H., 2011. Impacts of climate change and urban development on water resources in the Tualatin River Basin. Oregon. Ann. Assoc. Am. Geogr. 101, 249–271. https://doi.org/10.1080/00045608.2010.544934. Pratomoatmojo, N.A., 2018. LanduseSim algorithm: land use change modelling by means of cellular automata and geographic information system. IOP Conf. Ser. Earth Environ. Sci. 202. https://doi.org/10.1088/1755-1315/202/1/012020. Rukundo, E., Doğan, A., 2019. Dominant influencing factors of groundwater recharge spatial patterns in Ergene river catchment, Turkey. Water (Switzerland) 11. https://doi.org/10.3390/w11040653. Saaty, T.L., 2004. Fundamentals of the analytic network process—multiple networks with benefits, costs, opportunities and risks. J. Syst. Eng. Syst. 13, 348–379. Shrestha, S., 2015. Assessment of Bagmati River Encroachment Through Application of GIS and Remote Sensing. Pokhara University. Shrestha, S., Bhatta, B., Shrestha, M., Shrestha, P.K., 2018. Integrated assessment of the climate and landuse change impact on hydrology and water quality in the Songkhram River Basin. Thailand. Sci. Total Environ. 643, 1610–1622. https://doi.org/10.1016/j.scitotenv.2018.06.306. Shrestha, S., Kafle, R., Pandey, V.P., 2017. Evaluation of index-overlay methods for groundwater vulnerability and risk assessment in Kathmandu Valley. Nepal. Sci.
15
Journal of Hydrology: Regional Studies 26 (2019) 100635
S. Lamichhane and N.M. Shakya
Total Environ. 575, 779–790. https://doi.org/10.1016/j.scitotenv.2016.09.141. Shrestha, S., Semkuyu, D.J., Pandey, V.P., 2016. Assessment of groundwater vulnerability and risk to pollution in Kathmandu Valley. Nepal. Sci. Total Environ. 556, 23–35. Singh, L.K., Jha, M.K., Chowdary, V.M., 2018. Assessing the accuracy of GIS-based Multi-Criteria Decision Analysis approaches for mapping groundwater potential. Ecol. Indic. 91, 24–37. Singh, R.P., Nachtnebel, H.P., 2016. Analytical hierarchy process (AHP) application for reinforcement of hydropower strategy in Nepal. Renewable Sustainable Energy Rev. 55, 43–58. https://doi.org/10.1016/j.rser.2015.10.138. Soares-Filho, B.S., Cerqueira, G.C., Pennachin, C.L., 2002. DINAMICA—a stochastic cellular automata model designed to simulate the landscape dynamics in an Amazonian colonization frontier. Ecol. Modell. 154, 217–235. Thapa, B., Ishidaira, H., Pandey, V., Bhandari, T., Shakya, N., 2018. Evaluation of water security in Kathmandu valley before and after water transfer from another basin. Water 10, 224. Thapa, B.R., Ishidaira, H., Pandey, V.P., Shakya, N.M., 2017. A multi-model approach for analyzing water balance dynamics in Kathmandu Valley. Nepal. J. Hydrol. Reg. Stud. 9, 149–162. https://doi.org/10.1016/j.ejrh.2016.12.080. Thapa, R.B., Murayama, Y., 2009. Examining spatiotemporal urbanization patterns in Kathmandu Valley, Nepal: remote sensing and spatial metrics approaches. Remote Sens. (Basel) 1, 534–556. https://doi.org/10.3390/rs1030534. Veldkamp, A., Fresco, L.O., 1996. CLUE-CR: an integrated multi-scale model to simulate land use change scenarios in Costa rica. Ecol. Modell. 91, 231–248. https:// doi.org/10.1016/0304-3800(95)00158-1. Verbist, K., Torfs, S., Cornelis, W.M., Oyarzún, R., Soto, G., Gabriels, D., 2010. Comparison of single-and double-ring infiltrometer methods on stony soils. Vadose Zone J. 9, 462–475. Verburg, P., 2010. The CLUE model. Hands-on Exercises. Course Material, in: The CLUE Model. Hands-on Exercises. Course Material. Institute for Environmental Studies, University of Amsterdam. Verburg, P., Soepboer, W., Veldkamp, A., Limpiada, R., Espaldon, V., Mastura, S.S.A., 2002a. Modeling the spatial dynamics of regional land use: the CLUE-S model. Environ. Manage. 30, 391–405. Verburg, P.H., Soepboer, W., Veldkamp, A., Limpiada, R., Espaldon, V., Mastura, S.S.A., 2002b. Modeling the spatial dynamics of regional land use: the CLUE-S model. Environ. Manage. 30, 391–405. https://doi.org/10.1007/s00267-002-2630-x. Wakode, H.B., Baier, K., Jha, R., Azzam, R., 2018. Impact of urbanization on groundwater recharge and urban water balance for the city of Hyderabad, India. Int. Soil Res. Water Conserv. 6, 51–62. Zheng, H.W., Shen, G.Q., Wang, H., Hong, J., 2015. Simulating land use change in urban renewal areas: a case study in Hong Kong. Habitat Int. 46, 23–34. Zhou, R., Zhang, H., Ye, X.-Y., Wang, X.-J., Su, H.-L., 2016. The delimitation of urban growth boundaries using the CLUE-S land-use change model: Study on Xinzhuang Town, Changshu City, China. Sustainability 8, 1182.
16
Contents lists available at ScienceDirect
Journal of Hydrology: Regional Studies journal homepage: www.elsevier.com/locate/ejrh
Alteration of groundwater recharge areas due to land use/cover change in Kathmandu Valley, Nepal
T
Suraj Lamichhane*, Narendra Man Shakya Department of Civil Engineering, Institute of Engineering, Pulchowk Campus, Lalitpur, Nepal
A R T IC LE I N F O
ABS TRA CT
Keywords: AHP CLUE-S model Groundwater recharge Kathmandu Valley Land use/cover change
Study region: Kathmandu Valley, Nepal. Study focus: The focus of this study is to project future LULC, delineate potential recharge areas, and evaluate encroachment in recharge areas due to future changes in LULC. New hydrological insights for this region: The consequences of urbanization in Kathmandu Valley (KV) have been observed in various forms such as change in runoff, groundwater recharge, water scarcity, and others. To sustainably utilize groundwater resources by ensuring adequate supply/ recharge to groundwater system, land use/cover (LULC) management is required. A set of models and tools such as the Conversion of Land Use and its Effects at Small regional extent (CLUE-S) model for future LULC projection; geographic information system (GIS) for spatial data management and analysis; analytical hierarchy process (AHP) to estimate appropriate weights for different layers that influence groundwater recharge; and in-situ field test and analysis for infiltration rate were used to achieve the objectives. Results showed that built-up area in the KV watershed is projected to change by +21.4%, agricultural land by -20.5%, and forest areas by -0.9%. between 2020 and 2050. In terms of recharge area, 6% of open land is projected to convert into impervious area every decade. The projected changes are expected to have implications in terms of depletion in groundwater levels and subsequent consequences in urban water environment, including base flows in rivers.
1. Introduction Subsurface water in various forms such as well, stone spout, and others are used by human being from history to the hitherto for various purposes such as domestic, agriculture, and industrial, among others (Gautam and Prajapati, 2014). As water is central to achieve the sustainable development goals (SDGs) and groundwater is a key source to supply water in many areas, it plays a vital role in achieving the SDGs. Groundwater use is governed by factors such as availability, accessibility, transportability, and cost-effectiveness. Reliable sources and affordable cost are the primary factors that prompt people to use the subsurface water (Shrestha et al., 2016). Groundwater, the renewable reserve of freshwater in many cases, is stored in aquifers. Groundwater is also the primary source of water in the river during lean seasons, particularly when rainfall is nominal (Belhassan, 2011). Groundwater impacts the socio-economic health of the urban areas. However, urbanization process affects the quality and quantity of the groundwater in the area (Wakode et al., 2018). Due to population growth, agricultural development and urban area extension means of LULC change, the volume, shape, and pattern of the groundwater recharge as well as surface runoff is affected (Ansari et al., 2016). Such a change can affect directly and/or indirectly the socio-economic conditions as well as eco-environment of
⁎
Corresponding author. E-mail addresses: [email protected] (S. Lamichhane), [email protected] (N.M. Shakya).
https://doi.org/10.1016/j.ejrh.2019.100635 Received 1 July 2019; Received in revised form 13 September 2019; Accepted 13 October 2019 Available online 01 November 2019 2214-5818/ © 2019 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/BY-NC-ND/4.0/).
Journal of Hydrology: Regional Studies 26 (2019) 100635
S. Lamichhane and N.M. Shakya
the surroundings in many ways (Han et al., 2017). This is true for the Kathmandu Valley (KV), one of the largest urban centers in Nepal, as well. Over four million people in this valley rely on groundwater as main source of water supply. Due to limited surface water resources, groundwater has been a natural choice as an alternate source for the water supply. Industries, hotels, and housing colonies in the valley are using groundwater as a safe, reliable, and cost effective source of water since the 1980s (Pandey et al., 2010). Kathmandu Upatyaka Khanepani Limited (KUKL), the major supplier of water throughout the valley, is fulfilling only 19% of the demand in the dry season and 31% in wet seasons (Thapa et al., 2018). This clearly depicts that subsurface water is being used by public to fulfil the deficit. Groundwater contribution to the river flow is estimated at the range of 30%–74% depending on the topography, soil type, climatic factor, among others (Pandey et al., 2010). However, rapid urbanization of the catchment is resulting increase in runoff and decrease in groundwater recharge (Thapa et al., 2018). The change in surface recharge will ultimately affects groundwater budget. Earlier studies have estimated groundwater recharge is the valley in the range of 4.6–14.6 Million-Cubic-Meters (MCM) per year (Pandey and Kazama, 2014) (Shrestha et al., 2017). However, groundwater level is depleting in the valley (Pandey et al., 2010) due to the combination of land use/cover (LULC) change and excessive pumping. To regulate either of the drivers, we need a scientific understanding of the extent of changes and associated impacts. Kaliraj et al. (2014) has reported that plain areas of the valley, especially, foothills, are the most suitable recharge areas, however, they are under rapid urbanization. Therefore, it is imperative to understand the extent of encroachment of recharge areas due to LULC change. It requires projection of future LULC change using a suitable model and translating the LULC change into potential groundwater recharge areas. Different types of models are available for LULC projection. They include, but not limited to, CLUE-S (Verburg et al., 2002a,b), Dinamica EGO (Soares-Filho et al., 2002), GEOMOD (Praskievicz and Chang, 2011), Landuse Sim (Pratomoatmojo, 2018) and Land Change Modeler (LCM) (Kim, 2010). The Conversion of Land Use and its Effects (CLUE) is a dynamic model which can simulate the LULC conversion with the reference of driving forces (Veldkamp and Fresco, 1996). In a small region or scale, CLUE output is not properly valid because in the local level there are many driving forces that are not built in CLUE, therefore, the Conversion of Land Use and its Effects at Small regional extent (CLUE-S) model was developed for the small extent analysis (Verburg et al., 2002a,b). CLUE-S model has been commonly used for the simulation of LULC in regional and small scale (Verburg et al., 2002a,b; Zheng et al., 2015; Zhou et al., 2016). To further translate LULC data for various scenarios into potential groundwater recharge areas, Geographic Information System (GIS) can be used as a tool, as it is capable to handle analysis of spatial data (Chowdhury et al.;, 2009; Kaliraj et al., 2014; Jhariya et al., 2016) Different thematic layers that influence groundwater recharge such as drainage density, LULC, geology, rainfall, and topography, among others, are generally used for delineating groundwater recharge areas (Chowdhury et al.;, 2009; Singh et al., 2018). The studies have integrated remote sensing and GIS techniques for the purpose. Depending upon relative importance of one layer against another, differential weights can be assigned to the thematic layers (Chaudhary et al., 2016). Analytical hierarchical process (AHP), a form of multi-criteria decision making (MCDM) techniques, can be used for deriving differential weights to the layers (Singh et al., 2018). Earlier studies have also attempted to delineate groundwater recharge areas in the valley. For example, JICA (1990) divided the KV into three zones (i.e., northern, central, and southern groundwater districts) for the purpose of groundwater development and management, and also delineated some pockets as potential recharge areas. Later on Thapa et al. (2017) also classified the five different area from the output of model. However, none of them attempted to analyze future LULC scenario and associated impact on recharge areas. This study therefore aims to evaluate the impact of LULC change in potential groundwater recharge areas in the KV by answering following research questions – i) what are projected changes in future LULC under different scenarios? ii) What are potential areas for groundwater recharge under current condition? iii) how and where the potential recharge areas are projected to encroach due to LULC change?
2. Study area Bagmati river watershed area is spread in 613 km2, and it is located between 27032'13″ and 27049'10″N and 85011'31″ and 85 31'38″E as shown in Fig. 1. It is bowl shaped valley with the elevation ranging from 1212 to 2762 m above the mean sea level. Bagmati is the major river that flows north to south within the basin originated from Shivapuri hill. Bisnumati, Hanumante, Dhobi Khola, Manohara, Balkhu, and Nakhhu are the major tributaries of the Bagmati river. Kathmandu Valley is enclosed by the Mahabharat mountain range surrounded by four hills namely Phulchowki (2762 m) in the southeast, Chandragiri/Champadevi in the southwest, Shivapuri (2722 m) in the northwest, and Nagarkot in the northeast, formerly known as the forts of the valley. Average mean monthly precipitation ranges from 4.2 mm in December to 402.1 mm in July. The total average annual rainfall of the basin is 1533 mm/yr and more than 80% of total rainfall is received during the monsoon period (June -September). The valley climate is semi-tropical and the monthly average maximum and minimum temperatures vary from 29.8 0C to 3.4 0C and average humidity of the basin is 75% (DHM, 2015). In dry season, there is minimum flow and high flow during monsoon. Within the basin, the average annual runoff of the Bagmati river at the Khokana station (Valley outlet point St. No 550.05) located at the southern tip of the valley is 456.01 million m3/year. Administratively, Kathmandu valley represents the three district Kathmandu, Lalitpur and Bhaktapur where five major cities Kathmandu, Lalitpur, Bhaktapur, Kritipur, and Thimi are situated. The outer periphery (hilly area) of KV is covered by mixed forest and the peri-urban areas are composed of a mixture of agricultural and built-up lands, and the central core of the valley is covered by built-up area (Thapa and Murayama, 2009). 0
2
Journal of Hydrology: Regional Studies 26 (2019) 100635
S. Lamichhane and N.M. Shakya
Fig. 1. Location of the study area.
3. Methodology Fig. 2 depicts the methodological framework adopted in this study. LULC map was generated for baseline case using remote sensing products. Future LULC was projected based on CLUE-S model under different projection scenarios. A suitable set of indicators/layers for delineating groundwater recharge areas were identified based on literature review, expert consultation, and data availability. Differential weights of the layers were estimated based on AHP method. Recharge areas delineated from the analysis of layers were evaluated based on infiltration rate map prepared as a part of this study. Based on these analyses, projected change in future LULC, current recharge areas, and encroachment in recharge areas were estimated and mapped. Detailed methods are elaborated in following sub-sections. 3.1. Watershed delineation The outlet of the Bagmati river basin was considered at Katuwal Daha (850 16, 41.96′' East and 270 34′ 54.6′' North) to incorporate the all urban and peri-urban areas of the KV. Watershed boundary was defined by considering the demographic area of the KV using the 30 m ASTER-GDEM (Global Digital Elevation Map). The watershed area was calculated as 613 km2.The entire basin land is divided into five categories: agricultural, built up, forest, water, and restricted area (like airport, park, governmental office area, cultural area etc.). 3.2. Data collection and processing Land sat image from 2010 to 2018 (USGS land sat 5 to 8) of the study area was downloaded for the analysis. By using GIS and land sat image, the LULC map from 2010 to 2018 were prepared. The 2010 LULC data awas downloaded from International Center for Integrated Mountain Development (ICIMOD) and the LULC was validated from ICIMOD (2010) data. Biophysical and socio-economic driving forces that influence the LULC change were selected for analysis. Topographical factors such as elevation, slope, aspect, etc. were taken from the 30 m ASTER global digital elevation map (GDEM) data. Human-geographic driving forces like river network, road network, location of settlement, etc. were taken from the GDEM and GIS database. Similarly, socio-economic data involving population was taken from CBS (2012). The geological factors which influence the LULC change like a geological formation, soil type, soil texture was taken from the geological map of Nepal. Hydro-meteorological data of the KV was taken from Department of Hydrology and Meteorology (DHM, 2015) to input in the analysis. 3
Journal of Hydrology: Regional Studies 26 (2019) 100635
S. Lamichhane and N.M. Shakya
Fig. 2. Overall analysis framework. Notes: GW:- Groundwater; AHP:- Analytical hierarchical process; MCDM:- multi-criteria decision making ; LULC:- Land use/cover; GIS:- Geographic Information System; O1:- Objective One
3.3. Past LULC change and validation of LULC map For the retrieved data as above, to classify the land category and prepare the LULC map, GIS was used as an analysis tool. Existing LULC map of the KV is shown in Fig. 3. After plotting the dataset of 2010–2018, the LULC trend of each land use type was outlined. ICIMOD also prepared the LULC map of Nepal in 2010. In this map, there are nine LULC type were identified, viz. Agricultural, Barren, Grass, River, Built Up, Needle-leaved close forest, Needle-leaved open forest, Broad-leaved close forest, Broad-leaved open forest. For the validation purpose, all the similar type of LULC was categorized in the similar group of ICIMOD data and the total LULC type was changed into five categories. The generated data of the study area was validated from the ICIMOD (2010) data by using Kappa statistical analysis (K). The Kappa values vary from 0 to 1. The maximum value (1) represents the most perfect agreement between two maps and zero means no level of agreement (Cohen, 1960). The K value can be represented by Eqn. (1) as follows:
Pr (a) − Pr (e ) ⎤ K =⎡ ⎢ 1 − pr (e ) ⎦ ⎥ ⎣
(1)
Where, Pr(a) is the observed relative agreement among all rasters and Pr(e) is the hypothetical probability of a chance of agreement. For the purpose of validation, five LULC categories were defined. Most of the forest area in the valley is situated in the peripheral 4
Journal of Hydrology: Regional Studies 26 (2019) 100635
S. Lamichhane and N.M. Shakya
Fig. 3. (a) 2010 LULC map of the study area; (b) ICIMOD, 2010 LULC map of the study area.
part and managed by community forest user groups (CFUGs). Therefore, from the past LULC data, no significant change in the forest coverage area was detected. This study therefore considered low conversion area. Areas such as Airport, Parks, Religious areas, etc were considered as a restricted area due to the governmental restriction not allowing their conversion to other LULC type. Due to the encroachment of the liner waterway of the river and water body of the KV cannot pass the peak runoff and suffer in numerous flood in the urban area (Shrestha, 2015). Currently, Bagmati Integrated River Basin project and River Road Corridor project is under implementation for the management of the waterway and river corridor. In this study, it is assumed that the LULC of the water body is not changing. LULC demand projection scenario was basically estimated from the increase in built up in map and rate of increase in population in the KV. For the sensitivity analysis of LULC demand and LULC change of the study area, five types of LULC demand scenario were analyzed based on built area and population increase rate: a) Normal land use change scenario b) Double land use change scenario c) Half of land use change scenario d) Population growth rate scenario e) Half of population growth rate scenario. 3.4. Future LULC projection The CLUE-S model was selected for future LULC projection. It simulates LULC based on an empirical analysis of location suitability combined with the dynamic simulation of competition and interactions between the spatial and temporal dynamics of LULC systems (Verburg et al., 2002a,b). The model is separated into two modules, namely, non-spatial demand module and spatially explicit allocation module. The demand module examines the demand for various LULC through past LULC trends or scenario-based LULC and then translates the demand for LULC for application by the spatial allocation module (Verburg, 2010). Ten driving forces (geology, soil type, elevation, slope, population density, distance from market, etc)were considered for the future LULC projection and the dependency of the each driving forces was analysis based on ROC value (Zhou et al., 2016). The CLUE-S uses the logistic regression model (Eq. 2) for the identification of driving factors and their probability of the LULC change. Ordinary logistic regression model, as used in other studies (e.g., Verburg et al., 2002a,b; Chen et al., 2010; Zhou et al., 2016; Shrestha et al., 2018), was used for the purpose. First, the study area was divided into grid units and the response variable (spatial distribution of land use type) was in a binary number (1 indicated the transition occur and 0 indicated that this grid did not exist and explained) (Verburg et al., 2002a,b). Then, the logistic regression model of the probability of the suitability of the LULC was analyzed based on the SPSS tools and goodness of fit of the analysis is checked by the receiver operating curve (ROC) value (Pontius and Schneider, 2001).
pi ⎞ log ⎛⎜ ⎟ = β + β x1, i -------+ β x n, i 0 1 n ⎝ 1 − pi ⎠
(2)
Where pi is the probability of specified LULC type i transition in the grid, and, βi is a coefficient to be estimated for each explanatory variable x n, i . The probability of each grid is calculated to generate the suitability map of different LULC. The model then combines the initial LULC map, suitability map and LULC transition rules to calculate total probability (Eq. 3). TPROPi,u = Pi,u + ELASu + ITERu
(3) 5
Journal of Hydrology: Regional Studies 26 (2019) 100635
S. Lamichhane and N.M. Shakya
Table 1 Indicators/Layers selected for delineating potential groundwater recharge areas. S.N.
Layer
Relation
Reference
1
Slope
2 3
River Distance Geology
Rukundo and Doğan (2019); Ansari et al. (2016) Rukundo and Doğan (2019) Rukundo and Doğan (2019)
4
Land Use/Cover (LULC)
Steeper slopes results in low contact period with the land surface and reduces infiltration and groundwater recharge Farther from the river corridor, lesser will be the infiltration and recharge Porous media (e.g., sand, gravel, etc.) has high recharge potential compared to non or less porous. Agriculture and forest areas have high recharge capacity compare the built-up areas
5
Precipitation
6
Aspect
7
Elevation
8
Population Density
9
Market Distance
10
Road Distance
Longer duration, more frequent, and less intense precipitation has higher recharge potential. East and south facing areas are drier and have less recharge potential compared to other aspects Flat/plan areas has may provide opportunity for pondage and therefore higher potential for recharge Densely populated areas are more impervious and therefore may have lesser recharge potentials Farther the distance from market centre, higher will be potential for recharge, as they are relatively less developed Farther the distance from road centers, higher will be the possibility for recharge
Rukundo and Doğan (2019); Ansari et al. (2016) Rukundo and Doğan (2019); Bashir et al. (2008) Bashir et al. (2008) Bashir et al. (2008) Bashir et al. (2008) Bashir et al. (2008) Bashir et al. (2008)
Where, TPROPi,u is the total probability of cell i is suitable for land use u. Pi,u is the spatial distribution probability obtained by logistic regression; ELASu is the transition elasticity of land use u and it varies from 0 (lesser stability of the corresponding land use) to 1 (higher stability of the corresponding land use), and ITERu is the iteration variable of land use u. From the CLUE-S model was generate the probability map and land use map. The result from the CLUE-S was found in the ASCII file and easily analysis from GIS tools. 3.5. Delineation of groundwater recharge area A set of 10 indicators (Table 1), covering three broad categories (Fig. 4), were selected based on literature review, their logical ink to recharge potential, expert consultation, data availability, and statistical test (e.g., ROC test). AHP method, as elaborated in other studies (e.g.,Chaudhary et al., 2016; Kaliraj et al., 2014; Singh and Nachtnebel, 2016), was adopted to estimate weights to the layers. The weights factor was obtained by pair-wise comparison of each layer against rest of the layers (Saaty, 2004). 3.5.1. Data standardization For using of each thematic layer of the study, each thematic layer grids have a different unique value and the standardization process converts them into the dimensionless same value of the layers. The standardization Eq. (4) and Eq. (5) are given below and obtained value of the study area was "larger the better" and "smaller the worst" of the alternative value of the grids in each thematic layer (Pei-Yue et al., 2010). Standardization used for larger the better:
yi =
xi − xi ( min) xi (max ) − xi (min)
(4)
Standardization used for larger the better:
yi =
x (max) − xi xi (max ) − xi (min)
(5)
Where, yi is the standardized value of the thematic grid and i is the index of thematic grid and xi , xi (max ) , xi (min) are the original, maximum and minimum values of the thematic grids respectively. 3.5.2. Calculation of the groundwater potential recharge area (GWPRA) index All the thematic layers were then aggregated in GIS to get a single score of potential recharge area (Eq. 6) (Malczewski, 1999). Sum of the product of each weights and its corresponding grid value gave the high potential recharge zone (Jhariya et al., 2016). All the assigned value in each layer grid and the weights factor of each layer computed the theoretical potential recharge areas. The higher value of the grid has higher potential and vice-versa. m
GWPRA =
n
∑ ∑ (wj × xi) (6)
j=1 i=1 th
Where, GWPRA is the groundwater potential recharge area of the study area; x i is the normalized weight of the i class of the thematic layer; wj is the weight derived from AHP of the jth thematic layer; and m is the total number of thematic layer and n is the total number of class in thematic layer. 6
Journal of Hydrology: Regional Studies 26 (2019) 100635
S. Lamichhane and N.M. Shakya
Fig. 4. Methodology for delineating groundwater potential.
The obtained recharge value from the analysis was categorized in the six different zone. The dimensionless output value was analysis by using the GIS application. GWPRA was classified as per Table 2 (Singh et al., 2018). 3.6. Evaluation of groundwater recharge map Theoretical recharge potential estimated as described in section 3.5 were further evaluated based on infiltration data collected from the field. In situ field test with double ring infiltration method was used to estimate information at different locations within the KV because double ring infiltrometer is likely to give better performance (Verbist et al., 2010). Infiltration rate at 86 locations as shown in Fig. 7 (c), with wide spatial coverage, was estimated from field test. By spatial interpolation, a raster map of infiltration rate was developed. Those two maps were compared using Kappa statistical analysis (K) for the evaluation. The Kappa values vary from 0 to 1, with one indicating the most perfect agreement between two maps (Cohen, 1960). After proper validation, baseline recharge area map was finalized.a Table 2 Groundwater recharge potential value. S.N.
Very Low
Low
Moderate low
Moderate high
High
Very High
Recharge potential value
2.1 - 3.0
3.1 - 4.0
4.1- 5.0
5.1 - 6.0
6.1 -7.0
7.1 - 8.0
7
Journal of Hydrology: Regional Studies 26 (2019) 100635
S. Lamichhane and N.M. Shakya
Fig. 5. (a) 2018 LULC map of the study area; (b) Comparative diagram of LULC.
3.7. Projection and analysis of encroachment of recharge area The validated pair-wise comparison of weights factor of each layer is used to generate the potential future groundwater recharge area by using the future LULC map generated by CLUE-S model. Then the future potential recharge areas were compared with baseline map to estimate encroachment in recharge areas by using GIS tools.
4. Results and discussion 4.1. Validation of satellite-driven LULC data (2010 data) LULC data was downloaded from the USGS website in Land sat image and processed by using Arc-GIS tools. The images in Figs. 3(a) and 5 (a) show LULC maps from 2010 and 2018. For the purpose of statistical data validation, the observed map of 2010 was overlapped into the ICIMOD, 2010 as shown in Fig. 3(b). Overall, the LULC map generated in this study was in good agreement with the one from ICIMOD (2010). With this validation, the Kappa coefficient (K) value 0.68 was obtained from the analysis. The result showed that the generated LULC of 2010 and the corresponding process had-considerable accuracy. The results presented in Table 3 show that the built-up area increased by 4.96% whereas the agriculture area decreased by 6.51% of the total land between 2010–2018. It can be noticed that due to the 2015 earthquake, 2016 built up area trend partially decreased but again rapidly increased during years 2017 and 2018. Due to the formation of community forest, user’s group encroachment of forest land is drastically decreased. As a result of rapid urbanization post 2015 earthquake, forest areas has been decreased partially, as reflected in the form of decrease in forest cover from 2016. Furthermore, increase in built-up area is reducing the open agriculture and barren land, thus, adversely affecting groundwater. For the future analysis, water body and the restricted area land use was assumed that there is no change in the future.
Table 3 Comparison matrix of LULC between.2010–2018. LULC
Land Use/Cover 2010 (%)
Land Use/Cover 2018 (%)
Agriculture Land Built-Up Forest Restricted Area Water Body Grand Total
Agriculture Land
Built-Up
Forest
Restricted Area
Water-Body
Grand Total
32.06 4.00 6.31 0.02 0.02 42.40
8.07 12.50 1.26 0.09 0.01 21.93
8.73 0.34 24.87 – – 33.94
– 0.11 0.00 1.31 – 1.42
0.04 0.02 0.02 – 0.21 0.30
48.91 16.97 32.46 1.42 0.25 100.00
8
Journal of Hydrology: Regional Studies 26 (2019) 100635
S. Lamichhane and N.M. Shakya
Table 4 ROC value of each LULC. LULC Type
Forest
Water Body
Built up
Agriculture
Restricted Area
ROC
0.92
0.87
0.91
0.83
0.94
4.2. Evaluation of CLUE-S model For the LULC change projection, ten different driving forces were considered as an input values which is already discussed in previous section. The ROC value validates the logistic regression model. If the ROC value is less than 0.7, the accuracy of the model is low and if the value is less or greater than 0.9 represents the accuracy is credible and high precision respectively (Pontius and Schneider, 2001). The ROC value of each LULC was greater than 0.83 and it shows that there was a strong relationship between the driving forces and the LULC type. These data indicated that the regression results are excellent and the regression coefficient from the driving forces was validated parameter for each LULC change is shown in Table 4. For the sensitivity analysis of the LULC demand of the study area was analyzed by five demand scenarios and it already describe in previous section. Normal land use demand scenario was quantified based on the LULC change from 2010 to 2018 and the population growth analysis was based on census population growth rate from 1971 to 2011 of the KV. Population growth rate demand projection is overestimated and non-realistic, and half of the normal land use change growth was underestimated. However, normal land use change is more practical or realistic compared to the other as shown in Fig. 6. According to the five different scenarios described in the previous section, the CLUE-S model was used to project the LULC change in the KV. Among these scenarios, normal land use change scenario was taken for further analysis. As water body and restricted areas are not expected to change, in all the scenarios, LULC change was considered in agriculture and built-up areas only. Simulated data from years 2010 to 2018 was used for the validation of the map and the forecasted data for year 2020 to 2050 is used for the analysis. The generated map from the CLUE-S model was validated through the Kappa (K) coefficient. The Kappa coefficient of years 2014, 2016 and 2018 was 0.75, 0.64 and 0.62 respectively. It showed that the analysis is of considerable accuracy. The model assumed that the urban development is concentrated into the existing market and dense village but the practically development is scattered due to high-cost of land and less open space in the existing market. From the simulation of CLUE-S, from 2020 to 2030, forest, water body and the restricted areas are not projected to change significantly. But the agricultural land is projected to decrease dramatically from 40.11% to 33.44% within that decade. Average annual encroachment of the fertile land is projected to be 0.67% (4.07 km2) by mostly built up area due to the construction of
Fig. 6. Future LULC change simulated by CLUE-S model (a) Normal LULC Growth rate scenario, (b) Population growth rate scenario. 9
Journal of Hydrology: Regional Studies 26 (2019) 100635
S. Lamichhane and N.M. Shakya
Table 5 Comparison matrix of LULC between 2020 and 2030, 2030 and 2040. LULC
Land Use/Cover 2030 (%) Agriculture Land
Land Use/Cover 2020 (%)
Agriculture Land Built-Up Forest Restricted Area Water Body Grand Total
33.4 – – – – 33.4 Land Use/Cover 2040 Agriculture Land 26.2 – – – – 26.2
LULC Land Use/Cover 2030 (%)
Built-Up
Agriculture Land Built-Up Forest Restricted Area Water Body Grand Total
6.2 23.4 0.7 – – 30.3 (%) Built-Up 7.2 30.3 0.5 – – 38.0
Forest
Restricted Area
Water-Body
Grand Total
0.5 – 34. – – 34.5
– – – 1.5 – 1.5
– – – – 0.3 0.3
40.1 23.4 34.7 1.5 0.3 100.0
Forest – – 34.0 – – 34.0
Restricted Area – – – 1.5 – 1.5
Water-Body – – – – 0.3 0.3
Grand Total 33.4 30.3 34.5 1.5 0.3 100.0
physical infrastructure. Meanwhile, the built area of the KV is projected to increase by 6.85% within a decade which can be clearly seen in the map and appended Table 5 and 6 as shown in Fig. 7. From the overall simulation model (2020 to 2050), agriculture and forest lands is projected to reduce by 20.4% (125.05 km2) and 0.9% (5.51 km2) of the total area respectively. The total built area is projected to increase by 21.3% (130.57 km2) of the total area. Within three decades, half of the agriculture land of the KV is projected to shift into the urban area. If the current level of urbanization continues as usual, more fertile land might convert into a constructed area. Also, the open space used as recreational area of the KV might decline. In upcoming days, the agricultural production will be decreased due to the loss of fertile land and increase the food deficiency. Meanwhile, the groundwater recharge by the encroachment of potential groundwater recharge area would lower the groundwater table. The encroachment of the open areas by the urbanization between 2010–2018 was 10.4% (63.54 km2) of the total basin area. The LULC map of the CLUE- S model (2020 to 2050) clearly shows that 21.4% (130.68 km2) of the open area is projected to change from fertile and permeable agriculture land to built-up area.
4.3. Delineation of potential recharge areas (Under current and future scenario) Ten thematic layers were used for the AHP to find the paired comparison weights factor. The grid value of each thematic layer was given by the pair-wise comparison from 1 to 9. The consistency ratio (CR) of each thematic layer was less than 8% which indicates the analysis is consistent and the generated data from each layer was suitable for the future AHP analysis. Overall ten thematic layer analysis, CR value 9.7% indicate that the pair-wise comparison was more consistent (Saaty, 2004) and the output weight age of each layer was more realistic. Comparison matrix of AHP and weights factor of each layer is as shown in Table 7. Matrix was generated from the experts’ view with the assumption that each layer influenced the rate of groundwater Table 6 Comparison matrix of LULC between 2040 and 2050, 2020 and 2050. LULC
Land Use/Cover 2040 (%)
Land Use/Cover 2050 (%)
Agriculture Land Built-Up Forest Restricted Area Water Body Grand Total
LULC
Land Use/Cover 2020 (%)
Agriculture Land
Built-Up
Forest
Restricted Area
Water-Body
Grand Total
19.7 – – – – 19.7
6.2 38.0 0.5 – – 44.7
0.3 – 33.5 – – 33.8
– – – 1.5 – 1.5
– – – – 0.3 0.3
26.2 38.0 34.0 1.5 0.3 100.0
Land Use/Cover 2050 (%)
Agriculture Land Built-Up Forest Restricted Area Water Body Grand Total
Agriculture Land
Built-Up
Forest
Restricted Area
Water-Body
Grand Total
19.7 – – – – 19.7
19.6 23.4 1.7 – – 44.7
0.8 – 33.00 – – 33.8
– – – 1.5 – 1.5
– – – – 0.3 0.3
40.1 23.4 34.7 1.5 0.3 100.0
10
Journal of Hydrology: Regional Studies 26 (2019) 100635
S. Lamichhane and N.M. Shakya
Fig. 7. (a) Land conversion map of Kathmandu Valley between 2010–2018; (b) Land conversion map of Kathmandu Valley between 2020–2050; (c) Infiltration test locations in Kathmandu Valley.
Table 7 Comparison matrix of AHP and layers weight factor. S.N.
Type
1
2
3
4
5
6
7
8
9
10
%
1 2 3 4 5 6 7 8 9 10
Slope River Distance Geology Land Use Precipitation Aspect Elevation Population Density Market Distance Road Distance
1.00 2.00 3.00 4.00 5.00 0.50 0.33 0.20 0.14 0.33
0.50 1.00 2.00 2.00 3.00 0.25 0.17 0.11 0.11 0.33
0.33 0.50 1.00 2.00 2.00 0.17 0.14 0.11 0.11 0.11
0.25 0.50 0.50 1.00 3.00 0.25 0.33 0.50 0.11 0.11
0.20 0.33 0.50 0.33 1.00 0.11 0.14 0.11 0.11 0.11
2.00 4.00 6.00 4.00 9.00 1.00 0.50 0.50 0.25 0.50
3.00 6.00 7.00 3.00 7.00 2.00 1.00 0.50 0.33 1.00
5.00 9.00 9.00 2.00 9.00 2.00 2.00 1.00 0.50 2.00
7.00 9.00 9.00 9.00 9.00 4.00 3.00 2.00 1.00 3.00
3.00 3.00 9.00 9.00 9.00 2.00 1.00 0.50 0.33 1.00
7.65 12.89 18.52 17.16 28.91 4.32 3.38 2.62 1.51 3.06
recharge. Form the analysis, precipitation (28.91%) is the main key factor for the groundwater recharge and the effect of geology (18.52%) and land use (17.16%) is relatively high. Bank of the river is mostly formed by the river deposition which has relatively high infiltration capacity and its weights factor 12.89% was also a major influencing parameter. Slope influences the contact time of runoff water and soil layer; higher slope has high velocity and less recharge and vice versa. By using GIS and AHP weights factor the potential groundwater recharge area was identified as shown in Fig. 8b. For validation of the AHP generated potential groundwater recharge area by the actual field condition, number of spouts (83 no) were identified as per the geological, hydro-meteorological, socio-economical changes location. Every variation of the geological formation, land use change in the location and change in the soil texture area is considered for the test point. The location point must be free from an anthropogenic disturbance and it must be an open space for the field test. Double ring infiltrometer (see for details: Lucke et al., 2014) test was conducted in 83 different locations of the study area. The overall location of the field test point is shown in the map (Fig. 7c). Actual field raster map was prepared by the field test infiltration rate (mm/hr) point and that was compared with the AHP generated map for the evaluation of AHP outputs. Due to lack of the information about litho-logical, moisture content, water table, soil density and so on, actual recharge rate was not accurately quantified. Due to the variable grid value unit (one infiltration value and other is unit less recharge factor), statistically, these maps could not be validated. However, from the inception of map and as shown in pervious literature (Pandey et al., 2013; Thapa et al., 2018; Adhikari et al., 2019) foothill of the mountain and the alluvial part of the river corridor were most potential area. Groundwater potential area was computed using Eq. (6) and the output of the study area was categorized into six parts as per quintile mapping in ArcGIS. Area coverage of the very low, low, moderate low, moderate high, high, and very high was 4.29%, 9.14% and 14.72%, 32.59%, 32.89%, and 7.37% respectively of the total area coverage as shown in Fig. 8. Currently, high and very high areas are still in open area in the foot hill of northern and eastern part of mountain of the valley. Central part of the valley also has the potential recharge area but now it is mostly in concrete jungle. Northern and eastern part of the valley were more groundwater recharge potential zone and the peripheral boundary of the valley are less recharge potential due to high hill and geological formation. From the comparison of the observation of field test map, AHP method map and previous literature and maps (Pandey et al., 2013; Thapa et al., 2018; Adhikari et al., 2019) potential recharge location may be a similar area. The southern and northern foothill 11
Journal of Hydrology: Regional Studies 26 (2019) 100635
S. Lamichhane and N.M. Shakya
Fig. 8. (a) Potential recharge area mapping by field test; (b) Potential recharge area mapping by AHP method.
of mountain areas are the main possible recharge location of the KV. If that potential area could be preserved and managed, and artificial recharge system will be adapted which can overcome the water security situation of the KV in the future. Form the analysis, the theoretically generated map from the different layers using MCDM techniques is realistic and useable for the location which has less groundwater recharge information like Nepal. 4.4. Change in recharge areas under projected land use/cover Rapid urbanization mostly influences the groundwater recharge. Future encroachment of recharge area was identified by using AHP weights factor and the LULC map generated from the CLUE-S model. The projected future potential recharge zones are shown in Fig. 9. From the maps shown in Fig. 9, it can be observed that every decade the open area (agriculture and forest land) of the valley is projected to encroach by the urbanization. From the data, every decade on an average of 6% of open land is projected to decrease and each year 3.66 km2 of open land will change into the urban area. The likely reason for not consistent trend across future time periods in Table 8 is due to consideration of only one factor (i.e., LULC) whereas recharge areas are affected by various other factors as well. Such decrease will lead in the reduction of groundwater recharge. Reduction in groundwater leads to the lowering of groundwater table and base flow contributions in river basin. In addition, the rate of groundwater recharge will reduce over the time primarily due to the increment of the urban population and increases the volume of the water demand. Such reduction also effects in quantity and quality of water and other urban environments. 5. Conclusion A methodological framework has been developed and implemented in this study that helps to classify the LULC type, identification of the potential recharge area, possible protection of the agriculture land for the groundwater recharge, and future possible encroachment in the groundwater recharge areas. Identification of the LULC, the future projection of LULC and encroachment of recharge area generated in this study will be useful in effective urban planning, optimal utilization of water supply demand, and water resources planning and management. The sum of the present study could be summarized in the following conclusions: a) The trend of past and present LULC analysis highlights that the agricultural land is projected to decrease by 6.51% leading to an increment of 4.96% in built-up areas. It is interesting to note that between 2010–2018, the forest area is found to be increased by 1.48% due to community forestry programs. Conversion of pervious recharge areas into impervious constructed areas would lead into the reduction of groundwater recharge as the increment of forest area does not necessarily contributes significantly into recharge. b) The predicted scenario between 2020–2050 using CLUE-S model highlights that the built-up area is projected to increase by 21.4% in the next 30 years and the agricultural and forest areas is projected to decrease by 20.5% and 0.9% respectively. The 12
Journal of Hydrology: Regional Studies 26 (2019) 100635
S. Lamichhane and N.M. Shakya
Fig. 9. (a) Potential recharge map for 2020 (b) Potential recharge map for 2050 (c) Recharge encroachment area by potential urbanization from 2020 to 2050 (d) Recharge encroachment scenario by potential urbanization from 2020 to 2050.
results outline that the fertile lands would be converted into urban areas leading to the reduction in agricultural products and productivity. Lack of irrigation water may also be a crucial case which imparts additional pressure over the groundwater. c) Precipitation, geology, land use/cover, distance from the river, slope, among others were identified to be the significant factors to contribute in delineation of potential groundwater recharge area for the AHP/MCDM technique. Notably, 7.37% of the total area was obtained to be the high recharge area, 32.89% under high, 32.59% moderately high, 14.72% moderately low, 9.14% low, and 4.29% very low recharge areas. The high and moderately high recharge areas are situated in the urban areas with greater tendency of expansion. Thus, there might be the case that due to the LULC change; more high to moderate recharge areas would be converted to urban areas leading to additional pressure on groundwater. 13
Year
2020 2030 2040 2020
S.N.
1. 2. 3. 4.
to to to to
2030 2040 2050 2050
4.14 3.15 3.42 10.71
Increase in recharge area (km2)
Table 8 Decadal encroachment area and percentage of GWPRA.
0.68 0.52 0.56 1.75
Change in percentage 40.41 43.83 39.78 124.02
Decrease in recharge area (km2)
6.62 7.81 6.51 20.31
Change in percentage
5.94 6.66 5.95 18.55
Actual decrease percentage
S. Lamichhane and N.M. Shakya
Journal of Hydrology: Regional Studies 26 (2019) 100635
14
Journal of Hydrology: Regional Studies 26 (2019) 100635
S. Lamichhane and N.M. Shakya
d) Future LULC change analysis highlighted that 6% (3.99 km2) of rechargeable land area are projected to convert to impervious land in a decade. As Kathmandu Valley is already facing water scarcity, the urban expansion would put additional pressure leading to the graver water scarcity. The application of CLUE-S, GIS and AHP/MCDM techniques are more realistic, cost-effective, simple, and reliable approaches, and scientifically sound methods for the analysis. The methodology used in the study area can also be used into other urban areas in the country as well as worldwide. The findings from the study provide valuable information for decision-makers, planners, and managers for the sustainable management of the land and water resources. The impact of LULC on groundwater and surface water would ultimately affect the river runoff pattern. To address this issue, further analyses would be carried out in the future. Declaration of Competing Interest The authors declare no conflict of interest of any kind apply to the reported work. Acknowledgments The first author acknowledges the funding provided by Institute of Engineering, Tribhuvan University for the doctoral research grant. Authors acknowledge the support provided by University of Yamanashi, Japan, NTNU, Norway, Hiroshi Ishidaira, Bhesh Raj Thapa, Dipendra Gautam, Sarita Dawadi, Saraswati Thapa, Hari Poudel, and other people who supported in various ways. References Adhikari, B.R., Shrestha, S.Das, Shakya, N.M., 2019. Future Urban Water Crisis in Mountain Regions: Example of Kathmandu Valley, Nepal. Springer, Urban Drought, pp. 169–182. Ansari, T.A., Katpatal, Y.B., Vasudeo, A.D., 2016. Spatial evaluation of impacts of increase in impervious surface area on SCS-CN and runoff in Nagpur urban watersheds. India. Arab. J. Geosci 9. https://doi.org/10.1007/s12517-016-2702-5. Bashir, A.M., Fakhry, A.A., Philip, E.L., James, W.L., Mostafa, M.S., 2008. Environmental Hydrogeology. Taylor & Francis Group, USA. Belhassan, K., 2011. Relationship between River Flow, rainfall and groundwater pumpage in Mikkes Basin (Morocco). Iran. J. Earth Sci. 3, 98–107. CBS, N., 2012. National Population and Housing Census 2011. National Report. . Chaudhary, P., Chhetri, S.K., Joshi, K.M., Shrestha, B.M., Kayastha, P., 2016. Application of an Analytic Hierarchy Process (AHP) in the GIS interface for suitable fire site selection: A case study from Kathmandu Metropolitan City. Nepal. J. Socio-Economic Plan. Sci. 53, 60–71. Chen, F., Hu, Y., Peng, X., Wang, L., 2010. Simulation of land use/cover change based on the CLUE-S model. 2010 18th International Conference on Geoinformatics. IEEE 1–5. Chowdhury, A., Jha, M.K., Chowdary, V.M., Mal, B.C., 2009. Integrated remote sensing and GIS based approach for assessing groundwater potential in West Medinipur district, West Bengal. India. Int. J. Remote Sens. 30, 231–250. Cohen, J., 1960. A coefficient of agreement for nominal scales. Educ. Psychol. Meas. 20, 37–46. DHM, 2015. Hydrological Records of Nepal, streamflow Summary, Updated Version. Department of Hydrology and Meteorology: Goverment of Nepal, Ministry of Water Resources, Kathmandu. Gautam, D., Prajapati, R.N., 2014. Drawdown and dynamics of groundwater table in Kathmandu Valley. Nepal. Open Hydrol. J. 8. Han, D., Currell, M.J., Cao, G., Hall, B., 2017. Alterations to groundwater recharge due to anthropogenic landscape change. J. Hydrol. (Amst) 554, 545–557. https:// doi.org/10.1016/j.jhydrol.2017.09.018. ICIMOD, 2010. Land Cover of Nepal 2010 [Dataset]. Available online at:. International Center for Integrated Mountain Development (ICIMOD), Kathmandu, Nepal. http://rds.icimod.org/Home/DataDetail?metadataId=9224. Jhariya, D.C., Kumar, T., Gobinath, M., Diwan, P., Kishore, N., 2016. Assessment of groundwater potential zone using remote sensing, GIS and multi criteria decision analysis techniques. J. Geol. Soc. India 88, 481–492. JICA, 1990. Groundwater Management Project in Kathmandu Valley. Nepal Water Supply Corporation, Kathmandu, Nepal. Kaliraj, S., Chandrasekar, N., Magesh, N.S., 2014. Identification of potential groundwater recharge zones in Vaigai upper basin, Tamil Nadu, using GIS-based analytical hierarchical process (AHP) technique. Arab. J. Geosci. 7, 1385–1401. Kim, O., 2010. Comparison of Two GIS based Land change modules for constructing REDD baselines in Bolivia. AAG Annu. Meet. Lucke, T., Boogaard, F., van de Ven, F., 2014. Evaluation of a new experimental test procedure to more accurately determine the surface infiltration rate of permeable pavement systems. J. Urban Plan. Transp. Res. 2, 22–35. Malczewski, J., 1999. GIS And Multicriteria Decision Analysis. John Wiley & Sons. Pandey, V.P., Chapagain, S.K., Kazama, F., 2010. Evaluation of groundwater environment of Kathmandu Valley. J. Geogr. Environ. Earth Sci. Int. 60, 1329–1342. Pandey, V.P., Kazama, F., 2014. From an open-access to a state-controlled resource: the case of groundwater in the Kathmandu Valley, Nepal. Water Int. 39, 97–112. https://doi.org/10.1080/02508060.2014.863687. Pandey, V.P., Shrestha, S., Kazama, F., 2013. A GIS-based methodology to delineate potential areas for groundwater development: a case study from Kathmandu Valley. Nepal. Appl. Water Sci. 3, 453–465. Pei-Yue, L., Hui, Q., Jian-Hua, W.U., 2010. Groundwater quality assessment based on improved water quality index in Pengyang County, Ningxia, Northwest China. J. Chem. 7, S209–S216. Pontius Jr, R.G., Schneider, L.C., 2001. Land-cover change model validation by an ROC method for the Ipswich watershed, Massachusetts, USA. Agric. Environ. Ecosyst. 85, 239–248. Praskievicz, S., Chang, H., 2011. Impacts of climate change and urban development on water resources in the Tualatin River Basin. Oregon. Ann. Assoc. Am. Geogr. 101, 249–271. https://doi.org/10.1080/00045608.2010.544934. Pratomoatmojo, N.A., 2018. LanduseSim algorithm: land use change modelling by means of cellular automata and geographic information system. IOP Conf. Ser. Earth Environ. Sci. 202. https://doi.org/10.1088/1755-1315/202/1/012020. Rukundo, E., Doğan, A., 2019. Dominant influencing factors of groundwater recharge spatial patterns in Ergene river catchment, Turkey. Water (Switzerland) 11. https://doi.org/10.3390/w11040653. Saaty, T.L., 2004. Fundamentals of the analytic network process—multiple networks with benefits, costs, opportunities and risks. J. Syst. Eng. Syst. 13, 348–379. Shrestha, S., 2015. Assessment of Bagmati River Encroachment Through Application of GIS and Remote Sensing. Pokhara University. Shrestha, S., Bhatta, B., Shrestha, M., Shrestha, P.K., 2018. Integrated assessment of the climate and landuse change impact on hydrology and water quality in the Songkhram River Basin. Thailand. Sci. Total Environ. 643, 1610–1622. https://doi.org/10.1016/j.scitotenv.2018.06.306. Shrestha, S., Kafle, R., Pandey, V.P., 2017. Evaluation of index-overlay methods for groundwater vulnerability and risk assessment in Kathmandu Valley. Nepal. Sci.
15
Journal of Hydrology: Regional Studies 26 (2019) 100635
S. Lamichhane and N.M. Shakya
Total Environ. 575, 779–790. https://doi.org/10.1016/j.scitotenv.2016.09.141. Shrestha, S., Semkuyu, D.J., Pandey, V.P., 2016. Assessment of groundwater vulnerability and risk to pollution in Kathmandu Valley. Nepal. Sci. Total Environ. 556, 23–35. Singh, L.K., Jha, M.K., Chowdary, V.M., 2018. Assessing the accuracy of GIS-based Multi-Criteria Decision Analysis approaches for mapping groundwater potential. Ecol. Indic. 91, 24–37. Singh, R.P., Nachtnebel, H.P., 2016. Analytical hierarchy process (AHP) application for reinforcement of hydropower strategy in Nepal. Renewable Sustainable Energy Rev. 55, 43–58. https://doi.org/10.1016/j.rser.2015.10.138. Soares-Filho, B.S., Cerqueira, G.C., Pennachin, C.L., 2002. DINAMICA—a stochastic cellular automata model designed to simulate the landscape dynamics in an Amazonian colonization frontier. Ecol. Modell. 154, 217–235. Thapa, B., Ishidaira, H., Pandey, V., Bhandari, T., Shakya, N., 2018. Evaluation of water security in Kathmandu valley before and after water transfer from another basin. Water 10, 224. Thapa, B.R., Ishidaira, H., Pandey, V.P., Shakya, N.M., 2017. A multi-model approach for analyzing water balance dynamics in Kathmandu Valley. Nepal. J. Hydrol. Reg. Stud. 9, 149–162. https://doi.org/10.1016/j.ejrh.2016.12.080. Thapa, R.B., Murayama, Y., 2009. Examining spatiotemporal urbanization patterns in Kathmandu Valley, Nepal: remote sensing and spatial metrics approaches. Remote Sens. (Basel) 1, 534–556. https://doi.org/10.3390/rs1030534. Veldkamp, A., Fresco, L.O., 1996. CLUE-CR: an integrated multi-scale model to simulate land use change scenarios in Costa rica. Ecol. Modell. 91, 231–248. https:// doi.org/10.1016/0304-3800(95)00158-1. Verbist, K., Torfs, S., Cornelis, W.M., Oyarzún, R., Soto, G., Gabriels, D., 2010. Comparison of single-and double-ring infiltrometer methods on stony soils. Vadose Zone J. 9, 462–475. Verburg, P., 2010. The CLUE model. Hands-on Exercises. Course Material, in: The CLUE Model. Hands-on Exercises. Course Material. Institute for Environmental Studies, University of Amsterdam. Verburg, P., Soepboer, W., Veldkamp, A., Limpiada, R., Espaldon, V., Mastura, S.S.A., 2002a. Modeling the spatial dynamics of regional land use: the CLUE-S model. Environ. Manage. 30, 391–405. Verburg, P.H., Soepboer, W., Veldkamp, A., Limpiada, R., Espaldon, V., Mastura, S.S.A., 2002b. Modeling the spatial dynamics of regional land use: the CLUE-S model. Environ. Manage. 30, 391–405. https://doi.org/10.1007/s00267-002-2630-x. Wakode, H.B., Baier, K., Jha, R., Azzam, R., 2018. Impact of urbanization on groundwater recharge and urban water balance for the city of Hyderabad, India. Int. Soil Res. Water Conserv. 6, 51–62. Zheng, H.W., Shen, G.Q., Wang, H., Hong, J., 2015. Simulating land use change in urban renewal areas: a case study in Hong Kong. Habitat Int. 46, 23–34. Zhou, R., Zhang, H., Ye, X.-Y., Wang, X.-J., Su, H.-L., 2016. The delimitation of urban growth boundaries using the CLUE-S land-use change model: Study on Xinzhuang Town, Changshu City, China. Sustainability 8, 1182.
16