Spatial Prediction of Soil Total Phosphorus in a Karst Area: Comparing GWR and Residual-Centered Kriging

Abstract

Accurate soil total phosphorus (TP) prediction is essential to support sustainable agricultural practices and formulate ecological conservation protection policies, particularly in complex karst landscapes with high spatial variability and high phosphorus and cadmium content and interactions, complicating nutrient management. This study uses GIS and geostatistical methods to analyze the spatial distribution, influencing factors, and predictive modeling of soil TP in the karst region of northern Mashan County, Guangxi, China. Using 427 surface soil samples, we developed five predictive models: ordinary kriging (OK), regression kriging (RK) and geographically weighted regression kriging (GWRK) combined with environmental variables such as land uses, soil types, and topographic factors; residual mean-centered kriging (MM_OK), and residual median-centered kriging (MC_OK). Our results indicate that higher TP levels were observed in agricultural lands (paddy fields and dry land, at 766 and 913 mg·kg⁻¹, respectively) may due to fertilization, while forests and shrublands showed lower TP levels (383 and 686 mg·kg⁻¹, respectively), reflecting natural phosphorus cycling. The high-value areas of soil TP concentration are in the karst areas in the west and east of the study area, and the low-value area is in the Hongshui River valley in the north of Mashan. The spatial distribution of soil TP is affected by land use, soil type, and topography. The GWRK model exhibited superior accuracy (80.6%), with predicted concentration of TP closely aligning with observed TP values, effectively capturing fine spatial variations, and showing the lowest mean standardized error, average standard error, and mean absolute error. GWRK also achieved the highest R² (0.67), demonstrating robust predictive capability. MM_OK and MC_OK models performed well and showed smoother spatial transitions, while the OK model displayed the lowest predictive accuracy (62%). By utilizing spatially adaptive weighting, GWRK and its residual-centered kriging method improve soil TP’s prediction accuracy and smoothness in karst areas, providing a reference for targeted soil conservation and sustainable agricultural practices in spatially complex karst environments.

1. Introduction

Phosphorus is a critical indicator of soil quality, profoundly influencing soil fertility and plant growth [1,2]. As an essential nutrient, soil phosphorus supports photosynthesis and root development [3], which are vital for maintaining ecological stability and enhancing agricultural productivity. The phosphorus cycle is highly dynamic, and land use practices and soil types interact to shape phosphorus distribution, availability, and cycling efficiency within ecosystems [4,5,6]. In karst landscapes with shallow soils and severe erosion, the phosphorus cycle is likely to be disrupted, as shallow soils have low phosphorus retention, erosion can cause loss of topsoil, and the calcium in the landscape may bind phosphorus, reducing its availability [7,8,9,10]. Applying Cd-rich phosphate fertilizers to cultivated land will also increase the risk of safe use of soil in karst areas with high Cd geological backgrounds. Accurately mapping total phosphorus (TP) is crucial for understanding its distribution, enabling targeted fertilization, enabling soil conservation, and formulating ecological conservation protection policies [11].

As a widely applied regression technique, Ordinary Least Squares (OLS) is used to identify linear relationships between TP and auxiliary variables, such as soil properties, topography, and climate. Geographically Weighted Regression (GWR) extends OLS by allowing regression coefficients to vary spatially, thereby accounting for local variations in the relationships between TP and predictor variables. This approach has demonstrated strong potential for capturing spatial heterogeneity and reducing autocorrelation within the model, thereby enhancing prediction accuracy [12,13]. However, in complex karst terrains, where topographical and spatial variability are high, both OLS and GWR may generate inconsistent residual patterns, revealing limitations in these approaches.

Residual analysis is essential in regression modeling, as residuals reveal unexplained spatial variability and help assess model robustness [14]. In spatial prediction, residual kriging methods are frequently employed to enhance predictive accuracy. Co-kriging, which integrates additional environmental variables, has shown promise in improving prediction accuracy for soil properties, such as organic carbon and nitrogen, within karst landscapes [15]. However, conventional kriging assumes spatial autocorrelation, which contradicts the expectation that regression residuals should be randomly distributed. This inconsistency can decrease the predictive accuracy of kriging for residuals, particularly if residuals are not centered around zero.

To mitigate these limitations, centering residuals—by adjusting them around mean or median values—has emerged as a promising approach to improve model stability, reduce multicollinearity, and enhance spatial compatibility for kriging. Centering has been shown to mitigate distortions in spatial prediction and improve parameter estimation by addressing multicollinearity issues [16]. Despite advancements in TP spatial mapping, accurately capturing spatial variability in complex karst regions remains challenging, underscoring the need for refined methods tailored to these unique environments.

While considerable research has explored TP distribution, formation processes, and ecological impacts in karst areas [17,18,19], few studies have focused on regression kriging (RK) and geographically weighted regression kriging (GWRK) in these regions, mainly due to data acquisition challenges associated with fragmented karst terrain [20]. This study hypothesizes that centering residuals before kriging, using both mean and median values, can enhance TP spatial prediction accuracy in karst regions by improving spatial autocorrelation alignment and reducing regression bias. The specific objectives are: (1) to investigate the relationship between TP and environmental variables (e.g., land use, soil type) through regression modeling, evaluating how these factors influence TP spatial prediction accuracy; and (2) to assess the predictive accuracy of various models—including ordinary kriging (OK), RK and GWRK combined with environmental variables such as land uses, soil types, and topographic factors; residual mean-centered kriging (MM_OK) and residual median-centered kriging (MC_OK)—in mapping TP within a representative karst landscape using high-density sampling data (approximately 8 samples/km²) from a geochemical land survey.

This study focuses on a karst terrain in southwestern China, characterized by Carboniferous carbonate formations and diverse land covers including forest, shrubland, and farmland, representative of the broader karst regions in southwest China. By comparing spatial prediction approaches, this research aims to advance TP spatial mapping techniques, supporting sustainable soil management and conservation practices in ecologically sensitive karst landscapes. The study’s findings will provide valuable insights into model selection for soil nutrient prediction in karst areas, with implications for guiding sustainable fertilization practices for farmers, supporting soil quality monitoring by scientists, and informing policy decisions for regulators.

2. Materials and Methods

2.1. Overview of the Study Area

The study area is in northern Mashan County, Guangxi, China, within a typical karst landscape (Figure 1). Positioned at coordinates 108.09–108.2 °E and 23.73–23.85 °N, the region experiences a subtropical monsoon climate, with an average annual temperature of 21.3 °C and rainfall ranging from 1480 to 1667 mm [17]. The topography includes peak-cluster depressions and valleys typical of southern China’s karst regions, with formations of dolomitic limestone, Carboniferous limestone layers, and Permian to Triassic shales, silty shales, and mudstones in the north-central region [17,21]. This geological diversity represents a mix of carbonate and non-carbonate rocks commonly found in Southwest China.

The study area spans 60 km² and encompasses diverse land uses: paddy fields (10.55%), dry land (29.39%), shrublands (29.49%), and forests (12.86%). Dry land, which includes ferralsols, lixisols, brown cambisols, and fluvisols, is found in karst depressions and valley slopes, whereas paddy fields, which are made up of anthrosols, are found in karst valleys, according to the WRB soil classification system. Shrublands occupy karst hills with predominant brown cambisols, and forests are found on clastic rock hills with Lixisols and ferralsols [22].

2.2. Sample Collection and Testing

The data on the spatial variability of soil total phosphorus in this study are consistent with the data from the soil ecological stoichiometry study conducted by Ma et al. [17]. Surface soil samples were collected using a 500 m × 500 m grid from May to June 2017, with adjustments based on land use and soil type as informed by a 1:50,000 land-use map featuring contour lines and grids. Three to five soil columns of equal volume and depth were collected within a range of 20–50 meters and combined into a single sample, resulting in a density of 4–12 samples per km², depending on land use. Litter was removed before sampling soil columns from the 0–20 cm depth, with uniformity maintained by excluding roots, stones, and other debris. Samples were processed in an anti-pollution environment, which included grinding and screening, to prepare them for subsequent analysis.

The concentration of TP was measured via X-ray fluorescence spectrometry (Nippon Rigaku, ZSX Primus II, Tokyo, Japan), while available phosphorus (AvP) was determined using ammonium fluoride-hydrochloric acid or sodium bicarbonate solutions followed by ICP-OES ( iCAP™ 7600, Thermo Fisher Scientific, Waltham, MA, USA). Soil pH was measured using a glass electrode (Shanghai Lei Ci, PHS-3C) with a soil-to-water ratio of 1:2.5. Soil organic matter (SOM) was quantified through potassium dichromate-sulfuric acid oxidation, titrated with ammonium ferrous sulfate, following the “Land Quality Geochemical Evaluation Specification” (DZ/T 0295-2016).

2.3. Data Sources and Research Methods

2.3.1. Data Sources and Processing

The study utilized 30-meter resolution digital elevation model (DEM) data from the Geospatial Data Cloud (https://www.gscloud.cn). The slope, aspect, and elevation of each sampling location were extracted by ArcGIS 10.6 to reveal the contribution of these topographic factors in the soil TP regression analysis. These plots can be found in [22]. Land use data were derived from the 2014 land-use survey, while soil-type data were obtained from the Second National Soil Survey [17]. The 2014 land-use survey utilized high-resolution satellite imagery, supplemented by field validation to ensure accuracy, with an overall classification accuracy exceeding 85%. Soil-type data from the Second National Soil Survey were derived from extensive soil sampling and laboratory analyses conducted nationwide, offering a high degree of detail and reliability. These datasets were carefully evaluated for their applicability to the study area, ensuring consistency with local conditions. Qualitative variables, such as soil type and land use, were encoded as dummy variables for statistical analysis, with forest land serving as the reference category.

Point-biserial correlation in SPSS 20 assessed the relationships between land use, soil type, and TP. Following outlier removal using the mean ± 3 × standard deviation, 427 samples were selected for TP spatial prediction. ArcGIS 10.6 statistical modules divided samples into a training set (75%) and a validation set (25%), with descriptive and significance testing conducted in SPSS 20.

The phosphorus activation coefficient (PAC), calculated as AvP/TP, serves as an indicator of phosphorus availability, with high PAC values denoting high availability and low values indicating strong phosphorus fixation [23].

2.3.2. Spatial Prediction Models: OK, RK, GWRK, MM_OK and MC_OK

This study employed five spatial prediction models—OK, RK, GWRK, MM_OK, and MC_OK combined with environmental variables such as land use, soil type, and topographic factors to predict TP. A detailed description of these models and interpolation can be found in [22]. The RK model was used to decompose a predicted value and a residual component and perform ordinary kriging interpolation (OK) on the residual component. The final prediction result was generated by summing the above two components [24]. We used the mean or median of the GWR residuals for centering and then interpolated and added the predicted values to get the results of MM_OK and MC_OK.

Model accuracy was primarily verified using a training set for internal verification and a validation set for external verification [22]. The best results of internal verification include the following: the standard mean error (MSE) is closest to 0, the standard root mean square (RMSS) prediction error is closest to 1, the root mean square (RMS) prediction error is the smallest, and the average standard error (ASE) and the root mean square (RMS) prediction error are closest. For better external verification, results are as follows: the smaller the root mean square error (RMSE), the more stable the model, and the smaller the mean absolute error (MAE), the higher the model accuracy. Accuracy can evaluate the accuracy of the prediction. The value range is 0~1, and the closer to 1, the more accurate the prediction. A methodology flow chart is summarized in Figure 2.

3. Results

3.1. Geochemical Characteristics of TP in Soil

3.1.1. Descriptive Statistics and Correlation Analysis

The concentration of TP in soil within the study area ranges from 168 to 1791 mg·kg⁻¹, with an average of 796 mg·kg⁻¹, significantly exceeding the national average of 500 mg·kg⁻¹ for agricultural topsoil reported in China’s second national soil survey (

Table 1. Average Soil TP Concentrations and Basic Physical-Chemical Properties by Sample Type.

Parameters		Sample Size	pH	TP/mg·kg−1	AvP/mg·kg−1	PAC (%)	SOM/g·kg−1
Land uses	Paddy fields	83	6.16 b	766 b	17.50 b	2.25 b	20.7 a
	Dry land	247	5.61 c	913 a	28.86 a	3.20 a	11.1 c
	Forest	51	4.77 d	383 c	9.45 c	2.21 b	14.0 b
	Shrubland	46	6.94 a	686 b	2.32 d	0.40 c	15.1 b
Soil types	Anthrosols	98	6.14 a	852 b	22.52 ab	2.59 b	18.5 a
	Ferralsols	116	5.43 bc	697 c	21.24 b	2.92 b	12.5 b
	Lixisols	64	5.47 bc	666 c	20.67 b	2.69 b	12.1 b
	Cambisols	95	6.33 a	986 a	15.89 c	1.38 c	13.2 b
	Fluvisols	54	5.10 c	733 b	31.09 a	3.97 a	10.2 c
Total samples		427	5.76	796	21.50	2.60	13.7

Note: PAC = AvP/TP × 100%. Lowercase letters indicate statistically significant differences (p < 0.05) between land use or soil types.

). Data distribution analysis indicates skewness and kurtosis values of 0.74 and 3.14, respectively. The kurtosis value is slightly greater than the ideal value for a normal distribution (3). This suggests that the distribution may have slightly heavier tails or a sharper peak than a perfectly normal distribution. While the deviation is small, it supports the conclusion of non-normality in combination with other measures. Skewness, with a value of 0.74, greater than 0, indicates right skew (longer tail on the right). Following logarithmic transformation, the data distribution approximates normality (p > 0.05).

Mean TP concentrations vary significantly across different land uses, following the descending order of dry land > paddy fields > shrubland > forest (Table 1 and

Table 2. Descriptive Statistics of TP Concentration by Sample Type.

Parameters		Samples	Minimum	Maximum	Mean	SD	CV(%)
Land uses	Paddy fields	83	290	1550	766 b	261	34.13
	Dry land	247	376	1791	913 a	334	36.57
	Forest	51	168	848	383 c	138	35.98
	Shrubland	46	169	1326	686 b	281	41.04
Soil types	Anthrosols	98	304	1791	852 b	347	40.71
	Ferralsols	116	168	1704	697 c	363	52.15
	Lixisols	64	205	1269	666 c	220	32.96
	Cambisols	95	242	1753	986 a	368	37.33
	Fluvisols	54	230	1086	733 b	165	22.55
Total samples		427	168	1791	796	343	43.1

Note: Minimum, maximum, mean and SD are in mg·kg⁻¹; SD, standard deviation; CV, coefficient of variation; lowercase letters in the upper right corner of the mean value indicate significant differences (p < 0.05).

). This trend underscores the strong influence of land use on soil TP concentration. One-way ANOVA results indicate significant differences (p < 0.05) in mean TP concentrations across these land uses, suggesting that land use is a key factor affecting TP distribution, which is critical for spatial mapping. Among soil types, cambisols exhibit the highest mean TP levels.

Correlation analysis revealed significant relationships between TP concentration and both land uses and soil types, showing a positive correlation in dry land and a negative correlation in shrubland (

Table 3. Pearson correlation coefficient analysis of TP and sample types.

Land Uses		Soil Types
Paddy fields	−0.047	Cambisols	0.296 **
Dry land	0.400 **	Anthrosols	0.085
Forest	−0.041	Ferralsols	−0.160 **
Shrubland	−0.122 *	Fluvisols	0.458 **
		Lixisols	0.176 **

** and * Significantly correlated at the 0.01 and 0.05 level (bilateral), respectively.

). Both dry land and fluvisols show significant positive correlations with TP levels, suggesting these conditions are associated with higher TP in soils, possibly because these areas receive nutrient inputs like fertilizers or are subject to minimal phosphorus loss. This aligns with studies indicating that TP variability is governed by anthropogenic activities and land management practices, as supported by similar findings in agricultural soils in China. Shrubland and ferralsols show significant negative correlations with TP, indicating lower TP concentrations potentially due to phosphorus loss from slopes caused by soil erosion [9]. Cambisols also have a significant positive correlation with TP, likely due to high phosphorus retention capacity in calcareous soils. While the correlation between TP and paddy fields is weak and not significant, this suggests that paddy fields may experience phosphorus depletion due to regular waterlogging and leaching.

3.1.2. AvP and PAC Characteristics

Soil AvP levels range from 0.57 to 170.08 mg·kg⁻¹, with an average of 21.50 mg·kg⁻¹ (Table 1). The distribution of AvP is heterogeneous, as 54% of samples have concentrations between 10 and 40 mg·kg⁻¹, while 39.58% are below 10 mg·kg⁻¹. These values align with the average AvP concentrations reported for major rice-producing regions in China. AvP concentrations are highest in dry land, followed by paddy fields, forest, and shrubland. The positive correlation between TP and AvP (Figure 3,

Table 4. Partial correlation results between soil phosphorus and basic physical and chemical properties.

	TP	AvP	PAC	SOM	pH
TP	1	0.417 **	0.094	0.011	0.346 **
AvP	0.417 **	1	0.891 **	−0.202 **	−0.176 **
PAC	0.094	0.891 **	1	−0.204 **	−0.314 **
SOM	0.011	−0.202 **	−0.204 **	1	0.323 **
pH	0.346 **	−0.176 **	−0.314 **	0.323 **	1

Note: ** represents an extremely significant (p < 0.01) correlation.

) suggests that TP influences the availability of phosphorus in soils across land uses. Higher PAC values in cultivated lands suggest that anthropogenic activities, particularly phosphorus fertilization, enhance phosphorus availability.

3.1.3. Soil pH and Phosphorus Availability

The average soil pH is 5.76, ranging from 4.14 to 8.24, indicating mildly acidic to slightly alkaline soil profiles across different land uses (Table 1). Soil pH varies across land uses, decreasing from shrubland to forest, and across soil types, with cambisols exhibiting the highest pH. pH has a significant positive correlation with TP and a significant negative correlation with AvP (Table 4). High pH in the karst area promotes the formation of calcium–phosphate compounds [25,26,27], which may be moderately soluble depending on specific conditions. Low pH promotes phosphorus fixation by aluminum and iron oxides [28], reducing AvP levels even when TP concentrations are high.

3.2. Spatial Variation of Soil TP

The coefficient of variation (CV) for TP is 43.1%, indicating moderate spatial variability (Table 2). Among soil types, ferralsols have the highest CV (52.15%), suggesting considerable heterogeneity, while fluvisols have the lowest CV, indicating a uniform TP distribution (Table 2). Semi-variance analysis using GS+ 9.0 software suggested that Gaussian models fit cultivated land and the total dataset very well, while the spherical model for non-cultivated land is less effective (R² = 0.883). (

Table 5. Semi-variogram function model for soil TP and its corresponding parameters.

Sources	Models	R2	C0	C0 + C	C	C0/C0 + C (%)	A (m)
Cultivated land	Gaussian	0.987	0.0243	0.1216	0.098	20.0	2024
Non-cultivated land	Spherical	0.883	0.0001	0.1112	0.1111	0.1	3880
Total	Gaussian	0.995	0.0290	0.1340	0.1050	21.6	1850

Note, C₀, nugget effect; C, The structural variance; C₀ + C, finite sill value; C₀/C₀ + C, Nugget effect; A, range parameter.

). Determination coefficients (R²) above 0.88 across models indicate a strong fit, with spatial correlation classified as high, particularly in cultivated land due to factors like fertilizer input. Cultivated land and the total dataset have moderate spatial correlation, with substantial structural variance contributions. Non-cultivated land exhibits extremely high spatial correlation (C₀/(C₀ +C) = 0.1%), indicating minimal random variability or error. Non-cultivated land has a larger spatial correlation range (A = 3880 m), implying broader or more homogenous spatial structures compared to cultivated land.

3.3. Soil TP Prediction and Mapping

3.3.1. OLS Model

The OLS model yielded an adjusted R² of 0.40 and an Akaike Information Criterion (AICc) of 6918, indicating that the model accounts for approximately 40% of the variance in total phosphorus (TP) concentration. This suggests that land use, soil type, and topography contribute moderately to the explanation of TP variability at a global scale, though the model’s ability to capture local spatial patterns is limited.

3.3.2. GWR Model

The GWR model substantially improved upon the OLS fit, yielding an AICc of 6861 and a reduction of 57, which signifies a better balance between model fit and complexity. The GWR model exhibited an adjusted R² of 0.45, reflecting its enhanced capacity to account for an additional 5% of the variance in TP concentration compared to OLS. This improvement underscores the importance of incorporating spatial heterogeneity into the modeling process.

Moreover, the global R² value of 0.47 for the GWR model highlights its strong predictive performance for spatial variation in TP concentration. By allowing model parameters to vary locally rather than assuming global uniformity, the GWR model effectively captures localized influences of environmental variables, such as spatial variations in land management practices, soil composition, and microtopographic effects, which are critical drivers of TP distribution.

3.3.3. Comparison of Prediction Models

Table 6. Global Moran’s I, Z-statistic, and p-value statistical analysis results of the four models.

Models	OLS	GWR	MM	MC
Moran’s I	0.552	0.568	0.365	0.365
Z-scores	3.742	3.854	2.479	2.479
p-value	0.00018	0.00012	0.01318	0.01318

shows that the MM and MC models do not fully eliminate the spatial autocorrelation in the model residuals, and they can still be used for OK interpolation. Additionally, as the relationship between a predicted value and the average of a set of predicted values (Z-score) decreases, the residuals become more concentrated, which may indicate smoother data.

Five prediction models were used to evaluate TP spatial interpolation (Figure 4), with mean TP predictions across models closely aligning with observed values (796 mg·kg⁻¹). The GWRK model exhibited superior performance, with a prediction range from 189 to 1593 mg·kg⁻¹ and an average of 794 mg·kg⁻¹, closely matching observed TP levels. The GWRK map displays higher phosphorus concentrations on both sides of the lower middle part, indicated by darker shades (higher TP ranges) (Figure 4c). The presence of various color bands suggests a strong sensitivity to local spatial variations. GWRK appears to capture finer spatial details, reflecting small-scale variations in TP concentrations. This method’s strength is in accommodating spatial heterogeneity, which may explain why areas with slight concentration differences are more visible. The map generated by MM_OK also shows variation across the region but appears smoother than GWRK (Figure 4d). The transitions between different TP levels are less distinct, implying a generalized approach that may not capture very localized variations. MC_OK seems to provide results comparable to GWRK in most areas, particularly with the prominent yellow and orange zones (mid-level TP concentrations) (Figure 4e). However, it also has smoother transitions like MM_OK, with less precise distinctions between small spatial units.

3.3.4. Model Accuracy Assessment

The internal validation results of the five models are shown in

Table 7. Accuracy assessment of predicted TP concentrations using different models.

Models		Internal Verification			External Verification
	MSE	RMSS	RMS	ASE	MAE	RMSE	Accuracy (%)	r	R2
OK	0.59	196.12	3279.88	66.00	153.3	211.7	62.0	0.71	0.62
RK	−0.0008	1.08	219.11	200.14	156.9	208.7	79.4	0.68	0.63
GWRK	−0.0013	1.08	214.21	189.68	145.9	197.4	80.6	0.80	0.67
MC_OK	0.0088	1.09	207.84	194.30	149.2	198.9	79.7	0.70	0.66
MM_OK	−0.0209	1.06	208.25	193.55	153.4	205.1	80.0	0.80	0.64

. RK’s MSE is closest to 0, suggesting it is highly accurate internally. GWRK also performs well with an MSE of −0.0013, close to RK. MC_OK (0.0088) and MM_OK (−0.0209) have slightly higher MSEs, while OK has the highest MSE (0.59), indicating poorer accuracy in internal predictions. RMSS values close to 1 indicate well-calibrated models. MM_OK exhibits the RMSS closest to 1 (1.06), followed by RK and GWRK (both at 1.08) and MC_OK (1.09). OK has the lowest RMSS (196.12), suggesting it deviates more from ideal standardization compared to the other models. Lower RMS values indicate lower error magnitudes in predictions. MC_OK achieves the lowest RMS (207.84), closely followed by MM_OK (208.25) and GWRK (214.21). RK has a higher RMS (219.11), while OK has a significantly higher RMS (3279.88), indicating much poorer internal predictive accuracy. GWRK has the lowest ASE (189.68), suggesting its errors are more consistent. ASE of MM_OK (193.55) and MC_OK (194.30) also performs well, while RK has the highest ASE (200.14), indicating greater variation in its prediction errors.

MAE reflects the average absolute difference between predictions and actual values. GWRK has the lowest MAE (145.9), meaning it has the smallest average error in external predictions. MC_OK (149.2) and MM_OK (153.4) follow closely, showing reasonable accuracy. OK and RK have higher MAEs, suggesting greater average error. Lower RMSE values indicate better predictive accuracy for external data. GWRK again performs best with the lowest RMSE (197.4). MC_OK has a slightly higher RMSE (198.9), while MM_OK is a bit higher still (205.1), followed by RK (208.7) and OK (211.7). GWRK and MM_OK both have the highest accuracy (80.0%), suggesting these models have the highest accuracy in external predictions. MC_OK is slightly lower at 79.7%, followed by RK (79.4%) and OK with a significantly lower accuracy (62.0%). GWRK has the highest R² (0.67), suggesting it explains a substantial proportion of variance in TP concentration. MC_OK follows with an R² of 0.66, and MM_OK has 0.64, both indicating good predictive power but slightly less than GWRK. RK and OK have the lowest R² values, indicating weaker predictive performance.

4. Discussions

4.1. Analysis of Influencing Factors on the Spatial Distribution of Soil TP in Karst Areas

The observed variability in TP across land uses and soil types highlights a complex interplay between human activities and soil properties. Elevated TP in agricultural lands, such as dry land and paddy fields, contrasts with lower levels in forests and shrublands, reflecting the impact of fertilization and land management versus natural cycling in undisturbed soils [29,30]. Notably, phosphorus concentrations in dry lands are higher than in paddy fields because irrigation water during rice growth increases phosphorus mobility [30]. The distribution of AvP across land uses, with higher levels in cultivated areas, indicates the influence of fertilization practices. Variability in AvP—54% of samples between 10 and 40 mg·kg⁻¹ and 39.58% below 10 mg·kg⁻¹—reflects differential fertilizer application and soil adsorption capacity. This uneven distribution is further supported by high PAC values in cultivated land, suggesting that repeated fertilization and soil management boost phosphorus desorption and reduce the fixation ratio. A strong positive correlation between TP and AvP concentration implies that rising TP generally increases AvP, though it also presents risks for runoff and eutrophication [31]. Thus, balanced phosphorus inputs and careful management are essential in agriculture to mitigate these environmental impacts.

Shrublands predominate on shallow carbonate soils, while forests grow in more fertile non-karst soils. TP, AvP, and PAC levels are higher in cultivated lands than in shrublands or forests [29,30]. A negative correlation between SOM and AvP also suggests phosphorus availability is impacted by organic matter mineralization [30], along with pH, mineral composition, and microbial factors [32]. Higher TP levels in cambisols than those in ferralsols suggest a significant role for soil parent material. Cambisols, with their higher calcium concentration, promote phosphorus retention due to calcium–phosphate binding [25,26,27]. The range in soil pH, from mildly acidic to slightly alkaline, significantly impacts phosphorus retention and availability. Higher TP levels in alkaline cambisols reflect enhanced retention, while phosphorus in acidic soils is more prone to fixation by aluminum and iron oxides, reducing its availability [27,33].

4.2. Spatial Variability and Prediction Model Performance

The observed moderate spatial variability of TP, indicated by a CV of 43.1%, reflects the inherent heterogeneity of karst landscapes. In these landscapes, differences in topography, soil composition, and parent material lead to significant variations in soil properties across small distances. Both natural and anthropogenic factors influence this variability, shaping the distribution and behavior of phosphorus [34]. The spatial patterns captured by gaussian and spherical models for cultivated and non-cultivated areas, respectively, suggest that these models effectively characterize TP heterogeneity across diverse land use types. Determination coefficients exceeding 0.88 further affirm that land use and soil type are crucial determinants of TP distribution, supporting previous findings on phosphorus variability in heterogeneous environments [35,36].

Differences in block base ratios between cultivated and non-cultivated areas highlight the role of agricultural practices, such as fertilizer application, in intensifying TP spatial heterogeneity. For non-cultivated regions, the data indicate that intrinsic factors like topography and soil type govern TP distribution, suggesting a need for conservation-oriented land management practices to maintain balanced phosphorus levels and safeguard environmental quality.

In karst areas, the intricate interaction between lithological variations and fragmented topography influences soil properties [22]. Among the evaluated prediction models, the GWRK model achieved the highest accuracy in estimating TP distribution in karst landscapes, with an R² of 0.67 and the lowest MSE at −0.0013. This robust accuracy, confirmed by both internal and external indicators (e.g., low MSE, RMSE, and MAE, alongside a high adjusted R²), demonstrates GWRK’s capacity to capture the spatial heterogeneity of TP within the complex spatial patterns of karst landforms. This is consistent with previous studies that used GWRK combined with local environmental data to better predict fine-scale TP spatial changes [13,37]. The GWR model’s improved performance suggests that regional-scale management strategies tailored to local conditions could significantly enhance the effectiveness of interventions aimed at controlling TP levels.

While the MM_OK and MC_OK models also performed well in TP prediction, they each demonstrated specific strengths and limitations compared to GWRK. MM_OK and MC_OK yielded R² values of 0.64 and 0.66, respectively, presenting viable alternatives for TP mapping. MM_OK’s highest RMSS value (1.06) suggests its suitability for moderately variable landscapes, as its mean-based approach smooths microvariations, making it more applicable to large-scale TP mapping or cases where fine-scale accuracy is less critical. However, this smoothing limits MM_OK’s ability to capture local anomalies, reducing its effectiveness in landscapes with pronounced spatial variability [13]. Conversely, MC_OK exhibited lower RMSE and MSE, suggesting greater accuracy and precision in TP estimation. With its focus on the median rather than the mean, MC_OK is less sensitive to outliers—a key advantage in heterogeneous environments, enhancing prediction accuracy in areas with uneven TP distribution. Nevertheless, the smoothed predictions of MC_OK may lack the finer details needed for precision agriculture applications.

Unlike MM_OK and MC_OK, GWRK integrates local environmental variations, enabling predictions that more accurately reflect the spatial heterogeneity of TP in karst landscapes. The OK model, in contrast, performed poorly, with an R² of 0.62 and significant prediction errors (Table 7), indicating its limitations in capturing TP variability in highly heterogeneous karst regions. The limited range of TP values predicted by the OK model suggests that it applies a strong smoothing effect, thereby smoothing local variations at the expense of accuracy [13]. These findings underscore the strengths of models like GWRK and MM_OK in addressing the complex phosphorus distribution within environmentally diverse karst landscapes [12].

5. Conclusions

Land use and soil type significantly affect the spatial distribution of soil TP content in karst areas. Agricultural lands such as dry land and paddy fields have higher TP levels due to fertilization. In contrast, TP levels are lower in forests and shrublands in natural contexts. Also, soil type and pH have a strong impact on phosphorus availability. Cambisols in karst areas are rich in calcium and retain more TP, whereas alkaline soils tend to immobilize phosphorus, making them less available.

Among the models assessed, the GWRK model stands out for its superior predictive accuracy and stability, proving highly effective for capturing the spatial heterogeneity of TP in karst landscapes characterized by complex topography and varied land use. The model’s capacity to incorporate localized environmental variables allows it to provide detailed spatial predictions, making it particularly suitable for highly variable environments. MM_OK offers a broader distribution profile, which may be advantageous in applications where broader trends are more relevant than fine-scale details, while MC_OK balances the fine detail of GWRK with the smoothing effect of MM_OK, making it suitable for scenarios requiring moderate spatial precision, while the lower accuracy of RK and OK highlights the limitations of models that do not account sufficiently for spatial variability. Overall, these findings provide valuable insights for guiding agricultural productivity and ecological conservation in sensitive landscapes and reinforce the potential of spatially adaptive kriging models in developing sustainable phosphorus management strategies tailored to local environmental conditions, informing sustainable fertilization practices for farmers, monitoring soil quality for scientists, and guiding policy decisions for regulators in karst areas.