Optimizing Rainwater Harvesting Pond Locations in a Semi-Arid Region: A Geospatial Multi-Criteria Decision Analysis, Hammad Basin, Northeastern Jordan ()
1. Introduction
Jordan’s Water Scarcity
Jordan is one of the world’s most water-scarce countries. The increasing demand for water due to economic development, agriculture, population growth, and the expansion of human activities in local communities has placed significant strain on Jordan’s water resources. The nation has grappled with challenges in managing its surface water resources. To address this, Jordan has undertaken ambitious efforts to capture as much surface water as possible through the construction of concrete and earthen dams, as well as rainwater harvesting ponds, locally known as “hafayer”.
Rainwater Harvesting Projects and Their Role
Approximately 493 rainwater-harvesting ponds have been constructed across Jordan, primarily targeting semi-arid regions and low-gradient valleys. These ponds collect surface runoff during rainy seasons and serve multiple purposes, including fodder cultivation, the cultivation of native plants in their natural habitats, and providing drinking water for livestock during grazing seasons. Moreover, these projects contribute to combating desertification and mitigating the impacts of climate change on vegetation and biodiversity. Typically, these ponds are square, measuring approximately 150 meters on each side, with a total depth of 4 - 6 meters. The excavated material is used to create embankments around the pond, increasing its capacity to around 100,000 cubic meters (Figure 1).
For millennia, various civilizations have inhabited the lands that now constitute Jordan, establishing agricultural settlements and cities. These communities encountered a similar arid to semi-arid environment and developed innovative techniques to collect and store rainwater for use during dry periods (AbdelKhaleq & Alhaj, 2007).
With advancements in remote sensing and geographic information systems, numerous researchers have conducted extensive studies to identify the most suitable locations for constructing rainwater harvesting ponds and small dams. Notable studies include (AlAyyash et al., 2012; Al-Adamat et al., 2012; Al-Adamat, 2008; Albalawneh et al., 2015; Al-Adamat et al., 2010; Awawdeh, et al., 2010; Ziadat et al., 2006; Al-Sababhah, 2023).
Objective of the Current Study
This study aims to identify optimal locations for the construction of rainwater harvesting ponds with a capacity of approximately 100,000 cubic meters within the Hadalat Dam watershed. This capacity is adopted as the standard design capacity for such rainwater harvesting ponds (hafayer) implemented by the relevant governmental authorities in the region. By integrating a Multi-Criteria Decision Analysis (MCDA) model with geospatial data, this research seeks to develop a detailed suitability map for pond placement. The goal is to optimize rainwater capture, minimize negative environmental impacts, and support sustainable livelihoods, particularly in livestock grazing and agricultural practices.
Study Area
The study area is located in northeastern Jordan and covers an area of 13,202 square kilometers (Figure 2). Approximately 97.4% of the study area is situated in Jordan, while the remaining portion extends into Iraq, Syria, and Saudi Arabia. An international highway connecting Amman and Baghdad traverses the region. The area is home to several communities, most notably Ruweished and Manshiyat Al Ghayath, with a combined population of 7490. Additionally, the area attracts nomadic herders from across the kingdom during the grazing season.
Figure 1. Examples of rainwater harvesting ponds established in the study area.
Literature Review
Figure 2. Location map of the study area in Jordan.
Systematic literature reviews confirm that the integration of Geographic Information Systems (GIS) with Multi-Criteria Decision Analysis (MCDA) constitutes the most established methodology for pinpointing optimal rainwater harvesting (RWH) sites in arid and semi-arid regions globally (Wable et al., 2023; Sadushan & Neluwala, 2024). Most existing studies focus on incorporating biophysical, hydrological, and geological factors, such as slope, land use, and soil types (Adham et al., 2018; Gaikwad & Pawar, 2019). Methodologically, the vast majority of similar research in regions like Iraq and Egypt employs the Analytic Hierarchy Process (AHP) to determine criteria weights, often to ensure mathematical consistency (Hassan et al., 2025; Abdalla, 2025). However, AHP is often critiqued for its subjectivity and its lack of adaptability to the complex multi-objective requirements and severe local constraints characteristic of arid ecosystems, particularly where socio-environmental data are unevenly available (Ahmed et al., 2025). While some studies have attempted to integrate sustainability factors (Sadushan & Neluwala, 2024), few have successfully linked selection criteria to specific local objectives like vital rangeland management or livelihood support (Gharde et al., 2025). The current study, therefore, presents a novel methodology designed to overcome the functional and methodological shortcomings of previous research. Instead of relying on AHP-derived weights, an MCDA model was developed using an adaptive multi-objective weighting strategy that assigns the highest priority (45% weight) to maximizing the hydrological catchment volume a critical necessity in the Hammad Basin due to its exceptionally low rainfall. Furthermore, this work explicitly integrates crucial rangeland indicators (such as NDVI) and livestock distribution into the model, ensuring that the optimal sites are not only hydrologically feasible but also meet the pressing socio-economic needs of the local population. This constitutes a fundamental differentiation in selecting Hafair (pond) sites within the arid regions of Jordan.
2. Methodology
MCDM techniques have proven invaluable in the decision-making process for identifying suitable locations for rainwater harvesting ponds, particularly in arid and semi-arid regions. These regions often face water scarcity, making efficient water management a critical concern. By considering a multitude of often-conflicting criteria, such as hydrological characteristics, soil type, topography, and socio-economic factors, MCDM models can systematically evaluate potential sites. The methodologies employed in these studies typically involve weighting and combining various criteria to arrive at a comprehensive assessment. The flexibility of MCDM approaches allows for the inclusion of both quantitative and qualitative data, making them adaptable to diverse regional contexts. Ultimately, the goal of these studies is to provide decision-makers with a robust scientific basis for selecting sites that maximize the benefits of rainwater harvesting while minimizing negative environmental impacts.
The specific factors considered in MCDM studies for rainwater harvesting site selection, as well as the weights assigned to these factors, vary depending on the region, available data, and the researchers’ perspectives on the relative importance of each factor (Ammar et al., 2016). While there is general agreement on the most relevant factors, developing a standardized model with fixed weights is challenging. To foster sustainable development, it is recommended to implement rainwater-harvesting projects in phases, allowing for the evaluation and refinement of the MCDM model over time. This iterative approach involves constructing a limited number of ponds, assessing their performance, and adjusting the model by adding, removing, or re-weighting criteria. Additionally, parallel research on pond design and sedimentation rates can further optimize the process. Such a phased approach can help minimize errors, reduce costs, and maximize benefits.
2.1. Multi-Criteria Decision Analysis (MCDA) Methodology and Adaptive Weighting
The MCDA model utilized an Adaptive Weighting Strategy based on a Hierarchical Structuring of three primary, distinct objectives (goals) that guided the assignment of weights to all criteria. This approach was essential for achieving Multi-Objective Optimization tailored to the extreme environmental and socio-economic constraints of the arid Hammad Basin.
Goal I: Water Quantity and Accessibility (Maximum Runoff Potential)
The first and most critical objective focuses on water quantity and the ability of runoff to reach a proposed site, accounting for 45% of the total model weight. Given the extremely low precipitation rate of 87.3 mm/year and a low annual runoff depth of only 5.5 mm in this arid environment, maximizing the contributing catchment area is paramount. The primary goal of constructing harvesting ponds is to efficiently capture surface water before it is lost to evapotranspiration in wide mudflats or infiltration. To ensure the standard capacity of 100,000 m2 used in Jordanian water harvesting projects is met, the model requires a minimum effective catchment area of 100 km2 for each proposed pond.
This goal is primarily represented in the model by Flow Accumulation, derived from the Digital Elevation Model. However, due to the hydrological complexity of the area—including existing dams, small ponds (Sinks), and mudflats, the standard flow accumulation is potentially misleading. These existing features naturally or artificially intercept runoff, especially during low runoff depth events (which barely reach 1 mm). To address this, a detailed Tree Diagram of Nested Catchment Area analysis was conducted to produce two distinct flow accumulation rasters:
No-Sink Flow Accumulation: Represents runoff flow assuming all interception structures are inactive or full.
Sink Flow Accumulation: Represents runoff flow considering the upstream interception of Sinks.
These two flow accumulation rasters were processed using the Channel Density tool, yielding two final criteria: Sinks Channel Density and No Sinks Channel Density. The 45% weight for this goal was split equally between these two criteria (0.225 each) to cover the full spectrum of runoff depth events, from low intensity (where sinks capture all flow) to higher intensity (where sinks overflow and contribute to downstream areas).
Goal II: Justification of Need (Rangeland and Livestock Demand)
The second objective, representing 31% of the total weight, justifies the need for pond construction based on local activities, primarily the provision of drinking water for livestock during grazing periods. This goal is expressed through criteria that promote rangeland availability and sustain grazing activities:
Rainfall (0.10) and Runoff Depth (0.10): High weights were assigned as these are key factors enabling rangeland growth, particularly in wide wadis and at the edges of mudflats where seasonal rain alone is insufficient.
Actual Evapotranspiration (0.025): Sites with lower Actual Evapotranspiration are favored to reduce plant stress and preserve natural rangelands.
Potential Evapotranspiration (0.020): This weight is assigned to assess the suitability for establishing new managed rangelands, which are often supported by the Ministry of Agriculture using water from the harvesting ponds.
Soil Moisture (0.025): Higher soil moisture areas are favored as they allow vegetation to withstand drought for longer periods.
Land Cover/Land Use (0.04): This is a direct measure of demand, representing observed vegetation and livestock presence. This weight was subdivided into:
MSWI (0.02): Represents actual observed vegetation via satellite.
Livestock Density (0.02): Further subdivided based on local livestock composition: Sheep (0.015) and Goats (0.005), reflecting the higher ratio of sheep in the study area.
Goal III: Water Retention and Cost Efficiency
The final objective, accounting for 24% of the total weight, focuses on the construction and maintenance efficiency of the ponds. This goal seeks to minimize water loss (evaporation and seepage) and reduce construction/maintenance costs (excavation and sedimentation). The criteria included:
Potential Evapotranspiration (0.03): An additional weight is assigned here to evaluate direct water loss from the pond’s surface due to evaporation, favoring areas with lower Potential Evapotranspiration. (The total weight for this criterion is 0.020 + 0.03 = 0.05)
Total Insolation (0.045): The most important factor in this goal, selecting sites with the lowest Total Insolation to minimize evaporative losses due to direct solar radiation.
Curve Number (0.025): Although Curve Number is implicitly included in the runoff depth calculation (Goal I), an explicit weight is given here. Sites with a higher Curve Number are preferred to minimize deep percolation and seepage of collected water into the subsoil.
Soil Depth (0.025) and Estimated Ksat (0.025): Both are weighted to avoid sites with large soil depths or high estimated saturated hydraulic conductivity (Ksat), which lead to excessive seepage.
Topography (0.09): This collective weight is assigned to the topographic group of criteria (e.g., Slope, Valley Depth, etc.) to indirectly avoid areas with high excavation costs, excessive sedimentation rates, or structural collapse hazards.
Data Standardization and Aggregation
To ensure methodological coherence across the diverse datasets utilized:
Geospatial Data Standardization
Given the integration of geospatial datasets with varying resolutions (e.g., high-resolution DEM (12.5 m), soil data (250 m), and climate data (4 km), all spatial layers were rigorously standardized to a unified spatial resolution of 12.5 m using the appropriate Resampling technique (Nearest Neighbor for categorical data and Bilinear Interpolation for continuous data). Furthermore, all layers were transformed to a uniform Coordinate Reference System (JTM—Jordan Transverse Mercator) prior to the aggregation phase.
Criteria Aggregation
Figure 3. Main steps of method, the green color indicates the geospatial data, which has been assigned specific weights within the mcdm model, influencing the final suitability map.
The weights derived from the Multi-Objective Optimization process were aggregated using the Weighted Linear Combination (WLC) method, which is the standard methodology in MCDA-GIS research. The WLC formula utilized the three principal goal weights as follows:
where Wi: is the specific weight assigned to the criterion Ci based on the three objectives, and Ci: is the standardized value of criterion (i) after spatial classification.
In this study, an MCDM model was developed to identify suitable locations for rainwater harvesting ponds. The model primarily relied on six data groups: soil, topography, climate, land use, channel density, and insolation. Additionally, a hydrological network analysis was conducted to provide a more comprehensive understanding of the study area. Table 1 presents the final weights assigned to all criteria, while Figure 3 provides a summary of the methodology. A detailed breakdown of these data groups is as follows.
Table 1. Relative weights for criteria used in the MCDM model.
Factor’s Type |
Factor Type’s Weight |
Factor |
Factor’s Weight |
Channel
Density |
0.450 |
No Sinks Channel Density |
0.225 |
Sinks Channel Density |
0.225 |
Topography |
0.090 |
Relative Slope Position |
0.020 |
Valley Depth |
0.020 |
Topographic Wetness Index |
0.020 |
Profile Curvature |
0.005 |
Longitudinal Curvature |
0.005 |
Flow Line Curvature |
0.005 |
Plan Curvature |
0.005 |
Tangential Curvature |
0.005 |
Cross-Sectional Curvature |
0.005 |
Insolation |
0.045 |
Total Insolation |
0.045 |
Soil |
0.075 |
Curve Number |
0.025 |
Estimated Ksat |
0.025 |
Soil Depth |
0.025 |
Climate |
0.300 |
Rainfall |
0.100 |
Runoff |
0.100 |
Actual Evapotranspiration |
0.025 |
Potential Evapotranspiration |
0.050 |
Soil Moisture |
0.025 |
Land Cover Land Use |
0.040 |
MSWI |
0.020 |
Goats |
0.005 |
Sheep |
0.015 |
2.2. Hydrological Network Analysis
The study area was delineated within a single, point-defined catchment, centered on Hadalat Dam. Analysis revealed a bipartite catchment: a predominantly Jordanian northern section draining directly into Hadalat Dam, and a primarily Saudi southern section draining into the Qia’an basin, located southwest of the Jordanian Hammad Desert. Given the immense capacity of the Qia’an basin, potentially accommodating up to 493 million cubic meters, which significantly exceeds the total surface runoff generated by the upstream catchment, it was delineated as a separate hydrological unit. Consequently, the Qia’an basin was excluded from the study area and is represented in gray on some maps.
To identify water accumulation zones such as ponds, depressions, and sinks, including those associated with dams, a Normalized Difference Water Index (NDWI) analysis was conducted using Sentinel-2 satellite imagery collected during the rainy seasons of 2018-2022. Subsequently, all sub-watersheds within these zones were delineated. A tree diagram was employed to hierarchically classify the sub-catchments.
Two distinct flow-accumulation rasters were generated. The first, calibrated using a filled Digital Elevation Model (DEM) and referenced to the dam point, simulated a scenario of high surface runoff filling all water bodies and eventually discharging into the Hadalat Dam. The second raster was an aggregate of all sub-watersheds, with flow accumulation interrupted at each identified water body to simulate active channels under low flow conditions.
Channel Density
Based on the two previously mentioned flow accumulation versions, two versions of channel density were produced by dividing the flow accumulation raster by the highest value in the raster. These are channel sinks density based on the interrupted flow accumulation, and channel no sinks density based on the uninterrupted flow accumulation.
2.3. Topography
A high-resolution (12.5-meter) terrain-corrected Digital Elevation Model (DEM) served as the foundational dataset for our topographic analysis. Acquired from the Alaska Facility for Satellite Observing Systems (ASF) in 2008 using the ALOS PALSAR dataset, this DEM provided a detailed representation of the study area.
Valley Depth
To identify terrain depressions, valley depth was computed using SAGA GIS (Wood, 1996; Wood, 2009; Conrad et al., 2015; Wilson, 1988). Assigned a weight of 0.020 in the Multi-Criteria Decision Analysis (MCDA) model, this parameter indicated the suitability of areas for pond construction, with deeper valleys being more favorable.
Topographic Wetness Index (TWI)
Calculated using SAGA, the Topographic Wetness Index (TWI) quantifies a location’s predisposition to saturated conditions, considering both local slope and contributing upstream area. With a weight of 0.020 in the MCDA, TWI highlights areas with high water potential, making them suitable for pond construction.
Relative Slope Position
The relative slope position, computed using SAGA, determines a grid cell’s location relative to the highest and lowest points within a neighborhood. Assigned a weight of 0.020 in the MCDA, this parameter indicates the suitability of areas for rainwater harvesting. Areas with high relative slope positions, typically near ridgetops, were less suitable due to potential runoff and erosion.
Curvature
SAGA was used to calculate curvature, a composite parameter providing information on the land surface shape. With a weight of 0.06 in the MCDA, curvature was a significant factor in identifying suitable pond locations. Concave areas, indicated by negative curvature values, were preferred due to their water-concentrating properties.
2.4. Insolation
To assess the impact of solar radiation on site suitability, annual total insolation was integrated into the MCDA model with a weight of 0.045, and it was computed using SAGA.
2.5. Soil Data
Leveraging a comprehensive soil database provided by the Food and Agriculture Organization (FAO) and other sources covering the study area (Table 1), we developed a robust set of parameters to assess the soil’s hydrological behavior.
Estimated Saturated Hydraulic Conductivity (Ksat)
The Rosetta pedotransfer function (Zhang & Schaap, 2017) was utilized to predict Ksat values. This function, renowned for its accuracy and reliability, employs FAO soil data including texture (sand, silt, clay), bulk density, and water content at specific matric potentials. By selecting the most suitable Rosetta model based on the study area’s specific characteristics and available input data, we estimated Ksat for each grid cell.
To account for the influence of fine earthen material, we refined the Ksat estimates by multiplying them by the ratio of fine earthen material to the total material content (including coarse fragments). This adjustment enhanced the parameter’s sensitivity to the specific characteristics of the fine-grained soil fraction.
Curve Number
The Global Curve Number dataset at 250-meter resolution (GCN250) developed by Jaafar et al. (2019) was incorporated into the study. This dataset, derived from a combination of European Space Agency land cover data and hydrologic soil group data, provides a globally consistent representation of CN values. The GCN250 dataset was selected for its high spatial resolution, which is crucial for capturing the heterogeneity of hydrological processes within the study area.
Depth to Bedrock
To account for the influence of bedrock depth on the hydrological processes, the absolute depth to bedrock dataset from Hengl et al. (2017) was incorporated into the analysis. This dataset, derived from a combination of remote sensing data and soil profiles, provides global estimates of bedrock depth at a 250-meter resolution.
2.6. Climate Data
In this study, we employed the TerraClimate dataset (Abatzoglou et al., 2018) acquired from ClimateEngine to identify suitable locations for rainwater harvesting ponds in semi-arid environments. We opted for this dataset due to its high spatial resolution (~4 km) and extended time series (1958-2023), encompassing the period of our analysis (1970-2023). TerraClimate offers a comprehensive suite of climate variables, several of which were incorporated into our Multi-Criteria Decision Making (MCDM).
Rainfall
Rainfall data were obtained from the TerraClimate dataset at a spatial resolution of 1/24˚ (~4 km) and a monthly temporal resolution. To assess data accuracy, a comparison was conducted with the Global Historical Climatology Network (GHCN) database, demonstrating strong correlations (median Pearson’s correlation coefficient > 0.9) for annual rainfall.
Runoff
Runoff estimates were derived from a one-dimensional modified Thornthwaite-Mather climatic water-balance model (WBM) within the TerraClimate dataset at a spatial resolution of 1/24˚ (~4 km) and a temporal resolution of monthly. Validation against streamflow data from 587 pristine river basins indicated a reasonable capture of inter-annual variability (median Pearson’s correlation coefficient = 0.80).
Actual Evapotranspiration
Actual evapotranspiration data were calculated using the Penman-Monteith approach within the TerraClimate dataset at a spatial resolution of 1/24˚ (~4 km) and a temporal resolution of monthly. Comparison with observations from the FLUXNET network showed moderate correlations (median Pearson’s correlation coefficient = 0.77) for annual ET0.
Potential Evapotranspiration
Potential evapotranspiration was estimated using the Penman-Monteith approach within the TerraClimate dataset at a spatial resolution of 1/24˚ (~4 km) and a temporal resolution of monthly. While not explicitly validated, it was indirectly assessed through comparisons with actual evapotranspiration.
Soil Moisture
Soil moisture data were estimated within the TerraClimate dataset based on plant extractable soil water capacity data at a spatial resolution of 1/24˚ (~4 km) and a monthly temporal resolution.
2.7. Land Cover, Land Use, and Livestock
Modified Soil Adjusted Vegetation Index (MSAVI)
The Modified Soil Adjusted Vegetation Index (MSAVI) was calculated using an annual average to enhance the vegetation signal while minimizing soil background noise, providing a more accurate representation of vegetation cover, especially in areas with sparse vegetation (Qi et al., 1994).
Livestock
The Gridded Livestock of the World-2015 (GLW4) dataset, obtained from the Harvard Dataverse, was integrated into the MCDA model to assess the impact of goat and sheep densities on pond suitability (Gilbert et al., 2022a; Gilbert et al., 2022b).
3. Results
3.1. Hydrological Network Analysis
The study area encompasses seven dams, five of which are concentrated in the northwestern region within the Hammad Valley plain. The remaining two are situated in the northeastern part. Some of these dams, dating back to ancient times, have been restored and rehabilitated. Although characterized by large storage capacities due to the extensive upstream flatlands, these dams frequently experience complete drying during the summer months.
Figure 4. Active sinks, mudflats, dams, and rainwater harvesting ponds and their sub-watershed in the study area.
The area also contains 23 large, locally named mudflats that transform into temporary water bodies following precipitation events. These mudflats can persist for weeks or even months before drying out. Additionally, 41 rainwater-harvesting ponds are found distributed throughout the study area. While ponds within the western basalt plateau and central Hammad Valley are typically located within mudflats, those in eastern regions are situated along major valleys.
Besides the larger mudflats, 246 smaller mudflats, termed “Active Sinks” due to their efficient collection of surface runoff, have been identified using NDWI and satellite imagery. Many of these sinks merge with one or more adjacent depressions, resulting in shared sub-watersheds. Consequently, the total number of water-collecting depressions contributing to the surface runoff water balance is 245 (Figure 4).
Analysis of closed depressions using a 12.5-meter resolution enhanced DEM and NDWI revealed approximately 800,000 such depressions (inactive sinks) with depths exceeding 0.25 meters, distributed across the entire study area (Figure 6(A)).
Tree Diagram of Nested Catchment Area
Figure 5. Active channel at low runoff rates.
A tree diagram analysis of the nested catchment areas within the study region reveals a complex hierarchy of sub-watersheds, reaching up to the 18th order. This indicates that, in certain areas, surface runoff passes through as many as 18 natural or artificial ponds before ultimately reaching Hadalat Dam. The northern parts of the study area exhibit a maximum order of six, while the southern sections range from sixth to twelfth order. The longest sequence, consisting of 18 orders, is observed in the eastern regions (Figure 6(B)).
Flow Accumulation Analysis
Flow accumulation analysis, overlaid on the nested catchment diagram, quantified the drainage area in square meters for each grid cell (Figure 5). The findings indicate that the drainage area contributing directly to Hadalat Dam is approximately 700 square kilometers. Other ponds within the region exhibited drainage areas ranging from 50 to 400 square kilometers.
Figure 6. Active and inactive sinks starting at 0.25 m depth (A), tree diagram of nested catchment area (B), channel no sinks density (C), and channel sinks density (D).
A detailed analysis of rainwater harvesting ponds revealed that some ponds possess insufficient drainage areas relative to their storage capacity. For instance, a pond with a capacity of 100,000 cubic meters typically requires a drainage area of at least 100 square kilometers. Conversely, other ponds have significantly larger drainage areas than necessary, suggesting the potential to increase the number of rainwater harvesting ponds in these upstream areas.
Channel Density Analysis
The channel no sink density map, highlighting the major valleys and the spatial weights of uninterrupted flow accumulation used in the MCDM model, identifies three primary valleys. Additionally, the channel sinks density map depicts the spatial weights of interrupted flow accumulation (Figure 6(C) and Figure 6(D)).
3.2. Topographic Factor Analysis
Figure 7. Digital elevation model (A); relative slope position (B); slope (C); aspect (D); topographic wetness index (E); and valley depth (F).
Utilizing a Digital Elevation Model (DEM) as a foundational dataset (Figure 7(A)), spatial distribution maps of topographic factors integrated into the Multi-Criteria Decision Making (MCDM) model were generated. A detailed analysis of these factors follows:
Relative Slope Position
Despite the exceedingly low slope gradients within the study area, ranging from 0 to 7 degrees with a mean of 0.6 degrees (Figure 7(B)), the Relative Slope Position grid cell revealed spatially variable slope distributions. Notably, areas of relatively high slope inclination were identified, exhibiting linear extents predominantly oriented southeast-northwest in the eastern portion of the study area and southwest-northeast in the western part. Conversely, low-lying areas formed contiguous expanses across the entire study region, interspersed with isolated elevated areas (Figure 7).
Valley Depth
The Valley Depth grid cell highlighted the depth of valleys relative to adjacent elevations for each spatial unit on the map. Areas characterized by significant depths, optimal locations for establishing rainwater harvesting ponds, were prominently featured. Relatively deep valleys were concentrated in the eastern and western sections of the study area, while the central regions exhibited a dearth of such features. Notably, the Hammad Valley failed to manifest as a distinct depression on the map, instead appearing as a broad, flat expanse (Figure 7).
Topographic Wetness Index
The Topographic Wetness Index map indicated a moderate to high propensity for the region’s topography to retain moisture (values ranging from 6 to 15), signifying that the study area, on the whole, comprises terrain that inhibits surface runoff and promotes the dispersion of overland flow into scattered water bodies across the region (Figure 7). This finding substantiates the rationale for constructing small-scale rainwater harvesting ponds in dispersed locations, advocating for localized utilization rather than relying on large dams in lower-lying areas.
Curvature
Concave terrain is generally favored for locating rainwater harvesting ponds. This is because concave areas effectively collect and direct water towards these ponds. Different curvature analyses reveal the degree of concavity or convexity along specific lines or at specific points, influencing processes such as water flow acceleration or deceleration, water dispersion or convergence, sediment yield, erosion rates, excavation costs, and slope stability.
In detail, cross-sectional curvature measures curvature perpendicular to the steepest slope and is useful for identifying erosion and deposition zones (Figure 8(A)). Flow line curvature, measured along flow paths, also helps in identifying these zones (Figure 8(B)). Longitudinal curvature, measuring curvature along the steepest slope, is employed in this study to indicate water flow convergence and divergence, as well as acceleration and deceleration (Figure 8(C)). Plan curvature, which reflects the change in terrain aspect in the horizontal plane, is useful for determining water flow convergence or divergence (Figure 8(D)). Profile curvature measures curvature along the steepest slope and is used in this study as an indicator of slope stability around pond locations (Figure 8(E)). Tangential curvature, measured perpendicular to the gradient, is considered a factor in excavation costs (Figure 8(F)).
![]()
Figure 8. Cross-sectional curvature (A); flow line curvature (B); longitudinal curvature (C); plan curvature (D); profile curvature (E); and tangential curvature (F).
While all curvature measures provide insights into multiple parameters, their significance varies. For advanced analysis, the sequential arrangement of concave and convex areas along flow paths to proposed ponds can be studied. For example, a suitable concave location preceded upstream by a wide, shallow concave area on a flow line curvature map can reduce water velocity and promote sediment deposition before the water enters the pond.
3.3. Insolation
Figure 9. Total insolation.
The study revealed exceptionally high annual total insolation values, reaching up to 2420 kWh/m2 in certain areas. Even semi-shaded regions received substantial solar radiation, with a minimum of 1800 kWh/m2 annually and an average of 2222 kWh/m2 across the entire study area. Consequently, solar radiation plays a significant role in pond water evaporation. Selecting pond sites with lower total insolation could potentially reduce evaporation losses. However, this approach faces several challenges. For instance, spatial analysis of areas with total insolation below 2100 kWh/m2 indicated that 98% of semi-shaded zones had east-west extents less than 50 meters and north-south extents less than 21 meters. Moreover, finding suitable semi-shaded areas within the remaining shadow zones exceeding 40,000 square meters was challenging. Additionally, semi-shaded areas were concentrated on north-facing slopes, with their influence typically limited to the slope itself. Most valleys, being wide and shallow, were fully exposed to sunlight throughout the day (Figure 9).
3.4. Soil Data Analysis
Absolute Depth to Bedrock
The absolute depth to bedrock refers to the estimated depth from the land surface to the first impermeable or low-permeability rock layer. It represents the depth to the uppermost unconfined aquifer, although it may not necessarily consist entirely of soil in the strict sense. Rather, it typically denotes a hydrostratigraphic unit with high permeability extending down to the first confining layer. Analysis of absolute depth to bedrock data revealed significant spatial variability in the depth to the confining layer within the study area. Most of the western regions, covered by basaltic flows, appeared to lack any surficial permeable layer, soil, or permeable sediments, except for mudflats which exhibited relatively shallow depths. In contrast, the eastern regions showed permeable zones extending to depths of approximately 180 meters. This map can be instrumental in avoiding areas with thick permeable zones or in implementing measures such as lining the bottom and sides of rainwater harvesting ponds with impermeable materials to conserve the collected water.
Estimated Ksat
In this study, various soil properties, such as sand, silt, and clay content, as well as the water content at standard pressures and bulk densities, were integrated into the single parameter of saturated hydraulic conductivity (Ksat). Additionally, actual volumetric values were adjusted using the proportion of coarse fragments. The spatial distribution of Ksat values within the study area exhibited significant variability (Figure 10(C)). Higher Ksat values were concentrated in extensive areas in the central part of the study area, extending from south to north along a series of depressions and mudflats forming the so-called Hammad Wadi, situated between the basaltic flows in the west and the exposed Cenozoic carbonate rocks in the east. The eastern and western regions displayed lower Ksat values. This map can independently support decision-making regarding the selection of sites for rainwater harvesting ponds for various purposes, including on-surface direct use ponds or those designed to augment groundwater recharge. However, in this study, Ksat was considered an inverse property since the focus was on identifying suitable locations for on-surface direct use rainwater harvesting ponds. The Ksat values ranged from 6 cm/day to 37 cm/day, with an overall average of 21.8 cm/day.
Curve Number
Based on the analysis of global gridded curve number data (Jaafar et al., 2019) for the study area, the region exhibited a high propensity for surface runoff generation. The majority of the study area, as depicted in Figure 10(B), demonstrated a curve number of 80 under ARC1 conditions, indicating a significant potential for surface runoff. In particular, valleys with very low slopes exhibited even higher curve numbers, reaching up to 85, likely due to the accumulation of fine sediments. Conversely, relatively lower curve numbers were observed in scattered locations, especially in the southernmost part of the study area. Statistical analysis revealed a mean curve number of 80, with a minimum of 57 and a maximum of 85, resulting in a range of 28. Given the generally arid nature of the region, ARC2 and ARC3 conditions were deemed negligible for this analysis.
Figure 10. Absolute depth to bedrock (A); curve number (B); and estimated ksat (C).
3.5. Climate Factor Analysis
Annual Rainfall
Analysis of rainfall data for the region from 1970 to 2023 reveals extremely low annual averages. The southeastern highlands exhibited the highest annual rainfall, reaching 113 mm. Generally, eastern regions recorded slightly more than 100 mm of annual rainfall. Conversely, the western basaltic highlands had annual averages approaching 100 mm. Extensive areas in the central parts of the study area recorded comparatively lower rates, ranging between 70 and 80 mm annually. Overall, the spatial average of annual rainfall was 87.3 mm (Figure 11(A)).
Figure 11. Annual rainfall (A); annual surface runoff (B); annual soil moisture (C); annual actual evapotranspiration (D); annual potential evapotranspiration (E); and evapotranspiration ratio (F).
Surface Runoff
Annual surface runoff in the study area, as modeled by TerraClimate, was also exceptionally low. The spatial distribution of surface runoff showed a maximum of 5.5 mm per year and a minimum of only 3 mm. The overall average for the entire study area was 4 mm per year. Higher values were concentrated in the eastern and western parts of the study area, while the lower central regions recorded the lowest values (Figure 11(B)).
Soil Moisture
Annual soil moisture levels indicated widespread aridity in the region. Soil moisture content, as depicted in the monthly data series, consistently showed minimal retention. While most soil moisture evaporates within a few days of a rainfall event, the monthly cumulative soil moisture often approached 0 mm. Overall, the annual average soil moisture ranged between 0 and 0.2 mm (Figure 11(C)).
Annual Actual Evapotranspiration
Due to the region’s sparse vegetation, evapotranspiration was largely governed by abiotic factors such as temperature, wind speed, and relative humidity. Analysis of evapotranspiration data revealed that the ratio of evapotranspiration to annual precipitation was approximately 95%, with minimal spatial variation across the study area (Figure 11(F)). Actual annual evapotranspiration values ranged from a maximum of 106 mm to a minimum of 62 mm, correlating with the spatial distribution of average annual precipitation (Figure 11(D)).
Annual Potential Evapotranspiration
To align with the study’s objectives and climatic data, the target period for livestock grazing was defined from the beginning of October to the end of June. Analysis of the data indicated that the annual (October-June) potential evapotranspiration ranged between 640 and 720 mm. The central regions of the study area had the highest potential evapotranspiration values, while the lowest values were observed in the eastern parts and some limited areas in the western part of the study area (Figure 11(E)).
3.6. Land Cover, Land Use, and Livestock
Land Cover and Land Use
The Land Cover, Land Use (LCLU) map (Esri, 2024) provides a standardized classification of land uses within the study area. The majority of the study area is classified as rangeland, making it suitable for agricultural development, grazing, and vegetation enhancement. Conversely, the central region is predominantly characterized by bare ground, indicating less favorable conditions for these activities. While limited built-up areas are also present, LCLU was not considered in the MCDM model for this study. The map serves to provide a foundational understanding of the area (Figure 12(A)).
Modified Soil Adjusted Vegetation Index (MSAVI)
The Modified Soil Adjusted Vegetation Index (MSAVI) reveals that a significant portion of the study area exhibits sparse or seasonal vegetation, with values ranging from 0.05 to 0.13. Nevertheless, the potential for rangeland and agricultural development persists. MSAVI values exceeding 0.1 indicate suitable grazing conditions, provided that adequate water sources are available for livestock. The eastern and central regions offer promising prospects for livestock grazing, while the western regions, characterized by basaltic outcrops, generally lack natural pastures. Consequently, investing in naturally suitable areas in the eastern part of the study area is preferable to investing in less prepared western regions (Figure 12(B)). Areas with an MSAVI of 0.13 were assigned the maximum weight in the MCDM model for this factor.
Distribution of Sheep
An analysis of sheep distribution revealed approximately 77,000 sheep within the study area, relying on natural pastures for fattening cycles and primarily raised for meat production. The highest sheep density, reaching 48 sheep per square kilometer, was concentrated in extensive areas of the eastern and western parts of the study area. Grazing activities were virtually absent in the central region. Areas with the highest sheep density were assigned the maximum weight of 0.015 in the MCDM model (Figure 12(C)).
Figure 12. Land cover land use (A); modified soil adjusted vegetation index (B); distribution of sheep (C) sheep; and distribution of goats (D).
Distribution of Goats
The analysis of goat distribution indicated approximately 7600 goats, primarily raised for milk production. The highest density reached 5 goats per square kilometer, exhibiting a spatial distribution pattern similar to that of sheep. Areas with the highest goat density were assigned the maximum weight of 0.005 in the MCDM model (Figure 12(D)).
4. Discussion
4.1. MCDM Model
Figure 13. Suitability index map.
It is imperative to note that following the construction of the initial MCDM model, a thorough evaluation of the results is necessary. This involves recalibrating the criteria weights to ensure that the resulting suitability map exhibits discernible variations. Subsequently, a random sampling of points with high suitability scores should be conducted. The criteria contributing to the prominence of these points on the map should be meticulously examined. If a particular point is found to be primarily influenced by a limited number of criteria, it suggests that the weight assigned to that criterion may be excessive and requires adjustment. It is crucial to emphasize that the nature of this calibration is inherently constrained by the specific study area. The most suitable areas are determined relative to their immediate surroundings, and their suitability is not absolute. This approach aligns with the principles of sustainable development, which prioritize the utilization of local resources.
![]()
Figure 14. Insolation and topographical index.
Based on the final weights assigned to the criteria in the MCDM model, as discussed previously, a Suitability Index map was generated. This map visually highlights the most suitable locations for constructing additional rainwater harvesting ponds within the study area. Areas exhibiting a Suitability Index value exceeding 0.3 are considered moderately suitable, with suitability increasing as the index value rises. The maximum Suitability Index value attained in the study area is 0.55 (Figure 13).
Insolation and Topographical Index
To pinpoint the most precise locations for establishing rainwater harvesting ponds, the insolation and topographical criteria were integrated into a single map termed the “Insolation and Topographical Index”. This map, with a spatial resolution of 12.5 meters, aims to delineate fine-scale spatial variations and mitigate the impact of criteria exhibiting minor spatial variability, such as climatic factors (Figure 14). The blue-shaded area in Figure 14 exemplifies a region identified as highly suitable.
Proposed Rainwater Harvesting Pond Locations
Figure 15. Existing and proposed rainwater harvesting ponds with their suitability index.
Through a meticulous analysis of the Suitability Index map, the Insolation and Topographical Index map, and the active channel map at low runoff rates, 52 suitable locations for Proposed Rainwater Harvesting Ponds were identified. This represents the maximum number of ponds that can be feasibly constructed within the study area based on the current assessment. Considering the existing 41 Rainwater Harvesting Ponds, the total number of ponds would amount to 93 (Figure 15). Exceeding this threshold could potentially compromise the hydrological feasibility of surface water collection and evaporation mitigation, as well as the financial viability of the project.
Figure 16. Histogram of existing and proposed rainwater harvesting ponds via suitability index.
The Suitability Index values of the Proposed Rainwater Harvesting Ponds indicate a relatively higher suitability compared to the Existing Rainwater Harvesting Ponds. The average suitability of the Proposed Rainwater Harvesting Ponds is 0.37, whereas the Existing Rainwater Harvesting Ponds have an average suitability of 0.32 (Figure 16).
4.2. Critical Discussion and Comparison
The findings of this study in the Hammad Basin gain significant context and strength when critically compared to similar regional and global MCDA-based Rainwater Harvesting (RWH) site identification studies in arid and semi-arid zones. The uniqueness and comparative advantage of our model are primarily demonstrated in two areas: Methodological Innovation in Weighting and Criteria Inclusivity.
Methodological Superiority and Functional Efficiency
Our model’s approach fundamentally deviates from the prevailing methodology. Studies across the Middle East, such as those by Al-Sababhah (2023) in Jordan and Hassan et al. (2025) in Iraq, rely heavily on the Analytic Hierarchy Process (AHP) to assign criteria weights. While AHP ensures mathematical consistency, it is often criticized for injecting subjectivity and yielding generic results. In contrast, our study employs an adaptive multi-objective weighting strategy deliberately tailored to the extreme hydrological constraints of the Hammad Basin. Assigning an overwhelmingly high weight to Channel Density (0.450) ensures absolute priority is given to maximizing potential surface runoff, an essential factor for pond success in a hyper-arid environment. This data-driven, function-first weighting represents a crucial shift toward functional efficiency and robustness, surpassing the generic biophysical prioritization typical of many AHP-based RWH models.
Integration of Socio-Ecological Sustainability
A second critical differentiation lies in our model’s criteria comprehensiveness. Many regional RWH models, like those by Khan et al. (2022) in Pakistan, focus primarily on biophysical factors (e.g., slope, soil, rainfall, stream order) to determine technical feasibility. Our model expands on this by explicitly incorporating rangeland vitality and livestock distribution as a weighted criterion. This is consistent with the urgent call for ecology-inclusive hybrid frameworks for RWH site selection noted in recent systematic reviews (Ahmed et al., 2025). By prioritizing sites that are not only hydrologically viable but also strategically located to support Bedouin livelihoods and efficient rangeland management, the framework elevates the findings from purely technical suitability to addressing socio-economic sustainability. This makes the 52 identified sites more viable and sustainable for local implementation in the Hammad context.
Contextualization and Limitations
The output of 52 highly suitable RWH sites with a standardized design capacity is a significant practical result, aligning with global trends that advocate for integrating geospatial evaluation with risk assessments (Suni et al., 2023). However, a common methodological constraint encountered in regional studies, including our own, is the spatial resolution mismatch (4 km for climate data versus 12.5 m for topography) (Wable et al., 2023). While local calibration was implemented to mitigate this, the findings are best viewed as a robust first-tier selection, requiring detailed field-level hydraulic validation to confirm the precise runoff volume for each of the proposed locations.
5. Conclusion
This study successfully employed a Multi-Criteria Decision Analysis (MCDA) model to identify optimal locations for constructing rainwater harvesting ponds within the Hadalat Dam watershed in Jordan. By integrating diverse geospatial and climatic data, the research effectively assessed the region’s potential for rainwater harvesting, which was found to be substantial due to low soil moisture, high evapotranspiration rates, and limited surface water resources.
The MCDA model pinpointed 52 suitable sites for new ponds, focusing on areas with high suitability indices, adequate drainage, and favorable topographic conditions. These proposed ponds are expected to significantly contribute to enhancing water security, improving rangeland management, and supporting local livelihoods, particularly in the livestock and agricultural sectors.
The study’s findings underscore the importance of a comprehensive, multi-criteria approach for identifying suitable locations for rainwater harvesting ponds in arid and semi-arid regions. The integration of hydrological network analysis, topographic analysis, soil data, climate data, and land cover/land use information provides a robust foundation for decision-making.
It is crucial to emphasize that the identified suitable locations represent the maximum feasible number of ponds based on the current assessment. Implementing these ponds should be approached systematically, with careful monitoring and evaluation to assess their impact on hydrological processes and socio-economic conditions.
Future research could focus on refining the MCDA model by incorporating additional factors, such as socio-economic indicators (for example, distance to settlements, nomadic grazing routes, or land tenure) and water demand assessments. Additionally, detailed hydrological modeling and economic feasibility studies can provide further insights into the long-term sustainability and benefits of rainwater harvesting in the region.
By adopting a phased approach to pond construction and continuously refining the decision-making framework, it is possible to maximize the effectiveness of rainwater harvesting initiatives in addressing water scarcity challenges in the Hadalat Dam watershed and similar regions.