Assessment of Water Erosion in the Houran Valley Using the Gavrilovic Erosion Potential Method and Geomatics Techniques

Considering the vast area of the valley, which cannot contain all its boundaries within a single image, more than one satellite image with spatial accuracy (30 m) was used to cover the boundaries of the study area. The images were merged, and the boundaries of the study area were cut. These boundaries were obtained from the USGS website for the Landsat 8 satellite for the year 2024. All operations were carried out on them using the ArcGIS program.

The research is based on combining data derived from remote sensing with the application of a model. Gavrilovic for estimating the amount of water erosion of soils, through the following.

2.1 Data collection

Obtained Climate data (rainfall, temperature) from the NASA website for the year 2024 to assess their impact on erosion. Soil distribution maps and soil types were obtained from the Iraqi Ministry of Agriculture to assess their erosion susceptibility. Evaluate the impact of human activity, including agricultural expansion, overgrazing, and resource exploitation, on increasing erosion rates in the study area.

2.2 Application of the EPM

This model was developed during the 1950s and modified in subsequent years, and was applied to eroded areas in Italy, Switzerland, and Greece. It demonstrated great accuracy in estimating soil erosion through canals and cover erosion. After comparing the results extracted from the model with field measurements, a high degree of convergence was found in the model, indicating its potential for adoption in estimating soil loss in different regions. This, in turn, reduces effort and economic costs [7].

The EPM relies on remote sensing and GIS techniques to extract spatial values for the parameters affecting water erosion: slope factor, vegetation cover coefficient, soil erosion susceptibility coefficient, and rainfall factor. The process begins by collecting spectral, climatic, and geomorphological data from space and ground sources. Satellite images were then processed to derive slope and vegetation cover maps using spectral indices. In a GIS environment, each parameter is represented as a spatial layer with uniform spatial resolution. Spatial overlays are then performed to integrate these parameters into the model's main equation.

This process results in a quantitative estimate of eroded soil rates (m³/km³/year) for each spatial unit, with erosion severity classified into spatial categories that can be displayed on maps and analytical tables. This enables understanding of the spatial distribution of erosion risks and the development of appropriate management plans to mitigate them. The main equation is applied as follows [7, 8]:

$W=H\times T\times \pi \times \sqrt{{{Z}^{3}}}$          (1)

Since W is the annual rate of erosion (m³/km₂/year), H is the annual rainfall rate (mm), T is the temperature coefficient, $\pi $ is a constant value (3.1415), and $Z$ is the erosion coefficient potential. It is extracted using the following equation [9]:

$T=\sqrt{\frac{C}{10}+0.1}$        (2)

where, C is the annual average temperature.

It is considered that the Z-factor is one of the most important variables on which the Gavrilovic index depends, as it contains indicators that show the continuous and occurring changes in the study area, so it gives realistic results. It allows tracking and monitoring spatial changes in erosion levels over time and land use patterns. The Z-factor is extracted by the following [10]:

$Z=Y\times Xa\times \left( \varphi +\sqrt{Ja} \right)$        (3)

Since Y is the soil erodibility coefficient, Xa is the soil protection factor, $\varphi $ is the erosion development coefficient, and $Ja$ is the slope rate (%).

Erosion levels were classified according to the Z-factor values in Table 1. Gavrilovic highlighted some important indicators that illustrate the most significant factors revealing the changes occurring in the valley as a result of the changing values of these factors and indicators. Satellites were relied upon to estimate the final results of assessing the amount of erosion in the valley. Soil data were used to classify soils into three categories based on the physical and chemical properties and characteristics of the studied area. These factors were extracted from remote sensing data and the ArcGIS program. Then, the areas were classified according to the degree of erosion risk, and water erosion was estimated using the ArcGIS program's outputs.

Table 1. Potential erosion levels for Z-factor values

2.3 Extracting the potential erosion coefficient indicators

2.3.1 Soil erodibility index (Y)

It depends on the physical and chemical properties of the soil (texture, consistency, clay-sand ratio, organic matter content, and permeability). The higher the clay and organic matter content, the higher the soil's ability to resist erosion, while sandy or loose soils with weak cohesion are more susceptible to erosion. Y values are classified according to specific categories that link the degree of susceptibility to texture and cohesion properties.

The values of this index were calculated based on the soil map classified by the Iraqi FAO Company, which is a primary reference for determining the morphological characteristics of soil types in the study area. Laboratory analyses of field samples, including the chemical and physical properties of the soil, were also employed to determine the factors affecting its susceptibility to erosion. Studies indicate that the best apparent density for Iraqi soils is 1.20-1.40 g/cm³ [11, 12]. The soils of the area were classified and divided into three groups (Table 2).

Table 2. Soil chemical and physical properties of the study area

The following equation was applied to obtain the index values [10]:

$Y=\binom{0.00021 \times(12-O M) M 1.14+3.25}{(S-2)+2.5(P-3) / 100} \times 1.58$           (4)

Since $Y$ is the erosion coefficient, $OM$ is the percentage of organic matter, $M$ is the texture = (percentage of silt + fine sand) × (100 - percentage of clay), $S$ is the structure symbol, and $P$ is the permeability coefficient.

The obtained results are shown in Figure 3 and Table 3.

tu_pian_46.png

Figure 3. Soil types in the study area

Table 3. Erosion susceptibility coefficient y for the study area soil

2.3.2 Calcareous loamy soil group

The soils in this group are classified as (Je-10-2/3b) within the soil mapping unit, according to the FAO classification system (FAO, 2015).

Since:

Je = These are calcareous loamy soils that collect in valley bottoms and contain calcareous materials of no less than 50-805,

10 = the presence of linked soil units,

2 = represents the medium soil texture class,

3 = represents the fine soil texture class,

b = represents the simple undulating slope class with a slope degree of (0-8%).

The group of advanced calcareous desert soils and desert Recosols soils (Table 2 and Figure 3).

2.3.3 Soil protection index (Xa)

Xa represents the vegetation cover and its density. The coefficient values were extracted according to the vegetation cover values in the study area, which play a role in the soil erosion process. The plant works to limit the speed of water, which leads to a decrease in the amount of eroded soil if the vegetation is dense and vice versa. It is the rate at which rain falls on the Earth's surface.

Then, the collaboration was carried out in designing the cover based on the data of the Landsat 8 satellite imagery, where the bundle was effectively employed to derive the interactive factor of the robotic researcher. Then, its results were streamed into the EPM to estimate the soil protection index. This study combined remote sensing techniques and model analysis capabilities. It divided the Hauran Valley into five main categories of vegetation orientation ranges (Table 4), which enabled a quantitative spatial description of the extent of the contribution to reducing the percentage of variability in the desert and semi-arid environments.

Table 4. Standards for the soil protection factor (Xa) to Gavrilovic

To accurately determine vegetation cover values, the NDVI equation was applied to extract the Xa index (soil protection index), focusing only on the winter season due to the region's dry continental climate, where rainfall is limited during this season and nonexistent during the summer. The absence of rain in the summer practically eliminates rain erosion, making the analysis of vegetation cover in the winter more realistic for assessing the role of plants in soil protection. linking the NDVI results to the EPM model outputs, showing that areas with high vegetation density recorded low erosion rates due to reduced runoff velocity and the mitigation of raindrops.

Meanwhile, areas with low vegetation cover exhibited high erosion rates, reflecting the close correlation between remote sensing data and the erosion indices calculated by the model. The results of which we can use to extract the Xa index values using the following equation [13]:

$NDVI=\frac{NIR-RED}{NIR+RED}$          (5)

The near-infrared (NIR), the fifth spectral band of the Landsat-8 OLI sensor, is one of the most important spectral indices used in vegetation analysis. The red (R) and near infrared (NIR) bands, typically represented as Band 5 and Band 4 in certain sensors (e.g., Landsat 8), are used to calculate the Normalized Difference Vegetation Index (NDVI). This is one of the most widely used spectral indices for monitoring and assessing vegetation health and density. NDVI values range from +1 to -1; positive values indicate the presence of vegetation cover, and vegetation density increases as the value approaches +1, while. Negative values are often associated with aquatic areas, as water absorbs near-infrared radiation at a greater rate than visible infrared radiation, resulting in negative NDVI values [14, 15]:

$Xa=\left( Xa\,NDVI-0.61 \right)\times \left( -1.25 \right)$         (6)

Xa is the soil protection index, and $XaNDVI$ is the vegetation coverage index adjusted for the soil protection index.

The vegetation cover index was applied to satellite images for the year 2024. Positive values of the index represent dense vegetation cover, while negative values represent barren lands.

The study area was classified into five categories according to the Gavrilovic classification (Table 4). Noticing that the soil protection index Xa is high in areas with dense vegetation cover, which are mostly distributed in the eastern parts (Figure 4), and this is a result of the flatness of the surface. The deposition of soil and its continuous renewal, in addition to the high moisture content of the soil, provide a great opportunity for the germination process.

The coefficient recorded low values in areas devoid of vegetation cover or barren lands, which are usually found in the western parts as a result of the ruggedness of the region, the presence of rocky edges, and the thinness of the soil resulting from severe erosion. From the above, it is clear that the importance of vegetation cover in limiting soil erosion or reducing its severity is evident, as the area has little vegetation cover and weak resistance to erosion processes (Figure 4).

It is clear from the results of Figure 4 and Table 5 that the soil protection index categories varied across the study area. The weak protection category occupied the largest area, covering approximately 11,446.9 km² of the total area. This indicates the prevalence of areas highly susceptible to water erosion and their weak resistance to erosion.

tu_pian_48.png

Figure 4. Xa factor

Table 5. Xa, area, and percentage of each category

In contrast, the very weak protection category ranked lowest, with an area of only 4.2 km², reflecting its limited distribution.

These results indicate that most of the region's lands suffer from low levels of soil protection, underscoring the importance of vegetation cover as a key factor in reducing water erosion and maintaining soil stability.

2.3.4 Slope index (Ja)

The Wadi Hawran basin is characterized by a general slope extending from west to east, with its highest elevation reaching approximately 949 meters above sea level in the western part. In comparison, it gradually decreases to reach its lowest point at 80 meters above sea level in the eastern part, which represents the valley floor and its final outlet.

The topographic characteristics of the basin were analyzed using DEM data, which allowed for accurate slope values to be extracted within the study area. Elevation is a factor that directly affects water erosion processes, as steep slopes accelerate surface water flow, leading to soil fragmentation and its transport to lower areas via floods. The analysis results showed a direct relationship between the severity of the slope and the rate of erosion, with water erosion rates increasing in areas with high slopes, which are concentrated in the western parts, represented by isolated hills and rocky edges. Moreover, a decrease in flat and plain lands, which are distributed in the central and eastern parts, is represented by the oases, alluvial, and fan plains areas. Figure 5 illustrates the spatial distribution of the basin's slopes and their impact on water and soil movement, helping us understand the dynamics of hydrological processes that affect the sustainability of natural resources in the region.

tu_pian_51.jpg

Figure 5. Slope index

From Table 6 and Figure 5, it is clear that the largest area of the valley is controlled by the weak and very weak categories, respectively. These areas represent the valley bottoms and plains, with areas of 8324.18 and 6060.07 km^2, respectively, and a percentage of 46.3% and 33.7%, respectively. The third place was occupied by the medium category, with an area estimated at 2772.43 km² and a percentage of 15.4%, as they represent the low hilly areas that are distributed in most areas of the valley. The fourth and fifth categories represent the areas of intense and very intense highlands, and they occupy an area estimated at 655.61 and 132.67 km², respectively. The percentages of 3.65% and 0.73% for each of them, respectively, as these highlands represent areas of extreme erosion and drift due to height and lack of vegetation. In addition to the soil's ability to quickly disintegrate due to the original materials from which it was formed, there are the dry climatic factors in the region.

Table 6. Regression index area and percentage of each category

2.3.5 Current erosion index ($\varphi $)

This indicator illustrates the effect of land use (agriculture, grazing, urbanization) on soil erosion susceptibility. Exposed agricultural land or overgrazing reduces soil protection, while forests, natural pastures, and farms enhance soil resistance.

The index coefficients are extracted based on the satellite images of the Landsat 8 satellite through the attached table for the image by dividing the square root of the third band by the maximum radiation value (QMAX) as the relationship is directly proportional between radiation and erosion intensity (the radiation ratio increases with the increase in erosion intensity) as in the following equation [14]:

$Q=\frac{\sqrt{TM3}}{QMAX}$          (7)

where, $TM3$ is the third band in space visuals and $QMAX$ is the maximum radiation value.

tu_pian_52.png

Figure 6. Land use

The values of these parameters are obtained from the metadata of the satellite imagery (Figure 6). The Houran Valley region is a remote area in the western plateau of Iraq. Its distance from urban centers contributes to the weakness of organized agricultural activity and the absence of effective land-use management. This emptiness allows for overgrazing, leading to the degradation of natural vegetation cover and increased soil vulnerability to exposure and erosion. Data on the erosion index ($\varphi $), extracted from Landsat 8 imagery, show high levels of radiation in the third band, reflecting the severity of erosion due to surface exposure and poor vegetation cover.

2.3.6 Potential erosion coefficient (Z)

It is considered the most important element in the Gavrilovic model, which was used to estimate the amount of soil eroded by water, as a value can be extracted from a group by combining all the previous indicators into a single, comprehensive mathematical formula that captures water erosion within the basin, the Z values were classified into five main categories to determine the spatial risk level. The results in Table 7 reveal the risk levels in the study area.

Table 7. Z-level categories, area, and percentage of each category

tu_pian_53.png

Figure 7. Potential erosion of the Z factor

The temporal and spatial variation of rainfall has a direct impact on hydrological activity [14, 15]. After the necessary transactions were combined to extract the values of the Z factor, Figure 7 was derived for the potential erosion in the Houran valley, where the smallest area was occupied by the category (very intense), with an area of 871.9 km² and a percentage of 4.8%. The largest area was occupied by the medium category, with an area of 6441.8 km² and a percentage of 35.8%. From the above, we note that all the extracted categories are distributed over all parts of the basin except for the steep category, which is confined to areas with a steep slope, and is represented in the middle parts of the valley, which are devoid of vegetation cover and extend from west to east. These results indicate that the Z-factor is an effective tool for identifying areas most vulnerable to erosion and providing accurate spatial indicators to guide soil management strategies and reduce environmental degradation in the basin [16, 17].

2.4 Application of the EPM model to assess water erosion

As the model requirements must be analyzed, the auxiliary indicators are as in Eq. (1) [7, 8].

2.4.1 Rainfall index (H)

Rainfall is one of the most prominent factors controlling the intensity of water erosion in the Wadi Hauran Basin. Its kinetic energy directly affects the transport of soil from high and sloping areas to lower parts. The scarcity of vegetation exacerbates this process, allowing more soil to be eroded by surface runoff [18]. Given the vast area of the basin and the lack of local climate stations covering all its parts, we relied on NASA data available at http://chrsdata.eng.uci.edu, where spatially distributed points were selected to provide accurate data on rainfall amounts [19]. The results showed that the annual rainfall average during the study period was approximately 125 mm. Figure 8 shows a relatively low average [20]. However, its danger lies in the form of short, high-intensity rainstorms, rather than the length of the rainy season. This rainfall pattern is one of the key inputs into the rainfall coefficient calculations in the EPM model, which is used to estimate erosion rates. The model shows that the intensity of rainstorms, more than the total amount of precipitation, is the decisive factor in determining soil loss rates, thus contributing to explaining the spatial variation in erosion levels in the basin.

tu_pian_54.png

Figure 8. Rain equals lines

2.4.2 Temperature indicator (T)

The temperature index is an important factor for applying Gavrilovic's data, as it affects soils. High temperatures lead to increased evaporation, which leads to a decline in water quantities, which in turn causes the disintegration of soil and its components [15, 16]. High and low temperatures have a significant impact on the disintegration and fracture of rocks through the expansion and contraction of the minerals from which they originated, which facilitates their transport to another location by running floods. To obtain the average temperature during the study period, some equations were followed, as follows:

The numeric values (DN) of each pixel in the thermal infrared were converted to radiance values (DN to radiance) using the following equation [17]:

${{L}_{\lambda }}={{M}_{L}}\times {{Q}_{cal}}+{{A}_{L}}$          (8)

where, ${{L}_{\lambda }}$ is the TOA spectral radiance (watts/(m²⋅sr⋅μm)), M_L is the radiance multiplicative band (No.), ${{A}_{L}}$ is the radiance add band (No.), and ${{Q}_{cal}}$ is the quantized and calibrated standard product pixel values (DN).

The following equation was used to calculate the surface temperature [17]:

$B T=\frac{K_2}{\ln \left(\frac{K_1}{L_\lambda}+1\right)}-273.15$            (9)

where,

$BT$ = Top of atmosphere brightness temperature (℃),

${{L}_{\lambda }}$ = TOA spectral radiance (watts/(m²⋅sr⋅μm)),

K₁ is the constant band (No. 1), which is related to the thermal radiation emitted from the surface, and K₂ is the constant band (No. 2), which is used to adjust the logarithmic relationship between radiation and temperature.

All this data is contained in the metadata of the satellite imagery.

tu_pian_58.png

Figure 9. Annual average of the T factor for the Houran valley

Through Figure 9, it is clear that there is a variation in temperatures from the west of the study area towards the east. This is due to the topography of the area and its lack of vegetation, which contributes to higher temperatures. Additionally, the distance of the western region from water bodies helps moderate the atmosphere and lower the temperature.

2.5 Calculating the internal water volume using the EPM

Model Zachar [20] identified five types of water erosion based on the amount of soil in each, Table 8, to prepare a geographic database for Wadi Houran.

The valley was divided into five categories according to the amount of erosion, as the lowest amount formed by the category (very severe erosion) with an amount of 1210-7300 m³/km²/year, and an area of 306.6 km² with a percentage of 1.7%, while the largest amount formed by the category (weak) with an amount of 573-50 m³/km²/year and an area of 7917.07 km² and a percentage of 44.11% (Figure 10 and Table 8).

Table 8. Volume of lost soil and its percentage according to the EPM model

tu_pian_57.png

Figure 10. Erosion levels according to the amount of drift in the Houran valley