The Spatiotemporal Variability of Soil Available Phosphorus and Potassium in Karst Region: The Crucial Role of Socio-Geographical Factors

: The contents of soil available phosphorus (AVP) and potassium (AVK) in karstic moun-tainous agricultural areas have changed rapidly in recent decades. This temporal variation displays strong spatial heterogeneity due to these areas’ complex topography and anthropogenic activities. Socio-geographical factors can reflect the changes in the natural environment caused by human beings, and our objective is to enhance understanding of their role in explaining the changes of AVP and AVK. In a typical karst region (611.5 km 2 ) with uniform soil parent material and low climatic variability, 255 topsoil samples (138 in 2012 and 117 in 2021) were collected to quantify the temporal AVP and AVK changes. Random forest (RF) and partial dependence plot analyses were conducted to investigate the responses of these changes to socio-geographical factors (distance from the nearest town center [DFT] and village density [VD]), topography, biology, and landscape pattern indexes. The mean values of AVP (48.25 mg kg −1 ) and AVK (357.67 mg kg −1 ) in 2021 were significantly ( p < 0.01) higher than those in 2012 (28.84 mg kg −1 and 131.67 mg kg −1 , respectively). Semi-variance analysis showed strong spatial autocorrelation for AVP and AVK, ranging from 7.29% to 10.95% and 13.31% to 10.33% from 2012 to 2021, respectively. Adding socio-geographical factors can greatly improve the explanatory power of RF modeling for AVP and AVK changes by 19% and 27%, respectively. DFT and VD emerged as the two most important variables affecting these changes, followed by elevation. These three variables all demonstrated clear nonlinear threshold effects on AVP and AVK changes. A strong accumulation of AVP and AVK was observed at DFT < 5 km and VD > 20. The AVP changes increased dramatically when the elevation ranged between 1298 m and 1390 m, while the AVK changes decreased rapidly when the elevation ranged between 1350 m and 1466 m. The interaction effects of DFT and VD with elevation on these changes were also demonstrated. Overall, this study examined the important role of socio-geographical factors and their nonlinear threshold and interaction effects on AVP and AVK changes. The findings help unravel the complex causes of these changes and thus contribute to the design of optimal soil phosphorus and potassium management strategies.


Introduction
Phosphorus (P) and potassium (K) are crucial elements of plant nutrition, and their deficiencies can contribute to reduced plant growth, yield, and health [1][2][3][4].Unlike nitrogen with an atmospheric source, terrestrial plants obtain most of their P and K only from soils [5].The portions of total P and K that are readily available for absorption by plant roots are referred to as soil available P (AVP) and soil available K (AVK) [6], which are very limited in natural soils [3,5,7].To prevent crop P and K deficiencies, farmers in major agricultural countries (e.g., China) have been intensively applying chemical P/K fertilizers over the past decades [8][9][10].Despite boosting crop yields, the excessive application of these fertilizers poses many ecological and economic concerns, such as resource wastage, low P/K efficiency, an unbalanced composition of soil nutrients, and P leaching [3,7,[11][12][13][14].The spatiotemporal variability of AVP and AVK in agroecosystems is determined by various environmental variables, such as natural and socio-geographical factors [1][2][3][4]7,[15][16][17].Therefore, understanding the dynamics of AVP and AVK and the factors influencing their variability in agricultural systems is important to achieve optimal soil management for crop production and environmental quality.
The effects of natural factors (e.g., topography) on soil P and K have been well documented in previous studies [1][2][3]7,15].For example, elevation can influence the weathering of bedrock and the transport and deposition of soil P and K [1,3,18].The spatiotemporal variability of soil P and K in agro-ecosystems can also be influenced by human beings through intentional actions (e.g., land use type and fertilization) or incidental effects (enrichment sources of P and K [e.g., rubbish dumps and manure piles]) [4].However, except for land use types and their derived variables (e.g., landscape pattern indices), other anthropogenic factors (e.g., fertilizer, enrichment sources) are generally difficult to investigate and quantify [19,20].By contrast, socio-geographical factors (e.g., distance from the nearest town center and village/population density), which reflect the transformation of the natural environment triggered by the productive activities of humans, can be accurately investigated and spatialized.Previous studies have also revealed the great potential of these factors in driving AVP and AVK variability [4,16,17,21,22].For example, Turner and Hiernaux [4] and Samaké [17] found higher levels of soil P and K in fields located close to villages and attributed such findings to the large amount of agricultural inputs and net nutrient transfers from outlying areas.Population density is also positively correlated with soil K and P content [16,21,22].Overall, socio-geographical factors have an important influence on AVP and AVK, but their relative importance compared with natural factors remains unclear.This outstanding issue might prevent an integrative understanding of the critical drivers of soil P and K dynamics in agricultural areas.
Many scholars suggest that the relationship between natural factors and soil nutrients (e.g., soil P) [2,23] and between socio-geographical factors and soil pollutants [24,25] may be nonlinear or has a threshold type.However, only a few studies have verified the nonlinear threshold effects of socio-geographical factors on soil P and K [17].Samaké et al. [17], who identified the potential nonlinear threshold effects of socio-geographical factors on soil P and K, found that soil P and K in fallow plots do not significantly differ when the distance from the village exceeds a specific threshold (500 m and 1000 m, respectively).However, they identified these thresholds by discretizing the independent variables [17], which cannot sufficiently reflect the sudden changes in the relationship between sociogeographical factors and soil P and K [25].Other studies suggest that the effect of one environmental factor on soil P and K is often driven by another factor [1,4,20,22].For example, after controlling for the locational factor, the soil P content remains negatively related to the distance to the village, although the significance of such a relationship has declined [4].Overall, socio-geographical factors may have nonlinear threshold and interaction effects on AVP and AVK.Exploring such potential effects may provide insights into the spatial variation mechanisms of soil P and K.
As the largest continuous karst region in the world, the southwest China karst region is also one of the least undeveloped areas in China [26].To escape poverty, tobacco is widely cropped by local farmers as an important source of income.In recent years, the rapid accumulation of P and K in tobacco soils has been witnessed [27,28].However, the influence of socio-geographical factors, as important factors characterizing anthropogenic activities, on soil P and K accumulation remains unclear.Therefore, the objectives of this study are to (1) investigate the temporal changes in AVP and AVK; (2) determine the impact of socio-geographical factors on these changes; and (3) explore the nonlinear threshold and interaction effects of socio-geographical factors on these changes.This study was conducted in a karstic sub-region (611.5 km 2 ) with uniform soil parent material and low climatic variability.Especially, random forest (RF) modeling and partial dependence plot (PDP) analysis were employed to examine the effects of socio-geographical factors on AVK and AVK changes.

Study Area
The study area is located in the karst region of southwest China (109°9′-109°32′ E, 30°29′-30°49′ N), which covers an area of 611.5 km 2 (Figure 1).This site has a typical subtropical monsoon climate with a mean temperature of 18 °C, mean annual precipitation of 1132 mm, mean annual total evaporation of 870.4 mm, relative humidity of 67%, and sunlight of 1242.6 h.The topography is mountainous, with an elevation varying between 244 m and 2072 m and a slope ranging between 0° and 75°.The Triassic system Daye formation carbonate rocks cover about 95% of the study area.The population density in the area is relatively low, at 182 person km −2 , while the distribution of settlements is scattered.Planting flue-cured tobacco (Nicotiana tabacum L.) is the main source of income for local farmers.Tobacco fields are mainly distributed on both sides of the roads.The growing season for tobacco lasts from May to September, and most tobacco fields are not cropped during the other months.The amount of nitrogen fertilizer consumed in this area reaches approximately 111.45 kg ha −1 every year, while those of P and K fertilizers are 90 kg ha −1 -120 kg ha −1 and 270 kg ha −1 -315 kg ha −1 , respectively.

Soil Sampling and Chemical Analysis
A total of 138 and 117 tobacco field topsoil (0 cm-20 cm) samples were collected from the study area in 2012 and 2021, respectively, before tobacco planting and fertilization (March).The sampling densities in 2012 and 2021 were 0.23 and 0.19 points per square kilometer, respectively.The average sampling distance in 2012 and 2021 was 0.59 km and 0.45 km, respectively.According to a highly influential review [29], the sampling density for studies within a spatial extent of 10 4 km 2 is usually between 0.01 and 0.1 points per square kilometer.According to the recent national-scale soil surveys in China (i.e., the third national soil census), the sampling distance of the topsoil cropland samples in this area is 1 km.Therefore, the sampling densities and average sampling distances in 2012 and 2021 fall within the commonly used ranges.
The sampling sites were selected from local unique organizational units, namely, basic planting units (BPU).Each BPU consists of multiple neighboring tobacco fields with a total area exceeding 8 hectares.These BPUs are randomly dispersed across the entire study area, and our sampling activities in 2012 and 2021 encompassed all of them.Moreover, in each BPU, at least one representative tobacco field with 600 m 2 to 1000 m 2 was selected for sampling.The spatial location of these sites was recorded using a handheld GPS receiver, and 3 to 5 subsamples were collected and well mixed to represent the AVP and AVK information for each site.About 2 kg of the well-mixed sample was air-dried in a laboratory and then gently crushed to pass through a nylon sieve.The AVP content was determined via HCl-H2SO4 extraction and colorimetry, while the AVK content was measured via ammonium acetate extraction and colorimetry [30].

Environmental Variables
Based on classical theories (i.e., Jenny's factorial model of soil formation and Scorpan model) [31,32] and existing studies [1][2][3][4]17,22,25,[33][34][35][36], a total of 13 environmental variables representing topography, biology, landscape pattern, and socio-geographical factors were used to assess the spatiotemporal variability of AVP and AVK (Table 1 and Figure 2).All samples were obtained from the same geological formation (i.e., Triassic system Daye formation) (Figure 1).The annual temperature of sampling sites is closely related to altitude, and the soil samples received similar amounts of precipitation with a coefficient of variation (CV) of 0.05%.Consequently, geological and climatic variables were excluded from the environmental variable set.The annual normalized difference vegetation index (NDVI) of Landsat-8 images taken from 2013 to 2021 was computed from the Google Earth Engine platform using the maximum value composite procedure [37].In 2012, no available NDVI products from the Landsat missions were found.Specifically, the satellites from Landsat-1 to Landsat-5 ceased acquiring images in 2012.The Landsat-6 satellite launch failed.Landsat-7 products were not considered due to well-known striping issues (i.e., wedge-shaped scan-to-scan gaps).The Landsat-8 satellite had not yet been launched in 2012.The average NDVI values were calculated at a 30 m resolution.An elevation map (ASTER GDEM v3) with a resolution of 30 m was used in this paper [38].Six terrain indicators, namely, elevation (Ele), aspect (Asp), slope (Slp), plane curvature (PLC), profile curvature (PRC), and topographical wetness index (TWI), were calculated using SAGAGIS v7.2.0 [28].
Four classical landscape metrics, namely, mean patch size (MPS), landscape contagion index (CONTAG), aggregation index (AI), and landscape shape index (LSI), were employed to provide landscape pattern information [19,39].A higher MPS corresponds to a larger average area of patches in the landscape, higher CONTAG and AI correspond to higher connectivity, and larger LSI corresponds to more complex landscape shapes [39].The landscape metrics of each parcel were calculated at the landscape level using the moving window method and based on the land use map with a 30 m resolution and the FRAG-STATS 4.1 software [19,39].The land use map, surveyed by the Chongqing Provincial Department of Land and Resources, was not open to the public.The original data for the land use map is a 1:10,000 vector file, but we only had access to the 30 m resolution raster file derived from the vector file.The optimal sliding window of 840 m was determined through extensive testing following the procedures in Liu et al. [19].
The raster maps of village density (VD) and distance from the nearest town center (DFT) were generated using the point density and near tool of ArcGIS v10.3 [40], respectively.Village vectors were derived from the land use map taken in 2016, and the town centers were determined from the land use map and the Google Earth Pro platform.

Classical Statistics
Descriptive statistical analysis (i.e., mean, median, minimum, maximum, standard deviation [SD], coefficient of variation [CV], and skewness) was performed to characterize the AVP and AVK contents and environmental variables.The independent sample t-test was applied to assess the differences in the AVP and AVK contents between the two periods after natural log transformation.All statistical analyses were conducted in SPSS v22 [41].

Spatiotemporal Analysis of AVP and AVK
The spatial variability of AVP and AVK in 2012 and 2021 was identified through a semi-variance analysis in GS+ 7.0 software [42].The software automatically selects the optimal theoretical model, as indicated by the smallest residual sum of squares (RSS) and the largest coefficient of determination (R 2 ), and the corresponding parameters, namely, nugget variance (Co), sill variance (Co + C), and range.The ratio of Co and Co + C represents the nugget effect, which is an important parameter for describing spatial autocorrelation.Nugget effects of <25%, 25%-75%, and >75% represent strong, moderate, and weak spatial autocorrelations, respectively [42].Moreover, based on the basic parameters derived from the semi-variance analysis, the raster maps of AVP and AVK in 2012 and 2021 were produced via ordinary kriging interpolation in the geostatistical analyst tools of ArcGIS v.10.3 [40].Using geostatistical analyst tools, the directions and ranges of the major axis and minor axis were examined, and the anisotropy ratio (the ratio of the major axis range to the minor axis range) was calculated.If the anisotropy ratio is greater than 1, it means that soil properties exhibit greater spatial variability in the major axis direction, and vice versa [43].
The maps of AVP and AVK changes were then generated by performing raster subtraction on the predictions captured in 2012 and 2021.Afterward, information on the AVP and AVK changes was extracted from the 255 samples collected in 2 years using the extract multi values to points tool in ArcGIS v.10.3 [40].The spatial clustering of AVP and AVK changes was detected using the local indicators of spatial association (LISA) analysis developed by Anselin [38].The following clustering patterns were identified: (i) no significant spatial dependence; (ii) high-high, high value in a high-value neighborhood; (iii) high-low, high value in a low-value neighborhood; (iv) low-high, low value in a highvalue neighborhood; and (v) low-low, low value in a low-value neighborhood [28,44].

RF Modeling and Accuracy Comparison
The effects of topographical factors, landscape pattern factors, and socio-geographic factors on AVP changes and AVK changes were analyzed and compared by RF modeling.RF is a typical bagging ensemble learning model that combines the output of multiple decision trees to reach a single result [45].RF is developed based on two important ideas, namely, the bagging method and feature randomness.In bagging, a large number of trees are built based on bootstrap sampling of training data.All tree predictions are aggregated through averaging, and these averages are taken as the final predictions.Feature randomness means that only a subset of the original covariate set is considered in the tree-splitting process.RF has two main hyperparameters, which need to be set before training.The hyperparameter "n_estimators" (the number of trees in the forest) was set to 500 since our repeated experiments showed that this was sufficient to obtain stable results.The hyperparameter "max_features" (the number of maximum features provided to each tree) was set to the default value of the square root of the total number of covariates.
Five RF models were constructed based on different combinations of environmental variables.Model-A includes NDVI and the topographic variables; Model-B includes NDVI, the topographic variables, and landscape patterns; Model-C includes NDVI, the topographic variables, and socio-geographical factors; and Model-D includes all variables.Given that the reduction in the prediction performance of Model-D for AVK and AVP changes was due to information redundancy (Section 3.3) [46], the Boruta algorithm was used to select the optimal combination of variables for Model-E.The Boruta algorithm is a popular variable selection method used in machine learning, as proposed by Kursa and Rudnicki [47].The main steps of the Boruta algorithm are as follows: (1) enlarge the variable database and add >5 shadow attributes to each variable; (2) perform shuffle procedures on all attributes to remove correlations; (3) execute RF modeling based on the expanded data and calculate Z-scores (Equation ( 1)) for each attribute; (4) select variables based on Z-scores and generate a new dataset; (5) repeat the above steps until assigned criteria.
where OOBacc and OOBacc_shuffle are the out-of-bag accuracy of each decision tree with and without shuffle attributes, respectively.The mean and SD represent the mean and standard deviation of OOBacc-OOBacc_shuffle, respectively.A 10-times 10-fold cross-validation (100 replications) technique was used to evaluate the five RF models.The prediction performance of these models was assessed by root mean square error (RMSE) (Equation ( 2)) and R 2 (Equation (3)).
where xi and yi represent the measured and predicted values of the i-th sample, respectively, n is the number of soil samples, and ym is the mean prediction values.The RF model and Boruta algorithm were developed based on the Scikit-learn and BorutaPy libraries embedded in Python v3.8, respectively.

Variable Importance and Partial Dependence Plot Analysis
The effects of environmental variables on AVP changes and AVK changes were revealed by the nodes' impurity decrease method [45].In 100 repetitions of each RF model, the mean values of the relative contribution of environmental variables to sAVP changes or AVK changes can be identified.The contribution of the most important variable was set to 100%, and the relative importance of the other variables was scaled.
PDP analysis was performed to visualize the complex effects of environmental variables on AVP changes and AVK changes.This method can reveal the marginal effect of one or two features on the predicted outcome of the RF model [48,49], which is accomplished in this study by a one-way PDP analysis and a two-way PDP analysis, respectively.Considering the sample set size, all samples were used for PDP analysis.PDP analysis was executed in Python v3.8.

Descriptive Statistics
The AVP and AVK levels in 2012 and 2021 are shown in Table 2.The mean ± standard deviation of AVP in 2012 and AVP in 2021 were 28.84 ± 21.11 mg kg −1 and 48.25 ± 24.76 mg kg −1 , respectively, while their mean values were 131.67 mg kg −1 (32.12 mg kg −1 -371.44 mg kg −1 ) and 357.34 mg kg −1 (97.54 mg kg −1 -782.79 mg kg −1 ).The distributions of AVP and AVK in 2012 and 2021 are positively skewed (right-skewed).According to the independent samples t-test, the AVP and AVK contents increased significantly (p < 0.01) by 19.41 mg kg −1 and 225.67 mg kg −1 from 2012 to 2021, respectively.As indicated by their CV, the AVP and AVK contents in these two periods exhibited high variability (CV > 35%) [50], with the variabilities in 2021 being lower than those in 2012.The tobacco fields in the study area were distributed at high elevations (>1000 m) and in any slope direction and terrain curvature (concave or convex).The mean values of NDVI, Slp, TWI, MPS, AI, CONTAG, and LSI were 0.65, 7.55°, 10.2, 6.65 ha, 81.81%, 25.32%, and 1.73, respectively.The distances of tobacco fields from the town center ranged from 0.65 km to 7.77 km, with an average distance of 3.65 km.The village density was 29.44 ± 21.38.The distributions of NDVI, Ele, PLC, PRC, and AI are negatively skewed (left-skewed), while those of the other environmental variables are positively skewed (right-skewed).NDVI, Ele, AI, and LSI showed low variability (CV < 15%), and all other environmental variables demonstrated strong variability with CV values ranging from 37.33% to 128.78%.

Spatio-Temporal Variability of Soil Properties
AVP2012 was optimally fitted by the spherical semi-variogram model, while AVP2021, AVK2012, and AVK2021 were best fitted by the exponential model (Figure 3 and Table S1).Following the classification described in Cambardella et al. [42], a strong spatial autocorrelation was observed in AVP2012, AVP2021, AVK2012, and AVK2021 (nugget effect [Co/(Co+C)] < 25%).The spatial autocorrelation ranges of AVP2021 exceeded those of AVP2012, thereby suggesting that the spatial distribution of the former was more homogeneous than that of the latter [42].Similarly, compared with AVK2021, the spatial distribution of AVK2012 was more homogeneous.Besides, the anisotropy assessment results showed that the directions of the major axis were northwest-southeast, and the anisotropy ratios of AVP and AVK were all less than 1 (Table S1).This implied that AVP and AVK exhibited stronger variations in the minor axis direction (i.e., southwest-northeast), which was consistent with the topographical direction.The ranges for AVP and AVK in the two periods were greater than or close to the sampling distance, implying that our sampling system was adequate for detecting the spatial patterns of AVP and AVK [52][53][54].The optimal parameters identified in the semi-variance analysis were employed as input parameters for ordinary kriging interpolation.The spatial distributions of AVK and AVP differed between the two periods (Figure 4a,b,d,e).For example, the AVP content in the central part of the study area was higher than that in the eastern and western parts in 2021 yet was generally similar to that in the eastern and western parts in 2012.The spatial distribution of AVP and AVK changes (Figure 4c,f) was similar in the northern, eastern, and western parts of the study area.The southern part exhibited high AVP changes and low AVK changes, while the AVP and AVK changes in the eastern and western parts predominantly displayed low-low clustering.In the southern part, a high-high clustering of AVP changes and a slight low-low clustering of AVK changes were observed.

Model Performance
Statistics on the prediction performance of all five models for AVP and AVK changes are summarized in Table 3.Using only terrain variables, Model-A explained 14% and 15% of the variance in AVP and AVK changes, respectively.Model performance was slightly improved by the addition of landscape pattern indices and significantly improved by the addition of socio-geographical factors.Model-D employed all variables to predict AVP and AVK changes and obtained R 2 values of 0.36 ± 0.13 and 0.40 ± 0.11 and RMSE values of 7.34 ± 0.77 mg kg −1 and 58.91 ± 6.69 mg kg −1 , respectively.Model-E obtained the best prediction performance for AVP and AVK changes, with the highest R 2 of 0.41 ± 0.14 and 0.44 ± 0.13 and the lowest RMSE of 7.11 ± 0.79 mg kg −1 and 56.89 ± 6.57 mg kg −1 , respec-tively.Among the 13 environmental variables, the Boruta algorithm retained 8 and 9 covariates for AVP and AVK changes, respectively (Table S2).NDVI, DFT, VD, Ele, Slp, PLC, and AI were retained in the prediction of AVP and AVK changes.Other variables used for predicting AVP changes included TWI in Model-E, while those used for predicting AVK changes included CONATG and Asp.

Relative Importance of Covariates and PDP Analysis
Model-E was iterated 100 times, and the mean relative importance of the environmental variables in predicting AVP and AVK changes was recorded and presented in Figure 5. DFT and VD were ranked as the top two most important variables, followed by elevation.The relative importance of the other variables was relatively low, and their contribution was generally around 25% of the most important variable.

Nonlinear Threshold and Interaction Effects of Key Factors
Given that the relative importance of DFT, VD, and Ele far outweighed those of other factors, their associations with AVP and AVK changes were visualized via PDP analysis (Figure 6).Generally, AVP and AVK changes were negatively related to DFT and positively correlated with VD.AVP changes increased along with elevation, while AVK changes showed a decreasing trend.All three variables exhibited non-linear threshold effects on AVP changes and AVK changes.Specifically, the effects of DFT and VD on AVP and AVK changes were mutated at 5 km and 20, respectively.These changes flattened when DFT < 5 km and VD > 20, decreased rapidly when DFT > 5 km, and increased drastically when VD ranged between 0 and 20.The AVP changes also increased dramatically when Ele ranged between 1298 m and 1390 m, while the AVK changes decreased rapidly when Ele ranged between 1350 m and 1466 m.The bivariate interaction effects of the three most important environmental variables on AVP and AVK changes are visualized in Figure 7. Tobacco soils exhibited lower AVP and AVK accumulation in areas with VD < 20, and in areas with DFT > 5 km and VD ranging between 20 and 34.The effects of DFT and VD on AVP and AVK changes significantly varied across different elevations.Specifically, the relationship between DFT and AVP changes was roughly stable when Ele < 1298 m, while the effect of DFT on AVP changes exhibited a large abrupt change near DFT = 5 km in the area with Ele > 1390 m (Figure 7b).When Ele < 1466 m, the relationship between DFT and AVK changes fluctuated dramatically around DFT = 5 km, and when Ele > 1466 m, such a relationship was roughly stable (Figure 7e).In regions with Ele < 1390 m, the relationship between VD and AVP changes fluctuated dramatically when VD ranged between 20 and 33, but when Ele > 1390 m, such a relationship abruptly changed near VD = 20 (Figure 7c).At any altitude, the relationship between VD and AVK changes fluctuated at VD = 20, and this fluctuation was stronger at altitudes less than 1466 m (Figure 7f).

Temporal Changes of AVP and AVK from 2012 to 2021
From 2012 to 2021, the AVP and AVK contents in tobacco fields increased significantly (p < 0.01) by 67.3% and 171.39%, respectively (Table 2), which were much higher than those in Chinese arable land.Specifically, the AVK content of China's arable land increased by 15.1% between 1988-2007 and 2008-2018 [55].Such a rapid increase in AVP and AVK content was inextricably linked to long-term heavy fertilization [27,28].Given that tobacco is a high-value crop, farmers tend to apply higher amounts of chemical fertilizers in tobacco production than in staple crop production [27,28,55].For example, the annual application of chemical K fertilizers for tobacco production is 2.16 to 2.52 times the average level used in Chinese farmland [55].Moreover, the AVP and AVK contents in 2021 reached 48.25 mg kg −1 and 357.34 mg kg −1 , respectively, which far exceeded the average values for Chinese farmland (27 mg kg −1 and 139 mg kg −1 , respectively) [55,56].Overall, the growth rates and current levels of AVP and AVK in tobacco fields were extremely high, which might pose numerous risks for tobacco production and the ecosystem, including low P/K use efficiency, resource wastage, and P leaching [11,12].

Effect of Socio-Geographical Factors on AVP Changes and AVK Changes
The addition of VD and DFT to the RF models greatly improved the prediction of AVP and AVK changes (Table 3).Moreover, VD and DFT obtained higher relative importance than the other factors in the optimal model (Figure 5).These findings underscore the importance of socio-geographical factors in controlling the variation in AVP and AVK changes.The socio-geographical factors can reflect potential human impacts on soil P and K [4,17,21,22,57].Turner and Hiernaux [4] found that the distance to the nearest village exerts a greater influence on soil fertility parameters (e.g., soil P and K) than landscape position and soil physical properties.Similarly, Samaké et al. [17] revealed the close relationship of AVP and AVK contents with distance from villages in Sahel, Mali.The accumulation of P in soils in the European Union is closely related to the intensity of the livestock industry [57].Road density and population density exert a greater direct influence on AVK content than topography and annual precipitation [21,22].
Humans always input P and K into agricultural soil through intentional actions (e.g., fertilizer application) and incidental behaviors (e.g., migration and deposition of chemical elements from enrichment sources controlled by precipitation and wind) [4].In areas with larger VD and smaller DFT, such intentional actions are frequently observed, and the enrichment sources of P/K (e.g., rural livestock manure and domestic waste) are more abundant [4,17,21,22,58], thus explaining the negative and positive correlations of AVP and AVK changes with DFT (Figure 6a,d) and VD (Figure 6b,f), respectively.Similar findings have been reported in previous studies [4,17,21,22,59].For example, Samaké et al. [17] revealed higher levels of soil P and K in fields located closer to villages and attributed this finding to unknown quantities of household waste and organic manure.Turner and Hiernaux [4] suggested that those fields located closer to villages are more likely to receive more agricultural inputs and net nutrient transfers from outlying areas.Chi et al. [60] reported a negative correlation between distance-related factors (e.g., distance from the mainland) and soil P/K, which might be due to the fact that in uninhabited islands, human impacts on soils are predominantly destructive instead of replenishing, unlike in intensively agricultural areas [4,17].P and K elements from enrichment sources are also susceptible to transport and deposition to soils via precipitation and surface runoff [14], given that fields are typically located in sink areas, such as the foot of slopes or roadsides (Figure S1).Soils in areas located near town centers and in areas with high village density are also rich in black carbon, whose soil particle surface area is three times larger than that of normal soil, thereby reflecting its strong adsorption capacity to reduce the leaching loss of soil P and K [16,21].
In contrast to the linear response reported in previous studies [4,21,22], socio-geographical factors exerted nonlinear threshold effects on AVP and AVK changes in this study (Figure 6).The thresholds for the abrupt shift in the effect of DFT and VD on AVP and AVK changes were 5 km and 20, respectively.The AVP and AVK changes were flattened when DFT < 5 km and VD > 20, decreased rapidly when DFT > 5 km, and increased drastically when VD > 20.In other words, soil P and K accumulated more rapidly at locations near the town center (DFT < 5 km) and in areas with high village density (VD > 20).Fertilizer shops are mainly concentrated in towns.When tobacco fields are located too far from towns, the willingness of tobacco farmers to buy fertilizers may rapidly diminish due to the high transport costs and limited availability of transportation vehicles [58,61,62].Tan et al. [62] found that manure from livestock farms is mainly exported to crop farms within a 5 km distance.Moreover, when the village density exceeds a certain threshold, the P and K from anthropogenic activities may exceed the loading capacity of the soil and thus lead to more drastic increases in AVP and AVK contents.

Influence of Other Environmental Variables
Elevation was ranked as the third most important topographic factor (Figure 5), whose influence on soil P and K variations has been well documented in the literature [1][2][3][34][35][36].Elevation can control the (re)distribution of runoff and energy across landforms and therefore could explain the spatial variability in AVP and AVK changes [1,3].Elevation positively affected AVP changes yet negatively controlled AVK changes (Figure 6).Previous studies have also reported the positive and negative effects of elevation on soil P (positively: He et al. [18]; negatively: Cheng et al. [35]; Nimalka Sanjeewani et al. [15]; Wang et al. [36]) and K (positively: Li et [1]; Nimalka Sanjeewani et al. [15]; Blanchet et al. [3]; Wang et al. [3]; Winzeler et al. [34]).The increase in soil K and P along with elevation is generally due to the enhanced weathering of bedrock [1,18] or decreased microbial activity [15].The negative correlation of soil P and K with elevation is generally attributed to fine fractions (containing high P/K content) or surface-applied fertilizers, which are easily washed away to lower elevations by precipitation and surface runoff [34,36].The final driving direction of elevation and soil P/K may depend on which force predominates in the study area.As the most popular topographic variable, elevation has been widely reported for its nonlinear threshold effects on soil properties [2,63], and similar observations were made in this study.
The index of landscape aggregation/connectivity (e.g., CONTAG and AI) was more important than the shape (e.g., LSI) and area factors (e.g., MPS) (Table S2).These findings are consistent with those of previous studies focusing on the Iranian plateau [33] and the northeastern plain of China [19], respectively.A connected landscape is conducive to the movement of materials across the region [25].DFT and VD played a more important role in controlling the accumulation of soil P and K than topography and NDVI (Figure 5), thus supporting previous findings [4,21,22].Moreover, AI contributed more than any other natural factor except for elevation (Figure 5).These results altogether show that the anthropogenic impacts on soil nutrients (e.g., P and K) in agricultural soils tend to be more profound than the impacts of natural factors [20,64,65].

Interaction Effects of the Three Most Important Environmental Variables
The relationship among DFT, VD, and AVP and AVK changes varied across different elevations (Figure 7).The interaction effects of natural (e.g., elevation) and other factors on soil P and K were confirmed in previous studies [1,20,22].For example, topographic factors have significant effects on soil K, and these effects are enhanced by climate and soil properties [1].Turner and Hiernaux [4] and Wang et al. [22] revealed the interaction effect of socio-geographical factors on soil P and K.The negative correlation between soil P and distance to villages diminishes when locational factors are controlled [4].Population and road densities can affect AVK and K balance through complex interactions [22].Despite the preliminary evidence showing the interaction of socio-economic factors on soil P and K, to our knowledge, the variation in these interactions across different environments has been rarely quantified.Apart from revealing the interaction effects of socio-geographical factors and elevation on AVP and AVK changes, this study is also the first to quantify the differences in these interactions across different altitudes (Figure 7).These interaction effects should be considered to understand the underlying processes leading to AVP and AVK changes in karst regions and to reveal the spatial variability in these changes.

Implications and Limitations
The tobacco soils across the study area generally have high P/K fertility levels and increased magnitude (Table 2).The AVP content in 2021 far exceeded the agronomic soil P threshold (17.3 mg kg −1 ) and slightly exceeded the environmental soil P threshold (40.1 mg kg −1 ) for Chinese soils [13,66].Excessive soil P can severely inhibit the biological ability of crops to efficiently absorb and utilize P [66] and thus exacerbate the risk of non-pointsource pollution [14].Although excessive soil K does not imply an environmental problem, high surpluses of K may lead to an unbalanced composition of nutrients (e.g., low magnesium) in soils [3].Overall, large surpluses of P and K in soils are harmful to tobacco production and the environment.Given that P and K resources are very limited in China [5], local agricultural departments should strengthen soil P and K management, especially in areas located near town centers (DFT < 5 km) and with high village density (VD > 20).According to the PDP analyses (Figure 6), soil P and K were highly susceptible to accumulation in these areas.
The Chinese government launched the third national soil census in 2022 [67], one important task of which was to accurately map soil properties.Improving the accuracy of digital mapping has been an important challenge in areas with complex terrain, such as karst [68,69].Socio-geographical factors can improve the RF modeling accuracy for soil P and K (Table 3) and show a higher relative importance than the other factors (Figure 5).Therefore, socio-geographical factors can be useful indicators for improving the accuracy of soil property mapping.
This study focused on the effects of topography, socio-geographical factors, and landscape pattern indices on AVP and AVK changes in tobacco fields.The optimal model could explain 36% and 42% of the variance in AVP and AVK changes, respectively.The unexplained variations may be attributed to the following reasons: (1) soil erosion in mountainous areas can affect the accumulation of soil P and K [5], but high-precision soil erosion products are not accessible; and (2) similar to Minasny et al. [70], the spatial distribution of AVP and AVK in 2012 and 2021 was initially obtained by kriging interpolation, and then the spatial subtraction method was applied to calculate the temporal changes in AVP and AVK.However, kriging interpolation inevitably introduces uncertainty, thus influencing the modeling of AVP and AVK changes.

Conclusions
This study was conducted in a typical karstic mountainous agricultural region with a uniform geological type and low climatic variability.Results highlighted the temporal changes of AVP and AVK in tobacco fields, revealed the important influence of socio-geographical factors on these changes, and detected the nonlinear thresholds and interaction effects of these factors.The mean values of AVP (48.25 mg kg −1 ) and AVK in 2021 (357.67 mg kg −1 ) were significantly (p < 0.01) higher than those in 2012 (28.84 mg kg −1 and 131.67 mg kg −1 , respectively).Based on the optimal RF model, the socio-geographical factors DFT and VD were identified as the two most important factors affecting AVP and AVK changes, followed by elevation.The results of the PDP analysis revealed that these three socio-geographical factors affected AVP and AVK changes in a nonlinear threshold pattern.AVP and AVK showed stronger accumulation from 2012 to 2021 at locations near the town center (DFT < 5 km) and in areas with high village density (VD > 20).AVP changes increased dramatically when Ele ranged between 1298 m and 1390 m, while AVK changes decreased rapidly when Ele ranged between 1350 m and 1466 m.The interaction effect of socio-geographical factors with elevation on AVP changes and AVK changes was determined, indicating that the relationship among DFT, VD, and AVP and AVK changes varied along with increasing elevation.These findings highlight the important role and complex patterns of socio-geographical factors in controlling the spatiotemporal variations in AVP and AVK.Therefore, these factors should be fully considered in mapping studies of soil P

Figure 4 .
Figure 4. Spatial distribution maps for AVP in 2012 (a), AVP in 2021 (b), AVP changes and LISA cluster (c), AVK in 2012 (d), AVK in 2021 (e), and AVP changes and LISA cluster (f).AVP: soil available phosphorus; AVK: soil available potassium; high-high: high value in a high-value neighborhood; low-low: low value in a low-value neighborhood.

Figure 6 .
Figure 6.Partial dependence plots of the top three environmental variables for AVP changes (a-c) and AVK changes (d-f).The black bands along the x-axis visualize the distribution of data.Note: AVP: soil available phosphorus; AVK: soil available potassium; DFT: distance from the nearest town center; VD: village density; Ele: elevation.

Figure 7 .
Figure 7. Two-way interaction effects of the three most important environmental variables for AVP changes (a-c) and AVK changes (d-f).Note: DFT: distance from the nearest town center; VD: village density; Ele: elevation.

Table 2 .
[51]riptive statistics for soil available phosphorus (AVP), soil available potassium (AVK), and environmental variables.: /: PLC and PRC have both positive and negative values, which typically means that the CVs are meaningless[51].All sample points in 2012 (N = 138) and 2021 (N = 117) were used for the descriptive statistics. Note

Table 3 .
Performance of all five models for AVP changes and AVK changes using random forest models based on 100 repetitions (mean ± standard deviation).