Article

Long-Term Subsidence Assessment by LiCSBAS and Emerging

Hot Spot Analysis in Kathmandu Valley
Sagar Rawal 1, * and Guoquan Wang 2

1 Department of Earth and Atmospheric Sciences, Science and Research Building, University of Houston,
3507 Cullen Blvd, Room 231, Houston, TX 77204, USA
2 Department of Earth and Atmospheric Sciences, University of Houston, 3507 Cullen Blvd, Room 231,
Houston, TX 77204, USA; gwang@uh.edu
* Correspondence: srawal@cougarnet.uh.edu; Tel.: +1-832-744-5848

Abstract: Rapid urbanization in Kathmandu Valley has strained its aquifer system, caus-
ing significant land subsidence. This study employs LiCSBAS for InSAR processing of
Sentinel-1 data (2017–2024) to map subsidence-prone areas. The significant subsidence was
found in northwest (Baluwatar, Samakhusi, and Manmaiju), southern (Gwarko, Patan, and
Koteshwor), and northeast (Madhapur Thimi and Gathhaghar) regions with a maximum
subsidence rate ~21 cm/yr. Subsidence has also expanded towards the outskirts and open
areas in the eastern and southern parts of Lalitpur and Bhaktapur districts. Emerging hot
spot analysis reveals a slowing subsidence trend in high-risk zones, possibly linked to the
MWSP project reducing groundwater extraction from 58 MLD (2021) to 26 MLD (2024).
Many subsidence-affected areas are located over the Kalimati and Gokarna Formations in
highly urbanized areas. The key contributing factors to subsidence are soil compaction,
excessive groundwater use, and urban sprawl encroaching open areas and recharge zones.
These findings underscore the urgent need for sustainable groundwater management and
land-use planning to promote urban resilience.

Keywords: space time cube; emerging hot spot analysis; land subsidence

Academic Editors: Emre Havazli,

Talib Oliver-Cabrera and 1. Introduction
Francesca Cigna
Land subsidence (LS) is a pervasive geohazard that presents formidable challenges to
Received: 24 February 2025
urban environments, mainly in densely populated regions such as Asia. This phenomenon
Revised: 16 March 2025
is characterized by the gradual downward displacement of the land surface due to excessive
Accepted: 22 March 2025
Published: 26 March 2025
groundwater extraction resulting in the compaction of soil layers. The subsequent reduction
in pore pressure compromises the soil’s structural integrity, diminishing its capacity to
Citation: Rawal, S.; Wang, G.
Long-Term Subsidence Assessment by
support overlying loads. This elevates the vulnerability of edifices, thoroughfares, and vital
LiCSBAS and Emerging Hot Spot infrastructure to damage while amplifying the propensity for urban pluvial inundation,
Analysis in Kathmandu Valley. Land landslides, sinkholes, and an array of cascading geo disaster [1–5]. Many urban cities have
2025, 14, 700. https://doi.org/ suffered land subsidence due to the exploitation of groundwaters, which could potentially
10.3390/land14040700
impact 29% of total urban population [1,6–12].
Copyright: © 2025 by the authors. Several major cities, including Prato in Italy, Wuhan in China, Jakarta in Indonesia,
Licensee MDPI, Basel, Switzerland. and Metro Manila in the Philippines, have reported substantial subsidence due to rapid
This article is an open access article
urbanization and the over-extraction of groundwater resources [13]. Kathmandu Valley
distributed under the terms and
(KV), comprising three district areas including the country’s capital city, is a bowl-shaped
conditions of the Creative Commons
Attribution (CC BY) license
alluvial basin surrounded by mountains. As a nexus of socioeconomic vitality, political
(https://creativecommons.org/ authority, educational advancement, and industrial enterprise, the valley magnetizes sub-
licenses/by/4.0/). stantial rural-to-urban migration, intensifying its urban sprawl. This demographic surge

Land 2025, 14, 700 https://doi.org/10.3390/land14040700

Land 2025, 14, 700 2 of 25

exerts profound pressure on finite resources of groundwater, agriculture land, transporta-

tion, residential buildings, and other infrastructures. Surrounded by towering mountains
and national parks to the north, KV faces severe spatial constraints for urban expansion,
driving extensive encroachment onto riverbanks, open spaces, and cropland. The rapid
proliferation of residential structures within the valley is staggering; between the late
1990s and 2021, the number of newly built houses in the 50-square-kilometer Kathmandu
Metropolitan City (KMC) surged to 70,090 an alarming figure nearly seven times the area’s
sustainable capacity [14]. Lax regulatory oversight and enforcement deficiencies have
emboldened unchecked construction activities and widespread groundwater extraction
through tube wells, exacerbating environmental degradation. KV is the most densely
populated area in the country with a staggering population density of 5169 individuals per
square kilometer [15]. Between 2015 and 2024, the demand for water in KV escalated by
36% due to limited water sources and growing urban population. Kathmandu Upatyaka
Khanepani Limited (KUKL; the valley’s official water management authority) report shows
that groundwater pumping rates increased from 2.3 MLD in 1979 to 29.2 MLD in 1999,
peaking at 58.79 MLD in 2021. The rate of groundwater abstraction has surpassed the
replenishment rate since the mid-1980s, with this disparity continuing to widen over time.
The completion of the Melamchi Water Supply Project (MWSP) and additional sourcing
from the Bagmati River has enabled KUKL to increase its water supply to 47% of the total
demand (506 MLD) in 2024. This represents a notable advancement over the preceding
years where KUKL only managed to supply 28.8% in 2023 and a mere 21% in 2022 [16].
Despite these improvements, frequent landslides and floods disrupt water distribution
from MWSP and river sources, leading to prolonged supply shortages and exacerbating
Kathmandu Valley’s chronic water deficit. Due to a shortage in water supply, most residents
remain dependent on private vendors, third-party water suppliers, and unregulated tube
wells depleting the water level up to ~2.5 m annually [17]. Nepal ranks 11th in earthquake
vulnerability and 16th for multiple hazards globally as per the United Nation Development
program (UNDP) [18]. The most recent seismic disaster of 2015 Gorkha Earthquake (Mw
7.2) caused widespread devastation, inducing substantial liquefaction in KV and damaging
critical infrastructure [19]. If a future seismic event strikes with an epicenter near Kath-
mandu, the ongoing land subsidence will worsen the disaster. It will increase structural
collapses and heighten economic losses. Without strategic intervention, the combination
of subsidence and seismic hazards endangers the valley’s long-term resilience. This situ-
ation demands urgent mitigation measures. These measures must address groundwater
depletion and enforce sustainable urban planning policies.
Traditional methods like GNSS and leveling methods are most commonly used tech-
nique to precisely measure subsidence worldwide. HoustonNet operates and processes
over 250 permanent GNSS stations to monitor subsidence-related hazards with high pre-
cision in Houston, Texas [5]. But their prohibitive costs, labor-intensive nature, frequent
upkeep demands, and limited spatial coverage render them impractical for widespread
adoption in resource-scarce settings like Kathmandu Valley. Recent strides in remote sens-
ing have redefined large-scale geohazard surveillance, proffering cost-efficient, expansive
alternatives. Technologies including Unmanned Aerial Vehicles (UAVs), photogrammetry,
satellite imagery, Light Detection and Ranging (LiDAR), and Interferometric Synthetic
Aperture Radar (InSAR) have demonstrated versatility across domains from disaster relief
to ecological oversight and geohazard evaluation. Among these, InSAR has ascended as a
pre-eminent tool for discerning slow-evolving geohazards such as glacial shifts, landslides,
and subsidence with millimeter-scale fidelity.
Multiple studies have employed InSAR to investigate subsidence in Kathmandu
Valley. Bhattarai (2011) pioneered the application of Differential InSAR (DInSAR) using
Land 2025, 14, 700 3 of 25

ALOS PALSAR imagery (2007–2010), identifying seven distinct subsidence zones with
rates ranging from 1 to 9.9 cm/year [20]. Following the 2015 Gorkha Earthquake (Mw 7.2),
Krishna (2018) analyzed 16 ALOS-1 PALSAR and 26 Sentinel-1 A/B SAR data to measure
subsidence velocity before and after the event. Their result revealed an increase in sub-
sidence rates from 8.2 cm/year prior to the earthquake to 12 cm/year afterward. This
escalation underscores the seismic influence on subsidence dynamics [21]. Kumahara [22]
examined four active faults in a field visit on July 2025 (3 months prior to earthquake) and
created a detailed DEM from height of 10 m using Agisoft Photoscan Professional. His
findings found no evidence of surface ruptures, suggesting persistent, unrelieved tectonic
stress that could elevate the risk of future seismic events. Stallin (2017–2019) employed
Persistent Scatterer Interferometry (PSI) on descending-orbit SAR datasets and found a
maximum subsidence rate of 13 cm/year in the line-of-sight (LOS) direction [23]. More
recently, Huang et al. (2024) utilized Small Baseline Subset (SBAS) InSAR techniques on
Sentinel-1 data (2015–2021), reporting peak subsidence rates of 20 cm/year. He analyzed
the InSAR-derived subsidence rate based on the river channel’s DTM and found that
vertical velocity increased with elevation. Their findings showed a reduced subsidence
rate in areas near river streams, with seasonal fluctuation in time series providing crucial
evidence of connection between shallow aquifers and nearby streams [24].
Research on subsidence in Kathmandu Valley has primarily focused on measuring
rates with limited consideration of geological structures, lithological variations, or long-
term trends. Many InSAR studies rely on data from either ascending or descending orbits,
using line-of-sight displacement to estimate vertical motion. Combining both orbits allows
for more accurate vertical displacement measurements. Most existing studies cover less
than three years, failing to capture sustained patterns or underlying causes. Given the rapid
rate of subsidence and growing seismic risk, continuous monitoring and recalibration of
historical data are critical. Spatial variability across different geological formations and land-
use zones also requires further investigation. Previous validation efforts have depended on
a single GPS station, NAST (27.657, 85.328), which suffers from limited spatial coverage,
power outages, and maintenance issues. In this study, we applied Persistent Scatterer
InSAR (PS-InSAR), a histogram method, and GPS observations to assess the accuracy and
reliability of InSAR-derived results. We processed InSAR data from both ascending and
descending orbits from 2017 to 2024 to extract more precise vertical displacement. The
two years following the earthquake are excluded as post-seismic relaxation could distort
subsidence rates and introduce outliers into the analysis.
Emergency hot spot analysis (EHA) from ESRI is an effective tool for spatiotemporal
analysis, facilitating the identification of emerging patterns through comparison with histor-
ical data [25]. This approach offers a more efficient alternative to the time-intensive process
of examining the time series of thousands of pixels to detect trends. EHA conducts spa-
tiotemporal analysis by constructing a Space–Time Cube (STC) from input data, efficiently
isolating statistically significant regions based on their temporal evolution. Shuhab D. Khan
et al. (2021) fused EHA with InSAR data (2015–2019) and seismic records to probe the nexus
between seismicity and subsidence in Peshawar, Pakistan, pinpointing hot spots in densely
populated urban zones where groundwater overexploitation predominated [26]. Similarly,
Wei Zhang et al. (2023) integrated EHA with Persistent Scatterer InSAR (PS-InSAR) to
investigate settlement patterns in Beijing, uncovering persistent hot spot clusters that sig-
naled latent subsidence foci overlooked by conventional techniques [10]. The versatility
of EHA extends beyond subsidence, encompassing time-dependent phenomena such as
landslides, epidemiological outbreaks like COVID-19, environmental calamities, and stock
market fluctuations. This study pioneers the integration of EHA with LiCSBAS outputs,
Land 2025, 14, 700 4 of 25

aiming to discern shifts in subsidence trends within the study area, thereby advancing the
analytical framework for geohazard assessment in KV.
The main objective of this study is to (1) examine subsidence in KV and refine existing
records with updated findings; (2) find any changing pattern using EHA tool to assess
whether governmental interventions have effectively stabilized subsidence rates in specific
areas; (3) use the PS-InSAR and histogram methods along with GNSS to check their
reliability; (4) analyze subsidence pattern across the hydrogeological cross-sections by
integrating borehole lithology data along spatial profiles. This research will contribute
to the scientific understanding of land subsidence in Kathmandu Valley and provide
data-driven insights to help in future urban planning and hazard mitigation.

2. Materials and Methods

2.1. Study Area
The KV, located in the Lesser Himalayas, is characterized by a predominantly flat
terrain with some undulating features. The valley floor lies at an elevation ranging from
1300 to 1400 m above sea level. The KV consists of lacustrine and fluvial sediments deposits
of remnant ancient lakes over a basement rock [27]. The bedrock is composed of metasedi-
mentary rocks and exposed at the surface in few places, particularly around Pashupatinath
and Swayambhu. The sediment basin underlying the KV consists of lacustrine and fluvio-
deltaic deposits [28]. There is a spatial distribution of unconsolidated sandy layers in the
northern part, silty-clay in the central part, and sandy-gravel layers in the southern part of
the basin [27]. The sediment thickness gradually increases towards the south, reaching its
maximum in the central (650 m) and southern parts of the basin. Sakai proposed a new
classification of sediment stratigraphy such as Tarebhir, Lukundol, and Itaiti Formations in
the southern part; Bagmati, Kalimati, and Patan Formations in the central part; and Thimi
and Gokarna Formation in the northern part of the groundwater basin [27] (Table 1).

Table 1. Classification of principle of hydrogeological units. Stratigraphic units of valley-fill sediments

modified after ([27,29,30]).

Hydrogeological Geological
Region Age Sediment Types
Units Formations
Itaiti Southern 1 myr–present Gravel, sand, silty sand
Sand, gravel, sand clay,
Patan Central 19 kyr–11 kyr
Shallow Aquifer gravelly clay
Thimi (29 kyr–19 kyr),
Thimi and Gokarna Northern Gokarna Sand, sandy silt
(1 myr–29 kyr)
Clay, clayey sand,
Lukundol Southern 2.8 myr–1 myr
silty clay
Aquitard
Clay, clayey sand,
Kalimati Central 2.8 myr–30 kyr
lignite, silt, gravelly clay
Sand, sandy gravel,
Tarebhir Southern Older than 2.8 myr
sandy silt
Deep Aquifer
Sand, gravel, gravelly
Bagmati Central Older than 2.8 myr
clay, sand, gravelly silt

JICA categorized the valley into three distinct groundwater districts in 1990: north,
central, and south zones shown in Figure 1a. The north and foothills of the south region
mostly consist of a recharge zone due to minimum thickness or absence aquitard. Most of
Land 2025, 14, 700 5 of 25

Land 2025, 14, x FOR PEER REVIEW the urban areas lie in this central district zone and have a significant number of deep5 wells
of 25

used by industries and hotels mostly. This zone is also characterized by unconsolidated
coarse sediment at the bottom covered by thick black clay. This highly compressible and
impermeable
impermeable layer separates
separates the
the shallow
shallow and
and deep
deep confined
confined aquifer,
aquifer, with thickness
thickness range
range
from
from 5 m to more than 200 m. The annualannual average
average recharge
recharge rate
rate of
of acquirers
acquirers in
in valley
valley is
estimated
estimated to be 9.6 MCM/year
MCM/year[31].[31].

1. (a)
Figure 1.
Figure (a) Modified
Modifiedhydrogeology
hydrogeologyofofKathmandu
KathmanduValley [32].
Valley Fault
[32]. lines
Fault areare
lines represented in red
represented in line
red
modified from Asahi [33], and Green dashed line separates the three recharge zones; South,
line modified from Asahi [33], and Green dashed line separates the three recharge zones; South,Central
and Northern zones as defined by JICA [32]. (b) Schematic geological cross-section of profile CD
Central and Northern zones as defined by JICA [32]. (b) Schematic geological cross-section of profile
through Central Nepal modified from Cresswell [34].
CD through Central Nepal modified from Cresswell [34].
The tectonic activity has given rise to multiple faults in the region such as Thankot
The tectonic activity has given rise to multiple faults in the region such as Thankot
Fault, Chandragiri Fault, Thimi Fault, and Patan Bhaktapur Fault. Among these, only the
Fault, Chandragiri Fault, Thimi Fault, and Patan Bhaktapur Fault. Among these, only the
Thankot Fault remains tectonically active, with an observed vertical displacement rate of
Thankot Fault remains tectonically active, with an observed vertical displacement rate of
1 mm/yr [33]. Given the region’s high seismic activity, these faults have the potential to
1 mm/yr [33]. Given the region’s high seismic activity, these faults have the potential to
be reactivated by seismic events, leading to localized ground deformation [35]. Therefore,
be reactivated by seismic events, leading to localized ground deformation [35]. Therefore,
studying subsidence in the region is essential for ensuring infrastructure resilience and
studying subsidence in the region is essential for ensuring infrastructure resilience and
enhancing disaster preparedness to effectively mitigate earthquake hazards.
enhancing disaster preparedness to effectively mitigate earthquake hazards.
2.2. Interferometric Synthetic Aperture Radar
2.2. Interferometric Synthetic Aperture Radar
InSAR is a powerful technique for monitoring land deformation, utilizing satellite
InSAR
imagery to is a powerful
generate technique maps
displacement for monitoring
by analyzinglandphase
deformation, utilizing
differences acrosssatellite
multi-
imagery to generate displacement maps by analyzing phase differences
ple interferograms. It offers millimeter-level accuracy and is widely employed to detect across multiple
interferograms.
surface changes,It fault
offersmovements,
millimeter-level accuracyand
subsidence, anduplift
is widely
causedemployed
by bothtonatural
detect sur-
and
face changes, fault movements, subsidence, and uplift caused by both natural
human-induced factors. Persistent Scatter (PS-InSAR)and small baseline subset (SBAS) are and human-
induced
commonly factors. Persistent
used InSAR Scatter
method for(PS-InSAR)and
measuring slowsmall baseline subset (SBAS) are com-
deformation.
monly used InSAR method for measuring slow deformation.
For our study, we will use the SBAS InSAR technique to study the land subsidence
that For
usesour study, stacked
multiple we will use the SBAS
master slave InSAR
pairs oftechnique
Sentinel-1 to Satellite
study thedata.
land The
subsidence
results
that uses multiple stacked master slave pairs of Sentinel-1 Satellite
will be validated for consistency and accuracy using GNSS, PS-InSAR, and histogram data. The results will
be validated
methods. Thefor consistency
free and accuracy
InSAR processing using GNSS,
tool LiCSBAS PS-InSAR,
[36] is and histogram
used to download meth-
the LiCSAR
ods. The free InSAR processing tool LiCSBAS [36] is used to download
products (Unwrapped interferograms), which are publicly available to download through the LiCSAR prod-
ucts (Unwrapped
COMET-LiCS webinterferograms), which are publicly available to download
portal (https://comet.nerc.ac.uk/COMET-LiCS-portal/, through
accessed on
COMET-LiCS web portal (https://comet.nerc.ac.uk/COMET-LiCS-portal/,
22 October 2024). The datasets span for 7.6 years, from Jan 2017 to July 2024. The frame accessed on 22
October 2024). The datasets span for 7.6 years, from Jan 2017 to July 2024. The frame IDs
085A_06253_131313 and 019D_06217_131313 were used to download both ascending and
descending (Figure 2) interferograms. These interferograms cover 250 km swath at 5 m by
20 m spatial resolution, geocoded onto a 100 m grid.
Land 2025, 14, 700 6 of 25

Land 2025, 14, x FOR PEER REVIEW IDs085A_06253_131313 and 019D_06217_131313 were used to download both ascending 6 of 25
Land 2025, 14, x FOR PEER REVIEW 6 of 25
and descending (Figure 2) interferograms. These interferograms cover 250 km swath at 5 m
by 20 m spatial resolution, geocoded onto a 100 m grid.

Figure 2. Geographical location of study area. The area with white borders represents the clipped
Figure 2. Geographical location of study area. The area with white borders represents the clipped
Figure
area for2. Geographical
for our study. location and
study. The of study area.
boxThe area withthe
white bordersdiagram
representsof the
twoclipped
area our The red green
red and green box represents
represents schematic
the schematic diagram of two frames
frames
area for our study. The
019D_06217_131313 and red and green box in
085A_06253_131313 represents theand
descending schematic diagram
ascending orbits of COMET-LiCS
in two frames
019D_06217_131313 and 085A_06253_131313 in descending and ascending orbits in COMET-LiCS
019D_06217_131313
WebPortal
Web [37]. and 085A_06253_131313 in descending and ascending orbits in COMET-LiCS
Portal[37].
Web Portal [37].
During initialprocessing
During the initial processingsteps
stepsofof InSAR,
InSAR, 849849 datasets
datasets fromfrom descending
descending orbits
orbits and
and During
983 983 from
from the initial
ascending
ascending processing
orbits
orbits steps
were
were of InSAR,
retainedafter
retained 849 datasets
afterremoving from
removinginitially descending
initially faulty orbits and
faulty interferograms
interferograms
983 from ascending
identified
identified through aaorbits
through visualwere
visual retained
inspection
inspection of after
of the removing
the png
png file initially
file format
format of faulty
ofeach
each interferograms
interferogram.
interferogram. The
The
identified
network
network of through a
of available visual inspection
available interferometric of
interferometric pairs the
pairs alongpng
along withfile format
with their of each interferogram.
their temporal-baseline
temporal-baseline information The
information
network
of of available
of interferograms
interferograms interferometric
used
used for
for SBAS
SBAS and pairs
and along with
PS-InSAR
PS-InSAR theirin
isisshown
shown temporal-baseline
in Figure
Figure3.3. information
of interferograms used for SBAS and PS-InSAR is shown in Figure 3.

Figure
Figure 3. Spatialand
3. Spatial andtemporal
temporalbaseline
baselineofofInSAR
InSARanalysis(a)
analysis (a) SBAS;
SBAS; ascending
ascending orbits,
orbits, (b) (b) SBAS;
SBAS; de-
descending
Figure 3. orbits,
Spatial and and (c)
temporalPSI; ascending
baseline of
scending orbits, and (c) PSI; ascending orbits.orbits.
InSAR analysis(a) SBAS; ascending orbits, (b) SBAS; de-
scending orbits, and (c) PSI; ascending orbits.
The atmospheric correction was applied using GACOS external files. During loop
The atmospheric
closure and correction
quality checks, with was appliedthreshold
a coherence using GACOS
of 0.05,external files. During
63 interferograms loop
from the
closure and quality
descending checks,
frame and withthe
26 from a coherence
ascendingthreshold of 0.05,
frame were 63 interferograms
identified as faulty andfrom the
excluded
descending frame and 26 from the ascending frame were identified as faulty and excluded
Land 2025, 14, 700 7 of 25

The atmospheric correction was applied using GACOS external files. During loop
Land 2025, 14, x FOR PEER REVIEW 7 of 25
closure and quality checks, with a coherence threshold of 0.05, 63 interferograms from the
descending frame and 26 from the ascending frame were identified as faulty and excluded
from further processing. The remaining interferograms were processed using the small
from further
baseline processing.
inversion approachThetoremaining interferograms
generate velocity results,were
whileprocessed
time-seriesusing the small
analysis was
baseline inversion
performed using theapproach to generate
least squares method for velocity results,
each pixel, whileminimizing
further time-seriesresidual
analysisnoise.
was
performed using the least squares method for each pixel, further minimizing
LiCSBAS selects a stable, high-coherence location as the reference area. The same reference residual
noise.
area LiCSBAS
is utilized inselects
PS InSARa stable, high-coherence
to ensure location
consistent results asto
prior the referencebetween
validation area. The
thesame
two
reference area is utilized in PS InSAR to ensure consistent results prior
methods. The InSAR results are in a line-of-sight (LOS) direction and not direct verticalto validation be-
tween the
motion. two methods.
Equation Theto
(1) is used InSAR results ascending
decompose are in a line-of-sight (LOS) LOS
and descending direction andinto
velocity not
direct vertical motion. Equation (1)
horizontal and vertical components [38]. is used to decompose ascending and descending LOS
velocity into horizontal and vertical components [38].
! ! −1 !
Vu 𝑉 −cosθ
𝑐𝑜𝑠𝜃 a − cosα sinθ
𝑐𝑜𝑠𝛼 𝑠𝑖𝑛𝜃
a a V𝑉
asc
= (1)
Vh 𝑉 𝑐𝑜𝑠𝜃
cosθ d 𝑐𝑜𝑠𝛼
−cosα 𝑠𝑖𝑛𝜃
d sinθ d V𝑉
dsc

where 𝜃 is the incidence angle, 𝛼 is the azimuth angle. 𝑣 and 𝑣 are the LOS de-
where θ is the incidence angle, α is the azimuth angle. v asc and vdsc are the LOS deformation
formation in ascending and descending orbit direction. 𝑣 and 𝑣 are the vertical and
in ascending and descending orbit direction. vu and vh are the vertical and horizontal
horizontal displacement components.
displacement components.

2.3. Groundwater
2.3. Groundwater Level
Level Data
Data
Groundwater levels
Groundwater levels can
can influence
influence both
both elastic
elastic and
and inelastic
inelastic land
land deformation.
deformation. When
When
the water
the water level
level remains
remains above
above the
the pre-consolidation head, only
pre-consolidation head, only elastic
elastic deformation
deformation occurs.
occurs.
However, if the water level drops below the pre-consolidation head, permanent
However, if the water level drops below the pre-consolidation head, permanent deforma- defor-
mation
tion takes
takes place.
place. A hydrological
A hydrological cross-section
cross-section of Kathmandu
of Kathmandu ValleyValley reveals
reveals a three-
a three-tiered
tiered aquifer system comprised of an unconfined shallow aquifer, an aquitard,
aquifer system comprised of an unconfined shallow aquifer, an aquitard, and a confined and a con-
fined aquifer
aquifer (Figure(Figure
4b). 4b).

Figure 4.
Figure 4. (a)
(a) Distribution
Distribution of wells in Kathmandu Valley.
Valley. The
The red
red label
label represents
represents the
the borehole
borehole along
along
profileXY
profile XYwhose
whosecross-section
cross-sectionisisshown
showninin(b)
(b)The represents the
The bar graph represents the water
water level
level from
from 2014
2014 to
to
2020
2020 of
ofselected
selectedtubewell
tubewellA,A,B,B,C,C,DDand E. E.
and (b)(b)
Hydrogeological cross-section
Hydrogeological along
cross-section spatial
along profile
spatial XY
profile
modified
XY modifiedafterafter
Pandey [39].[39].
Pandey

The
Theshallow
shallowaquifer
aquifer(0–60
(0–60m),
m),primarily utilized
primarily forfor
utilized household purposes,
household is widely
purposes, dis-
is widely
tributed across
distributed the valley,
across while
the valley, deeper
while wellswells
deeper (>60 (>60
m), mainly used used
m), mainly for industrial purposes,
for industrial pur-
are concentrated in the central basin [39]. The two shallow and confined deep
poses, are concentrated in the central basin [39]. The two shallow and confined deep aquifers
aq-
uifers are separated by a silt-rich clay layer, with varying thicknesses (5 to 200 m) [32,40–
42]. The shallow aquifer gradually decreases in thickness from the northern valley to-
wards the center, while the deep aquifer is thickest in the central and southern regions
[41]. Spatial and temporal land-use and land-cover (LULC) analysis has revealed that
Land 2025, 14, 700 8 of 25

are separated by a silt-rich clay layer, with varying thicknesses (5 to 200 m) [32,40–42]. The
shallow aquifer gradually decreases in thickness from the northern valley towards the
center, while the deep aquifer is thickest in the central and southern regions [41]. Spatial
and temporal land-use and land-cover (LULC) analysis has revealed that Kathmandu
Valley has experienced the highest increase in built-up areas in the country. Over the last
30 years (1990–2020), built-up areas have expanded by 368.06% [43], significantly reducing
natural groundwater recharge zones. This urban expansion has also limited open spaces,
preventing rainwater from infiltrating the soil and replenishing groundwater reserves.
While monsoon rainfall recharges shallow aquifers, leading to a temporary minor rebound
in land elevation but prolonged droughts, soil compaction and unregulated groundwater
extraction without adequate replenishment can result in inelastic subsidence [44].
The total reported number of tube wells was 60 in 1989 [45], which increased to
218 according to a report by Metcalf and Eddy in 2000 [41]. By 2009, a KUKL report
recorded 379 tube wells within the valley. Currently, 104 deep tube wells extract approx-
imately 54 MLD (Million Liters per Day) of groundwater [16]. Understanding ground-
water dynamics is essential for evaluating subsidence risks. However, historical ground-
water level data are scarce. There is very limited groundwater level data available pri-
marily gathered by independent researchers, NGOs, and government agencies in frag-
mented efforts. The primary sources of water well data and locations were obtained from
(https://gw-nepal.com, accessed on 11 December 2024) and Smartphones4Water Nepal
(https://s4w-nepal.smartphones4water.org, accessed on 17 November 2024). Additional
undocumented well data were collected from independent researchers and local ward
offices. Figure 4a shows the distribution of wells in Kathmandu Valley, highlighting a
high concentration of deep wells in the northwestern and central regions. These areas are
also highly urbanized, with a significant presence of industries, commercial and education
institutional buildings. Despite gaps in water level data for many wells, five sites were
selected based on their continuous and long-term data availability. Figure 4a shows the
water level data for each well in February from 2014 to 2020, with the most significant
decline observed in Well D, located in the northwestern part of the valley. In this study,
InSAR time-series analysis will be conducted at these five sites to examine the correlation
between subsidence and groundwater levels. The recent well data will also be analyzed to
determine water level trends, further exploring the relationship between subsidence and
groundwater depletion.

2.4. Emerging Hot Spot Analysis

Emerging hot spot analysis (EHA) is a powerful tool in ArcGIS that utilizes the
Space–Time Cube (STC) to detect patterns in InSAR-derived subsidence data over a given
study period [26,46,47]. The STC is a three-dimensional framework for analyzing spatial
and temporal data, stored as discrete bins. In this structure, the XY plane represents
spatial dimensions, while the Z-axis represents time. The STC is stored in netCDF format,
organizing data in an (x, y, z, t) array, where x and y denote spatial coordinates, t represents
time, and z stores attributes such as subsidence magnitude, precipitation, or slope. To
analyze trends within each bin, the Mann–Kendall statistical test is applied, which evaluates
the direction and statistical significance of temporal changes using z-scores and p-values. A
negative z-score indicates a decreasing trend, while a positive z-score signifies an increasing
trend. A low p-value suggests statistical significance. The EHA utilizes the netCDF-based
STC as an input, applying Getis-Ord Gi* statistics to identify spatial and temporal patterns
in subsidence data. This method assesses whether a feature is a significant hot spot (clusters
of high values) or cold spot (clusters of low values) based on its neighboring bins. Figure 5
Land 2025, 14, 700 9 of 25

Land 2025, 14, x FOR PEER REVIEW 9 of 25

illustrates the EHA process, where individual bins with a temporal dimension along the
Z-axis are combined to form a STC.

Figure 5. Emerging
Figure 5. Emerging hot
hot spot
spot analysis
analysis workflow
workflow diagram
diagram modified
modified[48,49].
[48,49].

InSAR-derived displacement data were initially stored in .csv format, with latitude
InSAR-derived displacement data were initially stored in .csv format, with latitude
and longitude in the first two columns and cumulative displacement values recorded across
and longitude in the first two columns and cumulative displacement values recorded
multiple SAR acquisition dates. The dataset was transposed and pre-processed for STC
across multiple SAR acquisition dates. The dataset was transposed and pre-processed for
formation, using mean displacement as the summary field. A 200 m spatial resolution
STC formation, using mean displacement as the summary field. A 200 m spatial resolution
was selected, averaging all points within each grid cell to create temporal slices preserving
was selected, averaging all points within each grid cell to create temporal slices preserving
displacement trends. To analyze subsidence patterns, the 3D Space–Time Cube Explorer
displacement trends. To analyze subsidence patterns, the 3D Space–Time Cube Explorer
was used for visualizing individual bins and spatial trends. The generated STC served as
was used for visualizing individual bins and spatial trends. The generated STC served as
an input for EHA, where individual time steps were examined instead of the entire cube to
an input for EHA, where individual time steps were examined instead of the entire cube
avoid misinterpretations from broader regional trends. For this study, eight key patterns
to avoid misinterpretations from broader regional trends. For this study, eight key pat-
were selected to accurately represent subsidence trends (Table 2).
terns were selected to accurately represent subsidence trends (Table 2).
Table 2. Emergency hot analysis patterns.
Table 2. Emergency hot analysis patterns.
S.N EHA Pattern Observation
S.N EHA Pattern Observation
1. Intensifying hot spot 1. subsidence
High Intensifying hot spot overtime
rate, increasing High subsidence rate, increasing overtime
2. Persistent hot High subsidence rate, constant over time
2. Persistent hot High subsidence rate, constant over time
3. Diminishing hot High subsidence rate, decreasing over time
3. Diminishing hot High subsidence rate, decreasing over time
Alternates between high and low subsidence
4. Oscillating hot spot 4.
Alternates Oscillating
between highhot spot
and low subsidence phases with less than 90% high rate
phases with less than 90% high rate
5. Sporadic hot spot Random
5. subsidence
Sporadic rate with no fixedRandom
hot spot pattern subsidence rate with no fixed pattern
6. Diminishing hot spot Subsidence
6. rate decreasing
Diminishing hotover
spottime Subsidence rate decreasing over time
7. New hot spot 7.
Recently New hot
emerged as spot
a high subsidenceRecently
area emerged as a high subsidence area
8. Consecutive hot spot 8. Consecutive hot spot Continuous
Continuous presence of high subsidence rate presence of high subsidence rate

2.5.
2.5. Geological
Geological Composition
Composition and
and Land-Use
Land-Use Data
Data
Land
Land subsidence
subsidenceisissignificantly
significantlyinfluenced
influencedby the
by geological composition
the geological of the of
composition study
the
area.
studyGroundwater
area. Groundwaterexploitation, uncontrolled
exploitation, urban expansion,
uncontrolled large-scale
urban expansion, construction
large-scale con-
activities
struction and presence
activities and of unconsolidated
presence sedimentssediments
of unconsolidated are key factors
are keycontributing to subsi-
factors contributing
dence processesprocesses
to subsidence in urban environments. To investigate
in urban environments. these influences,
To investigate these we will overlay
influences, the
we will
InSAR-derived subsidence map
overlay the InSAR-derived with both
subsidence map thewith
geological
both the formation
geologicalmap and the map
formation land-use
and
map, generated
the land-use map,through GIS-based
generated through spatial analysis.
GIS-based Pandey
spatial and Kazama
analysis. Pandey conducted
and Kazama a
comprehensive hydrogeological study to map the spatial distribution
conducted a comprehensive hydrogeological study to map the spatial distribution of aq- of aquifer charac-
teristics, including the
uifer characteristics, thickness,
including thetransmissivity (T), hydraulic
thickness, transmissivity conductivity,
(T), and storage
hydraulic conductivity,
coefficient
and storage coefficient of both shallow and deep aquifers. The study involvedaconstruct-
of both shallow and deep aquifers. The study involved constructing hydroge-
ologic cross-section, differentiating
ing a hydrogeologic shallow aquifers,
cross-section, differentiating aquitards,
shallow and deep
aquifers, aquifers
aquitards, andbased
deep
on sediment stratigraphy from eight boreholes along profile line XY (Figure
aquifers based on sediment stratigraphy from eight boreholes along profile line XY (Fig- 4b) [39]. To
study
ure 4b)a [39].
correlation
To study between sediment
a correlation thickness
between and subsidence,
sediment thickness andwe will analyzewe
subsidence, InSAR
will
time-series trends at these borehole sites. For a broader regional assessment,
analyze InSAR time-series trends at these borehole sites. For a broader regional assess- we have
constructed
ment, we have columnar sections
constructed using additional
columnar borehole
sections using data from
additional GWDP
borehole andfrom
data JICAGWDP
along
profile line AB (south–north) and MN (west–east) (Figure 6a). This will
and JICA along profile line AB (south–north) and MN (west–east) (Figure 6a). This will allow us to examine
allow us to examine the spatial relationship between geological formations, land-use pat-
terns, and subsidence rates. Figure 6b illustrates the stratigraphic column along the west–
east (WE) section, revealing that clay thickness is greatest in the central region, a trend
Land 2025, 14, 700 10 of 25

Land 2025, 14, x FOR PEER REVIEW 10 of 25

the spatial relationship between geological formations, land-use patterns, and subsidence
rates. Figure 6b illustrates the stratigraphic column along the west–east (WE) section,
also observed
revealing that in thethickness
clay south-north (SN) cross-section
is greatest in the central(Figure
region,6c). The low-porosity,
a trend also observedhighly
in the
impermeable(SN)
south-north black clay layer within
cross-section the
(Figure aquitard
6c). zone becomes
The low-porosity, moreimpermeable
highly dominant toward
black
the center of the valley.
clay layer within the aquitard zone becomes more dominant toward the center of the valley.

Figure 6. (a) Hill Shed Map of study area, where blue represents the majormajor rivers.
rivers. Bagmati river is
the
the primary
primary drainage of valley
drainage of valley and
and main
main outlet
outlet for
for water
water flow
flow from
from valley.
valley. The
The sites
sites of
of boreholes
boreholes
used
used to construct columnar sections of sediment lithology are shown in yellow and red circles along
to construct columnar sections of sediment lithology are shown in yellow and red circles along
NS
NS and
and SW
SW direction
direction in
in map.
map. (b)
(b) The
The columnar
columnar section
section along
along MN.
MN. (c)
(c) Columnar
Columnar section
section along
along AB.
AB.

Esri has made the global land-use map available using Sentinel-2 data at a 10 m
Esri has made the global land-use map available using Sentinel-2 data at a 10 m res-
resolution, processed through Impact Observatory’s deep learning AI land classification
olution, processed through Impact Observatory’s deep learning AI land classification
model. The model was trained on five billion hand-labeled Sentinel-2 pixels collected
model. The model was trained on five billion hand-labeled Sentinel-2 pixels collected from
from 20,000 sites representing major global biomes. The maps are available from 2017
20,000 sites representing major global biomes. The maps are available from 2017 to 2023,
to 2023, classifying land cover into nine distinct categories. They are distributed under
classifying land cover into nine distinct categories. They are distributed under the Crea-
the Creative Commons Attribution (CC BY 4.0) license, requiring proper attribution. The
tive Commons Attribution (CC BY 4.0) license, requiring proper attribution. The web ap-
web applications can be accessed from (https://livingatlas.arcgis.com/landcoverexplorer/,
plications can be accessed from (https://livingatlas.arcgis.com/landcoverexplorer/, ac-
accessed on 20 January 2025). The webapp enables users to visualize global land-use
cessed on 20 January 2025). The webapp enables users to visualize global land-use cover,
cover, side-by-side comparisons across different years and animation features to observe
side-by-side comparisons across different years and animation features to observe tem-
temporal changes. For research and analysis, GeoTIFF files of the study region are available
poral changes. For research and analysis, GeoTIFF files of the study region are available
for download. This accessibility reduces computational time and resource consumption,
for download. This accessibility reduces computational time and resource consumption,
making it a valuable tool for large-scale environmental and land-use studies.
making it a valuable tool for large-scale environmental and land-use studies.
3. Results
3. Results
3.1. Land Subsidence Results
3.1. Land Subsidence Results
To analyze the spatial and temporal distribution of ground subsidence in the study
To ascending
region, analyze theand spatial and temporal
descending distribution
unwrapped of ground were
interferograms subsidence in the study
downloaded from
region, ascending and descending unwrapped interferograms were downloaded
LiCSAR for the period 2017 to 2024. The interferograms were processed using LiCSBAS from
LiCSAR
batch for the period
processing, 2017 to 2024.
incorporating TheGACOS
external interferograms were processed
files to generate velocity using LiCSBAS
maps and time-
batch processing, incorporating external GACOS files to generate velocity maps
series data. The results from the ascending and descending images were decomposed using and time-
series data.
Equation (1)The resultsvertical
to derive from the ascending and
displacements descending
(Figure 7a). images were decomposed us-
ing Equation (1) to derive vertical displacements (Figure 7a).
Land 2025, 14, x FOR PEER REVIEW 11 of 25
Land 2025, 14, 700 11 of 25

Figure 7. (a) Vertical Velocity InSAR Map from SBAS method using LiCSBAS (2017–2024). Red
dashed
Figure lineVertical
7. (a) denotesVelocity
the eightInSAR
areas selected
Map fromfor zonal
SBAS analysis
method based
using on significant
LiCSBAS subsidence
(2017–2024). rates
Red
(b) Velocity Map in Ascending Direction (2017–2024). (c) Velocity Map in Descending Direction
dashed line denotes the eight areas selected for zonal analysis based on significant subsidence rates
(2017–2024). These figures were made using Arc GIS Pro 3.4.0.
(b) Velocity Map in Ascending Direction (2017–2024). (c) Velocity Map in Descending Direction
(2017–2024).
The These figures
velocity mapswere madethat
reveal using Arc GIS Pro
subsidence is3.4.0.
unevenly distributed across the region.
Significant subsidence is concentrated in the northwest region, and most of the subsidence
The velocity maps reveal that subsidence is unevenly distributed across the region.
is distributed near the Ring Road. The maximum observed subsidence is ~21 cm/year,
Significant subsidence is concentrated in the northwest region, and most of the subsidence
with a vertical cumulative displacement of ~1.5 m over 91 months. Areas with mixed-use
is distributed near the Ring Road. The maximum observed subsidence is ∼21 cm/year,
zones such as Maharajgunj, Samakhusi, Baluwatar, Lazimpat, Dhumbarahi, Sukedhara,
with a vertical cumulative displacement of ∼1.5 m over 91 months. Areas with mixed-use
Bishalnagar, and Lainchaur, which are densely populated and highly urbanized with the
zones such as Maharajgunj, Samakhusi, Baluwatar, Lazimpat, Dhumbarahi, Sukedhara,
maximum number of hotels, industries, and the tallest business complexes, exhibited
Bishalnagar, and Lainchaur, which are densely populated and highly urbanized with the
the highest subsidence rates. Significant subsidence was also recorded in Koteshwor, a
maximum number of hotels, industries, and the tallest business complexes, exhibited the
mixed-use settlement area. Further, the old settlement areas like Patan and suburb areas
highest subsidence rates. Significant subsidence was also recorded in Koteshwor, a mixed-
like Gwarko, which also has many brick kilns, exhibited considerable subsidence. In the
use settlement area. Further, the old settlement areas like Patan and suburb areas like
outskirts of the valley near the northeastern region, in the area between Bode, Duwakot,
Gwarko, which also has many brick kilns, exhibited considerable subsidence. In the out-
and Bhatgaon, close to the Kathmandu Model College (KMC), a significant subsidence has
skirts of the valley near the northeastern region, in the area between Bode, Duwakot, and
been observed. To study the spatial extent of subsidence, eight Regions of Interest (ROIs)
Bhatgaon, close to the Kathmandu Model College (KMC), a significant subsidence has
were defined around areas showing significant subsidence, as indicated by dashed red
been observed. To study the spatial extent of subsidence, eight Regions of Interest (ROIs)
lines in Figure 7a. The Zonal Statistics tool in ArcGIS Pro was used to compute statistical
were defined around areas showing significant subsidence, as indicated by dashed red
summaries for these ROI zones (Table 3).
lines in Figure 7a. The Zonal Statistics tool in ArcGIS Pro was used to compute statistical
summaries for these ROI zones (Table 3).
Land 2025,
Land 2025, 14,
14, 700
x FOR PEER REVIEW 12 of 25
12 of 25

Table 3. Zonal
Table3. Zonal statistics
statistics summary
summaryof
ofsubsidence
subsidenceof
ofselected
selectedROI.
ROI.

AreaArea MaxMax MinMin Mean

Mean PCT90
PCT 90
S.NS.N Name
Name
(ha) (ha) (cm/yr)
(cm/yr) (cm/yr)
(cm/yr) (cm/yr)
(cm/yr) (cm/yr)
(cm/yr)
1 1Imadol Imadol 1142 1142 −13.4
−13.4 −1.8
−1.8 −9.2
−9.2 −−4.8
4.8
2 2Kaushaltar, Lokanthali
Kaushaltar, Lokanthali 494.7494.7 −12.6
−12.6 −2.3
−2.3 −8.6
−8.6 −−4.9
4.9
3 3Patan, Satdobato
Patan, Satdobato 606.92 −13.1
606.92 −13.1 −2.6
−2.6 −9.2
−9.2 −−5.7
5.7
4 4Kathmandu
Kathmandu
MedicalMedical
College College 473.42 −10.6
473.42 −10.6 −1.3
−1.3 −5.4
−5.4 −−2.4
2.4
5 5Baluwatar,
Baluwatar, Lazimpat
Lazimpat LainchaurLainchaur 1212.73 −21.1
1212.73 −21.1 −1.5
−1.5 −10.7
−10.7 −−3.6
3.6
6 Koteshwor 159.3 −10.8 −2.9 −7.4 −4.8
6 Koteshwor 159.3 −10.8 −2.9 −7.4 −4.8
7 Naya Baneshwor 326.46 −8.4 −3.1 −5.8 −4.1
7 8Naya Baneshwor
Dhungedhara, SanoBharyang 250.17326.46 −11.2−8.4 −2−3.1 −5.8
−5.2 −−4.1
3.2
8 Dhungedhara, SanoBharyang 250.17 −11.2 −2 −5.2 −3.2

Overextraction
Overextraction of of water
water has
has resulted
resulted in
in formation
formation of of bull’s-eye
bull’s-eye patten
patten formation
formation in
in
zone
zone 55 of
of subsidence
subsidence map,
map, which
which also
also exhibits
exhibits the
the highest
highest subsidence
subsidence rate
rate compared
compared to
to
rest
rest of
of the
the study
study area.
area. Figure
Figure 8a
8a shows
shows the
the time
time series
series of
of vertical
vertical displacement
displacement and
and mean
mean
subsidence
subsidence rates
rates at
at borehole
borehole sites
sites along
along the
the profile
profile XY
XY (Figure
(Figure4a).
4a).

Figure 8.
Figure 8. (a) InSAR-derived time
(a) InSAR-derived time series
series and mean vertical
and mean vertical velocity
velocity of
of eight
eight borehole
borehole sites
sites whose
whose
location isisshown
location shownininFigure
Figure
4a 4a across
across profile
profile lineline XY,columnar
XY, (b) (b) columnar section
section ofboreholes
of eight eight boreholes (c)
(c) Recent
Recent reported
reported damagesdamages seen in locations
seen in various various locations of our
of our study study
area sucharea such as landslide,
as landslide, fissures,
fissures, and and
sinkhole
formation [50].
sinkhole formation [50].

Equation
Equation (2)(2)isisused
usedtototransform
transform thethe
ascending
ascending LOSLOS
intointo
vertical displacement
vertical displacement of each
of
borehole site. H26
each borehole site.and
H26H02,and located in the in
H02, located central region,region,
the central recorded maximum
recorded subsidence
maximum sub-
rates of −
sidence 13 cm/year
rates and −17
of −13 cm/year andcm/year, respectively.
−17 cm/year, In contrast,
respectively. boreholes
In contrast, PH01,
boreholes B05,
PH01,
and B018, situated in the northern and southeastern regions where the
B05, and B018, situated in the northern and southeastern regions where the aquitard is aquitard is minimal
or absent,orexhibited
minimal negligible
absent, exhibited subsidence.
negligible This highlights
subsidence. that subsidence
This highlights in the valley
that subsidence in the
is primarily driven by sediment compaction, with the aquitard
valley is primarily driven by sediment compaction, with the aquitard playing playing a critical role.
a critical
The
role.aquitard, composed
The aquitard, composed of clay, actsacts
of clay, as aasbarrier
a barrier separating
separatingthetheshallow
shallowandand confined
confined
deep aquifers and plays a significant role in driving subsidence. During
deep aquifers and plays a significant role in driving subsidence. During monsoon season, monsoon season,
the
the shallow
shallow aquifer
aquifer isis recharged
recharged by by rainwater
rainwater and and nearby
nearby rivers,
rivers, which
which experience
experience peak peak
water
water flow, resulting in an elastic rebound of land surface. This is reflected by
flow, resulting in an elastic rebound of land surface. This is reflected by oscillatory
oscillatory
patterns
patterns inin the
thetime
timeseries
seriesfor
for regions
regions near
near thethe rivers.
rivers. Figures
Figures 4a 7a
4a and and 7a elucidate
elucidate that
that areas
areas experiencing maximum subsidence have large number of
experiencing maximum subsidence have large number of deep tubewells that extractdeep tubewells that extract
Land 2025, 14, x FOR PEER REVIEW 13 of 25
Land 2025, 14, 700 13 of 25

large
largeamount
amountofofgroundwater
groundwaterfromfromdeep
deepconfined
confinedaquifers.
aquifers.The
Thetime
timeseries
seriesofofthis
thisregion
region
has
has liner trend dominant with minimal seasonal fluctuations, indicating that theobserved
liner trend dominant with minimal seasonal fluctuations, indicating that the observed
extreme
extremesubsidence
subsidence is mainly associated
is mainly associatedwith
withdeep
deepgroundwater
groundwater pumping,
pumping, where
where sea-
seasonal
sonal rainfall has limited influence. The impermeable black clay layer prevents
rainfall has limited influence. The impermeable black clay layer prevents rainwater from rainwater
from penetrating
penetrating and recharging
and recharging the aquifer,
the deep deep aquifer, resulting
resulting in a persistent
in a persistent declinedecline
in water inlevels.
wa-
ter levels. This decrease in water levels is accompanied by soil compaction, further
This decrease in water levels is accompanied by soil compaction, further exacerbating the exac-
erbating
region’sthe region’s vulnerability
vulnerability to subsidence.
to subsidence.
𝑣
𝑣 vlos (2)
vu = 𝑐𝑜𝑠𝜃 (2)
cosθ
where 𝑣 is the vertical displacement, 𝑣 is the LOS displacement, and 𝜃 is the inci-
where
dent vu is the vertical displacement, vlos is the LOS displacement, and θ is the incident angle.
angle.
TheArcGIS
The ArcGISSurface
SurfaceProfile
Profiletool
toolwas
wasemployed
employedtotovisualize
visualizethe
thesubsidence
subsidencechange
change
across profile AB and MN (Figure 9). The primary objective was to analyze
across profile AB and MN (Figure 9). The primary objective was to analyze subsidence subsidence
ratesat
rates at each
each borehole
borehole location
locationand
andassess
assessthe correlation
the between
correlation sediment
between thickness,
sediment land
thickness,
use, lithology and vertical subsidence.
land use, lithology and vertical subsidence.

Figure
Figure9.9.(a)(a)
Sentinel 2 LULC
Sentinel 2 LULCMapMapfromfrom
ESRIESRI
[51]. [51].
(b) Geology of Kathmandu
(b) Geology Valley,Valley,
of Kathmandu modified from
modified
[52,53]. (c) Vertical
from [52,53]. MeanMean
(c) Vertical Velocity alongalong
Velocity profile line line
profile MN.MN.(d) (d)
Vertical Mean
Vertical MeanVelocity across
Velocity acrossprofile
profile
AB. The labeled borehole site’s location and columnar section are shown
AB. The labeled borehole site’s location and columnar section are shown in Figure 6.in Figure 6.

Figure9c,d
Figure 9c,d illustrate
illustrate thethe undulating
undulating pattern
pattern of subsidence
of subsidence alongalong both
both the the north-
north-south
south
(NS) and(NS) and east–west
east–west (EW)profiles.
(EW) spatial spatial profiles.
A comparison A comparison of the subsidence
of the subsidence profile
profile with the
with the land classification map (Figure 9a) and geological map (Figure
land classification map (Figure 9a) and geological map (Figure 9b) reveals that the highest9b) reveals that
the highestoccurs
subsidence subsidence
within occurs within
the Kalimati the Kalimati
Formation Formation
and Gokarna and Gokarna
Formation, Formation,
predominantly
predominantly in regions with a high density of urban development.
in regions with a high density of urban development. Both formations consist of highly Both formations
consist of highly
unconsolidated unconsolidated
sediments, sediments,
including including
sand, silt, sand,
clay, and silt, clay,
gravel [28].and
Thegravel [28].For-
Kalimati The
Kalimati Formation is primarily composed of black clay, silt, and fine sand,
mation is primarily composed of black clay, silt, and fine sand, indicative of a low-energy indicative
of a low-energy
depositional depositional
environment environment
associated with the associated with
prehistoric the prehistoric
Kathmandu Kathmandu
Paleo-Lake. The
Paleo-Lake. The substantial clay content of low permeability in this region
substantial clay content of low permeability in this region impedes the penetration of rain- impedes the
penetration
water in deepofconfined
rainwater in deepthereby
aquifers, confined aquifers, thereby
intensifying intensifying
subsidence. subsidence.
The expansion The
of built-
expansion of built-up areas in such region that absorbs minimal rainwater
up areas in such region that absorbs minimal rainwater increases surface imperviousness, increases surface
Land 2025, 14, x FOR PEER REVIEW 14 of 25
Land 2025, 14, 700 14 of 25

leading to heightened
imperviousness, waterlogging
leading to heightened increased runoff,increased
waterlogging and elevated flood
runoff, andrisk. The high
elevated flood
pumping
risk. The of
highgroundwater
pumping ofingroundwater
the GokarnainFormation,
the Gokarna consisting mainly
Formation, of fine-grained
consisting mainly of
sand and silt layers,
fine-grained sand andreduces pore water
silt layers, pressure,
reduces which
pore water exacerbates
pressure, theexacerbates
which risk of subsidence
the risk
in
ofthese areas. A
subsidence inbivariate analysis
these areas. of the vertical
A bivariate analysisvelocity (Figurevelocity
of the vertical 9c,d) and columnar
(Figure 9c,d)sec-
and
tion (Figure 6b,c) reveals that regions with the greatest thickness of sediment
columnar section (Figure 6b,c) reveals that regions with the greatest thickness of sediment deposits,
predominantly black clay,black
deposits, predominantly experience the highest
clay, experience thesubsidence. The borehole
highest subsidence. sites HS and
The borehole sites
NT along the AB profile across NS, characterized by a high clay content, exhibit
HS and NT along the AB profile across NS, characterized by a high clay content, exhibit the the high-
est subsidence
highest rates.rates.
subsidence Similarly, the borehole
Similarly, sitessites
the borehole NV,NV, NB,
NB,and
andSMSMalong
alongthe
theMN
MNprofile
profile
across EW, which are also predominantly composed of black clay, experience
across EW, which are also predominantly composed of black clay, experience significant significant
subsidence
subsidencerates.
rates.

3.2.
3.2.Cross-Validation
Cross-Validation
The
TheInSAR
InSARdeformation
deformationraterateisismeasured
measuredrelative
relativetotoaalocal
localreference
referencepoint,
point,requiring
requiring
quantitative
quantitativeassessment
assessmentfor foraccurate
accurateinterpretation.
interpretation.Our
Ourstudy
studyareaareahas
hasonly
onlyone
oneavailable
available
GPS
GPSstation
station(NAST),
(NAST), which
which includes a 26-month
26-monthdata datagap
gap(11(11December
December2019–27
2019–27November
Novem-
ber 2021),
2021), potentially
potentially impacting
impacting thereliability
the reliabilityofofsubsidence
subsidencerate rateanalysis.
analysis. Given
Given thethestudy
study
region’s
region’ssize
size (~310
(~310 sq. km) andand uneven
unevendeformation
deformationrates,
rates,a single
a single GPS
GPS station
station is insuffi-
is insufficient
for validating
cient InSAR-derived
for validating InSAR-derived surface motion.
surface To ensure
motion. To ensure accuracy, we we
accuracy, used PS-InSAR,
used PS-In-
differential
SAR, analysis
differential of ascending
analysis of ascending and and
descending deformation
descending deformation rates, andand
rates, GPS GPSdata for
data
thethe
for cross-validation
cross-validation of of
LOSLOSresults.
results.This multi-validation
This multi-validation approach
approach enhances
enhances reliability
reliabilityat
atGPS
GPS station
stationlocations,
locations,high-density
high-density PSPSareas, andand
areas, regions
regionscovered by both
covered ascending
by both and
ascending
descending
and descendingorbits.
orbits.

3.2.1.GNSS
3.2.1. GNSSand andHistogram
HistogramDifference
DifferenceTechnique
Technique
TheGPS
The GPSdata
datafor
forthe
theNAST
NAST station
station were
were obtained
obtained from
from thethe Nevada
Nevada Geodetic
Geodetic Lab-
Labor-
oratory (NGL) (http://geodesy.unr.edu/, accessed on 5 October 2024)
atory (NGL) (http://geodesy.unr.edu/, accessed on 5 October 2024) based on IGS14 as thebased on IGS14
as the reference
reference frame, covering
frame, covering the same theperiod
same period
as the as the InSAR
InSAR time series
time series [54]. While
[54]. While GPS
GPS measures deformation in NS, EW, and UD directions, InSAR records
measures deformation in NS, EW, and UD directions, InSAR records displacement in the displacement
in the satellite’s
satellite’s LOS direction
LOS direction (vl ), requiring
𝑣 , requiring conversion
conversion to vertical
to vertical displacement
displacement using
using Equa-
Equation (2) [55]. As shown in Figure 10a, both ascending and descending
tion (2) [55]. As shown in Figure 10a, both ascending and descending InSAR displacementInSAR displace-
ment trends
trends indicateindicate an increasing
an increasing subsidence
subsidence rate, aligning
rate, aligning with GPSwithobservations.
GPS observations.

Figure10.
Figure (a)Time-series
10.(a) Time-seriescomparisons
comparisonsofofGPS
GPSand
andascending
ascendingand
anddescending
descendingcumulative
cumulativedisplace-
displace-
mentwith
ment withprecipitation
precipitationfrom
from2017
2017to
to2024.
2024.(b)
(b)GPS
GPSstation
station(NAST)
(NAST)used
usedin
inthis
thisstudy
study[56].
[56].(c)
(c)Mean
Mean
velocity difference of ascending and descending orbit histogram.
velocity difference of ascending and descending orbit histogram.
 Land 2025, 14, x FOR PEER REVIEW 15 of 25
Land 2025, 14, 700 15 of 25

Atthe
At theGPSGPS site
site (85.328,
(85.328, 27.657),
27.657), the subsidence
the subsidence rate estimated
rate estimated by InSAR byisInSAR
−112 ±is29−112 ± 29
mm/yr,
mm/yr, compared
compared to −125.92 to −125.92
mm/yr mm/yr
from GPS.fromThis
GPS. This discrepancy
discrepancy may from
may result resulta from a 26-
26-month
month
GPS GPS
data gapdata
(11gap (11 December
December 2019–272019–27 November
November 2021),2021), differences
differences in datain data acquisi-
acquisition
tion frequency
frequency (InSAR:(InSAR: 12 GPS:
12 days, days,daily),
GPS: spatial
daily), resolution,
spatial resolution, satellite
satellite look look
angle, andangle, and
seasonal
seasonal variations.
variations. A differential
A differential analysis ofanalysis of ascending
ascending and descending
and descending deformationdeformation rates
rates provides
provides
an an alternative
alternative consistency consistency check, especially
check, especially when nowhen GNSSno GNSSor
station station
PS-InSARor PS-InSAR
data is
data is available
available [57,58].[57,58].
FigureFigure 10c presents
10c presents a histogram
a histogram (bin =size
(bin size 80)=and
80) and a normal
a normal curve
curve fit
fit for
for deformation
deformation raterate differences,
differences, with
with most
most values
values ranging
ranging between
between +40and
+40 and−−40 mm/yr.
40 mm/yr.
Themean
The mean(2.68
(2.68mm/yr)
mm/yr) and median
median (2.77
(2.77 mm/yr)
mm/yr) indicate
indicate minimal difference, confirming
difference, confirming
good consistency
good consistency between
between results.
results. AA skewness
skewness value
value ofof 0.09
0.09 suggests
suggests aa nearly
nearly symmetric
symmetric
distribution,while
distribution, whilea akurtosis
kurtosis of 2.26
of 2.26 (<3)(<3) reflects
reflects a platykurtic
a platykurtic distribution,
distribution, indicating
indicating fewer
fewer extreme
extreme values.values. Minor discrepancies
Minor discrepancies may bemay due be due to elevation
to elevation differences,
differences, incidence incidence
angles,
angles, heading
heading directions,
directions, or temporal
or temporal variations.
variations.

3.2.2.
3.2.2. PS-InSAR
PS-InSAR Method
Method
The
The PS-InSAR
PS-InSAR technique
technique [59][59] utilizes
utilizes permanent
permanent scatterers
scatterers (PS),
(PS), primarily
primarily man-made
man-made
structures,
structures, to analyze stable reflection signals for time-series generation using the
to analyze stable reflection signals for time-series generation using the SNAP-
SNAP-
StaMPS
StaMPSworkflow.
workflow.This technique
This technique is highly effective
is highly for monitoring
effective urbanurban
for monitoring land deformation,
land defor-
achieving millimeter-level
mation, achieving accuracyaccuracy
millimeter-level with a minimum of 20 images
with a minimum [60]. For
of 20 images thisFor
[60]. study,
this
31 Sentinel-1A
study, SLC images
31 Sentinel-1A SLC in ascending
images orbit (IWorbit
in ascending mode, (IW VVmode,
polarization) were acquired
VV polarization) were
from September
acquired 2023 to November
from September 2024 via2024
2023 to November ASF.via The
ASF.image from 23
The image fromMarch 2024 2024
23 March was
selected as theasmaster
was selected image,image,
the master with thewithremaining imagesimages
the remaining co-registered with master
co-registered image
with master
to generate interferograms. The details of the datasets are given in Table
image to generate interferograms. The details of the datasets are given in Table A1. The A1. The SNAP
toolbox was used
SNAP toolbox wasfor interferogram
used generation,
for interferogram followed
generation, by snap2stamps
followed by snap2stamps[61] to[61]
auto-to
mate processing for StaMPS [62]. The final LOS displacement was
automate processing for StaMPS [62]. The final LOS displacement was derived from derived from Matlab
(StaMPS) using seven
Matlab (StaMPS) processing
using steps. As shown
seven processing steps. in
AsFigure
shown 11a,
in the PS-InSAR
Figure 11a, theanalysis
PS-InSARalso
identifies the same regions, outlined by blue-bordered rectangle, experiencing
analysis also identifies the same regions, outlined by blue-bordered rectangle, experienc- significant
subsidence consistent
ing significant with SBAS
subsidence analysis
consistent (Figure
with SBAS7a). The findings
analysis (Figurefrom
7a). these two methods
The findings from
demonstrate a strong agreement.
these two methods demonstrate a strong agreement.

Figure 11.
Figure 11. (a)
(a) PS-InSAR
PS-InSAR deformation
deformation mapmap from
from ascending
ascending track.
track. The
The three
three large
large rectangles
rectangles outlined
outlined
in
in red shows region with highest subsidence. The small rectangle with red outline
shows region with highest subsidence. The small rectangle with red outline is common over-is common
overlapping
lapping areaarea between
between SBASSBAS
andand
PSIPSI pixel
pixel usedfor
used forcorrelation
correlationanalysis.
analysis. (b)
(b) Scatter
Scatter plot
plot between
between
LiCSBAS results and PS-InSAR velocity of region inside label A. (c) Scatter plot between
LiCSBAS results and PS-InSAR velocity of region inside label A. (c) Scatter plot between LiCSBAS LiCSBAS
and
and PS-InSAR
PS-InSAR velocity
velocityof
ofregion
regioninside
insidebox
boxlabel
labelB.
B.
Land 2025, 14, x FOR PEER REVIEW 16 of 25

Land 2025, 14, 700 16 of 25

However, subtle differences in subsidence rates may arise due to variations in pro-
cessing methods, temporal gaps, foreshortening, layover, shadow effects, and different
However, subtle differences in subsidence rates may arise due to variations in process-
orbital direction and incidence angle. The two regions shown by red rectangle in Figure
ing methods, temporal gaps, foreshortening, layover, shadow effects, and different orbital
11a labeled A and B were selected to compare results for an accuracy assessment of SBAS
direction and incidence angle. The two regions shown by red rectangle in Figure 11a labeled
and InSAR methods. The two regions were chosen based on their extensive coverage and
A and B were selected to compare results for an accuracy assessment of SBAS and InSAR
high pixel density, ensuring greater precision in the comparison. Region A, consisting of
methods. The two regions were chosen based on their extensive coverage and high pixel
1516-pixel points, demonstrated a correlation coefficient of 0.899 and a root mean square
density, ensuring greater precision in the comparison. Region A, consisting of 1516-pixel
error (RMSE) of 6.477 mm/year, indicating strong agreement. Similarly, Region B, which
points, demonstrated a correlation coefficient of 0.899 and a root mean square error (RMSE)
contains 671 PS points, exhibited a correlation of 0.822 with an RMSE of 5.518 mm/year,
of 6.477 mm/year, indicating strong agreement. Similarly, Region B, which contains 671 PS
reflecting a slightly
points, exhibited lower butofstill
a correlation substantial
0.822 with an RMSEconsistency
of 5.518inmm/year,
subsidence measurements
reflecting a slightly
(Figure 11b,c). These results confirm a strong correlation between PS-InSAR
lower but still substantial consistency in subsidence measurements (Figure 11b,c). and SBAS,
These
validating the subsidence patterns. Since this study primarily focuses on SBAS-InSAR
results confirm a strong correlation between PS-InSAR and SBAS, validating the subsidence
subsidence analysis,
patterns. Since PS-InSAR
this study wasfocuses
primarily only analyzed for 31 ascending
on SBAS-InSAR subsidenceorbit images
analysis, (2023–
PS-InSAR
2024) to check
was only the for
analyzed accuracy of SBAS
31 ascending results.
orbit images Further improvements
(2023–2024) in accuracy
to check the correlation and
of SBAS
RMSE may be achieved by increasing PS-InSAR temporal resolution to
results. Further improvements in correlation and RMSE may be achieved by increasingmatch the SBAS
time range.temporal resolution to match the SBAS time range.
PS-InSAR

3.3.
3.3.Temporal
Temporaland
andSpatial
Spatial Changes
Changes of
of Emerging
Emerging HotHot Spots
Spots (2017–2024)
(2017–2024)
The
The EHA tool was applied to examine the spatiotemporalpatterns
EHA tool was applied to examine the spatiotemporal patternsof
ofInSAR-derived
InSAR-derived
subsidence
subsidencetrends.
trends.This analysis
This integrates
analysis two
integrates keykey
two statistical methods:
statistical the Getis–Ord
methods: Gi*
the Getis–Ord
statistic and the Mann–Kendall trend test. Figure 12a presents the EHSA
Gi* statistic and the Mann–Kendall trend test. Figure 12a presents the EHSA results, results, high-
lighting subsidence
highlighting patterns
subsidence from from
patterns 2017 to 2024.
2017 to 2024.

Figure
Figure12.
12.(a)
(a)Emerging
Emerginghothotspot
spotanalysis
analysisof
ofInSAR
InSARvertical
verticaldeformation
deformationrates
rates(this
(thiswas
wasprepared
prepared
usingArcGIS
using ArcGISPro).
Pro).(b)
(b)Expansion
Expansionofofurban
urbanareas
areasacross
acrossvarious
variousareas
areas of
of Kathmandu
Kathmandu Valley
Valley [63,64].

Theintensifying
The intensifyinghothotspots
spotsare
arethe
themost
most prevalent
prevalent pattern,
pattern,accounting
accountingforfor56.17%
56.17%ofof
alldetected
all detectedhot hotspots.
spots.The
Theintensifying
intensifyinghot hotspot
spotareas
areasare
areprimarily
primarilyconcentrated
concentratedin inthe
the
Kathmandu and Lalitpur districts, where InSAR analysis indicates subsidence
Kathmandu and Lalitpur districts, where InSAR analysis indicates subsidence rates rang- rates ranging
from
ing 0.50.5
from cm/year
cm/yeartoto
1515cm/year.
cm/year.Persistent
Persistent hot
hot spots,
spots, covering 22.84% of
covering 22.84% ofthe
thevalley,
valley,are
are
aligned along the Ring Road and concentrated in high-population-density
aligned along the Ring Road and concentrated in high-population-density areas, includ-areas, including
Balaju
ing (west),
Balaju Madhyapur
(west), Madhyapur Thimi
Thimi (east),
(east),and
andLalitpur
Lalitpur(south).
(south). The clusters of
The clusters of persistent
persistent
hot spots in Gongabu and Madhyapur Thimi are interconnected along the Ring RoadRoad
hot spots in Gongabu and Madhyapur Thimi are interconnected along the Ring cor-
corridor,
ridor, suggesting
suggesting a relationship
a relationship betweenbetween subsidence
subsidence anddevelopment.
and urban urban development. The
The hot spot
hot spot clusters bifurcate near Madhyapur Thimi (east) and Gongabu (northwest),
clusters bifurcate near Madhyapur Thimi (east) and Gongabu (northwest), with branches with
branches extending
extending north andnorth
southandfromsouth from Madhyapur
Madhyapur Thimi,
Thimi, while while
those those
from from Gongabu
Gongabu extend
extend northward toward Pharping and westward toward Balaju. In southern Lalitpur,
 Land 2025, 14, x FOR PEER REVIEW 17 of 25

Land 2025, 14, 700 17 of 25

northward toward Pharping and westward toward Balaju. In southern Lalitpur, a cluster
aofcluster
persistent hot spotshot
of persistent suggest
spotsongoing
suggestsubsidence, likely driven
ongoing subsidence, by intensive
likely driven bygroundwa-
intensive
ter extractionextraction
groundwater in areas with a high
in areas withconcentration of hotelsofand
a high concentration a dense
hotels and apopulation. Mean-
dense population.
while, diminishing
Meanwhile, diminishing hot spots occupy
hot spots 15.87%
occupy of the
15.87% study
of the area,
study primarily
area, primarily located
locatedeast
eastof
ofthe
theTribhuvan
TribhuvanInternational
InternationalAirport
Airport(TIA)
(TIA) upup to Duwakot and northwest
northwestfrom fromBaneshwor
Baneshwor
totoTarkeshwor,
Tarkeshwor,indicating
indicatingthatthatsubsidence
subsidencehas hasslowed
slowedover overtime.
time.The
Thediminishing
diminishinghot hotspot
spot
cluster
clusterininthe
thenorthwest
northwest coincides
coincides with thethe
with highest
highestsubsidence
subsidenceratesrates
observed in the
observed inInSAR-
the In-
derived velocity
SAR-derived map, indicating
velocity a temporal
map, indicating shift inshift
a temporal subsidence intensity.
in subsidence This deceleration
intensity. This decel-
may be attributed
eration to the increased
may be attributed water supply
to the increased waterfrom
supply thefrom
Melamchi Water Supply
the Melamchi WaterProject
Supply
(MWSP), completed
Project (MWSP), in 2021, as
completed inwell
2021,asasfrom
wellthe
as Bagmati
from theRiver.BagmatiSince its implementation,
River. Since its imple-
the proportion
mentation, theofproportion
groundwater of extraction
groundwater in total water production
extraction has decreased
in total water production from
has43%
de-
increased
2021 tofrom
nearly
43% 8% inin 2024
2021 to [16]. However,
nearly 8% in 2024 further analysis isfurther
[16]. However, necessary to gain
analysis a deeper
is necessary
understanding
to gain a deeper of understanding
the underlying of mechanisms
the underlyingdriving this trend.driving this trend.
mechanisms

3.4.
3.4.Groundwater
GroundwaterStatus
Status
The
The rapid urbanizationand
rapid urbanization andpopulation
populationgrowth
growthhavehaveled
ledtotoananexponential
exponentialincrease
increase
ininwater
water demand within KV, rising from an estimated 214.6 in 2001 to 506 MLD in2024.
demand within KV, rising from an estimated 214.6 in 2001 to 506 MLD in 2024.
During
Duringthe thesame
sameperiod,
period,the
thepopulation
populationgrew grewfrom
from1.5
1.5million
milliontotoapproximately
approximately33million.
million.
However,
However,a alimited
limitedgovernment
government initiative was
initiative wasundertaken
undertaken to enhance
to enhance surface water
surface watersupply,
sup-
placing
ply, placing significant stress on groundwater resources. To address this growing crisis,
significant stress on groundwater resources. To address this growing water water
the MWSP
crisis, aimed aimed
the MWSP at delivering 170 MLD
at delivering of water
170 MLD through
of water a 26.4
through km diversion
a 26.4 km diversiontunnel,
tun-
which was initiated in 1998 and completed in 2021. Following the completion
nel, which was initiated in 1998 and completed in 2021. Following the completion of of MWSP
inMWSP
2021, water
in 2021,production has improved
water production significantly.
has improved However,However,
significantly. a deficit of 40% remains
a deficit of 40%
relative to total demand. Figure 13c shows a steady increase in water
remains relative to total demand. Figure 13c shows a steady increase in water supply supply each year,
each
but without timely interventions, climate change and prolonged droughts
year, but without timely interventions, climate change and prolonged droughts are likely are likely to
exacerbate water shortages in the future, increasing the deficit percentage
to exacerbate water shortages in the future, increasing the deficit percentage further. further.

Figure 13. (a) Water level data of five well locations from 2012 to 2014 (left) and InSAR time series
Figure 13. (a) Water level data of five well locations from 2012 to 2014 (left) and InSAR time series
from 2017 to 2024 (right). (b) IDW interpolation of groundwater level rate of change of well from
from 2017 to 2024 (right). (b) IDW interpolation of groundwater level rate of change of well from
S4W from 2018 to 2022 across 46 wells; black boundary represents the zonal ROI based on Figure 7a;
S4W
(c) from
KUKL 2018supply
water to 2022statistics
across 46 wells; black boundary represents the zonal ROI based on Figure
(2029–2024).
7a; (c) KUKL water supply statistics (2029–2024).
Land 2025, 14, 700 18 of 25

The line of sight (LOS) deformation from ascending orbit was transformed into vertical
displacement for five selected well monitored by GWRDP between 2014 and 2020 using
Equation (2) (Figure 13a). These sites were selected based on data frequency and time span.
Among the observed wells, Well E exhibited the highest subsidence, followed by Wells C
and D, all located in highly urbanized and industrialized areas. In contrast, Well A showed
minimal subsidence, while Well B exhibited a slight uplift of 0.15 cm/year. These wells are
situated on the outskirts of the valley, where shallow bedrock and minimal sediment thick-
ness contribute to greater resistance against subsidence. Since these wells primarily extract
water from shallow aquifers, which are naturally replenished by rainfall and nearby rivers,
their water level data may not sufficiently correlate with InSAR-derived subsidence trends,
particularly if subsidence is driven by deeper aquifer depletion. This limitation results in
a weak correlation between observed water level fluctuations and ground deformation.
Currently, data from deep boreholes is scarce, restricting comprehensive analysis. Future
studies incorporating water level records from deeper tube wells will be essential for im-
proving the understanding of groundwater extraction’s role in subsidence and enhancing
the accuracy of InSAR-based assessments.
The continuous subsidence trend with a steep negative slope, without significant
seasonal oscillation, suggests that deep aquifer extraction, which lacks direct rainfall
recharge, is the primary driver of land deformation in the study area. Despite limited
historical groundwater data, recent records from S4W-Nepal (S4W-Nepal) were obtained
from (https://smartphones4water.org, accessed on 12 September 2024). S4W-Nepal, a
non-profit research organization founded in 2017, leverages mobile technology to collect
hydro-meteorological data in Nepal by mobilizing young researchers to enhance water
management. To ensure consistency, a filter was applied to retain only wells with at
least 70% data coverage, identifying a period where the maximum number of wells had
continuous observations. The results showed that 44 wells met this criterion between 2018
and 2022, the highest number available. Inverse Distance Weighting (IDW) interpolation
was applied based on the linear trend of groundwater levels, allowing for the visualization
of spatial patterns across the region, as illustrated in Figure 13b. The groundwater level
is declining at a rate of up to 0.7 m/year in areas with high density of deep boreholes
(Figure 4a). These areas are also experiencing accelerated subsidence, as indicated by InSAR
results, highlighting that groundwater extraction from confined deep aquifer is primary
driver of subsidence in KV.

4. Discussion
Vigilant monitoring of land subsidence constitutes a critical endeavor to forestall
future geohazards and avert irreversible consequences. In high-density urban areas with
extensive built environments, subsidence control is both technically challenging and eco-
nomically prohibitive if early warning signals are ignored, leading to substantial financial
and infrastructural losses. The situation is particularly alarming in seismically vulnera-
ble cities like Kathmandu, where active fault lines lie in proximity. This study employs
LiCSBAS to process InSAR data for a comprehensive analysis of subsidence in the study
region. For the first time, we integrate LiCSBAS results into a Space–Time Cube framework
and utilize emerging hot spot analysis (EHA) for spatiotemporal assessment, identifying
significant subsidence hot spots and cold spots. The study reveals that the subsidence trend
is decreasing in area, which shows a maximum subsidence rate, probably due to a reduc-
tion in groundwater extraction by KUKL after the MWSP became operational. Also, the
northwestern region with high number of industries, hotels, and commercial buildings and
maximum concentration of deep tubewell over Gokarna Formation experiences subsidence
rates reaching up to ~21 cm per year.
Land 2025, 14, x FOR PEER REVIEW 19 of 25

Land 2025, 14, 700 19 of 25

Despite the severity of subsidence in Kathmandu Valley, Kathmandu has been con-
spicuously
Despiteabsent from global
the severity subsidence
of subsidence rankings. Bagheri-Gavkosh
in Kathmandu Valley, Kathmandu [65],
hasinbeen
theirconspic-
analy-
sis of 290 metropolitan cities based on the published literature, categorized
uously absent from global subsidence rankings. Bagheri-Gavkosh [65], in their analysis of subsidence
intometropolitan
290 five groups, with
citiesthe highest
based ratepublished
on the (20.1–39 cm/year)
literature,occurring
categorizedin 3.1% of casesinto
subsidence andfive
the
second highest (11.31–20 cm/year) in 9.8% of cases. Their study highlights
groups, with the highest rate (20.1–39 cm/year) occurring in 3.1% of cases and the second eastern China,
the Marand
highest Plain cm/year)
(11.31–20 in Iran, the in Texas
9.8% of Gulf Coast,
cases. northern
Their Jakarta, and
study highlights the Mexico
eastern China,Citythe
Metropolitan Area as having the most severe subsidence. Kathmandu, however,
Marand Plain in Iran, the Texas Gulf Coast, northern Jakarta, and the Mexico City Metropoli- is absent
from
tan these
Area assessments.
as having the mostSimilarly, a global land
severe subsidence. subsidence
Kathmandu, study by
however, Davydzenka
is absent [12]
from these
synthesizing 221 individual publications and open-source data places
assessments. Similarly, a global land subsidence study by Davydzenka [12] synthesizing Nepal’s maximum
recorded
221 subsidence
individual at merely
publications and 3.8 cm/year, positioning
open-source it near the
data places Nepal’s bottom recorded
maximum of global subsi-
rank-
ings. Another study by Tzampoglou [13] examined subsidence rates
dence at merely 3.8 cm/year, positioning it near the bottom of global rankings. Another across various cities
and found
study values ranging
by Tzampoglou [13]from 0.1 cm/year
examined to a maximum
subsidence of ~35various
rates across cm/year in Mexico
cities City.
and found
The study attributes approximately 40% of observed subsidence
values ranging from 0.1 cm/year to a maximum of ~35 cm/year in Mexico City. The to urbanization, particu-
larly in
study Mexico approximately
attributes City, Beijing, Semarang, and Tehran,
40% of observed with to
subsidence additional contributions
urbanization, particularlyfromin
industrial, agricultural, and anthropogenic impacts. Kathmandu remains
Mexico City, Beijing, Semarang, and Tehran, with additional contributions from industrial, excluded from
these analyses
agricultural, andasanthropogenic
well. This research aimsKathmandu
impacts. to bridge that gap by
remains providing
excluded from robust
theseevidence
analyses
that positions Kathmandu Valley among the world’s most subsidence-prone
as well. This research aims to bridge that gap by providing robust evidence that urban areas.
positions
Figure 14a delineates
Kathmandu Valley amongareas affected by subsidence
the world’s based on rate and
most subsidence-prone urbancumulative
areas. Figuredisplace-
14a
ment.
delineates areas affected by subsidence based on rate and cumulative displacement.

Figure 14.
Figure 14. (a)
(a) The
The first
first and
and second
second row
row images
images illustrate
illustrate regions
regions affected
affected by
by varying
varying subsidence rates
subsidence rates
and total displacement, as detected by InSAR from 2017 to 2024. (b) The bar diagram
displacement, as detected by InSAR from 2017 to 2024. (b) The bar diagram quantifies quantifies the
totaltotal
the areaarea
impacted
impacted by subsidence, categorizing
by subsidence, affected
categorizing regions
affected based
regions on the
based color-coded
on the classifi-
color-coded clas-
sifications in (a).
cations in (a). (c) The
(c) The trendtrend graph
graph depicts
depicts the water
the total total water
demand demand and groundwater
and groundwater extraction
extraction (meas-
(measured in million
ured in million liters liters perMLD)
per day, day, MLD)
from from 2019
2019 to to 2024.
2024.

A histogram (Figure 14b) visualization categorizes total affected areas into five classes,
A histogram (Figure 14b) visualization categorizes total affected areas into five clas-
revealing that ~174.7 square kilometers exhibit subsidence below 5 cm/year, while
ses, revealing that ~174.7 square kilometers exhibit subsidence below 5 cm/year, while
~155.4 square kilometers have undergone cumulative displacement of ~20 cm between
~155.4 square kilometers have undergone cumulative displacement of ~20 cm between
2017 and 2024, signifying a dire scenario. Less than 3 square kilometers, predominantly
2017 and 2024, signifying a dire scenario. Less than 3 square kilometers, predominantly in
in the northeastern region, experience subsidence exceeding 15 cm/year, while 3.3 square
the northeastern region, experience subsidence exceeding 15 cm/year, while 3.3 square
kilometers have suffered vertical displacement exceeding 1.6 m in just eight years. To
Land 2025, 14, 700 20 of 25

contextualize these findings, Bangkok has subsided 1.25 m over 113 years (2–3 cm/year),
New Orleans by 1.13 m over the same period (6 cm/year), Shanghai by 2–3 m since 1920
(5 mm/year), and Tianjin by 3–5 m (3–30 mm/year). Jakarta has recorded a subsidence of
2.5–4 m from 1970 to 2020 (25–30 cm/year), while Tehran has subsided 3–5 m (25 cm/year).
Ho Chi Minh City has recorded a 1–2 m subsidence over 50 years (20–25 mm/year) [13].
While these cities required decades to centuries to accrue over 2 m of subsidence, Kath-
mandu’s pace positions it among the world’s swiftest-subsiding urban centers. Our findings
indicate that deep aquifer depletion is a primary driver of land subsidence in Kathmandu
Valley. This conclusion aligns with the findings of Huang [24], who also identified ground-
water extraction as a significant factor contributing to subsidence in the region. Shallow
aquifers are primarily recharged by rainwater or nearby streams; however, this process con-
tributes minimally to the replenishment of deep aquifers. The expansion of built-up areas
into peripheral zones, encroaching upon open spaces near mountainous regions, is likely to
exacerbate the subsidence. Several studies on the interaction of shallow aquifer and nearby
rivers have confirmed the hydraulic connectivity between the two systems. Malla [66]
analyzed stable water isotopes in the Bagmati and Bishnumati Rivers, revealing significant
mixing with nearby wells. Prajapati [67] found that 12% of wells had higher water levels
than adjacent streams in pre-monsoon conditions, while 69% exhibited higher levels in
post-monsoon periods. Bajracharya [68] quantified stream–aquifer interactions and con-
firmed significant hydrological exchanges near river channels. The time-series analysis
from this study reveals oscillations in areas near shallow boreholes, suggesting elastic land
rebound following seasonal recharge. While research on the recharge mechanisms of deep
aquifers in our study area remains limited, Martinelli’s [69] study in Emilia-Romagna, Italy,
identified land subsidence resulting from the overexploitation of deep aquifers through
isotopic analysis. A significant study on the deep aquifers of Kathmandu Valley (KV)
by Shakya [70] revealed that recharge primarily occurs via lateral flow from the base of
surrounding mountains, augmented by the infiltration of high-frequency rainfall over pro-
longed periods. Therefore, additional research is essential to further elucidate the recharge
dynamics of deep aquifers in KV.
Shrestha [71] in his study about the recharge rate in the aquifer of KV determined that
infiltration rates are more influenced by land use than soil texture, with cultivated land
exhibiting higher infiltration due to root permeability. The rapid expansion of impervious
built-up areas, particularly in former agricultural zones, has significantly reduced natural
recharge zones, exacerbating groundwater depletion and intensifying subsidence. Similar
correlations between groundwater extraction, sediment consolidation, and subsidence have
been documented in Mexico City, Jakarta, and Tianjin [13].
KUKL reports that while water demand in KV has risen by ~17% from 2019 to 2024,
groundwater extraction has declined by 39.4% (Figure 14c). The additional water supply
from MWSP and the Bagmati River has significantly reduced dependency on groundwa-
ter. Thapa [72] projected a land rebound of up to ~5 m if MWSP reaches full operational
capacity, a prediction corroborated by Shrestha [73], who developed the first land subsi-
dence model using SHETRAN and MODFLOW. A similar trend was observed in Houston,
Texas, where a study of 95 GNSS stations between 2005 and 2012 recorded a slight land
rebound following enforced groundwater withdrawal regulations [1]. The MWSP and
Bagmati Water Supply Project have demonstrated a positive impact in curbing subsidence.
Migration trends indicate that major cities in Nepal’s Terai region are also experiencing
rapid urbanization, leading to increasing water scarcity same as KV. Reports have emerged
regarding declining tubewell water levels. Luitel [74] investigated subsidence in Pokhara
(another major metropolitan city) using Sentinel-1 SAR data and found an average rate
of 16 cm/year, indicating that subsidence is not confined to the Terai but also affects hilly
Land 2025, 14, 700 21 of 25

regions, with potential intensification in near future. Therefore, more research is crucial
to assess subsidence in urban regions across the country, in order to mitigate potential
geohazards and promote sustainable water and land management strategies. Our results
are accessible online via a web application (https://ktmsub.alwaysdata.net, accessed on
11 February 2025), developed by the authors to provide readers with interactive access to
our findings.

5. Conclusions
This study used LiCSBAS to process InSAR data from 2017 to 2024, enabling a spa-
tiotemporal analysis of land subsidence in the study area. The InSAR data were integrated
with the EHA tool to assess changes in subsidence patterns following the implementation
of the Melamchi Water Supply Project (MWSP) in 2021. Our findings indicate that signif-
icant subsidence is occurring across the region, with rates reaching up to ~21 cm/year
(2017–2024) and cumulative displacement exceeding ~1.5 m. The most affected areas are
densely populated urban centers with high concentrations of industries, hotels, educa-
tional institutions, and commercial buildings. The results align with previous research,
reaffirming that excessive groundwater extraction from deep aquifers is the primary driver
of subsidence. Similar patterns have been observed in other cities experiencing rapid
groundwater depletion, including Mexico City, Jakarta, and Houston, Texas [13].
The most impacted locations include Bagbazar, Lazimpat, Lainchaur, and Balaju (Kath-
mandu), Madhyapur Thimi and Lokanthali (Bhaktapur), and Gwarko and Patan (Lalitpur),
with subsidence extending into peripheral areas such as Gothatar, Bode, Imadol, Satungal,
and Dhapasi Tokha in the northeastern region. Previous studies predicted that subsidence
rates would decline with the increased availability of surface water from MWSP and Bag-
mati River. Our EHA analysis has shown decline in subsidence intensity in previously
high-risk zones. The EHA analysis of InSAR data shows a decline in subsidence in parts of
the northwestern KV. However, MWSP has been operational for only three years and has
faced temporary disruptions due to flooding. To ascertain this relationship with greater
precision, extended temporal data will be essential in the future, alongside comprehensive
ground-based measurements. This study has demonstrated the efficacy of three valida-
tion techniques that can facilitate accuracy assessments for various methodologies in the
absence of GNSS data.
While current mitigation strategies have shown promising results, water demand
continues to rise, leaving 40% of total demand unmet. Policy interventions must priori-
tize mandatory borehole registration, groundwater monitoring, and restrictions on new
infrastructure in vulnerable zones. Decentralized urban planning, infrastructure upgrades,
rainwater harvesting, and the establishment of artificial recharge zones should be integral
to future strategies aimed at preventing further subsidence. This research provides critical
insights for urban planners and policymakers, offering data-driven recommendations to
mitigate economic and infrastructural losses caused by land subsidence. By guiding long-
term urban resilience strategies, these findings can help ensure sustainable development
and safeguard the stability of rapidly expanding urban centers.

Author Contributions: Conceptualization, S.R. and G.W.; methodology, S.R.; software, S.R.; vali-
dation, S.R.; formal analysis, S.R.; investigation, S.R.; data curation, S.R.; writing—original draft
preparation, S.R.; writing—review and editing, S.R. and G.W.; visualization, S.R.; supervision, G.W.
All authors have read and agreed to the published version of the manuscript.

Funding: This research received no external funding.

Data Availability Statement: Sentinel data: SBAS InSAR time series was formed using unwrapped
interferograms obtained from https://comet.nerc.ac.uk/COMET-LiCS-portal/, accessed on 22 Oc-
Land 2025, 14, 700 22 of 25

tober 2024 using frame 019D_06217_131313 and 085A_06253_131313 for descending and ascend-
ing orbits from 2017 to 2024. PS-InSAR analysis was performed using data obtained from the
NASA’s Alaska Satellite Facility Distributed Active Archive Center (ASF DAAC) from 2023 to 2024
(https://search.asf.alaska.edu/, accessed on 22 May 2023). LULC maps were downloaded from
https://livingatlas.arcgis.com/landcoverexplorer (accessed on 25 December 2024). Metrological
Data: The precipitation data were acquired from the Department of Hydrology and Meteorol-
ogy, Nepal, upon request (https://www.dhm.gov.np/, accessed on 15 December 2024). GNSS
Data: The GNSS data for station “NAST” were obtained from https://geodesy.unr.edu/, accessed
11 December 2024. Water-well data: The water-well data were obtained from the Groundwater
Development Resource Project (https://gw-nepal.com/, accessed 20 December 2024) and S4W
(https://s4w-nepal.smartphones4water.org/, accessed 11 December 2024).

Acknowledgments: LiCSAR contains modified Copernicus Sentinel data (2017–2024) analyzed by

the Centre for the Observation and Modelling of Earthquakes, Volcanoes and Tectonics (COMET).
LiCSAR uses JASMIN, the UK’s collaborative data analysis environment (http://jasmin.ac.uk) (ac-
cessed on 22 October 2024). We are immensely grateful to the editor and anonymous reviewers for
their comments on the manuscript. I would like to thank Jack Thomson for a detailed review of
the manuscript.

Conflicts of Interest: The authors declare no conflicts of interest.

Abbreviations
The following abbreviations are used in this manuscript:

EHA Emergency Hot Spot Analysis

GACOS Generic Atmospheric Correction Online Service for InSAR
MWSP Melamchi Water Supply Project
KV Kathmandu Valley
DHM Department of Hydrology and Meteorology
KUKLLOS Kathmandu Upatyaka Khanepani LimitedLine of Sight

Appendix A

Table A1. Sentinel 1A Dataset Information for PS-InSAR Processing.

Orbit Path Period Frame Azimuth Master Total Incidence Sub Swath
September November 28 March
Asc 85 88 350.92◦ 31 39.46◦ IW2
2023 2024 2024

References
1. Kearns, T.J.; Wang, G.; Bao, Y.; Jiang, J.; Lee, D. Current Land Subsidence and Groundwater Level Changes in the Houston
Metropolitan Area (2005–2012). J. Surv. Eng. 2015, 141, 05015002.
2. Ge, L.; Ng, A.H.M.; Du, Z.; Chen, H.Y.; Li, X. Integrated Space Geodesy for Mapping Land Deformation over Choushui River
Fluvial Plain, Taiwan. Int. J. Remote Sens. 2017, 38, 6319–6345. [CrossRef]
3. Gao, M.; Gong, H.; Chen, B.; Li, X.; Zhou, C.; Shi, M.; Si, Y.; Chen, Z.; Duan, G. Regional Land Subsidence Analysis in Eastern
Beijing Plain by InSAR Time Series and Wavelet Transforms. Remote Sens. 2018, 10, 365. [CrossRef]
4. Yin, J.; Yu, D.; Wilby, R. Modelling the Impact of Land Subsidence on Urban Pluvial Flooding: A Case Study of Downtown
Shanghai, China. Sci. Total Environ. 2016, 544, 744–753. [CrossRef]
5. Wang, G.; Greuter, A.; Petersen, C.M.; Turco, M.J. Houston GNSS Network for Subsidence and Faulting Monitoring: Data
Analysis Methods and Products. J. Surv. Eng. 2022, 148, 04022008.
6. Li, Z.; Chen, Q.; Xue, Y.; Qiu, D.; Chen, H.; Kong, F.; Liu, Q. Numerical Investigation of Processes, Features, and Control of Land
Subsidence Caused by Groundwater Extraction and Coal Mining: A Case Study from Eastern China. Environ. Earth Sci. 2023,
82, 82.
Land 2025, 14, 700 23 of 25

7. Meng, D.; Chen, B.; Gong, H.; Zhang, S.; Ma, R.; Zhou, C.; Lei, K.; Xu, L.; Wang, X. Land Subsidence and Rebound Response to
Groundwater Recovery in the Beijing Plain: A New Hydrological Perspective. J. Hydrol. Reg. Stud. 2025, 57, 102127.
8. Tang, W.; Zhao, X.; Motagh, M.; Bi, G.; Li, J.; Chen, M.; Chen, H.; Liao, M. Land Subsidence and Rebound in the Taiyuan Basin,
Northern China, in the Context of Inter-Basin Water Transfer and Groundwater Management. Remote Sens. Environ. 2022,
269, 112792.
9. Wang, K.; Wang, G.; Bao, Y.; Su, G.; Wang, Y.; Shen, Q.; Zhang, Y.; Wang, H. Preventing Subsidence Reoccurrence in Tianjin: New
Preconsolidation Head and Safe Pumping Buffer. Groundwater 2024, 62, 778–794.
10. Zhang, W.; Zhang, T.; Fu, Z.; Ai, P.; Yao, G.; Qi, J. Emerging Spatio-Temporal Hot Spot Analysis of Beijing Subsidence Trend
Detection Based on PS-InSAR. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci. 2024, 48, 861–866.
11. Lees, M.; Knight, R. Quantification of Record-Breaking Subsidence in California’s San Joaquin Valley. Commun. Earth Environ.
2024, 5, 677. [CrossRef]
12. Davydzenka, T.; Tahmasebi, P.; Shokri, N. Unveiling the Global Extent of Land Subsidence: The Sinking Crisis. Geophys. Res. Lett.
2024, 51, e2023GL104497. [CrossRef]
13. Tzampoglou, P.; Ilia, I.; Karalis, K.; Tsangaratos, P.; Zhao, X.; Chen, W. Selected Worldwide Cases of Land Subsidence Due to
Groundwater Withdrawal. Water 2023, 15, 1094. [CrossRef]
14. Mushrooming Concrete Houses in Kathmandu: Harms Are Visible, but a Solution Is Nowhere in Sight | Aζ South Asia. Available
online: https://architexturez.net/pst/az-cf-219351-1618362841 (accessed on 14 March 2025).
15. Census Nepal. 2021. Available online: https://censusnepal.cbs.gov.np/Home/Index/EN (accessed on 19 August 2024).
16. Kathmandu Upatyaka Khanepani Limited (KUKL)—Water. Available online: https://kathmanduwater.org/ (accessed on
19 August 2024).
17. Dahal, A.; Khanal, R.; Mishra, B.K. Identification of Critical Location for Enhancing Groundwater Recharge in Kathmandu Valley,
Nepal. Groundw. Sustain. Dev. 2019, 9, 100253. [CrossRef]
18. United Nations Development Programme, Bureau for Crisis Prevention and Recovery. Reducing Disaster Risk: A Challenge for
Development: A Global Report. Disaster and Crisis Management; United Nations: New York, NY, USA, 2004; p. 143.
19. Acharya, P.; Sharma, K.; Acharya, I.P. Seismic Liquefaction Risk Assessment of Critical Facilities in Kathmandu Valley, Nepal.
GeoHazards 2021, 2, 153–171. [CrossRef]
20. Bhattarai, R.; Alifu, H.; Maitiniyazi, A.; Kondoh, A. Detection of Land Subsidence in Kathmandu Valley, Nepal, Using DInSAR
Technique. Land 2017, 6, 39. [CrossRef]
21. Krishnan, P.V.S.; Kim, D.; Jung, J. Land Subsidence Monitoring in the Kathmandu Basin, Before and After MW 7.8 Gorkha
Earthquake, Nepal by Sbas-Dinsar Technique. In Proceedings of the IGARSS 2018–2018 IEEE International Geoscience and
Remote Sensing Symposium, Valencia, Spain, 22–27 July 2018; IEEE: New York, NY, USA, 2018; pp. 525–528.
22. Kumahara, Y.; Chamlagain, D.; Upreti, B.N. Geomorphic Features of Active Faults around the Kathmandu Valley, Nepal, and No
Evidence of Surface Rupture Associated with the 2015 Gorkha Earthquake along the Faults. Earth Planets Space 2016, 68, 1–8.
[CrossRef]
23. Bhandari, S. Detecting Surface Displacement in Kathmandu Valley with Persistent Scatterer Interferometry. J. Geoinform. Nepal
2022, 21, 23–29. [CrossRef]
24. Huang, J.; Sinclair, H.; Pokhrel, P.; Watson, C.S. Rapid Subsidence in the Kathmandu Valley Recorded Using Sentinel-1 InSAR. Int.
J. Remote Sens. 2024, 45, 1–20. [CrossRef]
25. Privileges Documentation Esri Developer. Documentation; Esri Developer: Redlands, CA, USA, 2025.
26. Khan, S.D.; Faiz, M.I.; Gadea, O.C.A.; Ahmad, L. Study of Land Subsidence by Radar Interferometry and Hot Spot Analysis
Techniques in the Peshawar Basin, Pakistan. Egypt. J. Remote Sens. Space Sci. 2023, 26, 173–184. [CrossRef]
27. Sakai, H. Stratigraphic Division and Sedimentary Facies of the Kathmandu Basin Group, Central Nepal. J. Nepal Geol. Soc. 2001,
25, 19–32.
28. Stocklin, J.; Bhattarai, K.D. Geology of Kathmandu Area and Central Mahabharat Range, Nepal Himalaya. HMG/UNDP Mineral
Exploration Project; Technical Report; HMG/UNDP Mineral Exploration Project: Kathmandu, Nepal, 1977; 128p.
29. Paudel, M.R.; Sakai, H. Stratigraphy and Depositional Environments of Basin-Fill Sediments in Southern Kathmandu Valley,
Central Nepal. Bull. Dep. Geol. 2008, 11, 61–70.
30. Yoshida, M. Neogene to Quartenery Lacustrine Sediments In The Kathmandu Valley, Nepal. J. Nepal Geol. Soc. Spec. Issue 1984, 4,
73–100.
31. Pandey, V.P.; Chapagain, S.K.; Kazama, F. Evaluation of Groundwater Environment of Kathmandu Valley. Environ. Earth Sci. 2010,
60, 1329–1342. [CrossRef]
32. JICA (Japan International Cooperation Agency). Groundwater Management Project in the Kathmandu Valley, Final; Japan International
Cooperation Agency: Kathmandu, Nepal, 1990.
33. Asahi, K. Thankot Active Fault in the Kathmandu Valley, Nepal Himalaya. J. Nepal Geol. Soc. 2003, 28, 1–8.
Land 2025, 14, 700 24 of 25

34. Cresswell, R.G.; Bauld, J.; Jacobson, G.; Khadka, M.S.; Jha, M.G.; Shrestha, M.P.; Regmi, S. A First Estimate of Ground Water Ages
for the Deep Aquifer of the Kathmandu Basin, Nepal, Using the Radioisotope Chlorine-36. Groundwater 2001, 39, 449–457.
35. Ghimire, S.; Dwivedi, S.K.; Acharya, K.K. Pattern of Ground Deformation in Kathmandu Valley during 2015 Gorkha Earthquake,
Central Nepal. In Proceedings of the AGU Fall Meeting Abstracts; American Geophysical Union: Washington, DC, USA, 2016;
Volume 2016, p. S33D-2859.
36. Morishita, Y.; Lazecky, M.; Wright, T.J.; Weiss, J.R.; Elliott, J.R.; Hooper, A. LiCSBAS: An Open-Source InSAR Time Series Analysis
Package Integrated with the LiCSAR Automated Sentinel-1 InSAR Processor. Remote Sens. 2020, 12, 424.
37. Wright, T.J.; Gonzalez, P.J.; Walters, R.J.; Hatton, E.L.; Spaans, K.; Hooper, A.J. LiCSAR: Tools for Automated Generation of
Sentinel-1 Frame Interferograms. In Proceedings of the AGU Fall Meeting Abstracts; American Geophysical Union: Washington, DC,
USA, 2016; Volume 2016, p. G23A-1037.
38. Manzo, M.; Ricciardi, G.P.; Casu, F.; Ventura, G.; Zeni, G.; Borgström, S.; Berardino, P.; Del Gaudio, C.; Lanari, R. Surface
Deformation Analysis in the Ischia Island (Italy) Based on Spaceborne Radar Interferometry. J. Volcanol. Geotherm. Res. 2006, 151,
399–416. [CrossRef]
39. Pandey, V.P.; Kazama, F. Hydrogeologic Characteristics of Groundwater Aquifers in Kathmandu Valley, Nepal. Environ. Earth Sci.
2011, 62, 1723–1732. [CrossRef]
40. Gurung, J.K.; Ishiga, H.; Khadka, M.S.; Shrestha, N.R. The Geochemical Study of Fluvio-Lacustrine Aquifers in the Kathmandu
Basin (Nepal) and the Implications for the Mobilization of Arsenic. Environ. Geol. 2007, 52, 503–517.
41. Metcalf & Eddy, Inc.; CEMAT Consultants Ltd. Urban Water Supply Reforms in the Kathmandu Valley; Asian Development Bank:
Manila, Philippines, 2000.
42. Shakya, B.M.; Nakamura, T.; Kamei, T.; Shrestha, S.D.; Nishida, K. Seasonal Groundwater Quality Status and Nitrogen Contami-
nation in the Shallow Aquifer System of the Kathmandu Valley, Nepal. Water 2019, 11, 2184. [CrossRef]
43. Devkota, P.; Dhakal, S.; Shrestha, S.; Shrestha, U.B. Land Use Land Cover Changes in the Major Cities of Nepal from 1990 to 2020.
Environ. Sustain. Indic. 2023, 17, 100227. [CrossRef]
44. Galloway, D.L.; Jones, D.R.; Ingebritsen, S.E. Land Subsidence in the United States; Geological Survey (USGS): Reston, VA, USA,
1999; Volume 1182, ISBN 0607926961.
45. GWRDB (Groundwater Resource Development Board). Vulnerability Study in Budhanilkantha and Sorrounding Area Kathmandu
Nepal; Groundwater Resource Development Board: Kathmandu, Nepal, 2009.
46. Khan, S.D.; Gadea, O.C.A.; Tello Alvarado, A.; Tirmizi, O.A. Surface Deformation Analysis of the Houston Area Using Time
Series Interferometry and Emerging Hot Spot Analysis. Remote Sens. 2022, 14, 3831. [CrossRef]
47. Mohamadi, B.; Balz, T.; Younes, A. Towards a PS-InSAR Based Prediction Model for Building Collapse: Spatiotemporal Patterns
of Vertical Surface Motion in Collapsed Building Areas—Case Study of Alexandria, Egypt. Remote Sens. 2020, 12, 3307. [CrossRef]
48. Emerging Hot Spot Analysis (Space Time Pattern Mining)—ArcGIS Pro Documentation. Available online: https://pro.arcgis.
com/en/pro-app/latest/tool-reference/space-time-pattern-mining/emerginghotspots.htm (accessed on 13 December 2024).
49. Space-Time Cubes: Stack Time Like Lego—GIS Geography. Available online: https://gisgeography.com/space-time-cubes
(accessed on 15 December 2024).
50. Ekantipur. Available online: https://ekantipur.com (accessed on 21 December 2024).
51. Karra, K.; Kontgis, C.; Statman-Weil, Z.; Mazzariello, J.C.; Mathis, M.; Brumby, S.P. Global Land Use/Land Cover with Sentinel 2
and Deep Learning. In Proceedings of the 2021 IEEE International Geoscience and Remote Sensing Symposium IGARSS, Brussels,
Belgium, 11–16 July 2021; IEEE: New York, NY, USA, 2021; pp. 4704–4707.
52. Shallow Aquifer Mapping of Kathmandu Valley. Available online: https://www.academia.edu/27155659/Shallow_Aquifer_
Mapping_of_Kathmandu_Valley (accessed on 20 August 2024).
53. Gilder, C.E.L.; Pokhrel, R.M.; Vardanega, P.J.; De Luca, F.; De Risi, R.; Werner, M.J.; Asimaki, D.; Maskey, P.N.; Sextos, A. The
SAFER Geodatabase for the Kathmandu Valley: Geotechnical and Geological Variability. Earthq. Spectra 2020, 36, 1549–1569.
[CrossRef]
54. Blewitt, G.; Hammond, W.; Kreemer, C. Harnessing the GPS Data Explosion for Interdisciplinary Science. Eos 2018,
99, e2020943118.
55. Fuhrmann, T.; Garthwaite, M.C. Resolving Three-Dimensional Surface Motion with InSAR: Constraints from Multi-Geometry
Data Fusion. Remote Sens. 2019, 11, 241. [CrossRef]
56. UNAVCO. Available online: https://www.unavco.org/instrumentation/networks/status/other/overview/NAST (accessed on
22 January 2025).
57. Radman, A.; Akhoondzadeh, M.; Hosseiny, B. Integrating InSAR and Deep-Learning for Modeling and Predicting Subsidence
over the Adjacent Area of Lake Urmia, Iran. GIScience Remote Sens. 2021, 58, 1413–1433.
58. Wang, Y.; Guo, Y.; Hu, S.; Li, Y.; Wang, J.; Liu, X.; Wang, L. Ground Deformation Analysis Using InSAR and Backpropagation
Prediction with Influencing Factors in Erhai Region, China. Sustainability 2019, 11, 2853. [CrossRef]
Land 2025, 14, 700 25 of 25

59. Ferretti, A.; Prati, C.; Rocca, F. Nonlinear Subsidence Rate Estimation Using Permanent Scatterers in Differential SAR Interferome-
try. IEEE Trans. Geosci. Remote Sens. 2000, 38, 2202–2212.
60. Ciampalini, A.; Raspini, F.; Lagomarsino, D.; Catani, F.; Casagli, N. Landslide Susceptibility Map Refinement Using PSInSAR
Data. Remote Sens. Environ. 2016, 184, 302–315. [CrossRef]
61. Delgado Blasco, J.M.; Ziemer, J.; Foumelis, M.; Dubois, C. SNAP2StaMPS v2: Increasing Features and Supported Sensors in the Open
Source SNAP2StaMPS Processing Scheme; Zenodo: Stordalen Mire, Sweden, 2025. [CrossRef]
62. Hooper, A.; Spaans, K.; Bekaert, D.; Cuenca, M.C.; Arıkan, M.; Oyen, A. StaMPS/MTI Manual. Delft Institute of Earth Observation
and Space Systems; Delft University of Technology, Kluyverweg: Delft, The Netherlands, 2010; Volume 1, p. 2629.
63. ChangeKathmandu. Available online: https://x.com/ChangeKathmandu (accessed on 19 December 2024).
64. ICIMOD. Available online: https://geoapps.icimod.org/ktmbeforeafter (accessed on 10 January 2025).
65. Bagheri-Gavkosh, M.; Hosseini, S.M.; Ataie-Ashtiani, B.; Sohani, Y.; Ebrahimian, H.; Morovat, F.; Ashrafi, S. Land Subsidence: A
Global Challenge. Sci. Total Environ. 2021, 778, 146193. [CrossRef] [PubMed]
66. Malla, R.; Shrestha, S.; Chapagain, S.K.; Shakya, M.; Nakamura, T.; Malla, R.; Shrestha, S.; Chapagain, S.K.; Shakya, M.;
Nakamura, T. Physico-Chemical and Oxygen-Hydrogen Isotopic Assessment of Bagmati and Bishnumati Rivers and the Shallow
Groundwater along the River Corridors in Kathmandu Valley, Nepal. J. Water Resour. Prot. 2015, 7, 1435–1448. [CrossRef]
67. Prajapati, R.; Overkamp, N.N.; Moesker, N.; Happee, K.; van Bentem, R.; Danegulu, A.; Manandhar, B.; Devkota, N.; Thapa, A.B.;
Upadhyay, S.; et al. Streams, Sewage, and Shallow Groundwater: Stream-Aquifer Interactions in the Kathmandu Valley, Nepal.
Sustain. Water Resour. Manag. 2021, 7, 1–18.
68. Bajracharya, R.; Nakamura, T.; Shakya, B.M.; Nishida, K.; Das Shrestha, S.; Tamrakar, N.K. Identification of River Water and
Groundwater Interaction at Central Part of the Kathmandu Valley, Nepal Using Stable Isotope Tracers. Int. J. Adv. Sci. Tech. Res.
2018, 8, 29–41. [CrossRef]
69. Martinelli, G.; Chahoud, A.; Dadomo, A.; Fava, A. Isotopic Features of Emilia-Romagna Region (North Italy) Groundwaters:
Environmental and Climatological Implications. J. Hydrol. 2014, 519, 1928–1938. [CrossRef]
70. Shakya, B.M.; Nakamura, T.; Das Shrestha, S.; Nishida, K. Identifying the Deep Groundwater Recharge Processes in an
Intermountain Basin Using the Hydrogeochemical and Water Isotope Characteristics. Hydrol. Res. 2019, 50, 1216–1229. [CrossRef]
71. Shrestha, G.; Shakya, B.M.; Shrestha, M.B.; Khadka, U.R. Water Infiltration Rate in the Kathmandu Valley of Nepal amidst Present
Urbanization and Land-Use Change. H2Open J. 2023, 6, 1–14. [CrossRef]
72. Thapa, B.R.; Ishidaira, H.; Gusyev, M.; Pandey, V.P.; Udmale, P.; Hayashi, M.; Shakya, N.M. Implications of the Melamchi Water
Supply Project for the Kathmandu Valley Groundwater System. Water Policy 2019, 21, 120–137. [CrossRef]
73. Shrestha, P.K.; Shakya, N.M.; Pandey, V.P.; Birkinshaw, S.J.; Shrestha, S. Model-Based Estimation of Land Subsidence in
Kathmandu Valley, Nepal. Geomat. Nat. Hazards Risk 2017, 8, 974–996. [CrossRef]
74. Luitel, S.; Ghimire, S.; Chaulagain, S.; Sah, S.; Khatri Thapa, D.; Lamichhane, S.; Civil, B.E. Assessment of Land Subsidence
Pattern in Pokhara Valley Using Sentinel-1 InSAR Processing. J. Eng. Technol. Plan. 2024, 5, 27–37.

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual
author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to
people or property resulting from any ideas, methods, instructions or products referred to in the content.