Modeling Soil Erodibility Based on Silt Content: A Predictive Equation Proposal for Soils of Sorocaba (SP), Brazil ()
1. Introduction
Soil plays a vital role in the Earth's critical zone, serving as a dynamic interface between the atmospheric, hydrological, biological, and geological systems (Banwart et al., 2019). It supports biodiversity, regulates water and carbon cycles, and ensures structural stability in both urban and agricultural environments (Jackson et al., 2017). However, human activities, such as unplanned urban expansion, deforestation, and improper land use, have compromised these functions and intensified degradation processes, particularly erosion (Oliveira & Salles, 2020; Hernani, 2021).
Researchers define soil erodibility as the soil's inherent tendency to undergo particle detachment and transport. To quantify this property, they rely on attributes such as texture, organic matter content, structure, and permeability (Wischmeier & Smith, 1978; Wang et al., 2023). Among these, many studies have highlighted the silt fraction (%) due to its high mobility, low cohesion, and vulnerability to erosive forces, especially on sloped terrains with little vegetation (Keesstra et al., 2016; Bueno et al., 2022).
Several recent studies have demonstrated a strong correlation between silt content and erodibility, with coefficients exceeding 0.90. These findings suggest that soils with higher silt content are more vulnerable to degradation (Abiye & Dengiz, 2025). Li et al. (2024), for instance, showed that high silt levels reduce aggregate stability and promote fine particle loss during intense rainfall. Panagos et al. (2020) also observed that silt-dominated soils tend to form surface crusts and increase runoff, further intensifying erosion.
In urban areas, this issue becomes more critical due to soil sealing and unregulated land occupation. Vasconcelos et al. (2018) and Shojaeezadeh et al. (2024) identified higher erodibility values in peripheral urban zones, linking them to improper use of slopes and the lack of technical assessments. Moreover, when urban planners fail to integrate geotechnical and ecological factors, infrastructure becomes more fragile, and corrective projects become more expensive.
To address these challenges, researchers often combine physical soil characterization with predictive modeling techniques, including statistical regressions and pedotransfer functions (PTFs). These tools help estimate erodibility from easily obtainable variables, streamlining diagnostics and improving land-use planning (Veloso, Rodrigues, & Fernandes Filho, 2024). Recent approaches have also incorporated machine learning and advanced statistical techniques to improve predictions by considering land use and cover (Mallah et al., 2022).
Although methodological advances have emerged, studies that explore the relationship between silt (%) and erodibility in Brazilian urban settings remain limited, especially those that rely on field data and robust statistical models. This gap hinders the development of effective, region-specific soil conservation strategies.
In response, we aim to investigate how the silt fraction (%) influences soil erodibility by applying both simple and multiple regression models. Our objective is to construct pedotransfer functions that estimate the K factor based on textural variables. This approach provides a faster and more cost-effective alternative for diagnostics that support urban and environmental planning at the local level.
2. Study Area
The municipality of Sorocaba (São Paulo, Brazil), shown in Figure 1, covers 449.87 km2, with 77% of its territory classified as urban. The terrain is gently to moderately undulating, with elevations ranging from 550 to 750 m above sea level (Machado-Hess, 2020; IBGE, 2022a, 2022b).
The region is characterized by dynamic urban expansion, with land cover mosaic comprising urban, industrial, agricultural, and remnant native vegetation patches (MapBiomas, 2023). The climate is humid subtropical (Cfa), with an average annual temperature of 20.5˚C and mean precipitation of 1219 mm (Alvares et al., 2013; INMET, 2021). The original Atlantic Forest biome is heavily fragmented, currently restricted to riparian zones and legally protected areas (Kortz et al., 2014).
The dominant soils classes are Red Oxisols and Ultisols, formed from sedimentary rocks of the Paraná Basin and crystalline rocks of the Precambrian basement, resulting in notable physical and chemical variability (Sorocaba, 2024). However, the soils in Sorocaba are predominantly sandy-textured, with substantial silt fractions
Source: elaborated by the authors.
Figure 1. Geographic location of the study area in Sorocaba, São Paulo State, Brazil.
in some areas. This textural condition increases susceptibility to erosion, particularly in urban, peri-urban, and recently developed zones (Simonetti et al., 2018). The hydrographic network is structured by the Sorocaba River and major tributaries such as the Pirajibu, Ipanema, and Tatuí streams, which are critical for drainage and surface runoff (Nunes & Guandique, 2022).
We elected the municipality of Sorocaba to develop this study because Sorocaba is in an ecotone region, has a variety of soils compatible with the study, has a diversified land cover pattern and has little data about work (pedotransfer functions).
3. Materials and Methods
Soil erodibility was estimated using the EPIC model (Williams, Renard, & Dyke, 1983), which considers soil textural attributes and total organic carbon (TOC) content as input variables. Based on these data, correlation and regression statistical analyses were performed to develop pedotransfer functions capable of predicting erodibility from the silt fraction. The methodological steps included field sampling procedures, laboratory analyses, and statistical and geospatial data treatments.
We selected sixty-five sampling points across the municipality of Sorocaba (SP), distributing them based on land cover class stratification (MapBiomas, 2023) and performing random sampling within each stratum, following Zhao et al. (2023). The numbers of samples were: Natural Remnant Vegetation/Forest: 21 samples; pasture/grasslands: 24, bare soils: 9 samples; and urban settlements: 11 samples. We collected samples in a disturbed form from the soil surface layer (0 – 20 cm), placed them in plastic bags, and transported them to the laboratory for further analysis.
This study considers the influence of urbanization on the susceptibility of local soils to various forms of degradation, emphasizing that urban erosion is primarily driven by rainfall and wind in areas lacking vegetation and characterized by inadequate urban planning (Ferreira et al., 2018).
We conducted granulometric analysis of the surface samples using the adapted Bouyoucos method, which involves sieving and pipette sedimentation (Camargo et al., 2009). We defined the sand (0.05 - 2.00 mm), silt (0.002 - 0.05 mm), and clay (<0.002 mm) fractions according to the criteria established in the Soil Description and Field Sampling Manual (Santos et al., 2013). We classified the soils based on the IBGE textural triangle (IBGE, 2015).
We determined organic carbon content by loss on ignition at 440˚C for 3 hours, following the methodology described by ABNT (2017). We converted organic matter mass to TOC using the Van Bemmelen factor (1.724), assuming that organic matter contains approximately 58% carbon (Silva et al., 2023). We performed these analyses exclusively on surface layer samples.
We calculated soil erodibility of the surface layer using the EPIC model (Environmental Policy Integrated Climate), which is based on the Universal Soil Loss Equation (USLE), following Wischmeier & Smith (1978) and Williams, Renard, & Dyke (1983). Equation 1, used to calculate the K factor (KEPIC), directly considers the contents of silt, sand, clay, and total organic carbon (SOC).
(1)
Where:
KEPIC = soil erodibility (t·ha·h·ha−1·MJ−1·mm−1);
CLA = clay content (<0.002 mm) (%);
SIL = silt content (0.002 – 0.05 mm) (%);
SAN = sand content (0.05 – 2.0 mm) (%);
SOC = soil organic carbon content (%);
SN = 1 − (SAN/100);
0.1317 = conversion factor from metric to SI units.
Following Panagos et al. (2014), we included the very fine sand fraction (<0.05 mm) within the silt fraction.
We calculated the soil erodibility index (KEPIC) using data from the textural analysis and total organic carbon (TOC) content of surface soil samples. Next, we searched for statistically significant relationships between the physical and hydromechanical soil attributes to support the development of pedotransfer functions (PTFs) capable of estimating erodibility based on easily measurable variables, such as silt content (%) and TOC.
Due to the non-parametric nature of the data and the unequal sample sizes across categories, we chose for non-parametric statistical tests, one of which was Spearman’s rank correlation test, considering a significance level of 5% (p < 0.05). We performed statistical analyses with the assistance of BioEstat software version 5.3 (Ayres et al., 2007).
Additionally, we tested linear, logarithmic, exponential, and geometric regression models, using silt content (%) as the independent variable and the erodibility index (KEPIC) as the dependent variable. We assessed model fit using the coefficient of determination (R2), considering models with R² above 0.80 as adequate, due to their good predictive performance and applicability in erosion vulnerability studies.
We georeferenced the sampling points using high-precision GPS and integrated the data into a GIS environment. We conducted spatial analyses and generated maps using QGIS 3.36.0 (QGIS Development Team, 2024), applying Inverse Distance Weighting (IDW) interpolation to produce continuous representations of erodibility and silt content in the surface layer.
4. Results
We calculated the soil erodibility indices (KEPIC) at 65 sampling points distributed across the municipality of Sorocaba (SP) using data from the textural analyses and total organic carbon (TOC) contents of surface horizon samples. Based on these results, we conducted statistical analyses that revealed significant patterns of association between the physical soil attributes and its susceptibility to disaggregation.
Among the key findings, we highlight the strong correlation between the silt fraction (%) and the erodibility index, suggesting that this textural parameter acts as a sensitive predictor of the soil’s hydromechanical behavior in the region. Simple regression models based on the silt fraction demonstrated good statistical performance, with high coefficients of determination, underscoring the potential of these variables for developing pedotransfer functions (PTFs).
Among the models evaluated, the regression between silt content (%) and soil erodibility (K) demonstrated the greatest statistical robustness. Spearman’s analysis revealed a strong positive correlation (ρ = 0.8056; p < 0.0001; n = 65), indicating that an increase in the silt fraction is directly associated with higher susceptibility to erosion. Equation 2 expresses this relationship in a linear form.
(2)
The other regression models tested with variables from the surface horizon showed coefficients of determination below 80%, although some revealed statistical significance (p < 0.05). However, the low explanatory power of these equations limits their practical applicability. For this reason, these models were not considered suitable for developing pedotransfer functions (PTFs).
To visualize the spatial distribution of the analyzed properties, we generated maps using interpolation with the Inverse Distance Weighting (IDW) method at a 10-meter resolution. Figure 2 illustrates the spatial distribution of soil erodibility (K) in Sorocaba (SP), estimated by linear regression with the silt fraction (%) and the EPIC model.
Figure 2. Soil erodibility (K) maps in sorocaba: linear regression and EPIC model.
Figure 2 shows the spatial distribution of soil erodibility (K, in t ha h ha−1 MJ−1 mm−1) in the municipality of Sorocaba (SP), estimated using two different methods: on the left, the pedotransfer function based on the silt fraction (%) obtained by simple linear regression (K = 0.0138 × Silt + 0.0054); on the right, values simulated using the EPIC model. A general correspondence is observed between the spatial patterns, with areas of higher erodibility concentrated in the southern and eastern portions of the municipality. The regression estimates reflect the direct influence of soil texture, while the EPIC model incorporates additional factors such as land use and land cover.
The comparison between erodibility estimation methods allowed us to assess the accuracy and spatial behavior of the pedotransfer function formulated based on the silt fraction. For this purpose, we calculated the difference between the K values obtained from the linear regression equation and those generated by the EPIC model, shown in Figure 3. This approach aims to identify systematic patterns of underestimation or overestimation by the PTF in relation to the reference model, contributing to verifying its predictive robustness at the municipal scale.
Figure 3. Spatial difference between erodibility estimated by linear regression and the EPIC model (ΔK).
Figure 3 shows that the differences between erodibility values estimated by the two methods vary spatially across the territory of Sorocaba. The largest discrepancies are mainly concentrated in areas of the southeastern sector and isolated portions in the north, while central and western regions show closer values between the methods. The range of ΔK values indicates that, although there is a general correspondence between the models, some zones exhibit significant variations of the pedotransfer function compared to the EPIC model.
5. Discussion
Using textural analyses and total organic carbon (TOC) contents, we calculated the erodibility indices (KEPIC) at 65 surface horizon points in Sorocaba (SP). The correlations and regression models provide relevant insights into the mechanisms influencing soil susceptibility to erosion.
The strong Spearman correlation (ρ = 0.8056; p < 0.0001; n = 65) between silt content and the K factor highlights the importance of this fraction as an indicator of soil mobility and low cohesion, reinforcing results observed in different climatic contexts. The linear regression model between silt (%) and soil erodibility (K), proposed based on the sampled data, is represented by Equation 2:
(2)
The simple linear regression model (Equation 2) formulated from the silt content (%) showed a high explanatory power (R2 = 0.83), indicating that for each additional 1% of this fraction in the soil, erodibility (K) increases on average by 0.0138 t ha h ha−1 MJ−1 mm−1. This result reinforces the validity of the equation as an applicable pedotransfer function, especially in sandy soils where silt acts as a determining factor in structural instability. This strong association between silt and erodibility was also highlighted by Coelho et al. (2024), indicating that even small variations in the proportion of this fraction, known for its low cohesion and high susceptibility to mobilization, can result in significant increases in erosion vulnerability. Fernandes (2023) complements this understanding by demonstrating that the combination of fragile texture, reduced surface roughness, and lack of vegetation cover intensifies erosive processes, especially in areas with pronounced anthropogenic use or absence of conservation management.
The comparison between the spatial patterns of erodibility estimated by regression and by the EPIC model reveals significant agreement, especially in regions of higher susceptibility, indicating the applicability of the formulated pedotransfer function. The comparative maps show a general correspondence between the spatial distribution of erodibility estimated by the silt-based pedotransfer function and those obtained through the adapted EPIC model equation. However, localized discrepancies (ΔK) are observed, particularly in the southeastern and northern sectors of the municipality, where the values estimated by the two methods differ more markedly.
These variations do not arise from the inclusion of external factors in the EPIC model, such as slope, land use, or vegetation cover, since these variables were not considered in the present application. Both approaches were based on laboratory data, with the KEPIC equation incorporating textural fractions and total organic carbon content (TOC), while the pedotransfer function used only the silt fraction as a predictor.
These discrepancies largely reflect the inherent simplification of pedotransfer functions (PTFs), whose purpose is not to replace complex models but to offer a practical, accessible, and rapid alternative for preliminary erodibility estimates, especially in contexts constrained by time, resources, or analytical infrastructure. As highlighted by Weber et al. (2024), although statistical methods have advanced significantly, fundamental aspects such as applicability across diverse scenarios, scalability, and climatic suitability are often overlooked, underscoring the importance of strategies that balance simplicity and reliability.
As noted by Silva and Armindo (2016), the main premise of PTFs is that the parameters used in the equation should be easier to obtain than the attribute to be estimated, in this case, erodibility. Thus, PTFs serve as a technical resource enabling the extrapolation of point data based on available variables, optimizing the use of pedological information at operational scales.
This concept is supported by Van Looy et al. (2017), who state that the logic of PTFs lies in converting “data we have into data we need,” expanding possibilities for soil property characterization and modeling even in areas with low sampling density or lack of systematic monitoring.
The lower fidelity in specific areas is offset by operational agility, low implementation cost, and high replicability of the equation based solely on the silt fraction, an attribute easy to determine in the laboratory. In this context, the PTF proposed in this study can be considered a strategic tool for initial erosion vulnerability mapping, supporting territorial planning and decision-making in soil conservation projects, especially in areas dominated by sandy soils where the silt fraction plays a crucial role in superficial structural stability.
The results demonstrate that the pedotransfer functions (PTFs) formulated based on textural attributes and total organic carbon content (TOC) provide a practical and cost-effective alternative for estimating soil erodibility, especially in expanding urban areas with limited analytical data availability. However, these functions’ effectiveness depends on the scale of application and the soil’s surface layer, with reduced accuracy in environments exhibiting high heterogeneity in land use and cover.
In this regard, as argued by Weber et al. (2024), advancing PTFs requires incorporating spatial data, remote sensing, and machine learning techniques to overcome the limitations of simplified models and broaden their applicability across different edaphoclimatic contexts. Therefore, researchers recommend using PTFs as complementary tools to more comprehensive models like EPIC or integrating them into hybrid systems that combine laboratory variables with contextual attributes (Preetha & Joseph, 2025).
6. Final Considerations
This research demonstrated the applicability of pedotransfer functions (PTFs) as strategic tools for preliminary estimates of soil erodibility in expanding urban environments. The strong correlation between the silt fraction (%) and the KEPIC index (ρ = 0.8056; p < 0.0001) enabled the formulation of a linear regression model with high explanatory power (R2 = 0.83), validating the use of this variable as a predictor in simplified PTFs. The proposed model, based exclusively on silt content, offers technical and operational feasibility for initial diagnostics of erosion vulnerability, especially in sandy soils with low structural cohesion.
Comparison with estimates from the EPIC model showed overall spatial coherence, despite localized discrepancies. These differences arise from the simplified approach of the PTFs, which do not consider contextual variables such as land use and cover, but also reflect their main advantage: the possibility of application in scenarios with limited data and infrastructure. The philosophy of these functions, as emphasized by Silva & Armindo (2016) and Van Looy et al. (2017), relies on transforming easily obtainable variables into complex pedological information, expanding soil characterization potential even in low sampling density contexts.
Despite inherent formulation limitations, the PTF developed in this study exhibits high replicability, low cost, and practical applicability, making it useful for initial mapping, territorial planning, and conservation actions. The slight loss of precision in sectors with greater anthropogenic complexity is offset by methodological agility and efficiency, as also argued by Weber et al. (2024).
To overcome scalability limits and enhance model robustness, integrating PTFs with hybrid approaches that incorporate geospatial interpolation techniques, remote sensing data, and machine learning algorithms is recommended (Preetha & Joseph, 2025). This convergence of laboratory data and digital tools can boost the use of PTFs across diverse edaphoclimatic contexts, significantly contributing to environmental management and soil conservation.
Acknowledgements
The authors thank the Brazilian National Council for Scientific and Technological Development (CNPq) for financial support. Grant CNPq: 301955/2022.