The study area is the Charvak Reservoir basin in Tashkent Region, western Tien Shan, northeastern Uzbekistan (41.5–41.9° N, 69.8–70.5° E). The analysis window spans ~1,280 km² centered on the reservoir. Relief is rugged, rising from ~650 m at the shoreline to >3,500 m in the Chatkal, Ugam, and Pskem ranges; steep reservoir-facing slopes contribute to high susceptibility to slope instability (Figure 1). The Charvak basin occupies a tectonically active sector of the western Tien Shan where contrasting lithologies and relief combine with steep valley networks. Uplands are predominantly Paleozoic crystalline and sedimentary rocks, whereas foreland and intramontane zones include Neogene–Quaternary successions; extensive loess mantles are common and are known to lose strength upon wetting, contributing to slope instability. Regionally compiled geological mapping and active-fault datasets document major fault systems transecting the area and provide the structural context for elevated seismicity (Rosi et al., 2023).

The climate is strongly continental and semiarid; most precipitation falls in November–April, with sharp orographic gradients from the Syr Darya valley (~250–300 mm) to the Pskem valley (~800 mm). Hydrology follows two dominant phases—snow accumulation (October–March) and snowmelt-driven runoff (April–September) (Pachore et al., 2024). The Charvak Reservoir is a large multipurpose impoundment that regulates the Chirchik River system; it covers ~40 km² with storage on the order of ~2 km³, and shoreline fluctuations between drawdown and high-water stages locally affect nearshore moisture and bank stability (Juliev et al., 2023).

Recent landslide appraisals identified numerous failures in the lower Pskem valley and emphasized the presence of settlements around the artificial lake, underscoring coupled natural and human pressures (Tacconi Stefanelli et al., 2024). Downstream of Charvak, the Chirchik–Akhangaran system shows distinct geogenic and anthropogenic solute signatures, reflecting increasing water quality pressure along populated corridors (Fornasaro et al., 2024). Recent sediment surveys further indicate potentially toxic elements in recreational riverfront deposits, underscoring exposure risks where tourism intensifies (Fornasaro et al., 2025). Regional mapping also identifies the upper Pskem among Central Asian hotspots for landslide-damming susceptibility, relevant for reservoir inlets and valley corridors (Tacconi Stefanelli et al., 2024). At regional scale, transport corridors linking the Charvak basin appear as expected-loss hotspots, reinforcing the rationale for including infrastructure proximity in the indicator system (Caleca et al., 2024).

Human footprints cluster along reservoir margins and valley floors. In the wider Ugam–Chatkal context adjacent to Charvak, land use/land cover analyses over the last three decades showed measurable built-up expansion and produced scenario-based projections, indicating continued development pressure. Independently, regional risk mapping for Central Asia highlighted road–rail corridors in northeastern Uzbekistan that provide access to the Charvak basin as expected-loss hotspots, reinforcing the need to consider infrastructure proximity in vulnerability appraisals (Alikhanov et al., 2024).

The workflow produces one composite index (EVI) and three domain-specific subindices (TGRI, HERI, and CARI) corresponding to the topographic–geological, hydrological–ecological, and climatic–anthropogenic domains —computed as AHP-weighted sums of orientation-consistently normalized indicators within each domain to support diagnostic interpretation.

Ten indicators were selected to represent (i) topographic–geological predisposition (slope, lithological resistance, distance to faults, and SOC), (ii) hydrological–ecological stress (shoreline proximity and NDWI trend risk), and (iii) climatic–anthropogenic pressure (NDVI trend risk, NDBI trend risk, distance to infrastructure, and LST trend risk). All datasets were harmonized to Universal Transverse Mercator (UTM) Zone 42N and a 30-m analysis grid; SOC (250 m) was resampled to the common grid. Data sources and processing steps are summarized in Table 1.

For the retrospective assessment, warm season (May–August) Landsat Collection 2 Level-2 records (2000–2023) were used to derive trend indicators (NDVI, NDBI, NDWI, and LST), whereas terrain and ancillary layers (slope, lithology, distance to faults, shoreline proximity, and SOC) provided contextual controls. SOC was included as an integrative proxy of soil condition; additional soil properties (e.g., texture or depth) were omitted to limit redundancy with lithology/topography and because harmonized basin-scale layers of comparable reliability are often limited (Salgado et al., 2025). For the near-term component, observed conditions were mapped for 2018 and 2023 using year-specific seasonal composites, and a 2030 scenario was anchored to the 2023 baseline. NDVI and NDBI were extrapolated from 2023 using long-term (2000–2023) warm-season slopes, whereas hydrological pressure was represented by shoreline proximity and a shoreline-shift stress test; the NDWI trend risk layer was retained as a long-term indicator with exclusion of the nearshore band. Following common practice for long time-series analyses, indicators were orientation-normalized so that higher values consistently denote higher vulnerability (Yu et al., 2023).

An MCDA framework was applied using the AHP to weight indicators within three domains: topographic–geological, hydrological–ecological, and climatic–anthropogenic. The topographic–geological domain captures terrain and structural predisposition (slope gradient, lithological resistance, proximity to active faults) and soil susceptibility represented by SOC. The hydrological–ecological domain represents reservoir-related stress through shoreline proximity and moisture-change dynamics via NDWI trend. The climatic–anthropogenic domain characterizes vegetation and land use pressure (NDVI and NDBI trends), accessibility/development pressure (proximity to infrastructure), and thermal stress (LST trend).

Indicator weights were derived from expert-based pairwise comparisons using the 1–9 Saaty importance scale. For each domain, a pairwise comparison matrix was constructed and the priority vector (local weights) was computed; consistency was evaluated using the consistency ratio (CR), with CR ≤ 0.10 considered acceptable. Global weights were then obtained as the product of domain (group) weights and local weights and used in the integrated index. The AHP procedure follows the standard formulation and recent AHP–GIS applications in environmental and hazard assessment ((Aichi, 2024); (Taoukidou et al., 2025)). Local (within-domain) and global (Σ = 1.00) weights are reported in Table 2.

All geoprocessing and map algebra were implemented in Google Earth Engine (GEE). Landsat Collection 2 Level-2 scenes (TM/ETM+/OLI from Landsat 8/9) were scaled to surface reflectance using the USGS-provided factors, masked using QA_PIXEL cloud/shadow-related flags, and aggregated to annual warm season (May–August) median composites (Survey, 2024). All rasters were harmonized to UTM Zone 42N (EPSG:32642) on a 30-m grid and transformed to orientation-consistent vulnerability indicators in the [0,1] range using robust normalization where appropriate (e.g., p2–p98 scaling for trend and distance layers).

Directionality was defined so that higher values consistently indicate higher vulnerability (e.g., lower NDVI and more negative NDWI trends → higher vulnerability; higher NDBI and LST warming trends → higher vulnerability; shorter distances to faults, shoreline, and infrastructure → higher vulnerability via inversion of distance surfaces).

The composite EVI was computed as a weighted linear combination, as shown in Equation 1:

(1)

where 𝐼𝑖 are the normalized indicators and 𝑤𝑖 are the AHP weights.

For the 2000–2023 retrospective appraisal, EVI and thematic indices (TGRI, HERI, and CARI) are presented as continuous 0–1 rasters; where classed maps are shown, quantile-based breaks (q20, q40, q60, and q80) computed from the land-only distribution of the ranked surface are used to form five classes (1–5).

A near-term EVI projection to 2030 was generated using a baseline-anchored scenario design. The method combines an observed change component (2018–2023) with a forward scenario component (2023–2030) to produce consistent maps of recent dynamics and potential near-term shifts in geoecological vulnerability. Annual growing-season (May–August) Landsat composites were produced for each year using Collection 2 Level-2 surface reflectance with QA-based cloud/shadow masking and multisensor harmonization (Landsat 5/7/8/9).

For the baseline years (2018 and 2023), NDVI and NDBI were derived from seasonal composites, transformed from the native −1…1 domain to 0–1, and oriented to vulnerability (decreasing NDVI and increasing NDBI indicate higher vulnerability). For the 2030 projection, NDVI and NDBI were extrapolated from the 2023 baseline using pixel-wise linear slopes estimated from the long-term record (2000–2023; May–August), computed only where at least 12 valid years are available (Equation 2).

where 𝑋 is NDVI or NDBI and 𝛽 is the annual slope estimated from the long-term growing-season record (2000–2023; May–August) using pixel-wise linear regression, computed only where at least 12 valid years are available. Projected values were clipped to the valid −1…1 range and converted to vulnerability-consistent indicators prior to integration with the other layers.

Static and semistatic layers (slope, lithological resistance, distance to active faults, SOC, and baseline shoreline distance), together with long-term NDWI and LST trend layers (2000–2023), retained the same definitions across years unless modified by scenario assumptions. Hydrological–ecological pressure was represented by a renormalized HERI combining shoreline proximity and NDWI trend risk; the NDWI trend component was excluded within a 150 m nearshore band to suppress shoreline-driven artifacts while preserving the physical meaning of shoreline proximity.

A combined-stressor configuration was implemented through three adjustments applied prior to aggregation: (i) infrastructure expansion, represented by buffering roads by 250 m and recomputing proximity-to-infrastructure risk; (ii) high-water shoreline influence, implemented by shifting shoreline distance inward by 100 m before calculating shoreline proximity risk; and (iii) vegetation stress amplification, implemented by multiplying the NDVI-based vulnerability term by 1.2 and constraining it to 0–1.

Percentile-based normalization (p2–p98) for distance and trend layers was computed on land pixels only to prevent water surfaces from biasing scaling ranges. For comparability, domain-rank normalization and five-class thresholds were anchored to the 2018 reference distribution. Five classes were defined by 2018 quantile thresholds (q20, q40, q60, q80) and applied unchanged to 2023 and 2030. Outputs include ranked EVI maps for 2018 and 2023, the scenario-based EVI map for 2030, change maps for 2018–2023 and 2023–2030, class-area summaries, and class-to-class transition matrixes.

The three thematic indices highlight complementary drivers of ecological vulnerability in the Charvak Reservoir basin (Figure 5). The TGRI integrates slope, proximity to active faults, lithology, and soil susceptibility expressed by SOC. The spatial pattern (Figure 5a) shows high–very high values along the reservoir-facing escarpments of North Pskem and across the Kumbel–Kokand–Kenkol sector, where steep terrain coincides with weaker lithological units and lower SOC. Low TGRI dominates comparatively stable upland blocks and sectors under more resistant bedrock. Fault proximity contributes secondary, locally focused amplification of TGRI along mapped structures, consistent with structurally controlled slope weakening and drainage incision.

The HERI combines the NDWI trend risk and shoreline proximity (Figure 5b). High HERI values form a distinct nearshore belt driven primarily by short shoreline distances, whereas the NDWI trend component modulates the interior by elevating areas experiencing persistent drying. As a result, upland sectors away from the shoreline can still map as moderate–high where drying trends are sustained, whereas the immediate reservoir margin remains consistently high due to direct shoreline influence.

The CARI integrates NDVI trend risk, NDBI trend risk, proximity to infrastructure, and LST trend risk (Figure 5c). High CARI patches cluster along development and access corridors and in zones where built-up growth and infrastructure proximity overlap, with additional reinforcement where warming trends are pronounced. In contrast, more remote upland areas tend to remain low where development pressure is limited and vegetation trends are stable.

The five-class EVI map (Very low–Very high) integrates the three thematic domains and reveals coherent spatial structure across the Charvak basin (Figure 6). High–very high vulnerability (classes 4–5) is concentrated in the Kumbel–Kokand–Kenkol sector and along parts of the reservoir-facing North Pskem slopes, forming contiguous hotspots where predisposing conditions and multiple pressures overlap. In contrast, very low–low vulnerability (classes 1–2) is more prevalent across relatively stable uplands and remote interior sectors.

The spatial pattern is consistent with the domain indices by design: topographic–geological predisposition (TGRI) dominates where steep slopes and susceptible lithologies occur; HERI is highest in the reservoir-influenced belt and locally in drier interior areas; and CARI increases along accessible development corridors and thermal-stress patches. Global AHP weights indicate that slope and the built-up pressure component are the largest contributors, followed by shoreline proximity and infrastructure proximity, whereas fault proximity has the smallest global weight (Table 2). Class-area proportions for the long-term EVI map are reported in Table 3.

Figure 7a shows the 2023 EVI baseline map. For cross-year comparison, EVI was grouped into five classes (Very low to Very high) using fixed class thresholds derived from the 2018 reference distribution. High and very high classes are concentrated along reservoir-facing slopes and major valley corridors, whereas lower classes dominate the interior uplands and upper catchments. Open-water areas are masked. Within the nearshore band (≤150 m), the NDWI trend component is excluded and shoreline proximity governs the hydrological–ecological term, reducing shoreline-related artifacts in the change fields.

Figure 7b shows the observed EVI change between 2018 and 2023. Most of the basin exhibits low-to-moderate change magnitudes, with localized patches of increase and decrease. Class-based transition accounting for this period is summarized in Table 4.

Figure 7c presents the scenario-based EVI change from 2023 to 2030 under the combined-stressor assumptions. The strongest increases are spatially concentrated and follow areas most affected by the imposed shoreline and infrastructure stressors. Scenario transition outcomes are reported in Table 5.

Transition statistics were computed on the intersection of valid land pixels for each pair of years; therefore, total mapped land area in Tables 4–5 may differ from the single-date area totals reported in Table 3.

Domain decomposition improves transparency and interpretability in GIS–MCDA vulnerability mapping. A composite index is useful for zoning, but its compensatory nature can mask why a location attains a high overall score. Reporting the three domain subindices—TGRI, HERI, and CARI—addresses this by separating (i) terrain predisposition, (ii) hydrological–ecological pressure linked to reservoir influence and moisture-related change, and (iii) climatic–anthropogenic pressure associated with vegetation condition, built-up dynamics, accessibility, and thermal stress. This structure supports diagnostic interpretation of hotspots (identifying dominant drivers) and facilitates driver-specific recommendations without altering the composite workflow.

Dominant driver groups in the proposed framework are determined by the AHP weight structure rather than by a single indicator layer. The weighting assigns comparable overall influence to the topographic–geological and climatic–anthropogenic domains (0.40 each), with a smaller but meaningful contribution from the hydrological–ecological domain (0.20). At the indicator level, the model is most sensitive to slope- and accessibility-related pressures, including built-up and infrastructure proximity, as well as shoreline proximity. Lithology, SOC, and LST trend provide secondary modulation, whereas NDWI- and NDVI-related trends act as biophysical modifiers where persistent drying or vegetation decline is expressed. Fault proximity has the smallest global weight, reflecting its role as a contextual structural control that increases geological realism without dominating basin-wide patterns.

The topographic–geological component (TGRI) expresses terrain predisposition through slope steepness, lithological resistance, soil susceptibility (SOC), and structural setting. This configuration follows established process understanding for mountain environments, where steep gradients, weaker rock units, and structural discontinuities jointly condition slope instability and erosion-prone terrain. In the Charvak basin, elevated TGRI values along reservoir-facing escarpments and valley-side slopes therefore represent a plausible expression of background controls and are consistent with regional syntheses emphasizing lithological contrasts and tectonic structures ((Khakimov et al., 2024); (Liu et al., 2024); (Caleca et al., 2024); (Fazilova et al., 2025)). Given its comparatively small global weight (Table 2), fault proximity is best interpreted as a local modulator (e.g., fracture-controlled weakening and slope segmentation) rather than a basin-wide driver.

The HERI captures reservoir influence and moisture-related change by combining shoreline proximity with NDWI trend risk. High HERI values in the nearshore belt are expected because shoreline distance proxies bank undercutting, saturation, and drawdown effects typical for regulated reservoirs. Away from the shoreline, NDWI trend risk refines the signal by highlighting areas where persistent moisture change contributes to ecological stress. To suppress shoreline-driven artifacts in the trend term, NDWI trend risk is excluded in the immediate nearshore band, whereas shoreline proximity governs the hydrological–ecological pressure in this zone (Wang et al., 2024). This interpretation is consistent with reported effects of water level fluctuations on bank stability and riparian condition in mountain reservoirs ((Juliev et al., 2023); (Fazilova et al., 2026)).

The CARI integrates vegetation-condition change (NDVI trend risk), built-up dynamics (NDBI trend risk), infrastructure proximity, and thermal stress (LST trend risk). Its corridor-focused pattern is plausible because accessibility concentrates development pressure, land cover modification, and associated heat and fragmentation effects, particularly around tourism and recreation zones. Colocation of higher CARI values with built-up growth hotspots and proximity to infrastructure aligns with evidence that land use intensification around reservoir margins and along major access routes can amplify ecological stress ((Juliev et al., 2019); (Fazilova et al., 2026)). Including LST trend supports interpretation of CARI as a pressure index sensitive to warming-prone surfaces and land cover transitions.

Direct quantitative validation is constrained by limited inventories and field measurements; however, agreement between mapped thematic patterns and expected process behavior provides a first-order plausibility check. In addition, the relative influence implied by the AHP structure (Table 2) is consistent with the spatial signals: slope and built-up pressure have the largest global shares, with shoreline proximity and infrastructure proximity providing major additional contributions.

Recent ecological vulnerability studies provide the broader methodological context for such multifactor mapping, including AHP-based weighting frameworks (Hu et al., 2021), remote-sensing driven vulnerability indices at regional scales (Wu et al., 2021) hybrid modeling strategies (Dossou et al., 2021), and spatiotemporal evolution analyses across major plateaus and basins ((Xia et al., 2021); (Zhang et al., 2022); (Zhang et al., 2025)). Relative to these large-area applications—often dominated by broad climate–vegetation gradients and extensive land cover transitions—the Charvak basin exhibits a more localized signal shaped by reservoir shoreline processes and accessibility corridors. In this sense, the results remain consistent with corridor-related pressure patterns emphasized in vulnerability assessments that incorporate urbanization and infrastructure effects (Wu et al., 2022) while extending AHP-oriented mapping (Hu et al., 2021) through explicit reporting of domain subindices and a baseline-consistent near-term scenario stress test tailored to a reservoir-affected mountain basin.

The 2030 maps represent a scenario-based stress test rather than a deterministic forecast. The near-term horizon (2030 relative to the 2023 baseline) is chosen to support actionable planning and to remain aligned with the 2018–2023 reference window. The scenario combines baseline-anchored extrapolation of land surface indicators with explicit stressors that emulate plausible intensification of shoreline influence and corridor-related pressure; accordingly, the outputs should be interpreted as spatial sensitivity to these assumptions, not as exact future states.

The scenario indicates that vulnerability increases most where multiple pressures co-occur on predisposed terrain. Shoreline-related stress expands the nearshore belt of elevated vulnerability, whereas accessible valleys respond primarily to infrastructure- and land use-related drivers. This supports a two-zone planning logic: (i) shoreline belts where reservoir-level dynamics may amplify bank processes and ecological stress, and (ii) access corridors where development pressure may intensify disturbance. The domain indices provide a driver-specific basis for interventions—shoreline monitoring and protection where HERI dominates, land use regulation where CARI dominates, and hazard-aware restrictions where high TGRI indicates strong predisposition—supporting ecosystem-oriented mitigation and spatial planning in the Charvak basin ((Juliev et al., 2023); (Fazilova et al., 2026)).

From a water management perspective in Uzbekistan, the high-vulnerability shoreline belt is operationally relevant because seasonal drawdown and high-water phases can accelerate bank erosion and undercutting, affecting nearshore ecosystems and recreational assets; this motivates shoreline protection belts and targeted monitoring of erosion-prone segments. Hotspots along the main inflow valleys further indicate priority corridors where sediment delivery, slope processes, and infrastructure expansion may jointly increase management pressure, motivating coordinated monitoring and land use regulation (Juliev et al., 2023)(Fazilova et al., 2026).

Several limitations should be noted. First, the magnitude and ranking of vulnerability depend on indicator selection, orientation-consistent normalization, and the AHP weighting scheme; reporting both the composite EVI and the domain subindices (TGRI, HERI, and CARI) partially mitigates this by improving transparency. Second, results are constrained by the spatial resolution and temporal availability of the input datasets, including mixed-resolution ancillary layers (e.g., SOC), and by limited independent inventories for quantitative validation of mapped hotspots. Coupled NDVI–LST drought metrics, such as the temperature–vegetation dryness index (TVDI) and the normalized temperature drought index (NTDI), were not computed; long-term NDVI and LST trends were used instead as separate, interpretable indicators. Additional soil predictors (e.g., soil moisture/texture/depth) were not included because of limited harmonized basin-scale layers; this is a key limitation for predictive interpretation. Third, the 2030 output denotes a scenario-based stress test rather than a deterministic forecast, as it depends on assumptions about short-term trend persistence and prescribed stress adjustments anchored to the baseline. Finally, some predictors act mainly as contextual controls (e.g., fault proximity); thus, local structural effects may be underrepresented where detailed geotechnical data are unavailable.

Future work should test multiple scenario settings for land use/tourism development and reservoir operation, and apply alternative time-series approaches (e.g., nonlinear trends and breakpoint detection) where linear extrapolation is not appropriate. Robustness can be strengthened through uncertainty and sensitivity analyses of key modeling choices (weights, normalization, and fixed class thresholds). Validation could be expanded using independent evidence such as documented shoreline-instability segments, land use change inventories, and targeted field observations in representative shoreline and corridor hotspots.