Local calibration of bulk density models for agricultural soils in an inter-Andean valley of the Peruvian Central Highlands

Abstract

Soil bulk density (BD) is an essential indicator of soil quality, reflecting compaction, porosity, and water and nutrient retention capacity. However, its direct determination is costly and labor-intensive. This study developed and calibrated predictive models for BD in agricultural soils of an inter-Andean valley of the Peruvian Central Highlands using physicochemical properties (pH, electrical conductivity, organic matter, P, K, sand, silt and clay). Sixty-seven soil samples were analyzed using five modeling approaches: multiple linear regression (MLR), partial least squares regression (PLSR), artificial neural networks (ANN), classification and regression trees (CART), and random forest (RF). model performance was evaluated through fivefold cross-validation and goodness-of-fit metrics (coefficient of determination, R²; root mean square error, RMSE; mean absolute error, MAE; and performance to deviation, RPD). The RF model gave the best overall fit among the evaluated models (R² = 0.621, RMSE = 0.0862 g·cm⁻³, MAE = 0.0693 g·cm⁻³, RPD = 1.64), outperforming MLR (R² = 0.614), PLSR (R² = 0.608), ANN (R² = 0.584) and CART (R² = 0.514). Local calibrations significantly reduced errors compared to standard pedotransfer functions, which typically have standard errors of 0.14–0.20 g·cm^–3 and R² < 0.60. Results demonstrate the potential of local models, particularly machine learning approaches, to generate rapid, low-uncertainty BD estimates supporting soil fertility management and sustainable agricultural planning.

Graphical abstract

Keywords

Bulk density / local calibration / machine learning / pedotransfer functions / soil properties

Highlight

	● Local calibration markedly improves bulk density prediction accuracy.
	● Simple linear models match complex algorithms with limited sample sizes.
	● Global pedotransfer functions gave reduced accuracy outside their original domain.
	● Soil properties matter more than algorithm complexity for soil bulk density estimation.

1 Introduction

Soil bulk density (BD) is a key indicator in agricultural systems, as it reflects soil compaction and porosity, which directly determine productive capacity^[1,^2]. By describing soil particle packing, BD allows inference of physical properties associated with the movement and availability of water and air, both essential for crop growth and moisture retention^[3]. It also directly influences biological properties such as microbial populations, enzymatic activity, soil fauna and root development^[4]. Additionally, BD is essential for estimating nutrient stocks, including potassium and phosphorus, as it enables conversion of weight-based nutrient concentrations to volume- or area-based stock data^[5]. BD depends on texture, mineralogy, soil organic carbon and management practices, making it a central parameter for assessing soil quality, fertility and nutrient availability, thereby guiding sustainable management strategies in agroecosystems^[6,^7]. However, established direct measurements of BD are often slow, costly and labor-intensive, highlighting the need for more efficient predictive models^[4].

Given these constraints, numerous studies have developed prediction methods based on pedotransfer functions (PTFs)^[8,^9]. These approaches have evolved from simple linear or nonlinear regression models to more advanced data analysis techniques, driven by increased computational capacity and the development of machine learning algorithms^[6]. Several methods have been proposed to improve the accuracy and applicability of BD estimation under different soil conditions and agricultural systems^[4]. However, these models often perform inconsistently across different conditions, restricting their application in soil and crop management at regional scales, particularly regarding compaction assessment, hydrothermal regulation, carbon storage, hydraulic properties and soil erosion modeling^[10].

Recent studies have shown that the predictive performance of PTFs is strongly context-dependent, varying across soil types, regions, and calibration domains, which can limit their reliability when applied beyond the conditions under which they were developed^[4,⁷,¹¹,^12]. Although traditional empirical models have enabled BD prediction from soil physical and chemical properties, their validity is limited outside their calibration areas, leading to substantial over- or underestimations^[13]. This variability stems from the strong dependence of BD on soil structural factors and variables such as organic matter (OM) content, pH, electrical conductivity (EC), and textural fractions^[6]. This challenge is particularly relevant in inter-Andean agricultural valleys of the Peruvian Central Highlands, where soils are commonly represented by Fluventic and Ustorthent groups, which may have distinct BD responses compared to soils in other regions. Consequently, the development of locally calibrated models has become essential for improving prediction accuracy and regional applicability.

Machine learning methods have shown significant improvements in BD prediction by capturing complex nonlinear relationships among soil variables. Algorithms such as random forest (RF) and artificial neural networks (ANN) have outperformed traditional models in various contexts. For example, a study in Brazil found that PTFs developed with RF predicted BD more accurately than simple regression models^[14], while research in China demonstrated that a radial basis function Neural network improved BD estimation in agricultural basins^[11]. Incorporating region-specific data has proven essential for improving BD prediction. In Europe, machine-learning-based PTFs have been used to generate topsoil BD estimates from harmonized soil databases, highlighting the value of regional calibration and covariate integration^[3], while research in mountainous soils showed that a specific PTF combining ANN and stepwise regression markedly reduced errors compared to reference models^[4].

In South America, studies on BD have applied a wide range of approaches, from established techniques to advanced predictive modeling. Direct volumetric cylinder measurements have been used to relate BD with soil physical and chemical properties in Colombia^[15] and to estimate BD and calculate carbon storage in cacao systems in Ecuador^[16]. In Brazil, beyond direct measurements, PTFs based on simple and multiple linear regression have been developed, incorporating textural fractions and OM as explanatory variables^[17,^18]. More recently, nonparametric regressions and machine learning algorithms such as RF, support vector machines, classification and regression trees (CART) and gradient boosting have been implemented, providing improvements in both accuracy and adaptability to heterogeneous edaphoclimatic conditions^[6,¹⁹,^20]. These advances reflect a regional transition from traditional methods toward predictive modeling techniques that aim to overcome the limitations of empirical approaches in terms of representativeness and regional transferability.

Despite this regional progress, studies on BD in Peru remain scarce and have relied mainly on standard methods, limiting the precision and applicability of results under local conditions^[21,^22]. The unique edaphoclimatic conditions and diverse agricultural systems of the Central Highlands require locally adapted predictive tools. Addressing this gap requires combining standard statistical approaches such as linear regression, multiple linear regression (MLR) and partial least squares regression (PLSR) with machine learning techniques like decision trees (such as CART) and ANN, which have demonstrated the ability to capture nonlinear relationships and improve prediction accuracy when combined with cross-validation strategies^[19].

Therefore, this study aimed to develop locally calibrated predictive models for BD in agricultural soils from an inter-Andean valley of the Peruvian Central Highlands and to compare their performance with previously published models. Using a set of edaphic predictors including pH, EC, OM, available phosphorus, available potassium and textural fractions, both standard statistical approaches (including MLR and PLSR) and machine learning algorithms (including CART, ANN and RF) were developed and evaluated. The locally calibrated models were then compared with PTFs reported in the literature to assess whether local calibration improves prediction accuracy under the studied conditions. This research seeks to identify the most accurate and robust approaches for BD estimation in this local agricultural setting, providing practical tools adapted to the local edaphoclimatic reality to strengthen agricultural planning and sustainable soil management.

2 Materials and methods

2.1 Study area

The study was conducted at the Santa Ana Agricultural Experimental Station in the Mantaro Valley, El Tambo district, Huancayo Province, Junín Region (12°00′ S, 75°13′ W and 3245 m.a.s.l.). This inter-Andean basin is characterized by a temperate-cold climate with mean annual temperatures of 11–14 °C, a rainy season from October to April, and a prolonged dry season from May to September^[23].

This valley represents one of the main agricultural areas of the Peruvian Central Highlands, supporting intensive production of staple crops such as potato, maize, beans and peas, as well as high-value vegetables and forages including carrot, onion, broccoli, oats and alfalfa, which supply both regional and national markets^[24,^25]. Agricultural soils in the study area are predominantly associated with Fluventic and Ustorthent soil groups, reflecting depositional and geomorphological processes typical of inter-Andean valleys. Observed BD ranged from 1.15 to 1.71 g·cm^–3, indicating moderate spatial variability under comparable valley-scale edaphic conditions and intensive agricultural management.

2.2 Soil sampling, sample analysis and quality assurance

Sixty-seven undisturbed soil samples were collected at a fixed depth of 30 cm using galvanized steel cylinders (height 7 cm, diameter 5.3 cm and volume 154 cm³). Each sampling location was georeferenced using a D-RTK 2 GNSS system (DJI, Shenzhen, China). BD was determined at the Soil, Water and Foliar Laboratory of Santa Ana using the core method^[26]. Samples were oven-dried at 105 °C until constant weight and BD was calculated as the ratio of oven-dry soil mass to cylinder volume:

(1)BD=\frac{m_s}{V_c}

where, BD is the bulk density (g·cm^–3), m_s is the oven-dry soil mass at 105 °C (g) and V_c is the cylinder volume (cm³).

For physicochemical analysis, soil samples were air dried at room temperature, mechanical disaggregated and sieved to 2 mm. Soil pH and EC determination was performed using potentiometric and conductimetric techniques, respectively, at 1:2.5 soil:water ratio, according to US EPA Method 9045D and ISO 11265:1994 standard. OM was quantified by the Walkley-Black wet oxidation method. Available phosphorus was determined by extraction with 0.5 mol·L^–1 NaHCO₃ at pH 8.5 following the Olsen method. Exchangeable potassium was extracted using 1 mol·L^–1 ammonium acetate at pH 7 and subsequently quantified by flame photometry. Soil texture was characterized applying the Bouyoucos hydrometer method, using sodium hexametaphosphate as dispersing agent. Unless otherwise indicated, all analytical procedures were conducted following SEMARNAT (2002) guidelines^[27–^30].

2.3 Reference models for estimating apparent density

Locally calibrated predictive models were developed using both standard statistical methods MLR and PLSR and machine learning algorithms CART, ANN and RF. Model inputs included soil pH, EC, OM, available phosphorus and potassium and textural fractions (sand, silt and clay).

To establish a comprehensive benchmark, 66 published PTFs for BD estimation were compiled from the literature (Table S1). These reference models, based on OM, texture and depth, were developed under diverse pedoclimatic conditions spanning temperate, tropical, and arid regions. All models both locally calibrated and published PTFs were evaluated using root mean square error (RMSE), coefficient of determination (R²), and ratio of performance to deviation (RPD). This comparative framework enabled assessment of whether local calibration improves BD prediction accuracy relative to existing PTFs when applied to soils from the studied inter-Andean valley.

2.4 Model development and evaluation

All statistical analyses and predictive modeling were conducted in R version 4.3.2^[31]. Exploratory analysis included Pearson's correlation matrices visualized with corrplot^[32] and principal component analysis (PCA) performed using FactoMineR^[33] and visualized with factoextra^[34] to identify multivariate patterns among soil properties.

Five regression models were developed to predict soil BD using physicochemical properties as predictors. MLR established baseline linear relationships between soil properties and BD^[35], and PLSR addressed multicollinearity by maximizing variance through latent components^[36,^37]. Three machine learning algorithms captured nonlinear patterns: CART generated interpretable decision rules through recursive partitioning^[38–^40], ANN modeled complex interactions via interconnected neuron^[4], and RF used ensemble learning to improve prediction stability^[14].

All models were implemented using the caret package in R, with preprocessing based on centering and scaling of predictors. To ensure reproducibility of the resampling and training procedures, a fixed random seed was established in R using set.seed(123) prior to cross-validation and hyperparameter tuning, allowing the same data partitions, resampling process, and tuning results to be reproduced. Model calibration used fivefold cross-validation with hyperparameter optimization through grid search. For ANN, tested hyperparameters included hidden layer sizes (1, 3, 5 and 7 neurons) and weight decay penalties (0, 1×10^–4, 1×10^–3 and 1×10^–2). CART optimization evaluated complexity parameter values from 0.001 to 0.1 in increments of 0.001. PLSR automatically selected the optimal number of components between 1 and 15. For MLR, individual predictors were tested and the best performer based on RMSE was selected. RF models were trained using 500 trees, testing the number of predictors sampled at each split from 2 to the total number of predictors minus one, as defined by the tuning grid. Variable importance was computed to assess predictor relevance.

Model performance was evaluated using four metrics: coefficient of determination (R²) for variance explained, RMSE for prediction accuracy, mean absolute error (MAE) for average prediction deviation, and ratio of performance to deviation (RPD) for model reliability. The optimal model was selected based on simultaneously maximizing R² and RPD while minimizing RMSE.

3 Results

3.1 Statistical analysis

Table 1 gives the descriptive statistics for the 67 soil samples from agricultural plots. Soil pH (5.99 ± 0.50) and BD (1.46 ± 0.14 g·cm^–3) had low variability (CV < 10%), indicating uniform acidity and compaction. Available P and exchangeable K had high variability (CV = 76.1% and 47.2%, respectively), with strongly skewed distributions (skewness 1.44 and 0.60) reflecting heterogeneous nutrient management. OM (2.45 ± 1.10%, CV = 44.8%) and EC (2.73 ± 0.89 dS·m^–1, CV = 32.5%) had moderate to high variability. Particle size distribution averaged 41.8% sand, 32.1% silt, and 26.0% clay, corresponding to loam texture with moderate variability (CV = 23.4% to 25.8%).

3.2 Soil property interactions and influence on bulk density

Figure 1 shows the PCA biplot and soil textural distribution according to USDA classification. In Fig. 1(a), the first two principal components explain 60.1% of total variance, with PC1 (41.2%) representing the textural gradient and PC2 (18.9%) reflecting chemical variability. The BD vector approximately parallels that of sand content, indicating that soils with higher sand proportions have higher BD values due to lower porosity and greater compaction. Conversely, the clay, silt and OM vectors indicate negative correlations with BD consistent with finer-textured, organic-rich soils having lower BD.

Phosphorus has moderately positive association with BD, while potassium and pH have short vectors near the horizontal axis, indicating minimal contribution to variance. EC, positioned near the origin along PC2, indicates marginal influence on soil differentiation.

Figure 1(b) reveals samples predominantly classified as loam, sandy clay loam, and sandy loam, with some clay and clay loam. This distribution aligns with the PCA pattern, where sandier soils associate with higher BD and finer soils with elevated OM have lower compaction.

3.3 Predictive model performance for bulk density estimation

Figure 2 shows the comparison between observed and predicted BD values obtained through the five modeling approaches: ANN, CART, MLR, PLSR and RF. Each panel includes a red fitted regression line between observed and predicted values.

The results demonstrate acceptable and consistent performance across all approaches, with coefficients of determination (R²) ranging from 0.514 (CART) to 0.621 (RF). The RF, MLR and PLSR models provided the best fits, with R² values of 0.621, 0.614 and 0.608, respectively, and low prediction errors. Specifically, RF showed the lowest errors (MAE = 0.0693 g·cm⁻³; RMSE = 0.0862 g·cm⁻³), followed by MLR (MAE = 0.0710 g·cm⁻³; RMSE = 0.0873 g·cm⁻³) and PLSR (MAE = 0.0707 g·cm⁻³; RMSE = 0.0877 g·cm⁻³), demonstrating adequate predictive capacity. Notably, the RF model recorded the highest R² (0.621), indicating greater stability and precision in BD estimation.

In contrast, the CART model gave the lowest performance (R² = 0.514; MAE = 0.0779 g·cm⁻³; RMSE = 0.0987 g·cm⁻³), with a clear tendency to underestimate high BD values relative to the 1:1 reference line. The ANN model, although slightly outperforming CART (R² = 0.584; MAE = 0.0723 g·cm⁻³; RMSE = 0.0911 g·cm⁻³), had greater dispersion in predictions, particularly at the extremes of the observed range.

Table 2 shows the performance of the five modeling approaches applied to predict BD from physicochemical properties. Of these, the RF model gave the best performance, with a coefficient of determination (R²) of 0.621, an RMSE of 0.0862 g·cm⁻³, a MAE of 0.0693 g·cm⁻³, and an RPD of 1.64, reflecting the highest predictive capacity and good model stability against data variability.

The MLR model, which incorporated all textural variables, provided similar performance (R² = 0.614), although with slightly higher errors (RMSE = 0.0873 g·cm⁻³; MAE = 0.0710^** g·cm⁻³; RPD = 1.62^**). Likewise, the PLSR model, which also integrated all variables, reported comparable metrics (R² = 0.608; RMSE = 0.0877^** g·cm⁻³; MAE = 0.0707 g·cm⁻³; RPD = 1.61^**), indicating that the relationships between texture, OM and BD can be adequately captured through multivariate linear approaches.

Regarding machine learning models, the ANN approach gave an R² of 0.584, with only slightly higher errors (RMSE = 0.0911 g·cm⁻³; MAE = 0.0723 g·cm⁻³; RPD = 1.55), providing no substantial improvements over linear methods despite its greater structural complexity. In contrast, the CART model yielded the lowest overall performance (R² = 0.514; RMSE = 0.0987 g·cm⁻³; MAE = 0.0779 g·cm⁻³; RPD = 1.43), reflecting a limited generalization capacity for this dataset.

The RPD values ranged from 1.43 (CART) to 1.64 (RF). RMSE ranged from 0.0862 g·cm⁻³ (RF) to 0.0987 g·cm⁻³ (CART), indicating acceptable precision for practical applications.

3.4 Comparative performance of models in soil bulk density estimation

Figure 3 compares the performance of the five locally calibrated models developed in this study against 66 published PTFs for BD prediction. Although RF gave slightly better performance under internal cross-validation, with marginally lower RMSE and higher R² and RPD values, the differences from MLR and PLSR were small across all metrics and do not establish a clear superiority of ensemble methods at this sample size. MLR represents a more practical and interpretable alternative for operational applications, offering comparable predictive capacity with greater simplicity. ANN and CART gave slightly lower performance, although still within an acceptable range relative to many published functions.

Compared with the published models, the locally calibrated approaches developed here ranked among the best-performing models. While some literature-based equations gave comparable or even higher values for individual metrics, many published PTFs gave lower R² and RPD values than the locally calibrated models. Overall, this comparison indicates that local calibration can provide competitive predictive performance relative to models developed under broader and contrasting edaphoclimatic conditions. Model equations are provided in Table 2.

4 Discussion

4.1 Relationship of bulk density with soil properties

BD had a significant association with soil OM and with textural fractions, particularly sand, confirming patterns widely documented in the literature. A negative correlation was observed between BD and OM, reflecting the role of organic inputs in decreasing compaction and increasing soil porosity. Similar results have been reported in Amazonian soils^[18], in different biomes in Brazil^[17] and in soils from China^[41], where OM was identified as the main explanatory factor of BD. The logarithmic and exponential functions proposed in these studies further support this trend^[41,^42], indicating that OM acts as a dominant predictor across diverse soil environments. Nevertheless, the results also show that the magnitude of this effect depends on local conditions, which explains why global PTFs often underestimate or overestimate BD when applied outside their original calibration domain^[43,^44].

Textural properties also had a significant influence on the variability of BD, resulting a positive association with the sand fraction and a negative one with clay. This finding is consistent with the formulation of Adams^[45], who developed one of the first physically based equations linking BD with mineral and organic fractions. More recent studies have emphasized the interaction between clay and OM as key determinants of BD^[17,^18]. Our results support this interaction, as soils with higher clay and OM contents consistently had lower BD values. Likewise, contemporary approaches, such as the reference BD framework proposed by Keller and Håkansso^[46] and the textural entropy metric developed by Martín et al.^[47], reinforce the idea that both the proportion and the heterogeneity of textural fractions condition soil packing and, consequently, BD. In combination, these findings demonstrate that locally calibrated models integrating OM and texture provide more robust and reliable predictions than global functions, which is consistent with evaluations conducted in Mediterranean and alpine systems^[44,^48].

In our study, sand was identified as the main predictor of BD, having a positive relationship, whereas soil OM exerted a negative effect but of smaller magnitude. As shown in Fig. S1, variable importance across all models was highest for sand (100 on relative scale), whereas OM had a negligible contribution (0 on that scale) in all models evaluated. This pattern differs from several studies in which OM makes a stronger contribution to BD prediction^[17,¹⁸,^41], and can be attributed to the specific characteristics of Andean soils in our study area, which are characterized by low OM levels and high sand content, indicating that BD variability in this inter-Andean agricultural context is more strongly driven by textural composition than by OM. Although the negative direction of the OM effect remains consistent with the general principle that OM reduces compaction and increases soil porosity^[42], its relative importance is diminished in these sand-dominated, OM-poor soils, where coarser particle arrangements appear to govern packing efficiency and pore space distribution more decisively.

4.2 Performance of the developed models

The RF model achieved the best overall performance, followed closely by MLR and PLSR, which produced nearly equivalent results despite their differing approaches to dimensionality; full metric details are given in Table 2. These findings demonstrate that both multiple linear regression and dimension reduction techniques can provide robust BD predictions. ANN performed similarly to the linear models without significant improvements, which does not offset its greater complexity and dependence on larger datasets for adequate generalization^[10,^48]. Similarly, nonlinear approaches such as ANN and CART did not provide notable improvements, with CART recording the lowest fit among all models (Table 2), evidencing the limited generalization capacity of this algorithm in small datasets^[49].This behavior is consistent with international assessments highlighting that neural networks and tree-based algorithms require large and heterogeneous databases to fully express their predictive potential^[6,^19].

The efficiency of the RF model is explained by its ability to capture nonlinear and complex relationships among soil properties, automatically identifying the variables with the greatest influence on prediction in this case, OM and sand contents, which show opposite effects on BD^[18,⁴⁵,⁵⁰,^51]. This relationship pattern has been reported in multiple edaphoclimatic contexts, although the magnitude of its effects varies by region, reinforcing the need for local calibrations^[17,⁴²,^44]. Similarly, recent reviews emphasize that texture and OM remain the main axes in BD prediction and that the advantages of advanced models, such as boosting or deep learning, only materialize when large and diverse datasets are available^[6,^52]. As shown in Fig. S2, RF had consistently lower error metrics (RMSE and MAE) and higher RPD values than ANN across the fivefold cross-validation scheme, while R² had wider variability in both models, indicating that overall advantage the RF model was not driven by a single favorable partition but reflected a more stable predictive performance under the evaluated resampling conditions.

Several studies support the outstanding performance of RF in estimating BD and other soil physical properties. For example, Moreira et al.^[53] demonstrated that spectral-based models, such as RF, can accurately predict BD in tropical systems, while Bondi et al.^[38] and He et al.^[54] showed that RF improves the spatial representation of textural fractions and BD in regions with high edaphic heterogeneity. Likewise, Ali et al.^[55] reported that integrating environmental and textural covariates through RF significantly reduced estimation errors compared to linear methods, and Bouslihim et al.^[56] highlighted its stability when working with small datasets and its ability to handle nonlinear relationships among soil variables.

4.3 Comparison with previous models and pedotransfer function

The comparison of the developed models with established PTFs reveals that local calibrations offer clear advantages by reducing prediction errors and improving model reliability. Widely used equations, such as those proposed by Adams^[45], Bernoux et al.^[18] and Rawls et al.^[57], as well as the functions of Alexander^[58] and Curtis and Post^[59], established foundational relationships between BD, OC and texture. However, these global or regional functions often generate biases when applied outside the edaphoclimatic context in which they were developed. De Vos et al.^[10] confirmed that many PTFs tend to underestimate BD in surface horizons and recommended specific adjustments for each region or soil type, while Heuscher et al.^[51], working with more than 47,000 NRCS samples, demonstrated that the magnitude of relationships between organic carbon and BD varies notably across soil orders, limiting the applicability of universal equations.

The comparative analysis presented in Fig. 3 provides important insights into model performance beyond simple R² comparisons. While several published PTFs giving R² values within the range obtained in this study (0.50–0.62), many have notably higher RMSE values (0.14–0.20 g·cm^–3) and lower RPD coefficients. For example, models, such as Calhoun et al.^[60], Benites et al.^[61], and several functions derived from regional databases show acceptable R² values but fail to provide RPD values above 1.5, which is considered the minimum threshold for reliable predictions in soil science applications^[62]. This pattern underscores a critical limitation of relying solely on R² as a performance metric. Although this coefficient measures the proportion of variance explained, it does not adequately capture the magnitude of prediction errors or model stability across different data ranges^[63,^64].

The superiority of locally calibrated models becomes evident when considering multiple performance metrics jointly. A lower RMSE indicates smaller absolute deviations between predicted and observed values, essential for practical applications requiring precise BD estimates for carbon stock calculations or compaction assessments^[50,⁵⁹,^65]. Similarly, higher RPD values reflect greater model robustness and reduced sensitivity to data variability, making predictions more reliable under field conditions^[62]. The models developed in this study gave RPD values between 1.43 and 1.64, with RMSE values below 0.10 g·cm^–3, positioning them among the most accurate PTFs available for agricultural soils. In contrast, the models locally calibrated had lower errors (RMSE = 0.0862 g·cm^–3) compared to global functions, demonstrating that even simple equations based on sand and OM can outperform established PTFs when adapted to local conditions.

These findings align with recent observations emphasizing the limitations of global PTFs in specific agricultural contexts, while also reflecting the constraints inherent to locally calibrated models, such as the absence of mineralogical data and standardized management records, which in regions with volcanic and alluvial parent materials like the Mantaro Valley may represent relevant sources of BD variability not yet captured by the current predictors. Calhoun et al.^[60] and Botula et al.^[66] demonstrated that models based solely on texture and OM tend to overestimate BD in surface horizons and underestimate it in deeper horizons when applied across diverse environments. Bryk & Kołodziej^[6] further emphasized that global PTFs lose accuracy in region-specific applications. Therefore, the generation of calibrated functions in the Central Highlands not only improves the quality of estimates but also provides more reliable tools for agricultural planning and sustainable soil management tools that capture local edaphic particularities rather than depending on foreign models developed under different conditions. This reinforces the principle that local calibration, even with relatively simple models like MLR, can outperform global PTFs not only in explained variance but also in prediction precision and stability critical attributes for operational soil management^[67].

4.4 Practical implications and projection

The availability of locally calibrated models for estimating BD has direct implications for agricultural planning and soil conservation, since this parameter is indispensable for converting concentrations into carbon and nutrient stocks, as well as for assessing compaction and water retention capacity^[50,⁵⁹,^65]. Unlike standard PTFs based on universal approaches^[45,⁵⁷,^58], the functions calibrated significantly reduce errors, thereby increasing the reliability of carbon and nutrient stock estimation under different management scenarios. This allows for the generation of rapid, accessible and low-cost diagnostics for local farmers, while also contributing to the formulation of sustainable practices aimed at improving fertility and mitigating soil degradation processes. In this sense, the evidence aligns with recent studies that emphasize the need for PTFs adapted to regional conditions for effective fertilizer application strategies and land-use planning^[6,¹⁰,^51].

From an operational standpoint, it is important to emphasize that the MLR model, despite its slightly lower R² compared to RF, represents the most practical option for widespread adoption. The minimal difference in performance (ΔR² = 0.007; ΔRMSE = 0.001 g·cm^–3) does not justify the additional complexity of implementing ensemble machine learning algorithms, particularly in contexts where end-users lack access to specialized programming environments. The explicit mathematical formulation of MLR enables straightforward application using basic computational tools, facilitating its integration into routine soil assessments and extension programs. Also, the transparency of linear coefficients enhances user confidence and allows for direct interpretation of the influence of each soil property on BD predictions, which is essential for training and capacity-building initiatives. This aligns with recommendations from De Vos et al.^[10] and Botula et al.^[66], who advocate for simple, interpretable PTFs that can be readily adopted by diverse users without sacrificing prediction quality.

In line with previous studies on PTFs and model validation, the results confirm that local calibration is a key determinant of predictive performance, often more influential than algorithmic complexity. Recent studies indicate that advanced models, such as RF, can provide high accuracy when trained on large and representative datasets; however, under moderate sample sizes, more parsimonious approaches frequently give comparable performance and greater stability^[7,^68]. In addition, methodological studies emphasize that cross-validation constitutes an appropriate internal calibration procedure but may yield optimistic estimates if not interpreted with caution, particularly when spatial structure is present and independent validation datasets are unavailable^[62,^69]. Within this framework, model stability is understood as the consistency of predictive performance across cross-validation folds, as reflected by limited variability in RMSE, R² and RPD, rather than as evidence of external generalization^[63]. Accordingly, the marginal performance advantage of RF over linear models, such as MLR, observed in this study aligns with the well-documented trade-off between predictive accuracy, interpretability, and operational applicability. Finally, the comparison with a large set of published global PTFs confirms that extending models beyond their original edaphoclimatic domain tends to reduce predictive accuracy, highlighting the balance between spatial generalization and local precision.

Nevertheless, the scope of the models still presents limitations arising from the sample size and the spatial concentration within a single study area. Also, BD was assessed at a single soil depth, which precludes the analysis of its vertical stratification within the profile and the identification of potential subsurface compacted layers. Likewise, the moderate predictive performance obtained here indicates that BD variability may not be fully explained by the physicochemical properties included in the models alone, but also by additional factors not evaluated in this study, such as mineralogical composition and soil management history. Previous research has documented that the relationship between BD, texture and OM varies significantly across tropical, temperate and boreal environments^[18,⁴⁹,⁶⁵,^70], which underscores the importance of expanding the empirical basis to other agricultural systems and altitudinal zones of the Andes and the Amazon. In addition, contemporary studies have shown that integrating edaphic data with climatic, topographic and management covariates, together with the use of machine learning algorithms, can substantially improve the robustness of predictions^[3,⁴,^19]. Consequently, future research should prioritize depth-stratified sampling at 0–10, 10–20, and 20–30 cm increments to better characterize vertical BD variability and potential compacted layers, as well as replication across other inter-Andean agricultural settings and coupling with digital soil mapping techniques, which may support the development of more transferable functions and the consolidation of practical tools for soil resource management under comparable edaphoclimatic conditions.

5 Conclusions

This study developed and calibrated five predictive models for BD estimation in agricultural soils of an inter-Andean valley of the Peruvian Central Highlands, demonstrating that locally calibrated functions provide competitive predictive performance relative to global PTFs in terms of accuracy and stability. The key finding is that local calibration, rather than algorithmic complexity, determines predictive success, as simple linear models adapted to the studied local conditions can outperform more sophisticated global PTFs. Although the RF model provided marginally superior performance under internal cross-validation, the MLR model is recommended for operational applications due to its comparable predictive capacity, greater interpretability, and practical simplicity. The small performance differences observed across models do not justify additional implementation complexity. This makes MLR a practical alternative for operational use, particularly where model transparency and ease of implementation are required. The multi-metric evaluation further demonstrated that relying solely on R² as a performance criterion is insufficient. Several published PTFs give comparable explained variance but show substantially higher prediction errors and greater variability across performance metrics, indicating reduced robustness. In this study, model stability is understood as consistency of predictive performance across cross-validation folds, reflected by limited variability in RMSE, R², and RPD. Future research should extend model validation across contrasting soils, management conditions and elevation gradients to improve transferability. Integrating climatic, topographic, and management covariates, together with digital soil mapping approaches, may improve model generalizability while preserving local accuracy, thereby supporting the progressive development of regionally informed PTFs for sustainable soil resource management in Peru.

References

[1] Panagos P, De Rosa D, Liakos L, Labouyrie M, Borrelli P, Ballabio C . Soil bulk density assessment in Europe. Agriculture, Ecosystems & Environment, 2024, 364: 108907

[2] Zheng G H, Jiao C X, Xie X L, Cui X F, Shang G, Zhao C Y, Zeng R . Pedotransfer functions for predicting bulk density of coastal soils in East China. Pedosphere, 2023, 33(6): 849–856

[3] Chen S C, Chen Z X, Zhang X L, Luo Z K, Schillaci C, Arrouays D, Richer-de-Forges A C, Shi Z . European topsoil bulk density and organic carbon stock database (0–20 cm) using machine-learning-based pedotransfer functions. Earth System Science Data, 2024, 16(5): 2367–2383

[4] Bashir O, Ahmad Bangroo S, Amjid S, Nasta P, Shafai S, Palladino M, Senesi N, Rasool R, Choudhury B U, Jha P K, Abrol V, Spalevic V, Sestras P, Omidvar N, Ercişli S, Jaufer L, Benzougagh B, Bammou Y, Mohamed K, Lukić T, Kader S . Soil bulk density estimation by a novel pedotransfer function tailored on hilly terrains. Earth Systems and Environment, 2025, 9(3): 1615–1633

[5] Walter K, Don A, Tiemeyer B, Freibauer A . Determining soil bulk density for carbon stock calculations: a systematic method comparison. Soil Science Society of America Journal, 2016, 80(3): 579–591

[6] Bryk M, Kołodziej B . Pedotransfer functions for estimating soil bulk density using image analysis of soil structure. Sensors, 2023, 23(4): 1852

[7] Palladino M, Romano N, Pasolli E, Nasta P . Developing pedotransfer functions for predicting soil bulk density in Campania. Geoderma, 2022, 412: 115726

[8] Al-Qinna M I, Jaber S M . Predicting soil bulk density using advanced pedotransfer functions in an arid environment. Transactions of the ASABE, 2013, 56(3): 963–976

[9] de Souza E, Fernandes Filho E I, Schaefer C E G R, Batjes N H, dos Santos G R, Pontes L M . Pedotransfer functions to estimate bulk density from soil properties and environmental covariates: Rio Doce basin. Scientia Agricola, 2016, 73(6): 525–534

[10] De Vos B, Van Meirvenne M, Quataert P, Deckers J, Muys B . Predictive quality of pedotransfer functions for estimating bulk density of forest soils. Soil Science Society of America Journal, 2005, 69(2): 500–510

[11] Li A W, Cheng J L, Chen D, Li W D, Mao Y R, Chen X Y, Zhao B, Shi W J, Yue T X, Li Q Q . Spatial interpolation of cropland soil bulk density by increasing soil samples with filled missing values. Scientific Reports, 2025, 15(1): 8008

[12] Robinson D A, Thomas A, Reinsch S, Lebron I, Feeney C J, Maskell L C, Wood C M, Seaton F M, Emmett B A, Cosby B J . Analytical modelling of soil porosity and bulk density across the soil organic matter and land-use continuum. Scientific Reports, 2022, 12(1): 7085

[13] Makovníková J, Širáň M, Houšková B, Pálka B, Jones A . Comparison of different models for predicting soil bulk density. Case study–Slovakian agricultural soils. International Agrophysics, 2017, 31(4): 491–498

[14] dos Santos W P, Vaz C M P, Martin-Neto L, Anselmi A, Tomasella J, de Souza Costa F, Albuquerque J A, de Jong van Lier Q, Galbieri R, Perina F J . Predicting bulk density in Brazilian soils for carbon stocks calculation: a comparative study of multiple linear regression and Random Forest models using continuous and categorical variables. Discover Soil, 2025, 2(1): 7

[15] Salamanca Jiménez A, Sadeghian Khalajabadi S. Bulk density and its relationship with other properties in soils of the Colombian coffee-growing region. Cenicafé Journal, 2005, 56(4): 381–397 (in Spanish)

[16] Barrezueta-Unda S, Luna-Romero E, Barrera-León J. Carbon storage in several soils planted with cocoa in El Oro Province, Ecuador. Agroecosystems Scientific Journal, 2018, 6(1): 154–161 (in Spanish)

[17] Benites V M, Machado P L O A, Fidalgo E C C, Coelho M R, Madari B E . Pedotransfer functions for estimating soil bulk density from existing soil survey reports in Brazil. Geoderma, 2007, 139(1−2): 90–97

[18] Bernoux M, Cerri C, Arrouays D, Jolivet C, Volkoff B . Bulk densities of Brazilian Amazon soils related to other soil properties. Soil Science Society of America Journal, 1998, 62(3): 743–749

[19] Haddad D B, de Assis L S, Tarrataca L, da S. Gomes A, Ceddia M B, Oliveira R F, de P. Junior J R, Brandão D N. Brazilian soil bulk density prediction based on a committee of neural regressors. In: 2018 International Joint Conference on Neural Networks (IJCNN). Rio de Janeiro, Brazil: IEEE, 2018, 1–8

[20] Han J Q, Mao K B, Xu T R, Guo J P, Zuo Z Y, Gao C Y . A soil moisture estimation framework based on the CART algorithm and its application in China. Journal of Hydrology, 2018, 563: 65–75

[21] Carbajal C, Moya-Ambrosio F, Barja A, Ottos-Diaz E, Aguilar-Tito C, Advíncula-Zeballos O, Cruz-Luis J, Solórzano-Acosta R . Soil quality variation associated with land cover in the Peruvian jungle of the Junín region. Soil Security, 2025, 19: 100188

[22] Lin W R, Tadai O, Takahashi M, Sato D, Hirose T, Tanikawa W, Hamada Y, Hatakeda K . An experimental study on measurement methods of bulk density and porosity of rock samples. Journal of Geoscience and Environment Protection, 2015, 3(5): 72–79

[23] Correa K, Avalos G, Endara Huanca S M, Sosa J, Lavado-Casimiro W, Quevedo Caiña K, Tello C, Ortega M. Report on dry and wet conditions in Peru during the 2020–2021 hydrological year. Lima: National Service of Meteorology and Hydrology of Peru (SENAMHI), 2021. Available at SENAMHI Institutional Repository on May 14, 2026 (in Spanish)

[24] Pizarro S, Pricope N G, Figueroa D, Carbajal C, Quispe M, Vera J, Alejandro L, Achallma L, Gonzalez I, Salazar W, Loayza H, Cruz J, Arbizu C I . Implementing cloud computing for the digital mapping of agricultural soil properties from high resolution UAV multispectral imagery. Remote Sensing, 2023, 15(12): 3203

[25] Ccopi D, Ortega K, Castañeda I, Rios C, Enriquez L, Patricio S, Ore Z, Casanova D, Agurto A, Zuñiga N, Urquizo J . Using UAV images and phenotypic traits to predict potato morphology and yield in Peru. Agriculture, 2024, 14(11): 1876

[26] Blake G R, Hartge K H. Bulk density. In: Klute A, ed. Methods of Soil Analysis: Part 1 Physical and Mineralogical Methods. 2nd ed. Madison: American Society of Agronomy and Soil Science Society of America, 1986, 363–375

[27] Sarker A, Kashem M A, Osman K T, Hossain I, Ahmed F . Evaluation of available phosphorus by soil test methods in an acidic soil incubated with different levels of lime and phosphorus. Open Journal of Soil Science, 2014, 4(3): 103–108

[28] International Organization for Standardization (ISO). ISO 11265: 1994 Soil Quality: Determination of the Specific Electrical Conductivity. Geneva: ISO, 1994

[29] Secretaría de Medio Ambiente y Recursos Naturales (SEMARNAT). Mexican Official Standard NOM-021-RECNAT-2000 Which Establishes Specifications for Soil Fertility, Salinity and Classification: Studies, Sampling and Analysis. Mexico City: SEMARNAT, 2002. Available at Orden Jurídico Nacional website on May 14, 2026 (in Spanish)

[30] United States Environmental Protection Agency (U.S. EPA). SW-846 Test Method 9045D: Soil and waste pH. Washington: U.S. EPA, 2004. Available at U.S. EPA website on May 14, 2026

[31] The R Core Team. R: a language and environment for statistical computing. Vienna: R Foundation for Statistical Computing, 2024. Available at R Project website on May 14, 2026

[32] Wei T, Simko V. corrplot: Visualization of a Correlation Matrix. R Package Version 0.95. 2024. Available at CRAN website on May 14, 2026

[33] Lê S, Josse J, Husson F . FactoMineR: An R package for multivariate analysis. Journal of Statistical Software, 2008, 25(1): 1–18

[34] Kassambara A, Mundt F. factoextra: Extract and Visualize the Results of Multivariate Data Analyses. R package version 2.0.0. 2026. Available at CRAN website on May 14, 2026

[35] Weisberg S. Applied Linear Regression. 4th ed. Hoboken: John Wiley & Sons, 2014

[36] Abdi H . Partial least squares regression and projection on latent structure regression (PLS Regression). WIREs Computational Statistics, 2010, 2(1): 97–106

[37] Wold S, Sjöström M, Eriksson L . PLS-regression: a basic tool of chemometrics. Chemometrics and Intelligent Laboratory Systems, 2001, 58(2): 109–130

[38] Bondi G, Creamer R, Ferrari A, Fenton O, Wall D . Using machine learning to predict soil bulk density on the basis of visual parameters: tools for in-field and post-field evaluation. Geoderma, 2018, 318: 137–147

[39] Parra J S, de Souza Z M, de Medeiros Oliveira S R, Farhate C V V, Marques J, Siqueira D . Phosphorus adsorption prediction through decision tree algorithm under different topographic conditions in sugarcane fields. Catena, 2022, 213: 106114

[40] Souza V F C, Cicalese F, Laber E S, Molinaro M . Decision trees with short explainable rules. Theoretical Computer Science, 2025, 1047: 115344

[41] Wu H B, Guo Z T, Peng C H . Land use induced changes of organic carbon storage in soils of China. Global Change Biology, 2003, 9(3): 305–315

[42] Ruehlmann J, Körschens M . Calculating the effect of soil organic matter concentration on soil bulk density. Soil Science Society of America Journal, 2009, 73(3): 876–885

[43] Nanko K, Ugawa S, Hashimoto S, Imaya A, Kobayashi M, Sakai H, Ishizuka S, Miura S, Tanaka N, Takahashi M, Kaneko S . A pedotransfer function for estimating bulk density of forest soil in Japan affected by volcanic ash. Geoderma, 2014, 213: 36–45

[44] Sevastas S, Gasparatos D, Botsis D, Siarkos I, Diamantaras K I, Bilas G . Predicting bulk density using pedotransfer functions for soils in the Upper Anthemountas basin, Greece. Geoderma Regional, 2018, 14: e00169

[45] Adams W A . The effect of organic matter on the bulk and true densities of some uncultivated podzolic soils. Journal of Soil Science, 1973, 24(1): 10–17

[46] Keller T, Håkansson I . Estimation of reference bulk density from soil particle size distribution and soil organic matter content. Geoderma, 2010, 154(3−4): 398–406

[47] Martín M Á, Reyes M, Taguas F J . Estimating soil bulk density with information metrics of soil texture. Geoderma, 2017, 287: 66–70

[48] Yi X S, Li G S, Yin Y Y . Pedotransfer functions for estimating soil bulk density: a case study in the three-river headwater region of Qinghai Province, China. Pedosphere, 2016, 26(3): 362–373

[49] Tranter G, Minasny B, McBratney A B, Murphy B, McKenzie N J, Grundy M, Brough D . Building and testing conceptual and empirical models for predicting soil bulk density. Soil Use and Management, 2007, 23(4): 437–443

[50] Federer C A . Nitrogen mineralization and nitrification: depth variation in four new England forest soils. Soil Science Society of America Journal, 1983, 47(5): 1008–1014

[51] Heuscher S A, Brandt C C, Jardine P M . Using soil physical and chemical properties to estimate bulk density. Soil Science Society of America Journal, 2005, 69(1): 51–56

[52] Chen Y, Liao Q B, Chen T Y, Zhang Y C, Yuan W E, Xu J Q, Zhang X Y . Freeze-drying formulations increased the adenovirus and poxvirus vaccine storage times and antigen stabilities. Virologica Sinica, 2021, 36(3): 365–372

[53] Moreira C S, Brunet D, Verneyre L, Sá S M O, Galdos M V, Cerri C C, Bernoux M . Near infrared spectroscopy for soil bulk density assessment. European Journal of Soil Science, 2009, 60(5): 785–791

[54] He W J, Xiao Z W, Lu Q K, Wei L F, Liu X . Digital mapping of soil particle size fractions in the Loess Plateau, China, using environmental variables and multivariate random forest. Remote Sensing, 2024, 16(5): 785

[55] Ali A, Ismael I O, Mustafa H T, Krwanji D, Esmail A O . Enhancing soil texture and bulk density mapping using soil grids and machine learning: a comparative analysis with observed data. Asian Soil Research Journal, 2024, 8(4): 61–78

[56] Bouslihim Y, Rochdi A, El Amrani Paaza N . Machine learning approaches for the prediction of soil aggregate stability. Heliyon, 2021, 7(3): e06480

[57] Rawls W J, Nemes A, Pachepsky Y. . Effect of soil organic carbon on soil hydraulic properties.. Developments in Soil Science, 2004, 30: 95–114

[58] Alexander E B . Bulk densities of California soils in relation to other soil properties. Soil Science Society of America Journal, 1980, 44(4): 689–692

[59] Curtis R O, Post B W . Estimating bulk density from organic-matter content in some Vermont forest soils. Soil Science Society of America Journal, 1964, 28(2): 285–286

[60] Calhoun F G, Smeck N E, Slater B L, Bigham J M, Hall G F . Predicting bulk density of Ohio soils from morphology, genetic principles, and laboratory characterization data. Soil Science Society of America Journal, 2001, 65(3): 811–819

[61] Akpa S I C, Ugbaje S U, Bishop T F A, Odeh I O A . Enhancing pedotransfer functions with environmental data for estimating bulk density and effective cation exchange capacity in a data-sparse situation. Soil Use and Management, 2016, 32(4): 644–658

[62] Wilson O, Schoeman D, Bradley A, Clemente C . Practical guidelines for validation of supervised machine learning models in accelerometer-based animal behaviour classification. Journal of Animal Ecology, 2025, 94(7): 1322–1334

[63] Teodorescu V, Obreja Brașoveanu L . Assessing the validity of k-fold cross-validation for model selection: evidence from bankruptcy prediction using random forest and XGBoost. Computation, 2025, 13(5): 127

[64] Beutler S J, Pereira M G, de Souza Tassinari W, de Menezes M D, Valladares G S, dos Anjos L H C . Bulk density prediction for histosols and soil horizons with high organic matter content. Revista Brasileira de Ciência Do Solo, 2017, 41: e0160158

[65] Périé C, Ouimet R . Organic carbon, organic matter and bulk density relationships in boreal forest soils. Canadian Journal of Soil Science, 2008, 88(3): 315–325

[66] Botula Y D, Nemes A, Van Ranst E, Mafuka P, De Pue J, Cornelis W M . Hierarchical pedotransfer functions to predict bulk density of highly weathered soils in central Africa. Soil Science Society of America Journal, 2015, 79(2): 476–486

[67] Cienciala E, Exnerová Z, Macků J, Henžlík V . Forest topsoil organic carbon content in Southwest Bohemiaregion. Journal of Forest Science, 2006, 52(9): 387–398

[68] Szabó B, Szatmári G, Takács K, Laborczi A, Makó A, Rajkai K, Pásztor L . Mapping soil hydraulic properties using random-forest-based pedotransfer functions and geostatistics. Hydrology and Earth System Sciences, 2019, 23(6): 2615–2635

[69] Meyer H, Reudenbach C, Hengl T, Katurji M, Nauss T . Improving performance of spatio-temporal machine learning models using forward feature selection and target-oriented validation. Environmental Modelling & Software, 2018, 101: 1–9

[70] Schaetzl R J, Luehmann M D, Rothstein D . Pulses of podzolization: the relative importance of spring snowmelt, summer storms, and fall rains on spodosol development. Soil Science Society of America Journal, 2015, 79(1): 117–131

RIGHTS & PERMISSIONS

The Author(s) 2027. Published by Higher Education Press. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0)