Geodetic Landmarks Soil Loss Susceptibility Using RUSLE and Climate Change Scenarios

Geodetic landmarks (GLs) are essential for obtaining the precise height, horizontal coordinates, and the Earth's gravity field. Once physically implanted on the surface, they are susceptible to movement and displacement. This study aims to assess the soil susceptibility of GLs for past and future scenarios through the Revised Soil Loss Equation (RUSLE). So the soil loss estimations were made for the GLs in Brazil's southern Santa Catarina region. Our results showed average soil loss values, reaching 175915 t/ha/year, while the GLs were 2109 t/ha/year. There was an increase in GLs in the null class, mainly caused by urban infrastructure increase. At the same time, a decrease occurred in the low, very severe, severe, and moderate classes. In contrast, for future scenarios, an increase in the GLs average soil loss was found until 2100. However, it is essential to highlight that the most relevant increase occurred in the 2021-2040 period. After that, some scenarios as ssp126 remained stable, ssp245 and


Introduction
Geodetic Landmarks (GLs) are used to reference any survey and study performed in a region, which requires consistent and secure locations (AMIRI-SIMKOOEI et al., 2012). They made part of the national geodetic system in Brazil under the Brazilian Institute of Geography and Statistics (IBGE) responsibility.
They can be classified in altimetric (heights) and planimetric (horizontal coordinates).
Several GLs in Brazil are old, implanted since the mid-1940s, and may no longer exist (CARLOS et al., 2008). However, evaluating the susceptibility to erosion of the places where they were located is very important, especially regarding the possibility of restoring the network and replacing them in a safe area (FREITAS;FAGGION;MEDEIROS, 1998;SESTRAS et al., 2021). For this reason, the physical surface susceptibility needs to be estimated in order to assess areas for the GLs installation where they could be safe for a long time (YANG et al., 2020).
There are several methods for estimating soil loss, for example, the However, one of the most frequently used globally and with a wide worldwide bibliographical reference is the Revised Universal Soil Less Equation (RUSLE) (BIDDOCCU et al., 2020;BORRELLI et al., 2021;FAYAS et al., 2019;MAQSOOM et al., 2020).
The RUSLE formula is considered an update of the USLE equation to facilitate the implementation of computational solutions for soil loss estimation (GHARIBREZA et al., 2021;RENARD et al., 1997). The following was taken into account, there are several physical characteristics, such as (i) the energy with the soil raindrops impact, referred to as rainfall erosivity (R); (ii) the soil erodibility, which is a function of its composition, nature, and consolidation (K); (iii) the topographical factors of slope and slope length (LS); (iv) the land use and land cover (LULC) classes (C); and also (v) the factors that delay soil loss such as conservation practices (P) (MARKOSE; JAYAPPA, 2016).
One of the main factors behind soil loss is rainfall patterns and climate change (BEHERA et al., 2020). Global warming has contributed to the increase in extreme events that accelerate water erosion (UNITED NATIONS CLIMATE CHANGE, 2020). According to the latest report by the Intergovernmental Panel on Climate Change (IPCC) (2021), global temperatures will continue to rise over decades, primarily due to anthropogenic influences, and are expected to cause intense rainfalls in some regions associated with flooding.
The Coupled Model Intercomparison Project (CMIP) is used as the scientific basis of the IPCC and offers a wide range of numerical experiments that can simulate changes in rainfall for past, present, and future scenarios (EYRING et al., 2016;O'NEILL et al., 2016;PEDERSEN et al., 2021). Some studies have used CMIP scenarios to estimate soil erosion (PAL; CHAKRABORTTY, 2019; KUMAR et al., 2021) and found that soil erosion is closely related to simulated rainfall erosivity factors.
In this way, the climate models are powerful tools to evaluate future rainfall variability and soil loss, despite their inherent limitations. Therefore, modeling how climate change will affect soil susceptibility is considered a hot topic and lacks conclusive answers.
In this study, the use of RUSLE was chosen mainly because of its worldwide use. However, it was not found in the scientific literature their application to monitor GLs, nor using it associated with future climate scenarios to verify the soil erosion projected for the following decades, which is the main novelty of this research. So, this main study objective is to analyze the soil loss susceptibility for the GLs located in the southern region of Santa Catarina state in Brazil, considering a historical period of more than 30 years and future scenarios through projected climate simulation until the end of the century.

Study area
The Santa Catarina state southern region in Brazil has an area of approximately 9.7 thousand km². It is located between latitudes 27° 54' 18 "S and 29° 22' 43 "S and longitudes 48° 32' 00 "W and 50° 03' 26 "W, bordered to the north the Florianópolis city metropolitan region, and to the southwest the Rio Grande do Sul state, and to the east the Atlantic Ocean ( Figure 1).
The study area is part of the Atlantic forest biome characterized by dense rainforest and high anthropization. Some original forests remain on high slopes in the general mountain range escarpments (VIBRANS et al., 2013). Agriculture, pasture formation, forestry (eucalyptus), urbanization, and mining (coal) have occupied the less steep plains and slopes. The region has a humid mesothermal climate of the type Cfa (subtropical with warm summers), according to the Köppen classification (BECK et al., 2018), with an average annual temperature between 17 °C and 20 °C (EPAGRI, 2021) and well-distributed rainfall all year. The rainfall average is between 1220 mm and 1660 mm, with 98 to 150 rainy days (LOPES, 2015 (EMBRAPA, 1999).

Data acquisition
For the soil loss estimation using RUSLE, the following data were used for each component of the equation: 1. To calculate the rainfall erosivity (R) factor, data from 47 meteorological stations were used. They were operated by several institutions such as the National Water Agency (ANA), the Mineral Resources Research Company

RUSLE Model
The soil loss was calculated for the entire study area using the RUSLE model. However, before explaining it, we need to understand its predecessor, the USLE model, which was developed in the United States to provide a decision support tool for managing field studies in small water basins (WISCHMEIER; SMITH, 1965). Although these authors are often cited as the creators of the equation, the development of its concept began long before this with studies related to soil erosion rates (ALEWELL et al., 2019;SCHWERTMANN;VOGL;KAINZ, 1987).
In the 1990s, USLE evolved into RUSLE (RENARD et al., 1991), maintaining the USLE multiplicative form and implementing improvements in calculation processes such as the soil erodibility and rainfall erosivity (ALEWELL et al., 2019). Therefore, RUSLE combines the advantage of being based on the same USLE database with some process-based calculations for variable environmental effects on the erosive system (NEARING, 2013). In addition, the complete approach of all factors can be applied to a limited area (XU et al., 2008).
According to RENARD et al., (1996), the RUSLE is an erosion model developed to estimate the average annual soil loss by water on slopes and can be used for agricultural and non-agricultural purposes. However, it does not assess the sediments deposition (ZHANG; O'NEILL; LACEY, 1996) but estimates potential soil loss rates, which indicate the intensity of the erosive processes using the (Equation 1).

A=R×K×L×S×C×P
(1) Where A is the annual soil loss estimation (t/ha/year); R is the rainfall erosivity factor (Mj.mm/ha/hour/year); K is the soil erodibility factor (t.hour/Mj/mm); L is the slope length factor; S is the slope degree factor; C is the LULC factor, and P is the conservation practices all dimensionless factor.
The data used for each factor and the RUSLE process is shown in the flowchart ( Figure 3). The soil loss map was created with a spatial resolution of 30 meters, which is compatible with the digital elevation model and the LULC map's spatial resolutions. The erosivity map of the meteorological stations was created with the same spatial resolution. In the meantime, the soil constitution map, obtained from the SoilGrid data initially with a spatial resolution of 250 meters, was resampled to 30 meters.

Figure 3. Flowchart of the RUSLE soil loss calculation process
The soil loss estimates were made in the Santa Catarina southern region over a ten-year period from 1985 to 2019. In addition, we used a set of four different future climate scenarios (CMIP6-ssps) from 2015 to 2100, as described below.
Each of the steps in the RUSLE model equation are described in more detail below.

Rainfall erosivity (R) factor
The rainfall erosivity (R) factor quantifies the energy with which raindrops impact the soil. The erosive rainfall capacity on the earth is influenced by the raindrop's intensity and size (WISCHMEIER, SMITH, 1978). The calculation of the R-factor (equation 2) was calculated according to (LOMBARDI NETO;MOLDENHAUER, 1992).
Where: r is the monthly rainfall average (mm), and P is the annual rainfall average (mm).
The R factor was calculated for five years: 1985, 1995, 2005, 2015, and 2019 which were selected because they are part of the period covered by the LULC maps provided by the MapBiomas project. In addition, the year 1985, the beginning of our time series, is considered the peak of coal activity in the Santa Catarina state (MARTINS, 2005), reflecting several environmental impacts such as vegetation suppression and soil exposure.
To create the erosivity map, data from 47 meteorological stations were collected from the Hidroweb database (ANA, 2021) from 1985 to 2019. Of these, 31 stations were used after removing outliers ( Figure S1). The discard criteria followed the method of establishing an acceptability interval from the station average values and standard deviations.
The rainfall erosivity (R) maps were made using the meteorological station's data through the Inverse Distance Weighted (IDW) interpolation method.
The rainfall data for 1985,1995,2005, and 2015 were used to calculate erosivity for these specific years (Ryear). While rainfall data from 1985 to 2019 were used to calculate the average Erosivity average map (Raver) ( Figure S2).
The future rainfall soil erosivity (Rfut) map ( Figure  The ssps considered population growth in addition to the increased global mean temperature, as a response to anthropogenic radiative forcing (EYRING et al., 2016;O'NEILL et al., 2016;PEDERSEN et al., 2021).
The ssp126 represents the "best case" of the future from the sustainability point of view, with low radiative forcing 2.6 W/m² at the end of the century and CO2 emission decrease, reaching -8.62 Gt until 2100. In contrast, the ssp585 represents the upper limit of possible future paths, that is, the most "pessimist"

Soil erodibility (K) factor
The erodibility factor (K) is the relationship between soil loss per unit and rainfall erosivity (JIANG et al., 2020). Different soil types' propensity to erosion can be considered by two main factors: the infiltration capacity and the runoff process (WISCHMEIER; MANNERING, 1969).
The K-factor can be calculated from the sand, clay, silt, and organic compounds proportions present in the soil (Equations 3, 4, 5, 6 and 7) proposed by (WILLIAMS, 1995) and used for (ANACHE et al., 2015).
K=f csand ×f cl-si ×f orgc ×f hisand ×0,1317 Where fcsand is the factor for soils with high coarse sand content and soils with low sand content; fcl-si is the factor for soils with high clay and silt content; forgc is the factor for soils with high organic carbon content; while fhisand is the factor for soils with highly high sand content. The equations that calculate the factors are: Where ms, msilt, mc, and orgc are the sand, silt, clay, and organic carbon proportions in the soil, respectively, obtained from SoilGrids product (ISRIC, 2021).
Therefore, the soil erodibility map was generated at SoilGrid spatial resolution of approximately 237 m and resampled to 30 m for the soil loss calculation. So, the erodibility (K) map for the region was used for all future scenarios in this study ( Figure S4).

Slope length (L) and steepness (S) factors
The (L) The slope is measured in percentage; SlopeLenght is calculated by multiplying the accumulated flow by the pixel size; Constant is 22.1 for the metric system and 72.5 for the imperial system, and NN is a variable dependent on the slope accordingly (Table 1). For the topographic factor (LS) generation, the SRTM DEM with a spatial resolution of 30 m was used ( Figure S5). In flat regions, where no slopes are detected, (LS) values are equal to zero, indicating no soil loss at these locations.

Land use and land cover (C) and conservation practices (P) factors
The (C) factor is related to erosion processes related to LULC, with values ranging from zero for soils protected against erosion to one for soils without protection (GANASRI; RAMESH, 2016). At the same time, conservation practices (P) are related to actions to prevent or reduce soil erosion, which is an essential agent in soil loss (BAMUTAZE et al., 2021). Usually they are combined into one

Soil loss susceptibility analysis at the geodetic landmarks
The potential soil loss analyses at the GLs were performed in four steps. Finally, the future soil loss (Afut) for the following periods: 2021-2040, 2041-2060, 2061-2080, and 2081-2100 under ssp126, ssp245, ssp370, and ssp585 scenarios changing only the erosivity factor (Rfut) in the RUSLE model was estimated. In the third and final step, the GLs were classified as null, low, moderate, severe, and very severe soil loss values (Table 3) accordingly (ARNOLDUS, 1980;FAO, 1967). Very severe Source: (ARNOLDUS, 1980) 5 Results and discussion

Rainfall and erosivity
The minimum, maximum, and average rainfall and erosivity in the Santa  It was also observed that the higher (lower) rainfall average is associated with the higher (lower) erosivity average values.Therefore, the high erosivity average values shown in Figure 4 are response to high rainfall average values.
However, the maximum rainfall occurred in 2015, while the maximum erosivity in 1995. This probably happened because in 2015 the precipitation was better distributed over the whole year than in 1995, which had more concentrated extreme rainfall events in certain months, consequently increasing the maximum erosivity.

Soil erodibility factor
The soil erodibility (K) values in the Santa Catarina southern region ( Figure   S4) determined from the sand, clay, silt, and organic carbon proportions using the SoilGrids data (ISRIC, 2021), were compatible with other studies developed in Brazil, such as (GOMES et al., 2017) for similar soil classes. The soils of the region have maximum proportions of 59% clay, 51% silt, 47% sand, and 14% organic carbon.

Topographic Factor
The topographic (LS) factor, which is based on the slope length factor (L) and the steepness factor (S) together, showed a significant discrepancy due to the different slope classes present in the Santa Catarina southern region. The relief varies from the coastal plain to the general mountain range escarpments, from zero meters at sea level to approximately 1816 m of altitude.
The topographic factor values ranged from zero in flat areas, with no soil loss due to no sedimentary movement, to 23640 in steep areas, with severe soil loss, being classified up to LS factor of 20 and above ( Figure S5).

Land use and land cover (LULC)
The  It can be seen the GLs growth located in the urban infrastructure and its decline in the pasture areas over the years, while the other classes presented stability. A larger GLs quantity in the urban infrastructure class is associated with zero soil loss. However, it is known that this is not true because, in the urban perimeter of cities, towns, villages, and paved roads, there are regions that are not impermeable, often with exposed soil where erosive processes take place.
Furthermore, some GLs were also classified within water bodies, such as rivers, lakes and oceans, mainly located along the Port of Imbituba. That occurred mainly due to the LULC map's spatial resolution that generalizes the pixel value.
So this is considered an error and a problem due to the spatial data used.
Anyway, only about 1% of the GLs are classified in the water bodies class and consequently are not considered in our soil loss analysis. In addition, it is worth noting that of the 896 GLs listed by the IBGE for the Santa Catarina southern region, about 2% are in pixels with no data. Therefore they were also discarded from the analyses, and 877 GLs were assessed in the next stage.

Soil loss estimation
The average soil loss (A) values in the Santa Catarina southern region ranged from zero t/ha/year in flat areas and urban infrastructure to 175915 t/ha/year in sloping areas to the high runoff ( Figure 5). In contrast, the GLs soil loss values varied from zero t/ha/year, also found in urban and flat areas, to 6357 t/ha/year. Because this maximum (A) value is very different from the dataset, this GLs was checked, and its location was confirmed inside the water LULC class, which in theory is not possible. However, the spatial resolution of the data used to calculate the soil loss using RUSLE can cause this kind of inconsistencies.
Therefore, an acceptance interval was established using the mean and the standard deviation to discard GLs outliers. So, GLs with null and soil loss values above 2109 t/ha/year were discarded from our analysis. GLs remained for the final analysis. Although they are not spatially and uniformly distributed, if their number is divided by the study area size, a density of 0.03 GLs/km 2 exists. Therefore, descriptive statistical values of maximum, minimum, and average erosivity factors (Ryear) for these GLs are presented (Table 5). year. Otherwise, the soil loss estimation with the average erosivity factor (Raver) using the average rainfall from 1985 to 2019 (Table 6)

Geodetic landmarks soil loss susceptibility
All the GLs were classified as null, low, moderate, severe, and very severe soil loss values for all the years of the study (Table 7)    To calculate the GLs soil loss susceptibility two considered different approaches were taken into account. One uses the average rainfall for the entire period between 1985 and 2019 (Raver) ( Table 7), and the other uses the specific rain for each year studied (Ryear) (Table 8). However, changes between them were not found in the Null class, while in the other classes, the maximum differences were 1.9% for the Low (1985 and, 1% for the Moderate (2015), 1.5% for the Severe (1995), and 1.8% for the Very Severe (2015) classes. It can be seen that more than half of the GLs are or were in regions without soil loss and a large part in low to moderate locations, with only a small amount in severe and very severe susceptibility areas. This is mainly due to the care taken during the GLs construction, considering its preservation over the years and the city's growth, especially urban expansion.
The gravimetric and altimetric GLs showed the highest soil loss, with 31% and 16% located within moderate, severe, or very severe susceptibility classes.
Otherwise, the planimetric's triangulation and polygonation classical methods to obtain horizontal coordinates and the GPS and doppler methods that got planimetric information by satellite tended to have low or no soil loss susceptibility since only 0.5% of them are in such locations.
In relation to GLs classified as good, not found, destroyed, or not built status in 2019, it was confirmed that all the planimetric's are classified as good.
At the same time, the gravimetric's GLs, where the measurement is made directly on the ground without the need to build a physical structure since the gravimetric information is essential for the equipotential surface determination, are classified as good or not built.
However, 22% of the GLs were classified as good conditions, 4% as destroyed, 22% as not found, and 52% as not built, all of them gravimetric. Except for the not built GLs, between the null and low soil loss classes, about 27% of them were classified as good, 49% as not found, and 8% as destroyed. In the moderate, severe, and very severe soil loss classes, 3% were classified as good, 12% destroyed, and 1% not found.
It is worth noting that the percentage between GLs classified as good and destroyed in the different soil loss classes does not seem to have a close relationship between soil loss susceptibility and the GLs disappearance.
However, the difference between the 49% of GLs not found in soils with null and low loss class and 1% for the moderate, severe, and very severe classes may indicate the destabilization of the structure due to soil erosion or even its burial.
Even though a null values increase occurred over the years and stability regarding the soil loss classes was found.
Comparing the results using Ryear and Raver, there were no differences for the null class. However, the average rainfall results were superior for the low and moderate classes, while the opposite happened in severe and severe soil susceptibility. Hence, it can be said that there tend to be less severe and very severe soil losses in the region, with more moderate, low, and even zero losses due to soil sealing by urban infrastructure.
Only a few GLs were found in mining areas, so this class had no significant influence on the results. However, coal mining has severe consequences for land resources and exerts tremendous pressure on the environment (FENG et al., 2019). The surface open-pit mining can increase the potential for sedimentation even in areas in an abandonment state due to soil exposure to rainfall (PIETROŃ et al., 2017), modifying the topography, removing vegetation, and exposing previously buried material (HOPKINS et al., 2013).
Possibly more GLs may be located in mining zones in the study area, which is one of the most important economic activities in the region. However, some boundaries between classes may be ambiguous due to the spatial resolution used, resulting in uncertainties. Especially in the LULC classification, there are very similar spectral resolution features between categories such as mining and other non-vegetated areas, such as bare soil (POWELL et al., 2004;SUN et al., 2020).

Future geodetic landmarks soil loss scenarios
Future rainfall estimations considering the climate scenarios in the Santa Catarina southern region are presented (Figure 7). In the RUSLE model, only the rainfall erosivity (R) was used considering four possible climatic scenarios for the future according to IPCC, for 2021IPCC, for -2040IPCC, for , 2041IPCC, for -2060IPCC, for , 2061IPCC, for -2080IPCC, for and 2081IPCC, for -2100 The rainfall average has a decrease between the observed value in 2019 and the forecast for the period 2021-2040. Then, the values increase for all scenarios, except in ssp126 scenario, which remains stable. By 2100, the rainfall average is projected to be higher than the observed in 2019, except for the ssp126 scenario, which is projected to remain slightly lower. The minimum rainfall values follow the trend of the average values. In contrast, the projected maximum rainfall is projected evermore to be higher than the observed in 2019, with the ssp126 scnenario remaining stable from 2040 and onwards, while the other scenarios continue to increase.  year already had the erosivity value 2% higher than 1985, but 13% lower than 1995, 5% lower than 2005 and 25% lower than 2015.
From the Rfut values for the four periods and the K, L, S, and CP values from 2019, the soil loss (Afut) graphs for the four future scenarios were generated for the Santa Catarina southern region (Figure 7). and 61 t/ha/year, the maximum and the average, respectively, the projected maximum values are always higher, on average 685 t/ha/year, more than 50% which is considered as severe (Table 3), with value of 2364 t/ha/year until 2100.
If historical series were considered, the soil loss value for 1985, 1995, 2005, and 2015 only be surpassed from 2061-2080 in the scenarios ssp245, ssp370, and ssp585. There are no differences between the minimum values. In contrast, there are variations for the average values for all the scenarios, ssp126 (7.8 t/ha/yr), ssp245 (7.2 t/ha/yr), ssp370 (7.7 t/ha/yr) and ssp585 (6.2 t/ha/yr). It is interesting to note that, for the period 2021-2040, the ssp585 presents the lowest soil loss, reflecting the lower rainfall projected for the period and consequently erosivity. However, from 2041-2060, lead to have the highest Afut values as expected.
The GLs soil loss susceptibility in the Santa Catarina southern region for 2021-2040, 2041-2060, 2061-2080, and 2081-2100 using ssp126 and ssp245 ssp370, and ssp585 scenarios is shown in (Figure 8). While for the reference year 2019, it has already been presented in (Figure 6).
Between 2021 and 2040, GLs soil loss values shift from moderate to low class, evidenced by an increase in the proportion in the low class and a decrease in the moderate class, including a slight reduction in the severe and very severe classes. Despite that, there was considerable stability in almost all scenarios for the other periods, demonstrating that the rainfall erosivity alone would not be the most important factor for the GLS class's severity.
We remind that in the future scenarios, only rainfall erosivity (R) was taken into account, and the factors of LULC (C) and conservation practices (P) also contribute significantly to soil loss but were not considered in this study and can be used in future studies. Moreover, the soil composition for the erodibility calculation (K) and the topographic factor (LS), although important, have much lower variability over the years than (C) and (P), and their variability can even be neglected for soil loss assessment in short periods of (A).

Conclusion
In this study, we analyze the geodetic landmark's propensity of instability as a function of the water erosion susceptibility of the soils on which they are established. This study is essential for determining the GLs' useful life belonging since the GLs' high implementation cost requires the longevity of these structures. The soil loss quantification is a good indicator, and the RUSLE equation is considered a consecrated method for this purpose. The combined factors of LULC and conservation practices were the most significant contributor to GLs null soil loss. At first, it seems to be a desirable situation. However, with urban growth, there is the destruction risk of existing GLs. In addition, there is also exposed soil in regions classified as urban infrastructure, where soil loss values are neglected in LULC maps, often due to their spatial resolution.
For future soil loss estimation, a trend towards increased rainfall and consequently erosion and soil loss was observed, with only the most optimistic scenario, the ssp126, remaining relatively stable. However, most GLs' were found in moderate and low susceptibility regions, with few in severe and very severe soil loss susceptibility areas, demonstrating that rainfall erosivity alone would not be the most essential factor for the GLs class's severity. This situation is due to the fact that the future soil loss projections took into account only rainfall and erosivity. For more accurate results, future LULC forecasts and studies on soil composition and relief changes would be necessary even though they are more stable and do not change in a short period.