Model-based estimation of land subsidence in Kathmandu Valley, Nepal

ABSTRACT This study is the first to assess land subsidence in the Kathmandu Valley, Nepal. Land subsidence simulations were based on a fully calibrated groundwater (GW) flow model developed using a coupled surface–subsurface modelling system. Subsidence is predicted to occur as a result of deep aquifer compaction due to excessive GW abstraction. The north and north-east areas at the periphery of the GW basin are hotspots for this subsidence. The estimated subsidence is most sensitive to changes in land cover within the recharge areas. The model shows the Melamchi water supply project assists in the control of subsidence to some extent. In the absence of land subsidence measurements, this paper highlights the location and the potential levels of the subsidence hazard which will be useful for hazard prevention management. Additionally, this work provides a basis to design field investigations, monitoring networks for land subsidence and upgrading the present GW monitoring network. Although the study has presented a preliminary analysis, a more comprehensive model inclusive of clay subsidence is required to address the subsidence vulnerability of the central densely populated core of the valley, which reflects the need for a comprehensive database of the hydrogeology in the valley.


Introduction
Land subsidence is the movement of the Earth's surface as it shifts/ sinks/ settles downward relative to a datum. Land subsidence can have severe impacts in terms of the water sector. This can include reduction in groundwater storage capacity of a geologic strata, deterioration in groundwater quality, restrictions on groundwater pumping, intensification of the movement of ground fissure, destruction of natural drainages, thus, exposing the land to flooding, and river channel distortion (Balogun et al. 2011). Additionally, it can cause structural damage to highways, buildings, etc. Of possible causes, fluid withdrawal has been found to be the most important one (Poland et al. 1972, Larson et al. 2001, Hu et al. 2004, Xue et al. 2005). There are 150 major cities in the world where land subsidence is substantial (Krijnen & de Heus 1995;Dottridge & Abu Jaber 1999, Heath & Spruill 2003. In the U.S.A., withdrawal of groundwater has permanently lowered the elevation of about 26,000 km 2 (Holzer & Galloway 2005). In China, it has affected an area of 90,000 km 2 (Xue et al. 2008) resulting in huge economic losses. Land subsidence studies in the San Joaquin Valley, California (Larson et al. 2001); Antelope Valley, California (Galloway et al. 1998Leighton & Phillips 2003); Yangtze River Deltas, China (Xue et al. 2008, Balogun et al. 2011Mexico City (Osmano glu et al. 2011);Kerman, Iran (VaeziNejad et al. 2011); and Pingtung area, Taiwan (Hu et al. 2006), have concluded that a major portion of the subsidence resulted from excessive groundwater abstraction. Therefore, understanding the role of the groundwater environment in controlling land subsidence is an important aspect of groundwater management.
Land subsidence studies include measurements, modelling, or both measurements and modelling. Land subsidence can be indirectly measured by monitoring land surface displacements such as historic geodetic surveys, GPS records, InSAR (interferometric synthetic aperture radar), and T-LiDAR (tripod light detection and ranging). The latter two are remote sensing technologies. InSAR is an attractive option due to its historic records (from 1992), high spatio-temporal resolution (30 m, every 35 days), and potentially high resolution of displacement ( § 5 mm) (Galloway & Sneed 2013). However, it requires radiometrically favourable landscapes and lacks higher precision. On the other hand, T-LiDAR is ideal for regions that are actively deforming, but too labour intensive for scales greater than a few square kilometres. For sub-millimetre precision, greater accuracy and direct measurements of subsidence, borehole extensometers are used (Galloway & Sneed 2013).
Other researchers have used observed data or remote sensing or both for addressing the spatiotemporal variation of historical land subsidence (Balogun et al. 2011, Hung et al. 2011, Osmano glu et al. 2011). Paired historical groundwater-level and subsidence time series have been used to estimate the groundwater-level threshold corresponding to the onset of inelastic compaction in the aquifer system (Holzer 1981, Galloway & Sneed 2000. Some have used the records to validate the consistency of their coupled model of groundwater flow and aquifer system compaction (Galloway et al. 1998, Larson et al. 2001. Such calibrated models have been used to predict the future land subsidence for various groundwater management scenarios (Larson et al. 2001, Leighton & Phillips 2003. Xue et al. (2008) applied time-varying visco-elasto-plastic characteristics to the aquifer material.  estimated interbed storage parameters from land subsidence observations and calibrated a regional groundwater flow and subsidence model. In this study where there are no observed records of land subsidence, a model-based estimation cannot be validated. However, a model validated using other hydrological (surface run-off) and hydrogeological (groundwater flow) processes can be used to calculate land subsidence. This may yield a very good first estimate, although unvalidated, for the subsidence analysis especially in data-poor and research-lacking regions like the Kathmandu Valley, Nepal. In addition, numerical modelling prior to comprehensive hydrogeological investigation helps in conceptualizing the groundwater flow controls and highlights the field data requirement for comprehensive groundwater modelling (Rushton 1986(Rushton , 2003. The Kathmandu Valley (602 km 2 ) is located at the uppermost part of the Bagmati River Basin (BRB) (Figure 1). BRB is one of the most vulnerable basins in the country in terms of freshwater availability under environmental change . Kathmandu Upatyaka Khanepani Limited (KUKL) is the sole water supply operator providing water supply services in the urban and many rural areas of the valley (Shrestha 2012a). In total, the KUKL supply accounts for only 22.5 % of demand in the driest month, up to a maximum of 37.8 % in the wettest month (KUKL 2011). In search of more reliable sources, most of the residents depend entirely on shallow groundwater sources for their daily water needs. They do so through an unspecified number of dug wells, hand pumps, borings, and stone spouts (Shrestha 2012b). In addition, there are hundreds of private deep wells for commercial and institutional purposes. In recent times, groundwater from KUKL wells within the deep aquifers has become a major contributor to the city's water supply (Shrestha 2012b). It constitutes nearly half of the total water supply during wet season and 60%-70 % during dry season (Rana et al. 2007). Groundwater, therefore, has been and will continue to be an important source of water supply in the area . Excessive groundwater abstraction in conjunction with surface sealing from urbanization has appeared to have an impact on the groundwater reserves within the valley since the 1980s. The history of the groundwater status in the Kathmandu Valley as described in  is shown in Table 1.
The water level is gradually decreasing over the time in all the well fields ) of the valley which is caused by the heavy extraction of groundwater particularly from the deep aquifer and its inadequate recharge. The discharge of groundwater pumping wells is also declining Figure 1. Study area -Kathmandu valley; Insets: (A) Bagmati river basin and Nepal, (B) Bagmati river basin and its upper catchment of Kathmandu valley. Source: Prepared by authors based on GIS layers obtained from various sources. Country boundary and river networks from Department of Survey, the Government of Nepal; river basin boundary delineated from 90m resolution digital elevation model downloaded from http://gdem.ersdac.jspacesystems.or.jp/; groundwater basin extent from Department of Mines and Geology, the Government of Nepal; recharge areas delineated based on GIS layers from ; urbanization extent extracted from Google-Earth (2012); groundwater monitoring locations from GWRDB (2011). Haphazard pumping, considerable GW declination compared to their discharge at the time of installation . Bacterial contamination of the deep aquifer due to the improper or non-sealing of abandoned wells also indicates the depletion of the groundwater levels (Metcalf & Eddy 2000). On top of that, the population of the valley has increased by over threefold in the last three decades. Unfortunately, there is a continuing lack of implementation of the policy on sustainable water supply and groundwater resources. The effect of this is that the clay aquitard and deep aquifer are expected to get consolidated by the resulting reduction in pore water pressure. This raises concerns about the risks of land subsidence in areas with highly compressible clay and silt layers (Pradhanang et al. 2012).  have rated land subsidence monitoring of the valley as very poor or nil and  have suggested to initiate monitoring and study of land subsidence. Subsidence has been reported in the Shantibasti area, Hyumat tole, and Imadol, but detailed studies have not been carried out (Rana et al. 2007). Fortunately, the valley so far has not faced large amounts of local or differential subsidence. However, this does not imply that there is no land subsidence. Especially in regards to the fact that the valley has geological conditions similar to the locations of the world where maximum subsidence has been exhibited.  have recommended to develop a groundwater model of the Kathmandu Valley aquifers and thus suggest management scenarios incorporating the effect of climate change. Lack of monitoring and related studies make it clear that research on land subsidence in the valley is essential. The objective of this paper is to analyse the possibility of land subsidence over the catchment based on the dynamics of groundwater table (GWT) for future scenarios via hazard maps. Due to lack of land subsidence data for direct validation, land subsidence is predicted by the combined use of a validated groundwater model, resulting fluctuations in effective stresses and the subsurface information of the soil.

General methodological framework
This study uses the SHETRAN hydrological modelling system (http://research.ncl.ac.uk/shetran/) to conceptualize the catchment and its hydrological processes. The model is calibrated and validated using available data resources which is then subjected to future scenarios to get the land subsidence estimates. The framework of study is given in Figure 2.

Setting up SHETRAN hydrological model
SHETRAN is a physically based, distributed, deterministic, integrated surface and subsurface modelling system, designed to simulate water flow, sediment transport and contaminant transport at the catchment scale. It is designed to have the capability to predict the consequences of given changes in climate and land use (Ewen et al. 2000). SHETRAN has a modular design, in which each module or component is used to represent the different physical processes of the hydrological cycle in each part of the catchment. The water flow component of SHETRAN is an updated version of the Syst eme Hydrologique Europ een (SHE) (Abbott et al. 1986). Table 2 summarizes the methods used in SHETRAN to simulate the processes of its water flow component. For technical details, please refer document manuals (SHETRAN 2013a(SHETRAN , 2013b(SHETRAN ,2013c. The catchment is conceptualized as an ensemble of columns (refer Figure 3) and networks of stretches of channels, called river links. Each column is a stack of computational cells containing information of land cover/vegetation type at the top and sequential depth of soil/s horizons  underneath it. Each river link runs along the edge of a column (refer Figure 4(a)). SHETRAN works out the water balance in each column. The net outcome of hydrological processes at all the columns is the hydrological behaviour of the catchment. The Kathmandu Valley is modelled using cells of 700 m Â 700 m resolution in horizontal and varying (from few centimetres at the top to a few metres at the bottom) resolutions in vertical. The columns representing the deepest parts of the valley has a maximum of 100 cells. A no-flow boundary condition, i.e. assumption of catchment subsurface isolation is incorporated into the model. It means the groundwater basin has been modelled as impervious bed with no inter-basin groundwater movement. Apart from the natural processes, groundwater abstraction is also incorporated in the model. The deep aquifer abstraction is modelled using 260 abstraction points (refer Figure 4(b)) and their corresponding pumping capacity and screen levels. Whereas the shallow aquifer abstraction is represented by pumping wells, one at each of the columns within the urban land cover. The model considers 70% of daily abstraction to occur at peak hours while 30% of abstractions is returned back to the rivers (IDC 2009). A 2.4-times increase in abstraction in dry months (KUKL 2011) is also included.

Model calibration and validation
The flow hydrology for the Kathmandu Valley catchment has been calibrated using the hourly outlet discharge time series for 2011 and validated with the discharge for 2012. Whereas, the monthly groundwater level time series for 2011/2012 was used for the validation of groundwater hydrology. Due to requirement of hourly resolution of input hydrometeorological data of SHETRAN for the best results and their availability only since the establishment of the Flood Forecasting Project, a single year was used each for the calibration and validation of the model. Model calibration parameters are listed in Table 3. The initial values of these parameters for calibration are taken from literature and SHETRAN manuals (SHETRAN 2013a). The urban columns were provided with a few centimetres thin and almost impermeable clay layer with a minimal AET/PET ratio at the top to simulate quick drainage and lack of atmospheric exposure. A high conductivity sand layer was also used as the second layer in the urban and the first layer in non-urban areas to enhance the discharge calibration. During calibration, the simulated discharge hydrograph at the outlet of the catchment is assessed by the Nash Sutcliffe Efficiency (NSE) and discharge volume deviation (V) test defined as follows: where x i is the observed value, b x i is the value estimated by the model, x is the mean of the observed values, n is the number of data values, V s is the simulated discharge volume, and V o is the observed discharge volume. Nash-Sutcliffe criterion describes the ability of the model to explain the variability in the data while the volume deviation checks the mass balance of the model. A good model is one that has a good visual match to the observed data, with a value of NSE close to one and a low DV. The calibration and validation of model was efficiently achieved with the above-mentioned statistics employing SWAT Calibration Helper v1.0 (Shrestha 2016) which is a customized Visual BASIC for Applications (VBA)-coded Microsoft Excel file. During validation, the measure of success for the simulated GW elevation is measured with the correlation (r) between the measured and the simulated phreatic surface levels (PSL) supported by computation of probable error (PE) and standard error (SE): PE represents the reliability of the obtained value of r. A low probable error suggests that the soobtained value of r is reliable.

Population/land cover
The valley's population has increased by over threefold in the last three decades ( Figure 5) (CBS 2012). As a result, land cover has changed drastically; spreading the concrete jungle radially outward covering the valley floor and the foothills of the valley rim. It is assumed that the population of the valley will continue to increase with this same trend. The population forecast and the corresponding urban coverage in the model are as shown in Table 4. Here it is assumed that 40% of the population increase is accommodated by the existing urban cover while the remaining 60% of the growth is represented by increase in the urban cover. Since population growth is a continuous process, it is assumed that its effect can be averaged over time as shown in Table 5. Previous geological investigations show that the middle clay aquitard is at its thinnest at the northern extent of the valley floor and an interconnection between the shallow and the deep aquifer exists here. Hence this is the primary recharge area for the valley groundwater basin. If the urban cover seals this recharge area in future, further reduction in the recharge can be anticipated. Bearing this possibility in mind, three scenarios are created for the propagation of urban cover in the catchment: (1) Encroachment of northern recharge areas (2) Conservation of northern recharge areas (3) Uniform outward radial expansion The propagation of the urban land cover for the three scenarios is depicted in Figure 6.

GW abstraction
In a bid to meet escalating water demand as well as to reduce stress on groundwater resources, the Melamchi Water Supply Project (MWSP) is underway (Shrestha 2012a). MWSP is a comprehensive multi-donor water supply mega project that aims to alleviate the chronic shortage of potable water in the valley. The arrival of MWSP is sure to have positive impact on the groundwater abstraction trend in the valley. However, having an uncertain implementation history, people still doubt that the completion of MWSP will be achieved on time. Bearing this possibility in mind, two scenarios were produced regarding the GW abstraction pattern in the future: (1) Completion of MWSP on time or Melamchi (2) Failure to complete MWSP up to 2028 or No Melamchi   MWSP is expected to be completed in three phases: (1) 170 million litres per day (MLD) from Melamchi River; (2) 170 MLD from Yangri River; and (3) 170 MLD from Larke River (MWSDB 1998). Based on the current progress of MWSP, it will complete its first phase by 2016 and its second and third phases by 2019 and 2025, respectively (Shrestha 2012a). The average groundwater contribution in the water supplied by KUKL in the years 2011KUKL in the years , 2016KUKL in the years , 2019KUKL in the years , and 2025, and the corresponding total supply was taken from Shrestha (Table 6) (2012a) . For the deep aquifer water withdrawal, the same 260 abstraction points (Figure 4(b)) are used but with a variation in pumping rate. For No-Melamchi scenario, the pumping is determined by the population and considered as a function of population increase. Whereas for Melamchi scenario, the water supply eases the stress on groundwater, thus pumping is considered as a function of the decrease in GW % contribution in the KUKL water supply (Table 6). On the other hand, the water withdrawal from shallow aquifer is basically direct household groundwater consumption. Assumptions such as water requirement of 135 litres per capita per day (LPCD) and fulfilling of deficit from shallow (household) groundwater extraction are made. The shallow GW abstraction will increase with the difference between supply and demand, and with the proliferation of urban cover. The change in shallow GW abstraction rates from each urban cell for both the scenarios is shown in Figure 7. The figure also depicts the total GW abstractions made from the GW basin model throughout the simulation. The annual GW abstraction rises from 128 mcm/yr from the beginning to 245 mcm/yr at the end of the simulation for No-Melamchi scenario, while it dips down to 17 mcm/yr for Melamchi scenario. It can also be observed that the No-Melamchi and Melamchi scenarios are discretized into four and seven major simulation steps, respectively.

Future meteorology
The estimate of future rainfall was carried out from the analysis of historical meteorological data. Following the declining trend of the last decade (from 1710 to 1615 mm), the annual average rainfall was estimated such that the 10-year average of the estimate followed this declining trend (refer Figure 8). This annual precipitation was utilized along with the seasonal pattern of the available 2011À2012rainfall as input to all future scenarios. The annual PET was assumed to be constant at 1067 mm based on the annual PET of the last decade, previous literature, and meteorological database climate estimator software, LocClim. While, the seasonal pattern was again extracted from the 2011À2012PET data which was calculated following Monteith (1965).

Land subsidence estimation
After the completion of the scenario runs in SHETRAN, the simulated spatial data of the GWT is used in a GIS as a raster map. Cell-by-cell difference between two simulation times gives the spatial distribution of subsidence/ rise of GWT during the interval. The drop in pressure due to subsidence of GWT is expressed using the Terzaghis effective stress principle as dp ¼ Ds where s is the effective stress, g sand is the bulk specific weight of the bed containing phreatic surface, g sat is the saturated specific weight of the bed, and g w is the specific weight of water. In order to obtain the values of the specific weights, basic soil mechanics relationships of index properties are used: g ¼ ðG þ S r eÞ 1 þ e Â g w (8) where n is porosity, e is the void ratio, G is the specific gravity of solids, and S r is the degree of saturation. The values of porosity correspond to those from the calibrated groundwater model. The bulk specific weight corresponds to the soil that has become unsaturated due to falling phreatic level and lies almost immediately above the phreatic surface. There are several approaches for estimating land subsidence. The most popular being modelling using modular subroutines of MODFLOW (McDonald & Harbaugh 1988) such as Interbed Storage Package 1, 2, and 3 -IBS1, IBS2, IBS3, and the SUB package (Leake 1990, 1991, Leake & Prudic 1991. This involves coupling between the GW model and the subsidence model. However, this paper uses a non-coupled approach in which GW flow from SHETRAN is used in conjunction with a simplified form of a full analytical solution of land subsidence due to GW pumping. The following equation is the simplified form while the full analytical solution can be referred to in Bear and Cheng (2010): where d is the land subsidence corresponding to an increase in effective stress dp, n is the number of clay or aquifer inter-beds having compressibility a, and thickness z. This simplification is based on the assumption that only vertical compressibility prevails.
Using the available data points of storage coefficient (S) of the deep aquifer and the relationship between S » a, a raster of compressibility of deep aquifer is developed. However, such data corresponding to clay were unavailable. Thus the estimated subsidence corresponded to the compression of deep aquifer only which was heavily exploited for groundwater. Using these four raster (a, z, and two GWT rasters for dp), subsidence raster is estimated in GIS on a cell-by-cell basis employing the above-obtained simplified analytical expression for subsidence.

Data source and pre-processing
The data acquired and their sources are depicted in Table 7. The data are assessed for continuity and reliability. The continuity in data is achieved by infilling short in the data, especially for time series of rainfall. Assessing the reliability of data included, amongst others, scanning for erroneous high rainfall values, observing the lag times between the precipitation and its corresponding flood event.
The meteorological distribution array is prepared by Theissen polygon generated from stations locations. The spatial distribution of the vegetation and overland roughness are identical and both generated from the land cover map. Soil distribution array is formulated from interpolation of the three soil depths using Kriging method. In this study, 675 soil categories are used, each is a unique set of three soil depths. The GW abstraction locations were modelled as located at the centre of the pixel including it. All the processing involved in the preparation of GIS layers and images have been carried out using GIS software.

Results and discussions
We calibrated and validated a SHETRAN model to simulate groundwater flow in the study area; developed three future scenarios for land cover, two for groundwater abstraction pattern, and future meteorological conditions; and estimated land subsidence based on change in effective stress due to fluctuation in GWT for current and future conditions.

Model calibration and validation
The calibrated values of the model parameters are as shown in Table 8. The overland Strickler's coefficient varies from 0.25 for forest cover to 8 for densely populated urban cover. The values of 0.25 for forest and 1 for arable cover are comparable to the range of 1-10 applied for a grass and forest dominant cover using the same model in Bathurst et al. (1998). The calibrated value of saturated hydraulic conductivity, K (0.74 m/day) of the deep aquifer is in good agreement with values estimated in Pandey and Kazama (2011) (i.e. 0.32-8.78 m/day) and reported in Metcalf and Eddy (2000) (i.e. 0.51-8.16 m/day). The value for the shallow aquifer is 5.1 m/day which is slightly out of range as compared to 12.5-44.9 m/day as estimated by . The reason could be the high conductivity sand layer that is incorporated in the model for augmenting calibration. However, the almost 10 times larger value of K in horizontal as compared to that in vertical is in agreement with one of the previous attempts of GW modelling of the valley using MODFLOW (KC 2011). The initial values of the van Genutchen parameters a and b are taken from the SHETRAN manuals and are calibrated to adjust the soil moisture and conductivity characteristics. The value of specific storage (S s ) for the shallow aquifer is taken from the previous literature. Whereas, S s for deep aquifer is taken as mode of the values after interpolation of 12 data points over the valley . The values for the saturated (u sat ) and residual (u res ) moisture content are taken as provided in the table of soil properties (SHETRAN 2013a). For the calibration period, a good agreement is obtained between the simulated outlet discharge and the observations (refer Figure 9). This can be quantified by an NSE of 0.856 and volume deviation of 2.09 % for the calibration period. Application of the model to the validation period resulted in an NSE of 0.884 and 14.26 % volume deviation. All the calculations regarding discharge are made using hourly resolution of observed and simulated data. The volume deviation of discharge throughout the simulation period is found to be 6.67 % indicating that a robust model of catchment surface flow response has been obtained. The slightly greater volume deviation for the validation period comes from the pre-monsoon and initial monsoon period while the fit is remarkably good for the rest of the hydrograph. There seems to be some discrepancies between the observed and simulated discharge hydrograph in terms of timing of the peaks and the recession between the consecutive peaks. These are most probably due to routing errors borne from coarse resolution grid as in SHE-TRAN the river flows along the edges of the grid elements (Figure 4(a)). However, the result obtained can be considered as a good representation of the run-off response of the catchment under the diverse range of meteorological conditions of the Kathmandu Valley. For the validation of the model at the subsurface, the simulated phreatic levels are compared to the observed groundwater monitoring data at 19 GW monitoring stations (refer Figure 1). Comparative plots of the simulated phreatic levels versus the observed levels for selected months of each season with corresponding coefficient of correlation (r) are shown in Figure 10. The value of r for the overall simulation is 0.831 with a probable error of 0.047 and a standard error of 0.070. A low probable error suggests that the obtained values of r are reliable. This is an improvement over the previous attempt of distributed GW modelling of the valley (KC 2011). However, necessity of upgrading the existing GW monitoring system was sorely felt (number of stations, frequency of monitoring, and automation of monitoring).

Seasonal fluctuation in recharge
The modelled seasonal fluctuation of recharge is as shown in Figure 11. It gives us information on the movement of moisture in the GW basin. During the monsoon, almost entire catchment is under positive recharge due to the high amount of rainfall infiltration. In winter when there is low rainfall, the moisture movement is seen from the areas devoid of shallow aquifer to the areas possessing shallow aquifer. In the pre-monsoon, the model is showing a growth of positive recharge at southeast part of the catchment.

Land subsidence estimation for different scenarios
The estimated subsidence showed that out of 678 pixels of 700 m Â 700 m area of the GW basin almost 300 pixels did not show any subsidence at all. Statistical analysis of the subsidence exhibiting pixels for each scenario is done using box plots as illustrated in Figure 12. The average subsidence in the subsidence exhibiting pixels for the baseline scenario (III-nM) is about 1.6 mm/yr. Median values would better represent the central tendency of basin subsidence as the average values are higher even than the third quartile due to heavy influence from the outliers (red crosses). From comparison of median values (red horizontal lines inside the box), it is clear that scenario I is the worst one and that arrival of MWSP helps in controlling the subsidence to some extent. This is also verified by the increase in the number of outliers towards the upper whisker of the box plots after the introduction of Melamchi. The spatial distribution of estimated subsidence for No-Melamchi land-cover scenario III is presented in Figure 13(a). It can be seen that the peripheral regions are showing greater subsidence rather than the central urban core. There are two hotspots for subsidence vulnerability: north (Dhapasi) and north-east (Besigau). The latter shows subsidence rates up to 28 mm/yr. The occurrences of both the hotspots near the largest two clusters of KUKL deep abstraction wells (black dots) clearly show their direct impact on subsidence. It is evident that the arrival of the water supply project reduces the subsidence rates especially at the north-eastern hotspot and in the areas to the southwest of central urban core (Figure 13(b)). Comparative analysis of urban cover expansion patterns shows that scenario I in comparison to scenario III will result in greater subsidence rates at southern and south-eastern parts of the GW basin (refer Figure 13(c)). So sealing of recharge areas to the north is found detrimental to the southern parts of the valley. The complex geology and GW motion underneath the valley could be accredited for this spatial difference. Similar comparisons between scenarios II and III does not yield such a significant result (refer Figure 13(d)). This could be attributed to the fact that the motion of GW is most sensitive to the land cover of the recharge areas which is almost the same for these two scenarios.
With the estimated subsidence, Kathmandu valley is compared to other major subsidence events around the world in the past (Figure 14). The figure is a subsidence versus GW pumping graph where the latter is calculated using annual GW abstraction volume and the catchment area. It can be seen that with an average estimated deep aquifer subsidence of 1.6 mm/yr, Kathmandu valley is definitely not in the group of large subsidence locations of the world. However, it is the combination of groundwater pumping and vulnerable hydrogeology that produces subsidence, the latter of which has not been well incorporated in this study due to lack of comprehensive hydrogeological investigations.

Combined hazard map
Figure 15(a) shows recharge, groundwater subsidence, and land subsidence estimates in a single map for the complete picture of hazard induced by the excessive GW abstraction for the baseline scenario. The three locations at which subsidence have been reported in the past are also shown. There is spatial discrepancy between their location and model results. This indicates that it is either the clay compression behind these events or the model resolution is not suitable to simulate such localized effect. According to the subsidence expression, it is clear that high subsidence will occur with a larger declination of GWT, greater aquifer thickness (i.e. deep aquifer) and higher value of compressibility of the deep aquifer. Since a greater depth of clay aquitard is found at the central area of the valley, the GWT fluctuation here is minimum and so is the estimated subsidence. Moreover, from the recharge characteristics of the valley discussed previously, the central region of the valley is the most well-recharged portion of the GW basin. The greater compressibility of the deep aquifer (refer Figure 15(e)) and the GW subsidence are the reasons behind the two hotspots.

Summary and conclusions
This study is the first to predict land subsidence using a GW model of the Kathmandu Valley. This is a useful starting point for more detailed land subsidence modelling of the basin. The employed GW model is an improvement over the few previous attempts that have used distributed models of the study area, enabling better modelling of the subsidence vulnerability of the valley to be carried out. Pixel-wise analysis shows that most of the parts of GW basin are not subsiding. The average subsidence for the baseline scenario is found to be only 1.6 mm/yr with half of the values less than 0.38 mm/yr. Although this basin-level subsidence hazard seems to be low, the model shows high    deep aquifer subsidence vulnerability towards the peripheral areas of the groundwater basin. Two distinct hotspots are found in the northern areas exhibiting up to 28 mm of subsidence per annum. Analysis of various scenarios found that the completion of MWSP could control the subsidence rates to some extent. This was observed mainly in the north-eastern hotspot and in areas south-west of the central urban core. The estimated subsidence is most sensitive to the change in land cover of the recharge areas. Overall, the study shows that based on the estimations produced for future scenarios, the valley will not be in the list of the largest subsidence locations in the world. In spite of the lack of direct validation, the information from this model will be useful for planners and policy-makers in the future for achieving a sustainable GW environment and safety in the valley. In addition, the work can also be referred to for designing field investigations, establishing a subsidence monitoring network and upgrading the present groundwater monitoring network. In order to address the subsidence vulnerability of the central densely populated core of the valley, a more comprehensive model at higher resolution which also includes clay subsidence is needed. Such a detailed subsidence model could capture the localized effects missed by this model. For that a comprehensive database of the hydrogeology within the valley needs to be built. Another future area of research could be the analysis of the effect of ultra-urbanization on subsidence due to the increase Figure 15. (a) Combined hazard map for excessive groundwater abstraction; (b) shallow aquifer depth; (c) clay aquitard depth; (d) deep aquifer depth; (e) compressibility raster for deep aquifer with data locations for storage coefficient. Note that darker represents deep/high value. Source: Three locations of localized subsidence in the past (black dots in (a)) from Rana et al. (2007); aquifers-aquitard depths and data locations of storage coefficient from  in vertical load from multistory buildings. The present data scarce situation also favours a nonmodel, index-based approach that makes use of the existing information to address the subsidence vulnerability.

Disclosure statement
No potential conflict of interest was reported by the authors.