Abstract

Soil erosion is a major hydrogeological hazard in archaeological sites within the East Africa Rift System (EARS). Dating back to the 15^th century CE, and with its distinct stone-lined irrigation fields, Engaruka is one of the impressive abandoned Archaeological Islands of Intensive Agriculture in East Africa (AIIAEA). Aggressive arid climate coupled with socioeconomic factors e.g. overgrazing and farmland encroachments has resulted to extensive gullying that threatens to degrade the northern half of the site. In this study, we use the Unit Stream Power-based Erosion Deposition (USPED) model and UAV-derived (Digital Elevation Model) DEM to estimate the soil erosion-deposition spatial patterns in Engaruka North Terraced Area, ENTA. We then estimated the volume of sediment removed by erosion on an active section of Intermediate North Gorge (Gully Site), ENTA based on 2023 UAV data. We found the mean soil erosion-deposition rates (Mg∙ha⁻¹∙yr⁻¹) to be 33.9 and 31.1 respectively. Intensive erosion processes are dominant along the EARS foot slope while gullying is concentrated in the low-sloping alluvial areas. The eroded volume was measured to be 5936.4 m³ with a total length of 145.7 m of ancient terracing material eroded away on the surveyed gorge. These findings can be useful in planning for effective soil conservation measures within the site as well as future erosion pin measurement methods. This study emphasizes the influence of landscape changes on soil erosion in East Africa’s archaeological areas towards attaining Sustainable Development Goal (SDG) 11.4.1.

Keywords

Soil Erosion, USPED, UAV, ALOS-PALSAR-2, Engaruka, Archaeological Terraces

Share and Cite:

1. Introduction

Engaruka has a large reticulated network of stone-lined farm fields (c. 2 Km²) and agricultural hillside stone terraces (c. 9 Km²) making it arguably the most prominent site among the Archaeological Islands of Intensive Agriculture in East Africa (AIIAEA) (Davies, 2013; Sassoon, 1966, 1967; Stump, 2006; Sutton, 1978, 1986; Widgren & Sutton 2004). These AIIAEA sites are regarded as islands due to their remoteness and the labor-intensive agricultural techniques employed by the original inhabitants (Widgren & Sutton 2004). Several stratigraphic studies have demonstrated that the agricultural terraces on these sites are not terraces intrinsically but are instead sediment traps (Ferro-Vázquez et al., 2017; Lang & Stump, 2017; Stump, 2006). For instance, in Konso Ethiopia, Ferro-Vázquez et al. (Ferro-Vázquez et al., 2017) showed how the inhabitants would construct low-level dry-stone terraces in order to trap tillable top soil as well as mitigate soil erosion. They would then seasonally raise the walls of the terraces to improve efficiency. In Engaruka, after abandonment about 18th century CE followed by years of no or little maintenance efforts, most of these terraces have now been covered by soil (near-surface) due to denudation. This, combined with other factors like aridity, flash-floods, farmland encroachment and overgrazing has resulted to extensive gullying and severe land degradation (Figure 1 and Figure 2) particularly on the northern half of the site (Engaruka North). Westerberg and Hagberg (2018) demonstrated how a single high-magnitude process, the 1997-1998 El Niño denuded approximately 80 ha of the site (c. 6.6 mm denudation) and led to the rapid evolutions of gullies in Engaruka North (EN) (Westerberg & Hagberg, 2018), see Figure 2. Kabora et al. (2020) employed a medium resolution (30 m) global digital elevation model, DEM (ASTER DEM) in an agent-based model (Engaruka Sediment Transport and Trapping, ESTTraP) and estimated a denudation of approximately 3.5 mm for an area covering 3000 m² in Engaruka North, see Figure 2. These denudations translate to 108 & 42 Mg∙ha⁻¹∙yr¹ respectively assuming a 1.2 g/cm³ soil bulk density (Miller et al., 2021; Stocking, 1984). Stump (Stump, 2006) stratigraphically estimated soil accumulation of over 2 m deep in the southern end of Engaruka over a period of 300 years. Kabora and Katie (Kabora & Campbell, 2022) further utilized a Structure for Motion (SfM) approach on archival photos and UAV images over the central section of the site (Engaruka Central) and detailed its usefulness cultural heritage management e.g. reconstruction of earlier landscapes especially in such a data poor region. They acknowledged the low accuracy of the 3D products from the archival data due to the lack of camera calibration metadata and scan errors. However, these studies neither quantified the spatial pattern of lateral movement of soil (amount and distribution of soil erosion/deposition) nor the eroded volumes within the gully channels and gorges in this archaeological area.

Figure 1. (a) A 1958 historical aerial image (British Institute in Eastern Africa (BIEA)) (Stump, 2008; Sutton, 1978, 1986) showing the ancient stone-terraced landscape of Engaruka North Terraced Area (ENTA). The terraces and the farm plots appear as soil marks (light and dark linear variations); (b) A 2025 Google Earth image of the same extent as Figure 1(a) showing farmland encroachments onto the terraced site. Other visible landscape changes include settlement (the circular Maasai nomadic houses, Manyattas), gullies and transportation tracks and trails.

Figure 2. (a) A 1958 archive image (British Institute in Eastern Africa (BIEA)) (Stump, 2008; Sutton, 1978, 1986) depicting minimal gullying. The red ellipse indicates the location Figure 2(b) was taken. The blue polygon gives the ESTTraP site location Kabora et al. (2020) (b) A 2023 on-field picture by the author showing how land degradation is affecting a stone-lined irrigation channels (Stump, 2006; Westerberg & Hagberg, 2018). Notice the savannah landscape and the towering rift escarpment in the background.

A recent review of soil erosion studies in cultural areas indicates that quantifying erosion in archaeological areas is a complex task due to the complexity of erosion processes hence the very low number of publications in the past three decades (Papageorgiou et al., 2025). In recognition of this critical need, international bodies like the United Nations have laid out goals (UNSDG 11.4.1) to promote conservation (Papageorgiou et al., 2025). Among the methods summarized in (Papageorgiou et al., 2025) is Mitasova et al.’s (Mitasova et al., 1996; Mitasova & Mitas, 2001; Mitasova et al., 2001) Unit Stream Power-based Erosion Deposition (USPED) model which is a modification of the commonly used Revised Universal Soil Loss Equation (RUSLE). USPED parameterizes terrain curvature to map spatial pattern of erosion and deposition of a complex topography. It measures lateral movements of soil under sheet, rill and/or gully scenarios (Borrelli et al., 2015; Hadush et al., 2025; Honek et al., 2020). Notable works of this model in natural & cultural heritage sites include (Brandolini et al., 2023; Frisetti et al., 2026; Minervino Amodio et al., 2024; Minervino Amodio et al., 2023; Pagels et al., 2024). Some researchers e.g. (Gioia et al., 2021) have further refined the outputs of this model by experimenting with VHR UAV-derived (DEMs). Similar to the RUSLE model, USPED is limited in that it doesn’t estimate gully sediment volumes. Geomorphological studies (non-archaeological) by Frankl et al. and Zhuang et al. (Frankl et al., 2012; Frankl et al., 2015; Zhuang et al., 2024) not only showed the value of 3D UAV data in modelling gullies within East Africa Rift System (EARS) but also demonstrated how aerial photographs from the colonial period can be of use in such a data-poor region. Jorayev et al. (Jorayev et al., 2016) computed Olduvai’s high erosion volume (585 million m³) in the adjoining Ngorongoro basin using UAV and highlighted the value of gorge morphometric data in heritage conservation. The aims of the study are: 1) to determine the soil erosion-deposition rates and sediment routing in Engaruka North Terraced Area (ENTA) based on the USPED empirical model and VHR UAV-DEM and 2) to complement the USPED results with 3D volumetric sediment analysis on an active section of the Intermediate North Gorge, (Gully site), ENTA.

2. Materials and Methods

2.1. Study Area

Engaruka basin in which the Engaruka archaeological site is located, is located in the tropical savannah of EARS, North Tanzania (Ngorongoro-Lengai UNESCO Global Geopark (Pérez-Romero et al., 2023)) between latitudes 2˚53'S and 3˚05'S and longitudes 35˚49'E and 36˚05'E (Figure 3). We used the r.flow and r.watershed toolsets in GRASS (Geographic Resources Analysis Support System) GIS 7.8 software (Mitasova, 2014) to delineate the basin (Figure 4) focusing on the hydrology of the Engaruka site i.e. Loljoro, Makuyuni, Intermediate North Gorge, Engaruka and Olemelepo (Figure 4(b)) that supported agriculture until the 18^th CE. The steep western escarpment of the Rift Valley bisects the basin nearly into two regions; the towering western side (towards Ngorongoro Crater, c. 1000 m) and the eastern lowlands (towards Engaruka basin lake, c. 500 m). The Engaruka archaeological site is divided into three sections (Figure 4(a)) based on the type of dominant ancient remains and structures. Figure 5 depicts the category of dry-stone archaeological features found within the site (Sassoon, 1967; Sutton, 1978). Engaruka North (EN) (Figure 4(a)) has stone circles (Figure 5(c)), terraced fields (Figure 5(d)) and stone-lined irrigation channels (Figure 5(e)). Engaruka Central (EC) has tilled fields, cairns and rock piles (Figure 5(b)) while Engaruka South (ES) has an extensive reticulated network of stone-bordered farm plots (Figure 5(a)) and irrigation channels. In the pedological survey conducted by (Westerberg & Hagberg, 2018), the bottom of the Rift has predominantly torripsamment soils. Run-off from the undifferentiated basalt in the steep areas heightens surface erosion on the volcanic black and calcareous soils in the lowlands (Figure 4(a)). The area has a semi-arid, tropical climate with a bimodal precipitation pattern (March-May long rains and October-December short rains). The area records an annual rainfall of c. 460 mm (Tanzania Meteorological Agency). This study focused on the Engaruka North Terraced Area (ENTA), EN which is defined by the escarpment foot slopes, the southern EN-ES transition boundary and the extent of visible traces of terraces on the north and east directions, see Figure 4. ENTA (c. 6.82 Km²) was selected due to its susceptibility to gullying, socioeconomic land uses & level of soil erosion-deposition on the ancient terraces.

Figure 3. (a) Location of notable AIIAEA sites; (b) the geological setting of Engaruka archaeological site within the Engaruka basin, Ngorongoro-Lengai UNESCO Global Geopark.

Figure 4. The hydrogeological setting of Engaruka basin and location of Engaruka North Terraced Area (ENTA) within Engaruka North section. (a) Engaruka basin lithology (b) Sub-basins with Engaruka Basin.

Figure 5. Categorization of stone features within Engaruka archaeological site (Sassoon, 1967): (a) stone-bound farm plots (averagely 6 by 4 m) and irrigation furrows in Engaruka South (ES); (b) rock piles in Engaruka Central (EC); (c) stone-circles in EN; (d) irrigation terraces in EN; (e) eroded and exposed irrigation channels in EN. Notice the dry arid soils, susceptibility to both dry and wet bank erosions, the subsurface nature of the stone feature. (f) cairns in EC (both human figures in the picture are 1.78 m tall).

2.2. Data Sources

Table 1 shown below summarizes the specifications of the datasets applied in this study. Figure 6 to Figure 10 gives a pictorial representation of the data sources, acquisition equipment and preparation stages.

Table 1. Datasets in use for this study.

Figure 6. Examples of the photogrammetric datasets used in this study: (a) Oldonyo-Lengai Topographical sheet Map; (b) a 1958 aerial image. Notice the terraced EN landscape (covered terraces appearing as linear soil marks); (c) UAV (eBee Plus) orthomosaic for EN. Notice the intense evolution of the gully network between 1958 and 2023; (d) UAV (DJI P4P) orthomosaic for the Intermediate North Gorge (Gully site).

Figure 7. Example UAV data acquisition processes at Engaruka North Terraced Area (ENTA): (a) DJI P4P drone; (b) the same eBee Plus drone used by Tanzania Flying Labs in Engaruka, https://ag.dji.com/newsroom/pilot-story-leka-tingitana; (c) setting up the RTK GPS land survey equipment; (d) measurement of GCPs.

Figure 8. Evolution and characterization of the Intermediate North Gorge (Gully site), EN: (a) a zoom-in of the 1958 aerial image showing very minimal or no gullying; (b) a 2023 UAV (DJI) orthomosaic of the same extent as shown in (a) indicating dramatic growth of the gully; (c) an on-field picture indicating the approximately 1.5 m deep gully bank; (d) the U-shape floor of the active gully; (e) a measured GCP mark for drone survey (f) setting up GCP on the confluence of the main and side branches of the Intermediate Gorge.

Both the 1960 Topographical Sheet Map (Oldonyo-Lengai Sheet Number 39/4, 1:50,000) and the 1958 historical aerial image scans were provided by the British Institute in Eastern Africa (BIEA), Nairobi Kenya (Stump, 2008). The scan image scene which has a ground coverage of 4 km by 4 km at a scale of 1:14,000, were used as a visual reference for the Intermediate North Gorge (Gully site) area in the 1950s (to show gully evolution). The image scan was provided without important metadata like camera calibrations and orientations. Our check with the Tanzania Air Survey Department for the original copies was unfruitful. As and both Kabora and Katie (Kabora & Campbell, 2022) and observed, georeferencing these archival photos is a challenging task due to significant landscape changes. We used features like stone circles, cairns, visible terrace lines and hydrological features (e.g. points along the initial form of the Intermediate North Gorge (IG)) to georeferenced the image scan in ArcGIS Pro. Authors (Kabora & Campbell, 2022) and (Frankl et al., 2015) explained the complexity in attaining high height accuracies from such data using Agisoft’s workflow in the absence of camera calibrations. We contacted the original air service provider, National Collection of Aerial Photography, NCAP (NCAP_DOS_35/TN/0003_33. All rights reserved) who confirmed the unavailability of films plus the camera parameters. They however provided another 1960 1:20,000 scale (Figure 9) collaborative image that is among the oldest remote sensing products in Tanzania used in the production of the said 1958 Topographical Sheet Map. These historic data were used only as a visual interpretation of the Intermediate North Gorge (Gully site) evolution since the lack of camera metadata disqualified their use in a true multitemporal 3D volume change analysis.

The eBee Plus and DJI P4P drone data was acquired by Tanzania Flying Labs (TFL) in the May, 2023 field campaign, Figure 10. Table 2 gives the main technical specifications and accuracies of the UAV datasets obtained following conventional drone workflows in Pix4D software. Figure 11 shows an example of a 3D view product from the DJI P4P drone processing for section of the Intermediate North Gorge (Gully site) showing how transportation tracks influence development of rills on the site. The intensity of land degradation and erosion between this 1958-2023 timescale is very clearly noticeable (Figure 8). Due to prevalent water erosion, much of the ancient agricultural devices e.g. stone-terraces have been partially covered by topsoil and appear as linear soil marks on the black and white photographs, Figure 6(b).

Figure 9. (a) A 1960, 1:20,000 archive image (National Collection of Aerial Photography, NCAP_DOS_35/TN/0003_33. All rights reserved) (b) an inset indicating the Intermediate Gorge (preerosion). The red polygon indicates the Gully Site area, see Figure 6(b) & Figure 6(c).

Figure 10. Distribution of GCPs for the UAV surveys: (a) eBee Plus survey; (b) DJI P4P survey; (c) a zoom of the DJI P4P orthomosaic showing bank erosion (d) the U-shaped floor of the gully.

Figure 11. A 3D surface model of a section of the Intermediate Northern Gorge (Gully Site), see Figure 10(b) for this inset location. The model has a 1.5 vertical exaggeration. Notice the how the trail (mostly used by motor cycles) has impacted on the gully banks.

Table 2. Ground sample distance (GSD); ground control point (GCPs); root mean square error (RMSE) and other specifications for the UAVs in use.

2.3. USPED Model

The Unit Stream Power Erosion and Deposition (USPED) (Mitasova et al., 2001) model is an extended variation of the commonly used empirical Revised Universal Soil Loss Equation (RUSLE) (Renard, 1995; Wischmeier & Smith, 1978):

$T=R\ast K\ast LSt\ast C\ast P$ (1)

where, T (Mg∙ha⁻¹∙yr⁻¹) is the sediment flow at transport capacity, R (MJ∙mm∙h⁻¹∙ha⁻¹∙yr⁻¹) is the rainfall intensity factor, K (Mg∙h∙MJ⁻¹∙mm⁻¹) is the soil erodibility factor, LSt (dimensionless) is the topographic sediment transport factor, C (dimensionless) is the land cover factor and P (dimensionless) is the soil conservation factor. Unlike the 1-dimensional RUSLE model, the USPED model computes the average soil loss taking into account the lateral movement (i.e. both in the x and y directions) of top soil. The unique terrain factor LSt-factor in the USPED model parameterizes the profile as well as the tangential curvature of a complex terrain (Minervino Amodio et al., 2023; Pagels et al., 2024).

$LSt=\left(m+1\right){\left(Ac/\left(\alpha 0\right)\right)}^m{\left(\sin \beta 0/0.0896\right)}^n.$ (2)

where, LSt is the topographic factor, Ac is the upslope contributing area per unit of width (m²/m), α0 is the is the length (72.6 ft, equal to 22.13 m) of the standard terrain, β0 is the angle (9%, equal to a 5.16 slope degree) of the standard terrain, m and n (1.3 and 1.2) are empirical exponents. According to Mitasova et al. (Mitasova et al., 1996), m = 1 and n = 1 implies that sheet erosion is prevalent while m = 1.6 and n = 1.3 indicates that rill erosion is dominant. For this case we considered a combination of both sheet and rill erosion and adopted the model’s recommended m = 1.3 and n = 1.2 (Hadush et al., 2025). Total lateral movement of top soil (erosion and deposition, ED) is derived from the divergence of sediment flow (sediment transport capacity) based on the equation below (Mitasova et al., 1996):

$ED=\partial \left[\left(T\cos \alpha \right)/\partial x\right]+\partial \left[\left(T\sin \alpha \right)/\partial y\right]$ . (3)

where, ED is soil erosion or deposition (t/ha/year), α is the topographical aspect or the or the direction of flow measured in degrees. The terrain derivatives were obtained from the 0.5 m eBee Plus UAV DEM for the ENTA. Following after Moore (Moore, 1979) who observed a direct correlation between kinetic energy (KE) of storms and annual rainfall in East Africa, the Moore’s regression equation was applied to compute rainfall erosivity as follows (Moore, 1979; Wynants et al., 2018):

$KE=3.96\ast MAR+31.22$ . (4)

$R=17.02\left(0.029\right)\ast KE-26.0$ . (5)

where, MAR (mm) is the mean annual precipitation in and R (MJ∙mm∙h⁻¹∙ha⁻¹∙yr⁻¹) is the rainfall erosivity factor. As per the Tanzania Meteorological Authority, the only two weather stations (Karatu-Kilimo and Monduli) are averagely 10 Km away from the vicinity of the study site. The WAPOR (Hajirad et al., 2023) provides precipitation dataset at between 20 to 250 m spatial resolution for selected countries and hydro-basins. Unfortunately, Tanzania is excluded. Due to insufficient observed data, the long-term MAR was sourced from the kilometer-scale CHELSA-monthly climate dataset (Karger et al., 2017) which is among global datasets with the highest resolution. Monthly precipitation (January-December) was averaged over a period of 31 years (1990-2021) to account for interannual variability (Wynants et al., 2018).

$K=Fcsand\ast Fsicl\ast Forgc\ast Fhisand\ast 0.1317$ . (6)

$Fcsand=\left[0.2+0.3\exp \left(-0.0256SAN\left(1-SIL/100\right)\right)\right]$ .(7)

$Fsicl=\left[SIL/CIA+SIL\right]\wedge 0.3$ . (8)

$Forgc=\left[1-0.0256C/\left(C+EXP\left(3.72-2.95C\right)\right)\right]$ . (9)

$Fhisand=\left[1-0.70SN1/\left(SN1+\exp \left(-5.51+22.9SN1\right)\right)\right]$ . (10)

where SAN, SIL and CLA are the percentage of sand, silt and clay content, respectively; C is the organic carbon content and SN1 is the sand content subtracted from 1 and divided by 100. Fcsand gives a low soil erodibility factor for soil with coarse sand and a high value for soil with little sand content. Fsicl gives a low soil erodibility factor with a high clay to silt ratio; Forgc is the factor that reduces soil erodibility for soil with high organic contents. Fhisand is the factor that reduces soil erodibility for soil with extremely high sand content. WorldView-2 imagery of 0.5 m resolution was used to calculate the protective nature of vegetation cover (C) against erosion for Engaruka North catchment. The following equation was applied (Hategekimana et al., 2020):

$C=\exp \left(-2\ast NDVI/\left(\left(1-NDVI\right)\right)\right)$ . (11)

$NDVI=\left(\left(NIR-RED\right)\right)/\left(\left(NIR+RED\right)\right)$ . (12)

where, C is vegetation cover factor (dimensionless) and NDVI (dimensionless) is the Normalized Difference Vegetation Index derived from NIR and RED are the Near-Infrared and Red bands respectively.

A (P) factor of 1.0 (Wischmeier & Smith, 1978; Wynants et al., 2018) was applied as most of the farm support practices (ancient terraces) have since been degraded or covered by topsoil and debris (Kabora et al., 2020; Westerberg & Hagberg, 2018). Similar to other RUSLE/USPED works dealing with multiple datasets with varying spatial resolutions (Bagwan & Gavali, 2024; Balasubramani et al., 2015; Zakerinejad & Maerker, 2015), the Inverse Distance Weighted technique was used to bring all the above parameters to a uniform resolution of 5 m. All factors, see Figure 12 were georeferenced in the WGS84 ellipsoid/UTM zone 36 S projection system.

Figure 12. ENTA USPED paramaters (a) R factor (b) K factor (c) LS factor (d) C factor.

3. Results

3.1. Soil Erosion Deposition Rates for Engaruka North Terraced Area

The erosion-deposition (E-D) results were reclassified in the same number of classes as with (Gioia et al., 2021; Minervino Amodio et al., 2023). Figure 13 gives the spatial distribution of Soil Erosion Deposition (E-D) in ENTA. The mean rate of soil erosion was 33.9 Mg∙ha⁻¹∙yr⁻¹ while the average rate of soil deposition was also 31.1 Mg∙ha⁻¹∙yr⁻¹.

There is a pronounced pattern of erosion and deposition along the hydrological networks e.g. the gully systems in steeper relief areas of the foot slopes of the western EARS. The gently sloping lowland regions (alluvial plains) of the site have noticeably higher patterns of deposition due to low transport capacities. It’s noteworthy most of the stone circles in ENTA catchment area are in close proximity to gully networks, see Figures 14-19 as representative examples. Similar to (Gioia et al., 2021), we also noted increased gullying and E-D patterns on farmed area (land use), (Figure 17 & Figure 19). UAV gully volume analysis complemented USPED empirical model since it doesn’t predict volume of sediment eroded within gullies.

Figure 13. USPED model for ENTA. Map gives predicted annual soil loss in Mg∙ha⁻¹∙yr⁻¹ categorised as Erosion, Stable and Deposition regions. Intensive erosion processes are dominant along the EARS foot slopes while gullying is concentrated in the low-sloping alluvial areas.

Figure 14. Details of the inset E, see Figure 13; (a) a 1958 aerial photo showing minimal gully network in Engaruka North Terraced Area (before the 1998 El-Nino). The red ellipse indicates the same area the photo on Figure 2(b) was taken. (b) a 2023 eBee orthophoto for the same extent indicating intensity of gully evolution (c) soil erosion-deposition spatial pattern for the same area (yellow polygon) as mapped (overlaid) by USPED model. Notice high erosion pattern along the gully channels and human tracks.

Figure 15. Details of the inset F, see Figure 13; (a) a 1958 aerial photo showing minimal gully network in Engaruka North Terraced Area (before the 1998 El-Nino). The pink polygons indicate archaeological stone circles situated on a series of ancient stone-terraces. (b) a 2023 eBee orthophoto for the same extent indicating gully network evolution. The large circular visible features are non-archaeological features (Maasai Manyattas). (c) soil erosion-deposition spatial pattern for the same area as mapped (overlaid) by USPED model. Notice high erosion pattern along the gully channels and human tracks.

Figure 16. Details of the inset K, see Figure 13; (a) a 1958 aerial photo showing minimal gully network and terraced pattern in Engaruka North Terraced Area (before the 1998 El-Nino). The pink polygon indicates disintegrated stone circle (b) a 2023 eBee orthophoto for the same extent indicating gully network evolution (Intermediate North Gorge). Notice the proximity of the disintegrated stone circle to the gully network. The large circular visible features are non-archaeological features (Maasai Manyattas) (c) soil erosion-deposition spatial pattern for the same area as mapped (overlaid) by USPED model. The black polygon indicates the gully site for the DJI P4P drone survey. Notice high erosion pattern along the gully channels and human tracks. This figure also highlights the complementarity of UAV gully-volume analysis and USPED model.

Figure 17. Details of the inset G, see Figure 13; (a) a 1958 aerial photo showing minimal gully network and terraced pattern in Engaruka North Terraced Area (before the 1998 El-Nino). (b) a 2023 eBee orthophoto for the same extent indicating gully network evolution and encroaching farmlands (c) soil erosion-deposition spatial pattern for the same area as mapped (overlaid) by USPED model. Notice high erosion pattern along the gully channels, human tracks and the tilled areas.

Figure 18. Details of the inset L, see Figure 13; (a) a 1958 aerial photo showing minimal gully network and terraced pattern in Engaruka North area (before the 1998 El-Nino). (b) a 2023 eBee orthophoto for the same extent indicating gully network evolution and encroaching farmlands (c) soil erosion-deposition spatial pattern for the same area as mapped (overlaid) by USPED model. Notice high erosion pattern along the gully channels, human tracks and the tilled areas.

Figure 19. Details of the inset H, see Figure 13; (a) a 1958 aerial photo showing minimal gully network and cluster of stone circles along the EARS footslopes in Engaruka North Terraced Area (before the 1998 El-Nino). (b) a 2023 eBee orthophoto for the same extent indicating gully network development and trails (c) soil erosion-deposition spatial pattern for the same area as mapped (overlaid) by USPED model. Notice the increased erosion pattern due to the hillslope location (on the foot slopes) of cluster of stone circles.

3.2. Intermediate North Gorge (Gully Site) Volumetric Analysis

Figure 20. Topographic detail of the Gully Site area; (a) Elevation (b) and (c) vertical difference between UAV-DTM surface and RTK-measured surface for cross-sectional profiles M’-M and N-N’ respectively, see Figure 10, (d) the longitudinal profile P’-P of the Gully Site floor.

The P4P UAV survey covered an area of 42954.7 m² (Gully Site) of the Intermediate North Gorge. The elevation of this Gully Site area ranged between 915.18 m and 905.14 m above sea level (a.s.l.) (Figure 20) with an overall mean slope of 9.33⁰. RTK field survey points obtained along two profiles across the gully were used to check the errors associated with the UAV-derived bare earth model (DTM), see Figure 10 and Figure 20. The overall average vertical difference in the RMSE between the field points and the UAV-DTM was 0.11 m. The longitudinal profile for the main branch of the gully (the longest branch) had a total length of 303.4 m and an elevation range of 905.14 m to 914.64 m, with a mean 1:33.7 mean gradient. The gully area generally had low vegetation, although some the shadowed gully walls negatively impacted the photogrammetric process. The DTM was clipped along the gully rim so as to constrain the 3D volumetric analysis only on the gullied section only, see Figure 20. Gully morphometric analyses (slope, Terrain Wetness Index, depth and volume eroded) were performed using SAGA GIS, see Figure 21. In addition, as depicted on Figure 21(b), a total length of approximately 145.7 m of stone terracing material, affecting about 12 ancient terraces was cut (eroded away) by the gully between the years 1958 and 2023. The gully banks had the highest slopes of approximately 76.4⁰. Following after (Jorayev et al., 2016), the gully depth raster was formed by subtracting the elevation of the gully from a simulated flat surface. A volume map was then created by as a product of the depth values by cell area, giving the volume of sediment eroded per cell. As shown on Figure 21(a), the maximum depth of the gully was 2.64 m that translates to a maximum volume of 0.059 m³ per cell. The UAV model imagery estimated erosion of ap-proximately 5936.4 m³ of sediment.

Figure 21. (a) Slope, Terrain Wetness Index, TWI, Depth & Volume Eroded for the Gully Site area (b) a mapping of the affected stone terraces (red dash lines) and the eroded lengths (sections highlighted in blue color).

4. Discussion and Concluding Remarks

We conducted ground visits as well as interpretation of the UAV very high-resolution images to check the validity the resultant USPED erosion-deposition (see Figure 22). Due to lack of reliable and valid field gauging stations, we compared our findings with other soil erosion studies in and around the Engaruka basin (Zakerinejad & Maerker, 2015). The area’s erosion rate is slightly above the tolerable limits (25 t∙ha⁻¹∙yr⁻¹) for highland areas (Koirala et al., 2019). Wynant et al. (Wynants et al., 2018) utilized the RUSLE model in the neighboring Lake Manyara basin and mapped erosion hotspots dynamics without stating the erosion rates therein. Wynant et al. (Wynants et al., 2020) also used geochemical fingerprinting/source-tracing technique in the Lake Manyara basin and found sedimentation rates to be between the range of 80 - 181 t∙ha⁻¹∙yr⁻¹.

Figure 22. (a), (b), (c) and (d) on-ground validation pictures for the ENTA USPED model focusing on intensive erosion and deposition areas (Pagels et al., 2024).

Many other researchers have reported such erosion rates within the East Africa Rift Valley System (Alemu et al., 2025; Fenta et al., 2020; Kahsay et al., 2025; Watene et al., 2021; Wolka et al., 2015). Subhatu et al. (Subhatu et al., 2017) applied the RUSLE model in the terraced croplands of Minchet catchment, Ethiopia and found an annual soil loss of 32 - 37 t∙ha⁻¹. In the neighboring Meatu Distirct, Simiyu region, Stocking (Stocking, 1984) used low-level aerial photography and field survey with Universal Soil Loss Equation (USLE) and estimated erosion and deposition rates to 61.6 t∙ha⁻¹∙yr⁻¹ and 36.9 t∙ha⁻¹∙yr⁻¹ respectively. Defersha et al. (Defersha et al., 2012) reported erosion rates of 120 t∙ha⁻¹∙yr⁻¹ in the Mara River, Serengeti using Erosion-3D and Water Erosion Prediction Project (WEPP) models; Egbert (Egberts, 2020) reported sediment rates of 15.7 t∙ha⁻¹∙yr⁻¹ in Lake Manyara basin using Soil and Water Assessment Tool (SWAT); Egbert (Egberts, 2020) further used field erosion pins to study channel bed erosion in selected gullies within the Lake Manyarra national park and found maximum erosion pin measurements as 11.9 cm and 10.5 cm for erosion and deposition respectively. Frankl et al. (Frankl et al., 2012) in Ethiopia and Zhuang et al. (Zhuang et al., 2024) also reported rapid evolution of gully systems within the Kenyan EARS based on historical images and UAV data. Other soil accumulation studies within the Engaruka Basin include Kabora et al. (Kabora et al., 2020) (42 t∙ha⁻¹∙yr⁻¹) and Westerberg and Hagberg (Westerberg & Hagberg, 2018) (108 t∙ha⁻¹∙yr⁻¹). Our model results depicted high linear erosion pattern along linear features e.g. gullies, rills and road networks (see Figure 14 to Figure 19) as similar to other USPED studies (Gioia et al., 2021; Minervino Amodio et al., 2023, 2024; Pagels, et al., 2024). Moreover, Figure 17 and Figure 18 show increased pattern of erosion and deposition due to land tilling which is in line with work of (Gioia et al., 2021) that evaluated soil erosion and historical land use with USPED. In addition, Figure 14 indicates a more stable pattern in comparisons Figure 17 to Figure 18 due to the protective nature of vegetation i.e. trees. Gioia et al. (Gioia et al., 2021) further suggested that linear physical features like roads and trails can reduce surface erosion by allowing water to flow in a channeled manner. While this is true, such channeled flow, especially on arid, hilly areas like Engaruka, can result to intensive gullying if there are control measures e.g. gabions and check dams. Some of the limitations of the USPED model is use include its inability to estimate short-term soil losses e.g. amounts for a single El-Nino event plus inability to model complex physical process of soil erosion (Alewell et al., 2019). Therefore, there is need to plan for physical based models e.g. use or erosion pin measurements so as to improve such results and comprehensively understand erosion in the area. Use of equipment like terrestrial 3D laser scanners can further quantify detachment of stone features on the banks of gullies within the site (Alexiou et al., 2024) e.g. monitoring the exposed irrigation canal shown on Figure 2(b). Another recommendation is to incorporate field-scale techniques that can monitor the impact of trail erosion on the stone features in Engaruka (Jula & Voiculescu, 2025).

Westerberg and Hagberg (Westerberg & Hagberg, 2018) demonstrated how both high and low magnitude geomorphic processes are main the drivers of soil erosion at Engaruka. Figure 23 shows some of the anthropogenic factors (overgrazing, farms & settlements) that further exacerbate this hazard. It’s worth mentioning that as depicted by the Google Earth images on Figure 23, the local inhabitants have encroached on the ancient terraced area in Engaruka North by approximately 300m between the years 2017 and 2025. Poor farming practice (little soil conservation measures e.g. no re-terracing) is evident on the Google Earth imagery. Studies like that of Wolka et al. (Wolka et al., 2015) demonstrated how simple soil conservation measures e.g. terraces and bunds can reduce soil erosion rates by up 92%. Since gullying is now a severe problem that threatens to destroy the northern half of Engaruka site, it would be prudent for the concerned authorities (e.g. the Ngorongoro Conservation Authority), particularly in the interest of anticipated land use and climate change in the Serengeti region (Kariuki et al., 2022; Magang et al., 2024; Marchant et al., 2018; Mwabumba et al., 2022), to conduct thorough sensitization programs to educate the local community about heritage management.

Figure 23. Anthropogenic (socioeconominc) drivers of soil erosion in Engaruka archaeological site (a) overgrazing, Photo taken by Lars-Ove Westerberg (b) construction of Maasai settlements (Manyattas) resulting to structural damage of the stone bordered farm plots and terraces (c) 2017 Google Earth image showing farmland encroaching on to the terraced areas (notice the linear crop marks) (d) current (2025 Google Earth image) farmland boundary has shifted as indicated by the red arrows.

Soil erosion and deposition in Engaruka have led to some of the ancient remains being covered by topsoil and debris. Figure 24 shows more examples of some of surface and near-surface stone features on the site. Space radar has been shown to map subsurface archaeological features especially in homogeneous sand within desert areas (Chen et al., 2026; Chen et al., 2018; Guo et al., 2024). Although Engaruka doesn’t present such unique qualities favorable for radar penetration, for experimental purposes, we acquired two scenes of the second-generation, Phased Array type L-band Synthetic Aperture Radar (ALOS PALSAR-2) Spotlight (SL) mode data to test detectability of both surface and near-surface dry-stone features. The two PALSAR-2 SL datasets were acquired in the month of October which falls within the dry season for semi-arid Engaruka and in different specifications (incidence angle and direction). Similar to Balz et al. (Balz et al., 2016), Stewart et al. (Stewart et al., 2016), and Cigna and Tapete et al. (Cigna & Tapete, 2024), we hypothesized that detectability of eroded dry-stone archaeological features (e.g. >2 m diameter, the linear stone terrace and reticulated stone-bordered patterns) in the VHR radar SL data will be high given the low biomass or vegetation during the dry conditions. We used SNAP software to process the Single Look Complex (SLC) data by radiometrically calibrating and geocoding each scene to map geometry. Backscatter images were interpreted alongside optical imagery, including archival photos, WorldView-2, UAV orthomosaics and Google Earth as well ground survey data. We avoided multi-looking and also experimented with unfiltered images in order to preserve the SAR spatial resolution (Balz et al., 2016). Due to high backscattering (speckling effect), none of these stone circles were visible or detectable on the processed radar images. However, a keen eye can pick out the faint linear pattern of stone terraces in Engaruka North (see Figure 25). Figure 26 gives an on-ground view of these stone terraces as well their appearance on optical imagery (linear pattern of soil and crop marks). Similar to (Balz et al., 2016), the linear stone terraces pattern was more detectable on the SL scene with lower incidence angle (34.2⁰). Nevertheless, more radar scenes with varying radar geometry would be needed for conclusive experimentation (Chen et al., 2018).

Figure 24. Example of dry-stone surface features we expected to be detected on the VHR PALSAR-2 SL images: (a) low-level stone circles; (b) mound-like stone collection (rock-piles); (c) high-wall stone circles. The scale bar visible on the pictures is 1m long. The human figure in the picture is 1.78 m tall. Example of dry-stone features near-surface features we expected to be detected on the VHR PALSR-2 SL images: (a) exposed irrigation channel walls (Westerberg & Hagberg, 2018) (party excavated), photo by Rob Marchant (Marchant, 2023); Note the top left pedological profile inset from (Westerberg & Hagberg, 2018) indicting nature and stratigraphy of the Torripsamment soil type dominant in the area. In this inset, ‘A’ refers to Orchric epipedon layer, ‘C1’ for Sandy texture and ‘C?’ for Clayey or Silty texture (b) excavated irrigation walls from (Stump, 2006), the visible scales are 2 m by 1 m; (c) sediment traps and irrigation furrows from (Stump, 2006), scale is 1 m. These sediment traps cover an area of about 9 Km² in Engaruka North catchment (Kabora et al., 2020).

Figure 25. The faint detection of the linear stone terrace pattern in ENTA by PALSAR-2 SL data: (a) the linear stone terrace pattern as picked by the 34.2⁰ incidence angle data. Non-archaeological features (Maasai nomadic house structures, Manyattas) are marked with a yellow box; (b) The linear stone terrace pattern as detected by the image with 51.4⁰ incidence angle data. The pattern is more visible on the 2016 image with a lower incidence angle.

Figure 26. (a) On-ground photo of the same area shown in Figure 25. Notice the linear stone terraces (yellow arrows) and soil deposition in between the terraces; (b) a Google Earth image for the same extent showing crop marks; (c) a false-color WorldView-2 image for the same extent; (d) a true color WorldView-2 image for the same extent. The visible circular, non-archaeological features are Maasai nomadic house structures, Manyattas.

Engaruka North faces a severe gully erosion problem and direly needs urgent soil conservation measures. The USPED model indicates elevated (pronounced) erosion-deposition patterns along the gully networks as well as in areas with land use changes e.g. current farmlands. Vegetated areas e.g. those with trees had more stable erosion-deposition levels as compared to barren or gullied areas. Is it possible for the Ngorongoro Conservation Authority to reestablish the stone-terraces and raise them as in ancient times (Ferro-Vázquez et al., 2017) thereby regain control of surface runoff from the towering escarpment? Other possible measures may include erecting check dams, gabions as well as promoting tree planting in the lowland areas. This study highlights the importance of VHR DEM in understanding the distribution of sediment routing in an archaeological area. Our experimentation with radar data also forms a vital reference to other users intending to apply SAR data to study similar stone features in Sub Sahara region (Davies, 2013). The VHR datasets can also be incorporated in further studies in this data-poor area to monitor this crucial hydrogeological threat.

This study offers a valuable tool for the conservation and management of the da-ta-poor and semi-arid Engaruka archeological site. To the knowledge of the authors, this is the first attempt integrating different technologies (UAV, SAR and soil erosion/deposition models) in East African archaeological sites. However, it paramount to consider its limitations when applying the results for decision-making. It was difficult to quantitively validate the USPED empirical model results due to lack of measured/gauged stations. Future recommendation includes incorporation of complementary techniques such as physical erosion stakes, short-term sensors and more physical-based erosion models. The use of different datasets with different spatial resolutions can affect the results of empirical soil erosion models. Therefore, these results could be biased due to coarse spatial interpolation that compounds inherent limitations in the USPED model. USPED lacks a process-based structure thus fails to route sediment from its source to the sink (Hadush et al., 2025). Moreover, the data used is not real time but from different acquisition dates affecting the absoluteness of the results. Lack of current and precise data for R, K, C factors and general assumption of a unit (“1”) P factor reduced the accuracy of USPED and restricted its ability to predict short-term losses. The C parameter fails to consider tillage methods that affects erosion redistribution patterns. These USPED parameters can also be modelled in a future sense to account for anticipated land-use and climate changes (Assen et al., 2024). Presence of sparse vegetation cover, soil types and the stone-strewn landscape limited the detection capability of the SAR data. However, the comprehensiveness of this analysis can be improved by incorporating multiple SAR data from different satellite missions and with different viewing radar specifications as recommended by (Balz et al., 2016; Chen et al., 2016). The study estimated E-D rates for the Engaruka North Terraced Area and not the entire archaeological site. In addition, only a small section of the Intermediate Gorge was imaged and also didn’t sample other active gorges e.g. north of Makayuni river. Errors in the 3D gully volume estimates included gully rim boundary line demarcation, light and wind speed conditions during data take, interpolation and surface processing errors e.g. shadows from vegetation canopy cover (Yang et al., 2025), see Figure 10 and Figure 11. Similar to (Frankl et al., 2012; Lannoeye et al., 2016) multitemporal UAV or close-range imaging can provide better E-D change rates within the active gullies especially at the temporal scale of singular rainfall events. The archived aerial image scans weren’t applied in a true two-date 3D change analysis due to lack of sufficient cam-era calibration metadata. Future research could consider high-resolution data, both spatial and temporal e.g. terrestrial laser scanners so as to monitor gully formations, reveal in-depth erosion processes and their governing factors. Such techniques can quantify low-magnitude processes (dry-bank and aeolian erosion, debris flow and archaeological feature disintegration due to gullying see Figure 27). Socioeconomic factors e.g. overgrazing, nomadic lifestyle and improper cultivation play an important role on land degradation and soil erosion in Engaruka. Further studies can access these socioeconomic factors via detailed questionnaires about small-holder land-use practices and pastoralism.

Figure 27. Examples of loactions within ENTA recommended for very-high resolution erosion monitoring. (a) erosion due to overgrazing, settelments and trails (b) dry and wet bank erosion (red arrows) and effect on archeaologcal channel artefacts (c) an on-site photo showing wind and dust storms (d) archaoeological rock slides and movement due to rills and gullies (a disintergating stone circle in close proximity to a gully (yellow arrow)).

Despite its limitations, this study represents a critical step in the analysis of cultural heritage in data-scarce regions and is useful for conservation and mitigation planning. The study highlights the need to devise anti-erosive management particularly oriented towards the foot slopes of the western EARS and the gully networks. (Gashure & Wana, 2023) reported that smallholders in Konso, Ethiopia, a UNESCO AIIEA site responded to climate variability by adopting terracing, agroforestry, soil bunds, intercropping and reduced tillage. Other coping strategies they suggested include (1) reducing farm productivity and diversifying livelihoods into nonfarm/off-farm activities and (2) increasing community awareness on soil conservation so as to reduce encroaching onto the site (Assen et al., 2024). Similar measures can be adopted at Engaruka plus mechanical soil conservation measures for instance terrace reconstruction, check dams, gully plugs and minimizing the availability of rocks used to construct the Manyatta settlement wall fences contributing to the ruin of Engaruka’s dry-stone features.

Acknowledgements

The authors extend their gratitude to the University of Dar es Salaam and Department of Museums and Antiquities, Tanzania for their hospitality and assistance during field work. We also thank the Tanzania Commission for Science and Technology (COSTECH) for permitting us to obtain ground data within the Ngorongoro-Lengai UNESCO Global Geopark, Ngorongoro Conservation Area, COSTECH Permit No. 2019-682-NA-2019-386. Furthermore, the authors sincerely thank Professor Daryl Stump, his (Archaeology of Agricultural Resilience in Eastern Africa) AAREA team, Professor Amaury Frankl and Professor Lars-Ove Westerberg for their constructive insights. We appreciate the support we received from Mr. Christognas Ngivingivi, Ministry of Natural Re-sources and Tourism, Department of Antiquities, Tanzania. We finally would like to thank the reviewers who significantly enhanced this article.

Author Contributions

George Watene, Yueping Nie, Lijun Yu were responsible for the conceptual and research method design. George Watene is the principal researcher and wrote the manuscript. Thomas Biginagwa, Christognas Ngivingivi and Leka Tingitana were responsible for UAV flight permits, data acquisition and data analyses. Felix Mutua, Thomas Ngigi, Xianhu Wei played the project supervision and administration roles. Yves Hategekimana and Thomas Ngigi were responsible for reviewing, statistical analysis and discussion. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References