Spatial Variability of Soil Acidity and Lime Requirements for Potato Cultivation in the Huánuco Highlands

1. Introduction

Potato has historically played a vital role in global food security, with an annual worldwide production of approximately 3,740.7 million tonnes, making it the sixth most important crop after sugarcane, maize, wheat, rice, and oil palm [1]. Peru is the leading potato producer in the Andean region, with a total output exceeding 6 million tonnes [1]. This crop supports the livelihoods of over 700,000 families, most of whom are engaged in small-scale farming, contributing 10.1% to the country’s gross value of agricultural production (GVP) [2]. The Huánuco region stands out regarding cultivated area and production [2] and is one of the key centers of potato diversity [3,4]. However, Huánuco is among Peru’s poorest regions, with some of the highest extreme poverty levels [2]. It ranks as one of the most food-insecure areas, with a vulnerability index of 0.54, where values closer to 1 signify greater vulnerability [5]. This issue is further compounded by soil conditions, such as acidity, which significantly limit crop production in the region.

Soil acidity is one of the major constraints to crop production globally [6]. In South America, acidic soils account for 52% of the total agricultural area [7]. In Peru, pH values in the central Andes range between 4.2 and 4.6, classifying these soils as strongly acidic [8]. Acidic soils are typically characterized by severe deficiencies in essential nutrients like phosphorus, calcium, magnesium, and molybdenum while exhibiting toxic levels of aluminum, iron, and manganese [9]. These degraded, low-fertility soils exacerbate the vulnerability of native potato agroecosystems in the high Andean regions of Peru.

Several studies highlight that the integrated use of inorganic and organic fertilizers and liming amendments significantly enhances the physicochemical properties of acidic soils [10]. Specifically, liming raises soil pH to optimal levels, increases the availability of exchangeable calcium (Ca²⁺) and magnesium (Mg²⁺), and reduces phosphorus immobilization [11]. This practice also improves the soil’s physical and biological quality, increasing crop productivity [12].

However, soil properties and nutrient concentrations often fluctuate across crop fields due to the varying effects of soil-forming factors [13]. As a result, fertilizer applications can create areas of nutrient excess or deficiency [14]. In regions with excessive fertilizer levels, leaching can lead to water contamination, while insufficient fertilizer dosages can limit crop yields [15]. To address this, soil fertility maps, created using geostatistical methods, allow for precise scheduling of variable fertilizer and lime applications by predicting the spatial variation in soil properties [16].

The present study aims to predict the spatial variation of lime requirements for potato cultivation across five locations in the Pachitea province, Huánuco, Peru. The geostatistical interpolation method, Ordinary Kriging (OK), was employed to accurately assess the spatial variability of soil properties and generate characterization maps. Additionally, the study seeks to identify the physicochemical characteristics most strongly associated with soil acidity. Pearson correlation analysis and principal component analysis (PCA) were conducted on 23 soil physicochemical properties to achieve this.

2. Materials and Methods

2.1. Study Area

The study area is situated in the Peruvian highlands, specifically in the villages of Huarichaca, Chagragoto, Rumichaca, Plaza Punta, and Buenos Aires, located in the districts of Umari, Chaglla, Molino, and Panao, within the province of Pachitea, Huánuco region (Figure 1). The area experiences an average annual precipitation of 865.31 mm and a mean annual relative humidity of 78.46%. Minimum temperatures range from 9.36 °C in June to 5.4 °C in August, while maximum temperatures vary between 25.02 °C in September and 24.44 °C in November. These historical averages were calculated using data from the Chaglla meteorological station (9°50′49.9”S; 75°54′54.0”W; 2800 m.a.s.l.) managed by the National Service of Meteorology and Hydrology of Peru (SENAMHI) as shown in Figure 2.

2.2. Soil Sampling and Analysis

The R software [17] was used to determine the sampling points through the sample function from the sp package [18]. The R script generated the sampling coordinates by integrating a raster layer with high-resolution satellite imagery (e.g., Sentinel-2) as the base map, and stratified sampling was selected as the method. Ninety soil samples were collected at a depth of 30 cm from fields managed by native potato growers.

The soil samples were analyzed at INIA’s Soil, Water, and Foliar Laboratories network. The variables assessed were part of a comprehensive soil characterization using established reference methodologies. Soil texture (sand, silt, and clay percentages) was determined using the Bouyoucos hydrometer method [19]. pH was measured following EPA guidelines [20], while electrical conductivity (EC) was determined using the saturation extract method [21]. Organic matter content (OM) was analyzed using NOM-021-RECNAT-2000 [19] and total nitrogen (N) using ISO standards [22]. Available phosphorus for both neutral and acidic soils was measured using the Bray and Kurtz method [19], and available potassium (K), according to Bazán Tapia [23]. Exchangeable cations (H⁺, Al⁺³, Ca⁺², Mg⁺², K⁺, and Na⁺) were also quantified [19]. Bulk density (BD) values were obtained from the global digital soil mapping system, SoilGrids [24], with a spatial resolution of 250 meters, focusing on the 15-30 cm depth. The data were retrieved in TIFF format and processed using QGIS to extract the relevant values.

2.3. Estimation of the Liming Requirement for Soils

According to Havlin [25], calculations were made to estimate liming requirements.

$${\mathit{C}\mathit{E}\mathit{C}}_{\mathit{m}\mathit{E}\mathit{q}{\mathbf{100}\mathit{g}}^{-\mathbf{1}}}={({\mathit{C}\mathit{a}}^{+\mathbf{2}}+{\mathit{M}\mathit{g}}^{+\mathbf{2}}+{\mathit{K}}^{+}+{\mathit{N}\mathit{a}}^{+}+{\mathit{H}}^{+}+{\mathit{A}\mathit{l}}^{+\mathbf{3}})}_{\mathit{m}\mathit{E}\mathit{q}{\mathbf{100}\mathit{g}}^{-\mathbf{1}}}$$

(1)

$${\mathit{T}\mathit{S}\mathit{W}}_{\mathit{t}{\mathit{h}\mathit{a}}^{-\mathbf{1}}}={\mathit{B}\mathit{D}}_{\mathit{t}{\mathit{m}}^{-\mathbf{3}}}\times {\mathbf{0.3}}_{\mathit{m}}\times {\mathbf{10000}}_{{\mathit{m}}^{\mathbf{2}}}$$

(2)

$${\mathit{E}\mathit{A}\mathit{P}}_{\mathrm{\%}}={\displaystyle \frac{{({\mathit{H}}^{+}+{\mathit{A}\mathit{l}}^{+\mathbf{3}})}_{\mathit{m}\mathit{E}\mathit{q}{\mathbf{100}\mathit{g}}^{-\mathbf{1}}}}{{\mathit{C}\mathit{E}\mathit{C}}_{\mathit{m}\mathit{E}\mathit{q}{\mathbf{100}\mathit{g}}^{-\mathbf{1}}}}}$$

(3)

$${\mathit{L}\mathit{R}}_{\mathit{k}\mathit{g}{\mathit{h}\mathit{a}}^{-\mathbf{1}}}={\displaystyle \frac{({\mathit{E}\mathit{A}\mathit{P}}_{\mathrm{\%}}-{\mathit{P}\mathit{A}\mathit{S}}_{\mathrm{\%}})\times {\mathit{C}\mathit{E}\mathit{C}}_{\mathit{m}\mathit{E}\mathit{q}{\mathbf{100}\mathit{g}}^{-\mathbf{1}}}\times {\mathit{T}\mathit{S}\mathit{W}}_{\mathit{t}{\mathit{h}\mathit{a}}^{-\mathbf{1}}}\times \mathbf{2.8}}{\mathbf{550}}}$$

(4)

$${\mathit{D}\mathit{R}}_{\mathit{k}\mathit{g}{\mathit{h}\mathit{a}}^{-\mathbf{1}}}={\displaystyle \frac{({\mathit{E}\mathit{A}\mathit{P}}_{\mathrm{\%}}-{\mathit{P}\mathit{A}\mathit{S}}_{\mathrm{\%}})\times {\mathit{C}\mathit{E}\mathit{C}}_{\mathit{m}\mathit{E}\mathit{q}{\mathbf{100}\mathit{g}}^{-\mathbf{1}}}\times {\mathit{T}\mathit{S}\mathit{W}}_{\mathit{t}{\mathit{h}\mathit{a}}^{-\mathbf{1}}}\times \mathbf{2.8}}{\mathbf{300}}}$$

(5)

CEC (mEq∙100 g⁻¹) refers to the effective cation exchange capacity, determined by summing the exchangeable bases in acid soils. TSW (t∙ha⁻¹) indicates the topsoil weight at a depth of 30 cm that is suitable for potato cultivation. EAP (%) represents the exchangeable acidity percentage, calculated as the ratio of the soil’s exchangeable acidity to the effective CEC. LR (kg∙ha⁻¹) denotes the agricultural lime requirement, determined based on the mEq∙100 g⁻¹ of calcium necessary to neutralize the EAP up to the permissible acidity saturation (PAS) for potato crops. A threshold of 25% PAS is considered optimal for potato cultivation, and agricultural lime with a calcium oxide (CaO) content of 55% is used. For the agricultural dolomite requirement (DR), a CaO content of 30% is considered. Application of DR is recommended for plots where the percentage of exchangeable magnesium in the soil is less than 10%.

2.4. Geostatistical Interpolation

Understanding spatial variability through the appropriate geostatistical semivariogram model is essential for mapping and delineating the spatial variability of soil fertility, thereby optimizing fertilization programs [26]. Utilizing georeferencing, the soil analysis results from the province of Pachitea—on bulk density, percentages of sand, silt, and clay, pH, electrical conductivity (EC), organic matter (OM), available phosphorus (P) and potassium (K), effective cation exchange capacity (CEC), concentrations of exchangeable cations (H⁺, Al³⁺, Ca²⁺, Mg²⁺, K⁺, and Na⁺), and liming requirements—were employed for geostatistical analysis. Interpolations were performed using the Ordinary Kriging (OK) method, with the SAGA [27] and QGIS [28] software utilized throughout the process.

2.5. Model Validation

The nugget (C0), sill (C0 + C), range (R), and sill-nugget ratio (PSV) (6) are key parameters in semivariogram analysis that describe the spatial autocorrelation of the data. These parameters are derived from fitting the semivariogram model for geostatistical interpolation [29].

$$\mathit{P}\mathit{S}\mathit{V}={\displaystyle \frac{\mathit{S}\mathit{i}\mathit{l}\mathit{l}-\mathit{N}\mathit{u}\mathit{g}\mathit{g}\mathit{e}\mathit{t}}{\mathit{S}\mathit{i}\mathit{l}\mathit{l}}}={\displaystyle \frac{\mathit{C}}{\mathit{C}\mathbf{0}+\mathit{C}}}$$

(6)

The semivariogram model is selected based on criteria such as the root mean square error (RMSE) and the coefficient of determination (R²) [29]. Ideally, the chosen model should exhibit RMSE values close to zero (7) and R² values close to one, as these indicate the accuracy and quality of the model, respectively [30]. The kriging interpolations for each soil property were cross-validated using the leave-one-out method [31].

$$\mathit{R}\mathit{M}\mathit{S}\mathit{E}=\sqrt{{\displaystyle \frac{\mathbf{1}}{\mathit{n}}}{\sum }_{\mathit{I}=\mathbf{1}}^{\mathit{n}}[{\mathit{Z}}_{\mathbf{1}}({\mathit{x}}_{\mathit{i}})-{\mathit{Z}}_{\mathbf{2}}({\mathit{x}}_{\mathit{i}}){]}^{\mathbf{2}}}$$

(7)

Where n represents the number of samples, ${\mathit{Z}}_{\mathbf{1}}({\mathit{x}}_{\mathit{i}}$) and ${\mathit{Z}}_{\mathbf{2}}\left({\mathit{x}}_{\mathit{i}}\right)$ represent the predicted and observed values at site i, respectively.

3. Results

3.1. Descriptive Statistics for Soil Acidity

Table 1 presents the results of the soil samples’ physical-chemical analysis. The pH values ranged from 3.9 to 5.8, with a mean of 4.69 ± 0.39. The exchangeable acidity percentage (H⁺ + Al⁺³) varied between 0% and 72.75%, yielding a mean of 30.94 ± 24.87%. The variance for pH was 0.16, while the exchangeable acidity percentage (EAP) exhibited a variance of 618.73, indicating a high dispersion of potential acidity data within narrow pH ranges. All locations displayed mean pH values below the critical threshold (5.5) necessary for optimal potato crop growth. Plaza Punta and Buenos Aires were identified as the most degraded sites, characterized by pH values below the overall mean (Figure 4).

Regarding the potential acidity of the soils, the localities of Buenos Aires, Plaza Punta, Chagragoto, and Rumichaca exhibit EAP that exceeds the overall mean of 30.94 ± 24.87% (Figure 5). However, only the localities of Buenos Aires and Plaza Punta have averages that surpass both the general mean and the critical threshold (25%) for potato cultivation. Among these, the acid soils in Plaza Punta were the most degraded, showing a frequency of 100% of the data above the critical value. In contrast, the acid soils of Huarichaca were the least degraded, with 85% of the data falling below the critical threshold (Figure 5).

3.2. Principal Component Analysis of Soil Properties

Principal component analysis (PCA) was conducted on 23 soil physicochemical characteristics. As shown in Figure 5, the first two components account for 59% of the variance across the five studied locations. The biplot indicates that the variables contributing the most to the variance explained by Component I are exchangeable H⁺, exchangeable acidity percentage (EAP), and exchangeable calcium percentage (ECP). The soils in Plaza Punta and Buenos Aires show a stronger association with exchangeable acidity (H⁺ + Al⁺³), sand percentage, and organic matter (OM), which reflect a greater need for liming. Plaza Punta soils, in particular, exhibit high EAP and an inverse relationship with the exchangeable magnesium percentage (EMP) and cation exchange capacity (CEC), indicating a higher requirement for dolomite. Conversely, the variance related to acidity is not associated with the soils from Chagragoto and Huarichaca, which have higher pH, ECP, and bulk density (BD) values, indicating no liming requirement. However, the Rumichaca soil variance has no clear relationship with the studied physico-chemical variables.

3.3. Pearson’s Correlation Analysis of Soil Properties and Their Relationship with Soil Acidity

Pearson correlation analysis was conducted for the 23 evaluated variables (Figure 6). The most strongly correlated variable with soil pH is the exchangeable calcium percentage (ECP), with a positive correlation of 0.77 (Figure 7d). Additionally, there is a significant negative correlation between ECP and both the exchangeable acidity percentage (EAP) (Figure 7a) and exchangeable H⁺ (Figure 7b), with correlation values of -0.96 and -0.84, respectively. The exchangeable magnesium percentage (EMP) also shows a strong negative correlation with EAP, reaching a value of -0.74 (Figure 7c). Increasing ECP and EMP can reduce exchangeable H⁺ concentrations and raise soil pH. Furthermore, a noteworthy positive correlation exists between soil acidity with the sand percentage and organic matter (OM), with values of 0.63 and 0.71, respectively. Several variables also exhibit a high positive correlation with pH, including exchangeable Mg⁺², Ca⁺², K⁺, clay content, available K, and bulk density (BD), with respective values of 0.52, 0.65, 0.66, 0.59, 0.66, and 0.50. These relationships are critical, as soils with higher sand content and high OM levels require larger lime doses. In contrast, treatments that integrate different exchangeable bases and improve soil bulk density will enhance the effectiveness of liming on acidic Andean soils.

3.4. Spatial Variation of Soil Acidity in the Localities of Plaza Punta and Buenos Aires

Table 2 presents the parameters of the semivariogram models obtained using the Ordinary Kriging (OK) interpolation method and the cross-validation results based on RMSE and R² values for model evaluation. The Linear and Exponential semivariogram models provided the best parameterization for predicting unknown values of the soil parameters under study. In the locations of Plaza Punta and Buenos Aires, the soil physicochemical variables with the highest spatial predictability were Al⁺³ (mEq∙100 g⁻¹), Ca⁺² (mEq∙100 g⁻¹), Mg⁺² (mEq∙100 g⁻¹), exchangeable calcium percentage (ECP), and exchangeable acidity percentage (EAP), with R² values of 0.60, 0.76, 0.87, 0.57, and 0.53, respectively. These R² values indicate a strong correlation between observed and predicted values in unsampled areas. Notably, Ca⁺² (mEq∙100 g⁻¹) and Mg⁺² (mEq∙100 g⁻¹) emerged as significant variables for explaining the variance in the soils’ physicochemical properties in the PCA. Both variables show high negative correlations with EAP, at -0.71 and -0.66, respectively. Furthermore, they exhibit significant R² values in the prediction of spatial variability. Based on these findings, spatial variation maps were produced for EAP (Figure 8), Ca⁺² (Figure 9), and Mg⁺² (Figure 10), as these characteristics play a crucial role in explaining soil acidity and liming requirements.

3.5. Spatial Variation of Soil Acidity in the Localities of Huarichaca, Chagragoto and Rumichaca

Table 3 presents the parameters of the semivariogram models derived using the Ordinary Kriging (OK) interpolation method alongside the cross-validation results based on RMSE and R² values for evaluation. The semivariogram models with Linear and Gaussian parameterization proved to be the most suitable for predicting unknown values of the studied soil parameters. In the localities of Chagragoto, Huarichaca, and Rumichaca, the soil physicochemical variables with the highest spatial predictability were bulk density (g∙cm⁻³), exchangeable magnesium percentage (EMP), and exchangeable potassium percentage (EPP), with R² values of 0.79, 0.54, and 0.54, respectively. These values indicate a strong alignment between observed and predicted data in non-sampled areas, reflecting a high model fit. Notably, EMP (%) and EPP (%) emerged as significant variables in explaining the variance of the soil’s physicochemical characteristics in the principal component analysis (PCA). They showed strong negative correlations of -0.74 and -0.18 with exchangeable acidity percentage (EAP). Additionally, these variables exhibited substantial R² values for predicting spatial variability. Consequently, spatial distribution maps were generated for bulk density (Figure 11), magnesium (Figure 12), and potassium (Figure 13), which are crucial for understanding soil acidity and liming requirements (Figure 14) for these locations.

3.6. Liming Requirement

Figure 13 illustrates the spatial variation in liming requirements for the localities of Plaza Punta and Buenos Aires, where calcium (Ca²⁺) and magnesium (Mg²⁺) levels exhibited a significant correlation with observed and predicted values (R² values of 0.76 and 0.87, respectively) (Figure 14a). Both variables also showed strong positive correlations with soil pH, with correlation coefficients (R) of 0.64 and 0.51, respectively. Furthermore, linear regression analysis indicates that the observed values of Ca²⁺ closely align with the predicted liming requirements (Figure 14b).

The spatial variation map of soil acidity identifies five zones with differing requirements for lime and dolomite (Figure 15). Specifically, 34.07% of the areas necessitate liming with dolomite, while 18.68% require agricultural lime. In Plaza Punta, the liming requirement averages 5.1 ± 2.8 t·ha⁻¹, predominantly in regions with a higher demand for agricultural dolomite (Figure 15). In contrast, Buenos Aires has an average liming requirement of 3.4 ± 5.5 t·ha⁻¹, primarily in areas with moderate needs for agricultural lime. Rumichaca shows a liming requirement of 2.5 ± 4.4 t·ha⁻¹, more commonly found in regions with elevated agricultural lime needs (Figure 15). Chagragoto exhibits an average requirement of 1.6 ± 2.8 t·ha⁻¹, typically in areas with higher demands for agricultural dolomite. Finally, the locality of Huarichaca has liming values of 0.1 ± 0.3 t·ha⁻¹, which are more prevalent in areas with low requirements for agricultural lime.

4. Discussion

4.1. Soil Acidity in Potato Cultivation and Its Relationship with the Soil Physicochemical Properties

Potatoes are a crop with high nutritional demands, particularly for macronutrients. It requires 1.7 kg of nitrogen (N), 0.7 kg of phosphorus (P), and 3.15 kg of potassium (K) per ton of production [32]. However, when grown in acidic soils with pH levels below 5.5, potatoes experience severe nutrient deficiencies at the leaf level [33]. Furthermore, high exchangeable aluminum (Al³⁺) concentrations in acidic soils have been reported to be toxic to potato growth [34,35]. Elevated levels of potential soil acidity (Al³⁺ + H⁺), ranging from 20% to 25%, hinder root cell division and respiration, increase cell wall rigidness, and interfere with calcium and magnesium uptake and transport [36,37]. In this study, the soils from the five potato-producing localities exhibited average pH values of 4.69 ± 0.39 and an exchangeable acidity percentage of 30.94 ± 24.87%, creating limiting conditions for crop growth and development.

To address this issue, it is essential to prioritize liming in fertilization plans for potato crops to raise soil pH [38]. In this study, neutralizing the potential acidity method to a critical threshold of 25% for potato cultivation estimated an average liming requirement of 2.78 ± 3.77 t∙ha⁻¹. Furthermore, principal component analysis revealed that variables related to soil acidity, such as EAP and ECP, accounted for the most variance among the 23 physicochemical characteristics evaluated. These findings highlight significant spatial variability in the severity of acid soil degradation across Pachitea province.

The localities of Plaza Punta and Buenos Aires exhibit the lowest pH values (4.41 ± 0.2 and 4.48 ± 0.18, respectively), which can be attributed to high rainfall intensity, with mean annual precipitation of 1035.33 mm and 1159.37 mm, respectively. Elevated rainfall enhances the chemical weathering of soil clay minerals and introduces significant amounts of acid pollutants, such as SO₂ and nitrogenous compounds (NO, NO₂, and N₂O) [39]. However, the soils in Plaza Punta have lower cation exchange capacity (CEC) than those in Buenos Aires (4.8 ± 1.17 and 7.99 ± 1.98 mEq∙100 g⁻¹, respectively). This difference in CEC explains the variation in the severity of soil potential acidity between both localities. Plaza Punta has the highest exchangeable acidity percentage (EAP) (56.91 ± 13.49%) among the five localities and the highest liming requirement (5.1 ± 2.8 t ha⁻¹).

Differences in the severity of soil acidity and its negative impact on fertility potential between localities are likely due to mineralogical variations [40]. Soils rich in gibbsite (aluminum hydroxides) tend to have low CEC values, ranging from 1 to 4 mEq∙100 g⁻¹, while soils rich in kaolinite (1:1 phyllosilicates) typically have higher CEC values, between 4 and 15 mEq∙100 g⁻¹ [41]. This suggests that Plaza Punta soils may be high in gibbsite, whereas Huarichaca, Buenos Aires, and Chagragoto soils, with CEC values of 9.6, 8, and 6.2 mEq∙100 g⁻¹, respectively, are more kaolinite-rich. Gibbsite is a secondary mineral formed in soils experiencing intense 2:1 and 1:1 chemical weathering of aluminosilicates [42], which likely explains these mineralogical differences.

The results obtained in the Pearson correlation analysis suggest that the source of soil acidity is strongly linked to both clay and organic matter content. The clay percentage shows a strong negative correlation with exchangeable acidity percentage (EAP) (-0.71), while organic matter content is highly positively correlated with exchangeable soil acidity (0.71). This indicates that the dominant pedogenic processes responsible for soil acidity in these areas are clay illuviation and organic matter accumulation. One hypothesis posits that acidic kaolinitic ultisols, where these processes predominate, contain kaolinite minerals interbedded with minerals formed from hydroxylated vermiculite weathering [43]. Intercalation with hydroxylated minerals in such acidic soils significantly reduces their CEC, sometimes leading to misclassification of their mineralogy as kaolinitic [44]. This characteristic suggests that there may be greater potential to increase the CEC of these soils, as observed in the Chagragoto locality, where maximum CEC values exceeding 20 mEq∙100 g⁻¹ were found. This points to higher potential fertility in these soils due to the effect of the pH increase.

The positive effect of liming on potato crops is likely attributed to the soil’s increased CEC, resulting from a higher density of pH-dependent variable charges [39], which effectively enhances the soil’s nutrient storage capacity. Studies have confirmed that it is possible to increase potato yields on strongly acidic ultisols by raising the soil pH with dolomitic lime or by applying calcium and magnesium phosphate fertilizers without altering the pH [45]. This highlights the critical role of calcium and magnesium fertilization in managing acidic soils. In our study, Ca²⁺ (mEq∙100 g⁻¹) and Mg²⁺ (mEq∙100 g⁻¹) were the key variables in reducing soil acidity, showing a strong negative correlation with exchangeable soil acidity (-0.71 and -0.66, respectively) concerning the exchangeable acidity percentage (EAP). These findings suggest that the mass effect on the cation exchange of H⁺ and Al³⁺ ions for Ca²⁺ and Mg²⁺ is the most critical factor in acidic soil management [39].

In this context, recommended liming doses ranging from 2 to 14 t∙ha⁻¹ are expected to enhance the soil’s potential fertility. It is also likely that treatments aimed at increasing pH-dependent charges will improve the physicochemical properties of acidic high Andean soils. Several studies support the importance of these agricultural practices. Bordoloi et al. [46] highlight that applying 0.4 t∙ha⁻¹ of agrarian lime, supplemented with 1 t∙ha⁻¹ of vermicompost, increased soil pH from 4.54 to 5.12, reducing fertilizer usage by 50% and boosting potato yields by 36-54%.

4.2. Spatial Variability of Soil Acidity and Liming Requirement

Geostatistical analysis using semivariogram models and the Ordinary Kriging (OK) method enables precise identification of soil variables with the highest spatial predictability related to soil acidity [16]. Key indicators like the coefficient of determination (R²) and root mean square error (RMSE) are critical in evaluating the fit of these models between observed and predicted values in unsampled areas [47]. Models with R² values approaching 1 and low RMSE indicate a superior fit, making them effective tools for soil management in acidic soil regions and facilitating the optimization of practices such as liming and fertilization [48].

Several key variables have been identified in predicting soil acidity, with Ca²⁺ and Mg²⁺ being the most influential in determining the spatial variability of soil fertility, particularly in the localities of Buenos Aires and Plaza Punta. These exchangeable cations are strongly correlated with soil acidity, allowing for the generation of precise maps for agricultural management using geostatistical methods like Ordinary Kriging [49]. R² values ranging from 0.70 to 0.85, along with low RMSE in predictions for Ca²⁺ and Mg²⁺, are adequate for accurately forecasting soil fertility and acidity levels. Additionally, variables like bulk density and exchangeable potassium percentage (EPP) have proven to be critical in acidic soils found in grazing areas, with significant spatial model adjustments, particularly in the localities of Huarichaca, Chagragoto, and Rumichaca. This suggests that the relevance of predictive variables can vary depending on climatic conditions, land use, the extent of acid soil degradation, and other soil properties [50].

5. Conclusions

Using soil fertility maps enables the application of variable liming doses, leveraging geostatistical methods that accurately capture the spatial variation in soil acidity and potential fertility. This approach is relevant to food security and potato productivity among small farmers in the Pachitea province, Huánuco. The study identifies the key factors influencing soil acidity by conducting robust statistical analyses of soil chemical properties. Precision agriculture tools, such as these soil fertility maps, facilitate targeted liming interventions, enabling the calibration of optimal liming doses and the development of soil management strategies at a regional level.

The soil fertility analysis revealed that 47.25% of the total producers do not require liming, while 52.75% do need liming to enhance the chemical properties of their soil. Of the producers needing liming, 34.07% require dolomite, and 18.68% require agricultural lime. The average recommended application dose for lime and dolomite across all locations is 2.78 ± 3.77 t∙ha⁻¹.

It is essential first to restore the degraded acidic soils of the high Andean region, particularly those with pH values below 5.5 and exchangeable acidity percentages exceeding 25%. This can be achieved through liming treatments, complemented by strategies aimed at increasing the soil’s variable charges tailored to the specific needs of each locality. Such an approach will facilitate the development of optimal and sustainable fertilization plans for potato agroecosystems in the Peruvian Andean highlands. Furthermore, soil fertility maps generated from geostatistical data are invaluable for farmers, researchers, and policymakers, enabling improved soil management strategies, optimized liming applications, and enhanced potato crop productivity. The methodological approach applied in this study has effectively demonstrated the strong spatial dependence of the analyzed soils.