Quantitative soil erosion risk assessment due to rapid urbanization in the Cox’s Bazar district and Rohingya refugee camps in Bangladesh

Soil erosion due to anthropogenic interventions is an emerging threat to the south-eastern coastal districts of Bangladesh. In the recent decade, land degradation intensified this process in this region due to the settlement of the Forcibly Displaced Myanmar Nationals popularly known as ‘Rohingya’, who fled from Myanmar in 2017. To control soil erosion, it is necessary to identify the soil erosion locations and their intensity. The present study has attempted to quantify the regional soil loss (from 2015 to 2020) using the Revised Universal Soil Loss Equation model. In this method, soil erosion risk assessments are carried out by integration of different physiographic parameters in the spatial dimension. Results showed that the year 2015 and 2020 witnessed mean soil erosion at a rate of around 58.2 and 59 t ha− 1 yr− 1, respectively in the Cox’s Bazar district in which 20–21% of the study area was highly subjected to soil erosion. But in the Rohingya camps mean soil erosion for these periods is 59 and 78 t ha− 1 yr− 1, respectively with a cumulative increase of 32%. The sub-districts adjacent to the Rohingya refugee camps (Ramu, Ukhia, and Teknaf) have experienced intense erosion compared to the others. From 2015 to 2020, a large increase of about 4.54 (6%) t ha− 1 yr− 1 is observed in the Palongkhali union whereas a maximum decrease of about − 8.08 (23%) t ha− 1 yr− 1 is observed in the Pokkhali union. In Rohingya camps, the worst soil erosion-affected areas are the camps 8E, 10, 14, 15, 16, and 17 having more than 100 t ha− 1 yr− 1 mean soil loss. The geographic detector method is used to identify the influence of different factors on soil erosion. Land use change is found to be the major contributing factor to soil erosion increase from 2015 to 2020. The spatial variations are significantly controlled by slope and elevation factors. This work can help to understand the dynamics of soil erosion and the knowledge can be used for conservation planning in the Cox’s Bazar and Rohingya camp areas.


Introduction
Soil erosion is one of the major geo-environmental and agricultural problems in the present world (Kefi et al. 2011;Wu and Chen 2012;Borrelli et al. 2017). Although soil erosion was observed historically, its intensity in the recent time increases manifold due to large-scale anthropogenic interventions on the natural environment. Land use change through intensive agricultural practices and deforestation are the major causes of the intensification of soil erosion (Alkharabsheh et al. 2013;Li et al. 2014;Sun et al. 2020;Borrelli et al. 2020;Negese 2021 due to soil erosion is leading to loss of nutrient-rich fertile topsoil, increased surface runoff by exposing impermeable subsoil/bedrock, increase flooding/flash flooding extent by reservoir and valley siltation, and reduce water availability to plants (Ganasri and Ramesh 2016). During the last 40 years, one-third of the world's cultivable land resources have been lost by soil erosion. Every year 10 M ha cultivable lands will continue to be lost in the future (Pimentel et al. 1995). Therefore, estimation of soil erosion and identification of critical areas are the major concerns to combat the soil erosion problem. The soil erosion problem varies at different spatial scales due to large heterogeneity in the hydrological, geological, and climatological processes. Thus, the spatial distribution of soil erosion is a useful indicator to measure and locate the extent of erosion-affected regions and their actual causes. Spatiotemporal variation is a basic characteristic of soil erosion. Empirical models have been used to analyze the distribution of soil erosion, among which the universal soil loss equation (USLE) model or the revised universal soil loss equation (RUSLE) model is the most widely used (Amore et al. 2004). The modeling results can show the spatial variation in soil erosion, but the expression ability of the geographical intrinsic features must be further quantified (Ren et al. 2018). Apart from the spatial variation, the influencing factors both natural and anthropogenic causing these variations are also important to identify. Geographic detector methods and regression-based methods are widely used to identify the contribution of different factors to final soil erosion results (Gao and Wang 2019;Guo et al. 2021). The geographic detector method is the most suitable to identify spatial heterogeneity whereas the regression-based methods are mainly effective for fixed location-wise calculation. Therefore, in this study geographic detector method was used to evaluate the soil erosion influencing factors.
In Bangladesh, over 87% population is directly or indirectly engaged in agricultural activities (The World Bank 2016). It is mainly due to the high fertility and productivity of the topsoil deposits in the Bengal deltaic plains (Biswas et al. 2019). The increase in soil erosion threatens agricultural productivity and hence could impede the country's efforts to maintain the present food security in the future if soil erosion is not controlled within the limit. Siltation of the major river channel, reservoir, and hillslope valleys is causing flash floods and large irrigation problems. Which sometimes causes damage to planted crops and loss of property and lives. Lack of information on soil erosion severity and spatial variation, major challenges are faced by the government and other organizations to implement necessary steps to control soil erosion. Most of the watershed or river basins particularly in south-eastern Bangladesh are ungauged in terms of sediment dynamics and soil erosion. Therefore, spatial-based empirical models such as RUSLE are imperative to quantify the soil erosion dynamics in this region.
In this study, the spatial distribution of soil erosion for the periods of 2015 and 2020 was analyzed using the RUSLE method in the Cox's Bazar district and Rohingya camps. In this time frame, large-scale environmental degradation was observed here due to forced migration and settlement of Rohingya refugees in this area in 2017. The continuous development activities of the local government in Cox's Bazar district were also causing reduction of protected and hillslope forests and major land use change. These factors were the major support bases for choosing this area to study soil erosion distribution. No previous study was conducted in this region to quantify soil erosion. Hence, the main goals of this study were to (1) derive the soil erosion distribution for two temporal periods (2015 and 2020), (2) soil erosion in different administrative levels in two sub-domains: Cox's Bazar and Rohingya camps, and (3) determine the contribution of different influencing factors on soil erosion. The assessment of the soil erosion risk in the Cox's Bazar district and Rohingya camps is necessary because extensive soil erosion is adversely affecting the natural landscape and slope failure hazards. Determination of factors influencing soil erosion is also lacking in the study area, hence this study will give an idea of natural and anthropogenic factors that are contributing to intense soil erosion.

Site description
The study area is divided into two parts: The main domain is the Cox's Bazar district and the Rohingya camps are the subdomain. The study area is located in the southeastern portion of Bangladesh (Fig. 1). It is situated on the Chittagong-Tripura subduction belt and is composed of a series of anticlinal hills, hillocks, and synclinal valleys. The elevation of the Cox's Bazar district ranges from 0 to 261 m and Rohingya camps from 0 to 43 m (Fig. 1b, c). The average temperature of this area is about 25.6 °C. The predominant sub-tropical and coastal climate of Cox's Bazar district has more than 4000 mm annual rainfall. This area is mostly covered by clastic sand, silt, and clay-rich soil deposits. Agricultural activities on the hillslopes and plain land are extensively practiced in this area. The major cereal crop is rice and other different vegetables in Cox's Bazar (FAO 2020). After 2017, nearly a million Rohingya refugees migrated to this area from nearby Myanmar due to conflict and crime against humanity and settled in numerous makeshift camps. The Rohingya settlement causes large-scale environmental deterioration including extensive soil erosion (Hasan et al. 2020). Eroded hillslope soil mass caused extensive siltation of the valleys, which eventually resulted in an intensification of flash floods in the camp area (Quader et al. 2020).

Geological and geomorphological descriptions
The regional geology of the Cox's Bazar district consists of Tertiary clastic to Recent alluvium deposits (Reimann and Hiller 1993). The subduction of the Indian tectonic plate beneath the Burmese plate formed a series of alternating anticline and syncline geometry in the study area (Hossain et al. 2019). The anticlinal areas are elevated and form large hills, hillocks, steep cliffs, and small valley type geomorphological features. The average slope is also higher in the anticlinal area which is about 15° to more than 20° (Fig. 2a). The synclines are depressed low lying areas consisting of broad valleys, numerous stream channels and small hillocks

Data sources
Manifold data sets have been used for soil erosion assessment in the study area using the revised universal soil loss equation (RUSLE) method. The rainfall station data from Kutubdia, Cox's Bazar, and Teknaf having 24 h measuring frequency was collected from Bangladesh Meteorological Department (BMD: https://bmd.gov.bd) and used for erosivity estimation. The location of the rainfall stations is shown in Fig. 1. The Shuttle Radar Topography Mission (SRTM) elevation data sourced from https://earthexplorer.usgs.gov to prepare the slope and LS factor. USGS and GSB produced surface geological maps were used for soil erodibility estimation (Persits et al. 2001). For different parameter estimations of empirical erodibility equation, soil samples were collected from each of the lithological units during the fieldwork. The cloud-free Sentinel 2 satellite imagery type geomorphological features. In the valley area, average surface slope is less than 10°. In the elevated anticline older geological formation such as tertiary age Bhuban and Bokabil formations are exposed (Fig. 2b). Both of these formations are composed of alternating sandstone and shale type lithologies. In the syncline area, Tipam sandstone, Girujan clay and DupiTila sandstone formations of the Pliocene to Pleistocene age are exposed. The recent marsh clay, beach deposits and alluvium occurred as the thin surface of these geological units in the low-lying areas. Due to the undulating topography of the study area during and after the rainfall event infiltrated and surface runoff quickly moves to the valley or stream channel networks carrying a huge load of soil or sediments. estimation in this study. The equation for soil loss simulation using the RUSLE method is as follows: where SE is the annual soil erosion (t ha − 1 yr − 1 ), R is the rainfall erosivity factor (MJ mm ha − 1 h − 1 yr − 1 ), K is the soil erodibility factor (t h MJ − 1 mm − 1 ), LS is the slope length and steepness factor, C is the land cover and management factor, and P is the conservation support practice factor. The limitations of the RUSLE method in soil erosion prediction are that it is mainly sensitive to sheet and rill erosion and it does not account the gully erosion and dispersion of soils. The deposition of sediment before reaching the waterways is also not considered in this method. In the present study, soil erosion using the RUSLE method was estimated in a 10 m × 10 m GIS grid.

Rainfall-runoff erosivity factor (R)
The R factor is a major parameter for soil erosion calculation. The absence of detailed rainfall data in the study area particularly storm kinetic energy (E) and maximum 30 min storm intensity (I 30 ), makes the analytical calculation of R impossible. Therefore, indirect methods are suggested for R factor estimation using empirical equations (Renard and Freimund 1994). The spatial distribution of erosivity was estimated using Eq.
(2) of Singh et al. 1981: where P is annual rainfall in mm. The annual average rainfall of the Kutubdia, Cox's Bazar, and Teknaf stations were used for R estimation. The point-derived R values were then interpolated over the entire study region.

Soil erodibility factor (K)
The susceptibility of soil to be eroded by water is expressed by the K factor. It depends on the inherent properties of soil e.g., soil depth, texture, grain size, organic content, etc. (Wang et al. 2013;Auerswald et al. 2014). The nomographbased equation proposed by Wischmeier and Smith 1978 was used for soil erodibility factor estimation. The USGS and GSB joint produced surface geological maps of Bangladesh along with the field soil sample analysis was used to estimate different parameters of Eq. (3) for K estimation.
at 10 m spatial resolution sourced from European Space Agency (ESA: https://scihub.copernicus.eu) was used to prepare land use map of the study area. Digital shapefile of the different level administrative boundary of Cox's Bazar district was also used during the spatial analysis. The data type, derived product, and source of the datasets used for soil erosion estimation are given in Table 1.

RUSLE method
Basin-scale soil erosion assessment can be carried out by simple empirical models such as RUSLE and complex physical process-based models such as ANSWERS, EURO-SEM, SWAT, WEPP, GUEST (Beasley et al. 1980;Flanagan and Livingston 1995;Morgan et al. 1998;Rose 1999;Gassman et al. 2007;Sun and Wu 2022). The GIS based morphometric analysis of the drainage is also used for soil erosion estimation (Benzougagh et al. 2022). The RUSLE model (Renard et al. 1997), revised from the USLE model (Wischmeier and Smith 1978), has been widely used to estimate soil erosion worldwide (Ganasri and Ramesh 2016;Zerihun et al. 2018;Thomas et al. 2018;Teng et al. 2018;Rajbanshi and Bhattacharya 2020;Brahim et al. 2020). This model input datasets are easily available and less expensive to collect. The most important benefit of the RUSLE method is that it can be implemented in a grid-based GIS environment. It is more sensitive to the sheet and rill erosions but gully and channel erosions are not fully apprehended by this model. On the other hand, process-based models can evaluate detailed soil erosion from a different modes of occurrences but require high and sophisticated data and tools (Jetten et al. 2003;Li et al. 2017;Hajigholizadeh et al. 2018). To enable soil loss estimation in regions of scarce data and the use of these estimates by local stakeholders, it was vital to select a simple scheme or modeling method that is flexible for poor data availability and simulation of different scenarios based on local interventions on soil loss. Therefore, the RUSLE method was used for soil erosion where, α = 1 and β = 2. The NDVI value was calculated in spatial dimension from red and near-infrared bands of Sentinel 2 images using Eq. (6): where NIR = near-infrared band and Red = Red band.

Conservation practice factor (P)
The P factor represents the ratio of soil loss to the support practice on the up and down the slope of the agricultural (vegetation) conservation practices and is used to accommodate for the positive impacts of those support practices. The values closer to 0 indicate good conservation whereas the upper values (max. 1) are vice versa. In this study, the P factor value was assigned from previous studies based on land use type and is shown in Table 3.

Geographic detector method
The geographical detector method was used to detect the spatial heterogeneity of an event and assess the relationship between the event and potential risk factors (Wang et al. 2010). It is frequently used to measure the impact of the influencing factors on the derived final product (Gao and Wang 2019;Guo et al. 2021). In this study, it was used to quantify the impact of different natural and anthropogenic factors on controlling the spatial distribution of soil erosion. The basis of the method is that if the sum of the variance of the subareas in a region is less than the variance of the total region, spatial heterogeneity exists in the area (Wang and Xu 2017). With a geographical detector, the proportion of the spatial distribution of Y (the soil erosion in this study) (1) is very structured or particulate, (2) is fairly structured, (3) is slightly structured, and (4) is solid; c = profile permeability code in which (1) is rapid, (2) is moderate to rapid, (3) is moderate, (4) is moderate to slow, (5) is slow, and (6) very slow. The average value K estimated by Eq. (3) of the different soil types is shown in Table 2.

Slope length and steepness factor (LS)
The LS factor in the RUSLE method indicates the effect of topographic undulation on soil erosion (Suhua et al. 2013;Panagos et al. 2015). It is the combined impact of slope length (L) and slope steepness (S). The LS factor is dimensionless and obtained by using Eq. (4) of Wischmeier and Smith 1978. This equation is frequently used for LS factor calculation (Ganasri and Ramesh 2016). where, FA = flow accumulation; S = slope (percentage); and m = dimensionless exponent ranges from 0.2 to 0.5. The LS factor was directly calculated from SRTM 30 m DEM in the GIS grid and then resampled down to 10 m.

Cover management factor (C)
The C factor determines the impact of land management practice on soil erosion intensity. In the RUSLE method, this factor is estimated from land use and cropping pattern, canopy cover, surface cover, etc. The C factor can also be estimated from NDVI based equation that was applied for this study (Alexakis et al. 2013).  are higher in the downslope or coastal areas and lower in the highly vegetated hilly areas. A significant change in C-factor was observed in the Rohingya camps area from 2015 to 2020 (Fig. 3d, e). The soil erosion support practice is not regularly and properly implemented in the studied region. Therefore, the P-factor value is assigned based on the Land use pattern following Table 1.

Estimation of soil erosion in the Cox's Bazar district
The soil erosion calculated using the RUSLE method for the Cox's Bazar district is shown in Fig. 4. The soil erosion is heterogenous and distributed all over the study area. The heterogeneity remains constant for both temporal periods 2015 and 2020. The mean annual soil erosion rate for 2015 and 2020 are 58.2 (range: 0 to 1785) and 59 (range: 0 to 1876) t ha − 1 yr − 1 , respectively. The mean annual soil erosion from 2015 to 2020 is not increased significantly but has a slightly increasing trend. Soil erosion distribution shows a large range of uncertainty. On average 75-80% of the study area shows slight to moderate soil erosion rates. High and very high soil erosion region comprises 18-20%. Severe and very severe soil erosion area is only 1% of the study area (Table 4). Except for the slight and moderate soil erosion classes, all the soil erosion classes increase from 2015 to 2020. The severe, very high, and high soil erosion regions have coincided with the high elevation and sloping topography (Figs. 2a and 4). Low erosion regions are the flat and watersubmerged areas. The severe and very severe soil erosion zones occur in small patches that are mostly in the Teknaf sub-district regions. These patches are the hilltops of Nhila anticlines and are severely deforested. A significant increase in soil erosion from 2015 to 2020 is only observed in the Rohingya camps area.

Spatiotemporal comparison of soil erosion in the Cox's Bazar district
The mean annual soil erosion value for 2015 and 2020 both for the Upazila (sub-district) and union levels are shown in Fig. 6. At the Upazila level, the least soil erosion-affected area is the Pekua and the highest in the Ramu and Teknaf both in 2015 and 2020. Kutubdia, Maheshkhali, Pekua, and Chakaria Upazilas experience a relatively lower rate of soil erosion whereas Ramu, Ukhia, and Teknaf experience a high rate of soil erosion. Moderate soil erosion is observed in Cox's Bazar Sadar Upazila. From 2015 to 2020, soil erosion substantially increases in the Ramu, Ukhia, and Teknaf. A large increase is observed in the Ramu Upazila. Soil erosion decreases in the Kutubdia and Maheshkali Upazilas. that can be explained by X (influencing factors) is measured by a determinant (q value) as follows: where σ 2 is the variance of Y in the region, σ z 2 is the variation in zone Z divided by X, N is the number of sample units in zone Z, and L is the number of zones. σ 2 and σ z 2 can be expressed as Eqs. (8) and (9), where Y z,i, and Y j are the values of Y in the i-th sample unit of zone Z and the j-th sample unit of the entire region, respectively.
The input data for the geographical detector analysis include Y: soil erosion results, and X: influencing factors (slope, elevation, NDVI, rainfall, lithology, and land use) were categorical. The continuous factor maps were categorized using the natural break attributes of the ArcGIS software.

Estimation of RUSLE factors
The factors necessary for soil erosion estimation using the RUSLE equation are shown in Fig. 3. The R-factor calculated from rainfall data of 3 stations of Cox's Bazar district ranges from 1226 to 1608 MJ mm ha − 1 h − 1 yr − 1 . The R-factor value is higher in the southeast and decreases in the north direction (Fig. 3a). The R-factor is much higher compared to the other area of the world because the studied region receives a large amount of rainfall (> 4000 mm/yr). The soil erosion rate is highly sensitive to rainfall or R-factor (Wischmeier and Smith 1978;Dabral et al. 2008). The soil erodibility factor (K) value calculated from surface lithology ranges from 0 to 0.30 t h MJ − 1 mm − 1 . Soil lithology mainly consists of sandy loam, silt loam, loam, silty clay, and sandy clay loam ( Table 2). The majority of the area is covered by loam, sandy loam, and silt loam. The LS-factor determined from the SRTM elevation data ranges from 0 to 9.83. The LS-factor is also higher as the study area has undulated and hilly topography. The C-factor determined from NDVI using Eq. (5) ranges from 0 to 1. C-factor values Ramu increased at a significant rate. From all of the unions of Cox's Bazar district (n = 74) mean annual soil erosion rate increased in 44 unions, decrease in 28, and remains unchanged in 2 unions (Fig. 8). Soil erosion decreases in all the unions of Kutubdia. The majority of the. unions of Pekua, Chakaria, and Maheshkhali also show a decline. Whereas unions of Cox's Bazar Sadar, Ramu, The erosion rate remains unchanged in the Pekua, Chakaria, and Cox's Bazar Sadar Upazilas (Fig. 7). Some changes in union level soil erosion are also observed. Union level soil erosion remain the same or reduced in the northwest part of the study area between 2015 and 2020 (Fig. 6c, d). In the middle part of the study area, soil erosion from 2015 to 2020 remains constant but in the southern union and east of

Estimation of soil erosion in the Rohingya camps
Significant increases in soil erosion rate from 2015 to 2020 are observed in the selected Rohingya camps as shown in Fig. 9. The mean annual soil erosion rate in the Rohingya camps for 2015 and 2020 are observed at about 59 (range 0 to 375) and 78 (range 0 to 478) t ha − 1 yr − 1 , respectively. Large variations in the soil erosion rate are also observed in the Rohingya camps as in the Cox's Bazar district.
The standard deviation of the soil erosion rate is about 36 and 47 for 2015 and 2020, respectively (Fig. 5). In 2015, slight, moderate, and high soil erosion classes comprise about 42.06%, 45.42%, and 10.79% of the total area of the Rohingya camps. But in 2020, these classes change to 26.93%, 44.11%, and 21.61%, respectively ( Table 5). It is observed that significant reduction occurs in slight and moderate soil erosion classes while enlargement was observed in high erosion classes. Highly (100-150 t ha − 1 yr − 1 ) eroded soil area increased about 100% from 2015 to 2020. Very high, severe and very severe soil erosion-affected areas are less than 2% in 2015. But in 2020, nearly 8% area is affected by very high to very severe soil erosion. Changes in these Ukhia, and Teknaf show a rise in soil erosion. The highest amount of increase of erosion is observed in the Palongkhali union about 4.54 (6%) t ha − 1 yr − 1 where the Rohingya camps are located and the highest amount of decrease of erosion is observed in the Pokkhali union of Cox's Bazar Sadar about 8.08 (23%) t ha − 1 yr − 1 .  They combined have a 76% probability of controlling soil erosion distribution. In the Rohingya camps, the q-value of the slope and elevation are 0.51-0.53 and 0.01-0.04, respectively, and ~ 87% probability of controlling soil erosion. The rainfall factor influences about 5-7.8% (q = 0.033-0.086) in 2015 and 3.5-3.9% (q = 0.022-0.038) in 2020. The slowly varying lithology influences about 1-4% of soil erosion distribution. The anthropogenic influence is expressed by NDVI and Land use factors. The anthropogenic loss of vegetation is reflected by the NDVI signature and the loss of green space replaced by building expansion can be attributed to spatiotemporal land use change. So that, these two factors were considered to characterize the anthropogenic influence on soil erosion. Land use factor has 2.2% (q = 0.013) control in 2015 and 4.6% (q = 0.029) control in 2020 in the Rohingya camp area. The influence of land use on soil erosion increases by about 120% in the Rohingya camps between 2015 and 2020. But in Cox's Bazar area influence of the land use factor only change from 5 to 5.5%. The influence of NDVI from 2015 to 2020 increased in Cox's Bazar but decrease in the Rohingya camps (Fig. 12).
The combined impact of influencing factors expressed by the interaction q-value as shown in Fig. 13 is much higher than the solitary factors. The combined impacts of slope and NDVI (~ 0.64 and ~ 0.55), slope and rainfall (~ 0.58 classes from 2015 to 2020 range from 350 to 2000% due to Rohingya settlement.

Spatiotemporal comparison of soil erosion in the Rohingya camps
In the selected Rohingya camps (n = 26) soil erosion rate increased in all the camps from 2015 to 2020 and is shown in Fig. 10. In 2015, before the Rohingya settlement mean annual soil erosion rates were confined below 85 t ha − 1 yr − 1 . Hakimpur camps (14,15,16) are the relatively worst soil erosion-affected areas in 2015. In the Kutupalong -Balukhali expansion site camp 1 E, 1 W, 2 W, 8 E, 9, and 13 experienced the highest soil erosions in 2015 (Fig. 10a). But after Rohingya settlement in these camps soil erosion increases drastically to more than 100 t ha − 1 yr − 1 . At this time worst soil erosion-affected areas are the camps 8 E, 10, 17, 14, 15, and 16 (Fig. 10b). The graph of camp-level soil erosion changes from 2015 to 2020 is shown in Fig. 11. The mean annual soil erosion rate changes in the camp level range from 5 to 30 t ha − 1 yr − 1 with an average of 19 t ha − 1 yr − 1 . The highest change is observed at camp 4-EXT. A notable increase in erosion rate is also observed in camps 4, 5, 6, 17, 18, 20, and 20-EXT. The erosion rate remains constant in the KRC camp within the studied time interval. In the selected Rohingya camps, the overall soil erosion rate increases about 32% from 2015 to 2020.

Soil erosion influencing factors
The distribution of soil erosion in the Cox's Bazar and Rohingya camps area is influenced by a set of natural and of the world where soil erosion estimation was done using the RUSLE methods (Haregeweyn et al. 2017;Kijowska-Strugała et al. 2018;Fenta et al. 2021). The higher R factor value is due to the large-scale rainfall > 4000 mm/yr in Cox's Bazar district. The measured R-value was in good agreement with the value estimated for this area by Panagos et al. 2017 and nearby similar climatological regions (Ganasri and Ramesh 2016;Nakil and Khire 2016;Islam et al. 2020). The soil erodibility (K) factor was also high as the area is covered by weak and loose sand-rich surficial soil cover, which is prone to water-induced gully and rill erosion (Fig. 14).
In the Cox's Bazar district moderate and high soil erosion affected areas are predominant. Annual average rainfall for 2015 and 2020 was measured at about 58.2 and 59 t ha − 1 yr − 1 , respectively. These values show a good correlation For Rohingya camps slope and land use factor is the most important one that controls the temporal increase of soil erosion from 2015 to 2020.

Discussions
The RUSLE method was used to estimate water-laden soil erosion in the Cox's Bazar district and the heavily deforested Rohingya camps for the first time. Soil erosion rate remains almost the same from 2015 to 2020 in the Cox's Bazar district but largely increased in the Rohingya camps region. Measured soil erosion shows a higher degree of uncertainty both in the Cox's Bazar district and Rohingya camps (Figs. 5 and 11). In the studied region rainfall-runoff, erosivity (R) factor is much higher than in the other parts   Thomas et al. 2018). Ramu, Ukhia, and Teknaf sub-districts have the highest soil erosion due to their high topography and rugged terrain (Fig. 6a, b). The same reasons are applicable for higher soil erosion in the eastern and southern unions of Cox's Bazar (Fig. 6c, d). After 2017 Rohingya settlements in various camps in the Ukhia and Teknaf subdistricts caused large-scale gully and rill erosion that was observed during the fieldwork in these regions (Fig. 14a, c, f). About a 32% increase in the mean annual soil erosion rate was observed within a very short period from 2015 to 2020. This is mainly caused by large-scale deforestation and slope cutting for human settlement in the previously dense forested area. Weak and loose sandy lithology coupled with vegetation logging resulted in such extensive soil erosion in the Rohingya camps. All the Rohingya camps experienced an increase in soil erosion from 2015 to 2020. But camps 4,5,6,17,18,20, and 20-EXT had gone through serious soil erosion during this time interval. In response to extensive erosion of the slope and hill ridges' soil mass, extensive siltation has occurred on the valley floor. Which with the soil erosion assessed in the regions with similar geo-environmental conditions (Ganasri and Ramesh 2016;Rajbanshi and Bhattacharya 2020). In Cox's Bazar district, soil erosion is almost zero in the northwest and the far south. The coastal and low-lying swampy areas also have very small erosion rates sometimes negative. This was mainly due to the flat topography and deposition of the eroded soil mass from the hinterland (Ostovari et al. 2017;  the Cox's Bazar and Rohingya camps area using spatial statistical correlation. The slope and elevation are the dominant factors for higher soil erosion in Cox's Bazar district. Several authors have previously reported that in mountainous or eventually resulted in an increase in flash flood magnitude in recent times (Quader et al. 2020).
The geographic detector method successfully determined the influence of different physical factors on soil erosion in  Due to the lack of large-scale erosion experiments and gauged station data, the estimated soil erosion was compared with previous studies in regions with similar climatological and geological environments (Ganasri and Ramesh 2016;Nakil and Khire 2016). Significant uncertainties have existed in the RUSLE model-derived soil erosion which was common in the empirically-based models (Alewell et al. 2019). The different spatial and temporal resolutions of the datasets e.g., elevation, precipitation, land use, soil map, etc. were the important sources of uncertainty (Pappenberger and Beven 2006; Mondal et al. 2016). In the study area, gully and channel-cut erosion are prominent which are not properly accounted for in the RUSLE method (Renard et al. 1997). Therefore, further studies are required in the study area which should test and validate the RUSLE factors in the field. Sediment dynamics should be measured in the major watershed channels for validation of the final soil erosion results, in this case, Bakkhali river and for Rohingya camps nearby confluence of the Naf River could be the best place for monitoring. Gully and channel erosion should be assessed using field experiments and physical process-based modeling to quantify the overall soil erosion rate.

Conclusion
In this study, the distribution of soil erosion and influencing factors were estimated using the RUSLE and geographic detector methods for two temporal periods (2015 and 2020). Between these times major demographic and environmental changes have occurred in the study area. The soil erosion rate estimated using the RUSLE method for 2015 and 2020 showed that Cox's Bazar district and Rohingya camps have experienced about 58.2 and 59 t ha − 1 yr − 1 , and 59 and 78 t ha − 1 yr − 1 , respectively. In Cox's Bazar district, around 20-21% area was affected by high soil erosion, but temporal changes from 2015 to 2020 were not significant. The subdistrict and unions of Cox's Bazar Sadar, Ramu, Ukhia, and Teknaf were mainly affected by soil erosion than the others. But in the Rohingya refugee camps, a 32% increase in mean erosion was observed. Highly eroded area changes from 13 to 28% from 2015 to 2020. The major influencing factors of soil erosion were the slope and elevation, but the temporal increase was driven by mainly land use change both in Cox's Bazar and Rohingya refugee camps. The land use change impact on soil erosion was more prominent in the Rohingya refugee camps than in the whole Cox's Bazar district. This study provided sufficient quantitative information on the spatial and temporal distribution of soil erosion by rigorous scientific analysis. This information can be utilized to implement soil erosion controlling measures. hilly areas where an undulated topography exists slope and elevation are the predominant soil erosion controlling factors (Tamene et al. 2006;Guo et al. 2021). But the temporal variations from 2015 to 2020 were significantly controlled by land use change and associated vegetation (NDVI) response. Land use change is the major cause of soil erosion intensity change all over the world (Yang et al. 2003;Borrelli et al. 2017;Aneseyee et al. 2020). The influence of land use change on soil erosion increased from 2015 to 2020 in the Cox's Bazar district was not significant only 5.5%. But in the Rohingya camps, a 120% increase in land use change influence was quite notable. In terms of combined impact which was determined by the interaction q-value of geodetector analysis, slope ∩ land use factor dominantly controlled the soil erosion.

Data availability
The datasets used in this study are available for sharing on request to the corresponding author.

Declarations
Ethical approval Not Applicable.

Consent to participate Not Applicable.
Consent to publish All the authors consented to publish the paper.

Competing interests
The authors have no relevant financial or non-Acknowledgements The authors of this manuscript are very grateful to the Department of Geology, and the Department of Disaster Science and Climate Resilience, University of Dhaka for providing licensed software facilities. We have to convey our gratitude to the Prof Dr. Peter Sammonds, Institute for Risk and Disaster Reduction, University College London who initially checked our manuscript before final submission. We also appreciate the suggestions and recommendations from anonymous reviewers that essentially helped us to improve the manuscript throughout the review process.