Nutrient-Specific Spatial Dependence and Heterogeneity of Cultivated Soils in an Arid Irrigated Agroecosystem

1. Introduction

Agriculture has remained a foundational sector for China’s food security and rural development. The quality of cultivated land, not merely its area, has largely determined the stability of crop production and the sustainability of agricultural systems [1,2] Cultivated land has provided the material basis for food supply and socioeconomic development. Its quality has been closely linked to soil fertility, crop yield potential, and the resilience of agroecosystems under changing climatic and management pressures [3,4]. Accordingly, understanding and improving soil fertility at management-relevant scales has been a prerequisite for turning broad “cultivated-land protection” goals into measurable gains in productivity and environmental performance[5].

Since the beginning of the 21st century, China has placed increasing emphasis on protecting cultivated land and enhancing its quality. It has advanced policies that have upheld the cultivated-land “red line” while promoting a shift from quantity-oriented management to quality-oriented improvement. While ensuring the preservation of 1.8 billion mu (≈120 million hectares, not 1.2 billion hectares) of cultivated land, national strategies have increasingly stressed improving soil fertility to achieve the goal of “storing grain in the land.” Measures such as scientific fertilization, crop rotation, fallow practices, and pollution remediation support sustainable land use. Legal and accountability frameworks have further clarified protection responsibilities. Technical and ecological approaches have promoted soil improvement, erosion control, and remediation of polluted farmland. Farmland quality grading and cultivated-land quality evaluation frameworks have provided an institutional basis for enhancing and managing quality [6,7,8]. However, these frameworks have only guided effective intervention when the underlying soil indicators have been measured and interpreted with sufficient spatial detail to reflect within-region variability, rather than relying on coarse administrative averages. Systematic evaluation of cultivated-land quality has supported optimization of land-resource allocation and sustainable utilization. It has guided rational land use and conservation and advanced integrated management of cultivated land [9]. Yet the practical value of such evaluations has ultimately depended on whether key soil fertility indicators, especially nutrients, were quantified and mapped at scales relevant to farm management and local governance. Continuous improvement in cultivated-land quality has therefore remained central to China’s transition toward quality-oriented cultivated-land management [10].

Among the indicators used to characterize cultivated-land quality, soil nutrients have been particularly actionable because they directly constrain crop growth and mediate responses to fertilization and irrigation. Fertilizer inputs have been widely recognized to contribute substantially to yield gains in intensive cropping systems, but non-targeted fertilization has often lowered nutrient-use efficiency and increased environmental risks through nutrient accumulation, runoff, and leaching, as well as secondary soil degradation such as compaction or fertility decline. A key reason blanket fertilization has performed poorly is spatial heterogeneity: nutrient availability has varied strongly over short distances due to interactions among soil-forming factors (e.g., parent material, texture, and topography), water management, and field-level practices[11]. What has been well established is that geostatistical approaches can quantify nutrient variability and spatial dependence, and kriging-based interpolation can transform point observations into continuous nutrient surfaces for decision-making. What has remained insufficiently addressed in many regions, especially in arid and irrigated agroecosystems, is the combination of dense, recent sampling; transparent reporting of spatial dependence parameters; and management-oriented interpretation that identifies priority zones while acknowledging uncertainty. In particular, studies have sometimes reported descriptive statistics or model-fit indices while omitting independent validation of spatial predictions, or they have described policy motivations without clearly articulating the scientific gap and testable expectations linking soil processes to spatial patterns.

Wuwei City, Gansu Province, located in the eastern Hexi Corridor, has constituted a representative arid agricultural region that includes both irrigated and rainfed (dryland) systems. Water management and fertilization practices could intensify spatial contrasts in soil fertility. Despite the practical importance of soil nutrient management in Wuwei, cultivated-land assessment and fertilizer decision-making have remained constrained by coarse-scale information that fails to reveal nutrient-deficit zones, enrichment hotspots, or gradients associated with management and landscape context. Addressing this limitation has required spatially explicit quantification of key plough-layer nutrients across cultivated land, together with analysis of their variability, interrelationships, and spatial dependence.

In this study, we integrated 2022 foundational spatial datasets (soil type map, current land-use map, and administrative boundaries) with field-collected cultivated-soil test data to quantify and map the spatial heterogeneity of four key plough-layer nutrients: soil organic matter, total nitrogen, available phosphorus, and available potassium, across cultivated land in Wuwei City. Rather than treating “GIS processing” as the contribution, we used spatial databases and evaluation units as the analytical framework. This linked measured nutrient properties to spatial patterns that could support differentiated land-quality improvement. The study aimed to (i) characterize nutrient variability and distributional features across Wuwei’s cultivated land, (ii) examine nutrient correlations to identify coupled fertility indicators, and (iii) model nutrient spatial dependence and generate continuous nutrient surfaces to support targeted management. We hypothesized that nutrient contents in Wuwei exhibited substantial within-region heterogeneity and that soil organic matter and total nitrogen were strongly positively associated due to coupled nutrient dynamics and/or shared management controls. We also hypothesized that the four nutrients differed in spatial dependence (e.g., nugget-to-sill ratio and correlation range), reflecting different balances between structural controls (soil-forming factors) and field-level management influences. On this basis, we proposed spatially targeted measures to enhance cultivated-land quality, aligned with the observed nutrient patterns across Wuwei City.

2. Materials and Methods

2.1. Overview of the Pilot Zone

Wuwei City is situated in central Gansu Province, China, within the eastern section of the Hexi Corridor. Its geographical coordinates lie between 36°29′ and 39°27′ north latitude, and 101°49′ and 104°16′ east longitude. Positioned at the convergence of three major landform units, the Qinghai-Tibet Plateau, the Mongolian Plateau, and the Loess Plateau, it serves as a pivotal node city along the Silk Road Economic Belt. Administratively, Wuwei City is a prefecture-level administrative unit, overseeing one municipal district (Liangzhou District), two counties (Minqin County and Gulang County), and one autonomous county (Tianzhu Tibetan Autonomous County). Liangzhou District serves as the city’s political, economic, and cultural centre. The city borders Inner Mongolia Autonomous Region to the north, Qinghai Province and Baiyin City of Gansu Province to the south, Zhangye City to the west, and Ningxia Hui Autonomous Region and Lanzhou City of Gansu Province to the east, forming a crucial transitional corridor linking multiple northwestern provinces and regions. The region exhibits a typical temperate continental arid climate, with annual precipitation averaging 150–300 millimetres. Evaporation significantly exceeds rainfall, rendering agricultural production highly dependent on irrigation from Qilian Mountain snowmelt and groundwater resources. Soil types are diverse, yet commonly suffer from limiting factors such as low organic matter content and salinization. Land use within the region is dominated by irrigated agriculture, primarily cultivating maize, wheat, melons, and fruits, making it one of Gansu Province’s key commercial grain production bases. Figure 1 presents the administrative division map of the study area.

2.2. Data Sources

Soil sampling and laboratory analyses were carried out in accordance with the Grading of Cultivated Land Quality (GB/T 33469-2016), the Measures for the Survey, Monitoring and Evaluation of Cultivated Land Quality, and the Soil Testing standard series (NY/T 1121). The study area included cultivated land across four counties/districts of Wuwei City, Gansu Province. A total of 638 sampling sites were established at an approximate density of one site per 666.67 ha. At each site, plough-layer soil was collected using a five-point composite method. Sampling locations were recorded using GPS, and samples were sealed in plastic bags and transported to the laboratory.

In the laboratory, samples were air-dried at room temperature, gently crushed, and visible roots and gravel were removed. The soil was sieved through a 2-mm mesh and stored in sealed bags prior to analysis. Soil organic matter (SOM), total nitrogen (TN), available phosphorus (AP), and available potassium (AK) were determined following the corresponding NY/T 1121 procedures for each indicator. Analyses were conducted from 10 July to 21 September 2022 at the Laboratory of the College of Resources and Environment, Gansu Agricultural University. Instrument calibration was performed using standard solutions, and internal quality control included blank and duplicate samples. Spatial datasets comprised the Wuwei City soil type map (Second National Soil Survey), the 2022 administrative boundary map, and the 2022 land-use status map. All spatial data were managed in the China Geodetic Coordinate System 2000 (CGCS2000).

2.3. Research Methodology

This study integrated field-measured soil nutrient data with thematic spatial layers to analyse nutrient variability and spatial patterns across cultivated land in Wuwei City. Descriptive statistics and Pearson correlation analyses were performed in SPSS. For spatial modelling, nutrient data were examined for distributional skewness and transformed where necessary prior to semivariogram fitting. Experimental semivariograms were fitted in GS+ 9.0 using standard models (e.g., spherical, exponential, and Gaussian), and the best-performing model was selected based on goodness-of-fit criteria (e.g., RSS and R2). Nutrient distribution surfaces were generated in ArcGIS using geostatistical interpolation (ordinary kriging). To evaluate predictive performance, sampling sites were randomly split into a training set (n = 542) for model fitting and a validation set (n = 96) for independent assessment, and prediction accuracy was quantified using error statistics (e.g., ME and RMSE) by comparing observed and predicted values at validation locations.

3. Results

3.1. Basic Statistical Characteristics of Soil Nutrients

3.1.1. Basic Statistical Parameters of Soil Nutrients

To characterize nutrient status across cultivated land in Wuwei City and to support subsequent geostatistical modelling, we summarized soil organic matter (SOM), total nitrogen (TN), available phosphorus (AP), and available potassium (AK) using the minimum, maximum, mean, standard deviation, and coefficient of variation (CV) (Table 1). This approach quantified overall dispersion and enabled direct comparison of nutrient variability prior to semivariogram fitting and kriging. All sampling points were randomly divided into a training set (prediction points) for model construction and a validation set for accuracy assessment. The observed ranges of nutrient concentrations were broadly comparable between the two subsets, indicating that the random split did not create systematic differences in nutrient levels. For the validation points, SOM ranged from 6.40 to 39.60 g kg-1, TN from 0.26 to 1.68 g kg-1, AP from 4.70 to 118.60 mg kg-1, and AK from 16.50 to 804.00 mg kg-1. For the training points, SOM ranged from 1.28 to 78.40 g kg-1, TN from 0.03 to 4.27 g kg-1, AP from 2.90 to 125.20 mg kg-1, and AK from 11.80 to 870.00 mg kg-1. The wider ranges in the training set reflected the larger sample size and inclusion of extreme values, rather than substantive differences in underlying nutrient status.

Mean nutrient concentrations were also similar across the two subsets, supporting their representativeness for modelling and validation. Mean SOM was 17.32 g kg-1 (validation) and 18.85 g kg-1 (training), mean TN was 1.08 g kg-1 and 1.09 g kg-1, mean AP was 44.70 mg kg-1 and 46.85 mg kg-1, and mean AK was 151.35 mg kg-1 and 131.61 mg kg-1. Because subset means were close, subsequent prediction evaluation was unlikely to be confounded by biased sample partitioning, and differences in validation performance can be interpreted as model-related rather than sampling-related. CV values indicated substantial heterogeneity among nutrients across Wuwei’s cultivated land (Table 1). In the validation set, variability ranked as AK > TN > AP > SOM, with CVs of 97.05%, 60.19%, 54.90%, and 44.80%, respectively. In the training set, variability ranked as AK > AP > TN > SOM, with CVs of 97.08%, 80.38%, 59.63%, and 59.42%. Overall, AK showed the highest dispersion in both subsets (CV ≈ 97%), whereas SOM showed the lowest (CV ≈ 45–59%).

Importantly, CV describes overall dispersion but, by itself, does not measure spatial dependence; therefore, CV results were interpreted alongside semivariogram-based metrics (nugget, sill, range). Nevertheless, the relative rankings provide useful agronomic interpretation. The very high variability of AK likely reflected strong field-level management effects (fertilizer type and placement, irrigation scheduling, and residue return), combined with patchy soil texture and salinity conditions that influence potassium fixation and plant availability. In contrast, SOM typically varies more gradually because it integrates longer-term inputs and decomposition processes; thus, its lower dispersion is consistent with SOM functioning as a relatively conservative indicator shaped by cumulative management history rather than short-term fertilizer applications. TN showed intermediate dispersion, which is consistent with tight coupling between SOM and TN via organic matter accumulation and mineralization, while AP variability likely reflected heterogeneous phosphorus fertilizer inputs and strong sorption/precipitation reactions in calcareous or alkaline soils, which can create sharp contrasts in available P over short distances. From a management perspective, these results imply that AK (and, to a lesser extent, AP) may require more spatially differentiated recommendations than SOM and TN, because greater dispersion increases the likelihood of both deficiency and over-application within the same administrative unit.

3.1.2. Analysis of Correlations Among Soil Nutrients

The correlations among soil nutrients reflect the soil formation process and the activity of soil nutrients. Two factors significantly influence soil nutrients: soil organic matter and soil pH [12,13,14]. To provide a more intuitive analysis of the correlations among the four soil nutrients in the study area, SPSS 22.0 was used to conduct correlation analyses of organic matter, TN, AP, and AK. The results, as shown in Figure 2, indicate that organic matter exhibits a positive correlation with TN, AP, and AK. Its correlation with total nitrogen is most significant at the 0.01 level (two-tailed), with a correlation coefficient of 0.954. This was followed by AK at 0.378 and, finally, AP at 0.334. TN showed positive correlations with both AP and AK, with correlation coefficients of 0.368 and 0.395, respectively, both of which were significant. AP also showed a significant positive correlation with AK (r = 0.380). This correlation analysis indicated an interdependence among SOM, TN, AP, and AK. The correlation coefficients between indicators reflect the degree of their synergistic variation. When constructing soil quality evaluation systems, indicators that are strongly correlated can provide statistical justification for their complementary or substitutive use.

The relationships among soil nutrients can reflect shared controls from soil formation, organic matter turnover, irrigation regime, and fertilizer management. In many cultivated soils, SOM functions as a key integrator of fertility because it governs nutrient storage and mineralization; similarly, soil chemical conditions (e.g., pH) can modulate nutrient availability through sorption–precipitation reactions and cation exchange processes [15]. While soil pH was not measured in this study, an analysis focusing on the co-variation among the four analyzed indicators, SOM, TN, AP, and AK, can illuminate the interconnectedness of these nutrient pools. Across the study area, all nutrient pairs showed significant positive associations (p < 0.01), indicating that sites with higher fertility tended to exhibit elevated concentrations across multiple nutrients. SOM and TN were highly correlated (r = 0.954), indicating substantial overlap in the information they captured. This pattern is mechanistically consistent with the fact that a large fraction of soil N is stored in organic matter pools, and both SOM and TN respond to long-term residue inputs, manure application, and microbial mineralization–immobilization processes. In contrast, associations between SOM and AP (r = 0.334) and between SOM and AK (r = 0.378) were positive but modest, suggesting that AP and AK were shaped not only by SOM-related processes but also by management history (fertilizer type, placement, and frequency), irrigation-driven transport, and soil mineralogical constraints (P fixation/precipitation; K fixation on clay minerals). TN also exhibited modest correlations with AP (r = 0.368) and AK (r = 0.395), and AP was positively correlated with AK (r = 0.380), consistent with the co-application of nutrients in fertilizer regimes and the tendency for better-managed fields to simultaneously accumulate multiple fertility attributes.

These correlations represent statistical co-variation rather than causal “synergy.” Building on this, from the perspective of cultivated-land quality assessment, the near-collinearity between SOM and TN implies that including both indicators without adjustment may overweight a single underlying fertility gradient. In contrast, AP and AK provide complementary information that may better reflect management-driven differences and soil chemical constraints. Accordingly, evaluation systems can treat SOM and TN as partially redundant indicators (or use one as the primary metric and the other for verification), while retaining AP and AK to improve sensitivity to spatial differences in nutrient limitation and fertilizer management across Wuwei’s cultivated land.

3.2. Spatial Distribution Characteristics of Soil Nutrients

3.2.1. Normality Test

In geostatistics, semivariograms effectively analyze the structural characteristics and variability of soil nutrient spatial distributions, serving as the foundation for applying kriging interpolation to predict and simulate unknown values [16]. As the semivariogram calculation requires data to conform to a normal or approximately normal distribution (with skewness values between -1 and 1) [17], failure to do so may lead to proportional effects. This can distort the experimental variance function, inflate the base value and block gold value, increase estimation errors, and even obscure its inherent structure [18].

Proportional effects refer to the phenomenon in which, during spatial interpolation (such as Kriging), if the spatial distribution of data exhibits pronounced spatial trends or global structures, the statistical characteristics (particularly the variance function) lack spatial stability. This leads to local predictions being dominated by global trends, causing systematic biases in interpolation results. To mitigate the scaling effect in spatial interpolation, the original data are typically transformed to approximate a normal distribution, thereby satisfying the fundamental assumptions of geostatistical methods. This transformation process mathematically adjusts the frequency distribution pattern to reduce skewness and stabilize variance, thereby enhancing the robustness and accuracy of subsequent spatial prediction models. To further satisfy the fundamental prerequisites for statistical analysis, this study assessed the distribution characteristics of four soil nutrient data sets from Wuwei City. Results from the skewness-kurtosis test conducted in SPSS indicated that the raw data deviated from normality overall. As illustrated in Table 2, after comparing multiple transformation methods including logarithmic (Log), square root (SQRT), and Box-Cox transformations (Table 2), it was ultimately determined that the square root transformation would be applied to organic matter, while the logarithmic transformation would be applied to total nitrogen, available phosphorus, and readily available potassium. Following these transformations, the data distribution patterns improved significantly, becoming more normalized. This provides suitable conditions for subsequent spatial analysis.

The semivariogram is the key geostatistical tool for describing spatial dependence and provides the quantitative basis for kriging interpolation of soil nutrients [19]. In nutrient datasets, however, values often show right-skewed distributions and occasional extremes, particularly for fertilizer-responsive indicators. If these features are not addressed, a small number of high observations can disproportionately influence semivariogram fitting, inflate variance at short lags, and yield unstable estimates of the nugget, sill, and range, thereby reducing the interpretability and reliability of kriging outputs [20]. For this reason, nutrient distributions were evaluated and, where necessary, transformed to improve symmetry and stabilize variance prior to semivariogram modelling.

Soil nutrient concentrations may also exhibit a concentration-dependent variance (commonly referred to as a proportional effect), in which spatial variability increases as mean values increase. This pattern is plausible in cultivated land because nutrient levels are shaped by uneven fertilizer inputs, irrigation intensity, residue return, and soil-texture differences, which can generate high-value “patches” and stronger variance in more intensively managed areas[21]. Variance-stabilizing transformations reduce this effect by compressing the upper tail of the distribution, helping ensure that semivariogram estimation reflects spatial structure rather than being dominated by a few extreme values.

The distributional characteristics of SOM, TN, AP, and AK were assessed using skewness and kurtosis statistics, and the raw data deviated from approximate normality (Table 2; Figure 3). Several transformations were compared, including logarithmic (log), square-root (sqrt), and Box–Cox transformations. Based on the diagnostic statistics and visual inspection of frequency distributions (Table 2; Figure 3), a square-root transformation was applied to SOM, whereas logarithmic transformations were applied to TN, AP, and AK. After transformation, the nutrient distributions became more symmetric and closer to normal, providing a more stable foundation for semivariogram fitting and kriging interpolation [22]. To maintain interpretability, kriged predictions were subsequently back-transformed to the original units for mapping and reporting; the back-transformation approach (and any bias correction for log-transformed variables) should be stated explicitly to ensure reproducibility.

3.2.2. Trend Analysis

Regional soil properties commonly display large-scale spatial gradients (trends) and, in some cases, directional dependence (anisotropy) because they are jointly controlled by soil-forming factors (e.g., parent material and terrain) and spatially structured management (e.g., irrigation networks and fertilizer practices)[23]. If strong trends or anisotropies are ignored, semivariogram fitting and kriging parameters may be biased, leading to misrepresentation of the true spatial structure. Therefore, prior to semivariogram modelling, we evaluated global trends and directional structure using the trend-analysis and anisotropy tools implemented in the ArcGIS Geostatistical Analyst module. Trend analysis summarizes broad spatial patterns by fitting smooth surfaces to the data and visualizing how values change along the east–west and north–south directions [24]. Here, trend analysis was used as an exploratory diagnostic to identify whether nutrients exhibited systematic gradients that should be considered during variogram modelling and interpolation.

Figure 4 shows the trend analysis outputs for soil nutrients across Wuwei’s cultivated land. For SOM, the trend surface showed a clear spatial gradient: lower values in the north and higher values in the south, and higher values in the west than in the east, suggesting a broad-scale fertility pattern superimposed on local variability. AP exhibited a west-high, east-low tendency along the east–west direction, whereas variation along the north–south direction was comparatively weak, indicating a more spatially uniform distribution in that axis. TN showed relatively small-scale variation, but the overall pattern was consistent with SOM, with lower values toward the north/east and higher values toward the south/west, reflecting the close coupling between SOM and TN observed in the correlation analysis. AK showed limited directional variation along the east–west axis but displayed a north-low, south-high gradient along the north–south axis, suggesting that AK may be more strongly influenced by spatially structured management and/or soil–water conditions along that direction. Overall, the trend analysis results demonstrated that nutrients differed in the strength and direction of large-scale gradients, suggesting that spatial structure in Wuwei is shaped by a combination of regional controls and field-scale management. These diagnostics provided a basis for selecting appropriate semivariogram settings and for evaluating whether directional modelling (anisotropy) and/or trend handling were required in subsequent geostatistical interpolation.

3.2.3. Data Coordinate Transformation

For semivariogram modelling in GS+ 9.0, distance calculations must be based on planar (Euclidean) coordinates expressed in metres. Because the sampling locations were initially recorded as geographic coordinates (latitude–longitude) on an ellipsoidal surface, using these coordinates directly would yield spherical distances inconsistent with the Euclidean-distance framework used in standard variogram estimation and kriging parameterization[25]. Therefore, all sampling-point coordinates were projected to a planar coordinate system (CGCS2000) and converted to metric X–Y coordinates prior to semivariogram modelling, ensuring that lag distances, correlation ranges, and any anisotropy parameters were estimated on a consistent spatial scale. Latitude–longitude data for the 638 sampling points were converted using MapGIS, and the resulting projected X–Y coordinates were used throughout variogram fitting and spatial interpolation. Table 3 summarizes representative examples of the coordinate conversion. The complete coordinate dataset for all sampling points is retained in the study database and can be made available in an appendix or supplementary dataset to support reproducibility and secondary analyses.

3.3. Structure of the Semivariance Function

The spatial dependence of soil nutrients was using semivariograms, which describe how similarity between observations changes with separation distance and provide the parameter basis for kriging interpolation. For each nutrient (SOM, TN, AP, and AK), an experimental semivariogram was constructed and fitted with candidate theoretical models (e.g., linear, spherical, exponential, and Gaussian). The “best” model was selected by comparing goodness-of-fit statistics (R2 and RSS) and by ensuring that fitted parameters were geostatistically interpretable (i.e., plausible nugget, sill, and range values). In semivariogram terminology, C0 (nugget) represents microscale variability and measurement error, C0 + C (sill) represents total variance, and a (range) represents the distance over which values remain spatially correlated [26]. The strength of spatial dependence was evaluated using the nugget-to-sill ratio, C0/(C0 + C), which summarises the proportion of variance not captured by spatial structure [27,28]. Ratios < 0.25 are commonly interpreted as strong spatial dependence, 0.25–0.75 as moderate, and > 0.75 as weak [29,30]. This ratio should not be described as “randomness” in a literal sense; rather, high nugget proportions indicate that variability is dominated by short-range processes (including field-level management heterogeneity) or measurement noise at the current sampling scale.

For SOM (Table 4; Figure 5), the exponential model provided the best fit (R2 = 0.974; RSS = 4.262 × 10-4). The nugget-to-sill ratio (0.583) indicated moderate spatial dependence, implying that SOM patterns reflected both broader structural controls (soil-forming factors and regional gradients) and local management influences (e.g., residue return, manure inputs, and cultivation intensity). This is consistent with SOM acting as a relatively conservative fertility indicator that integrates longer-term inputs and decomposition, while still varying locally in response to management history.

For TN (Table 5; Figure 6), the spherical model showed the best fit (R2 = 0.988; RSS = 2.993 × 10-4), with a nugget-to-sill ratio of 0.697, indicating moderate-to-weak spatial dependence. Mechanistically, this is plausible because TN is closely linked to SOM through organic matter pools and mineralization–immobilization processes, yet it can still exhibit pronounced short-range variability due to heterogeneous residue incorporation, fertilizer management, and microsite differences in soil moisture and microbial activity. In practical terms, the higher nugget proportion suggests that TN may be harder to predict precisely at fine scales than SOM under the current sampling density.

For AP (Table 6; Figure 7), the Gaussian model provided the best fit (R2 = 0.992; RSS = 0.112), with a range of approximately 89.37 km and a nugget-to-sill ratio of 0.559 (moderate spatial dependence). The relatively broad range indicates that AP contained a substantial regional-scale component, which is consistent with AP being shaped by systematic differences in soil type and long-term fertilizer history. At the same time, the moderate nugget effect is expected because available P is strongly affected by local chemical constraints (e.g., sorption and precipitation in alkaline/calcareous soils) and spatially uneven fertilizer placement.

For AK (Table 7; Figure 8), the Gaussian model also yielded the best fit (R2 = 0.992; RSS = 2.70), but the nugget-to-sill ratio was 0.734, indicating weak-to-moderate spatial dependence and a strong short-range component. This pattern is mechanistically consistent with AK being highly responsive to localized management (application rate and placement, irrigation redistribution) and to soil mineralogical/texture controls that govern potassium fixation and release. Consequently, although the selected model fit the experimental semivariogram well, the high nugget proportion suggests that fine-scale AK mapping is likely to carry higher local uncertainty than SOM or AP, and interpretation should be supported by independent validation metrics and uncertainty outputs. Overall, the semivariogram results indicated that nutrient spatial structure differed by element: SOM and AP exhibited clearer regional-scale structure (moderate spatial dependence), whereas TN and, especially, AK exhibited greater short-range variability at the present sampling scale. This implies that differentiated management zones are more defensible for nutrients with stronger spatial structure, while nutrients with high nugget proportions may require either denser sampling, additional covariates (e.g., texture/irrigation indicators), or cautious interpretation at fine spatial resolution.

4. Discussion

This study employed GIS technology and geostatistical methods to systematically analyse the spatial distribution characteristics of four nutrients—soil organic matter, total nitrogen, available phosphorus, and readily available potassium—in cultivated soils of Wuwei City, Gansu Province[31]. It revealed their spatial variation patterns and structural features, providing scientific basis for precision management of regional soil fertility. Key findings are as follows:

(1) Soil nutrients exhibit moderate spatial variability, with readily available potassium showing the most pronounced variation: Descriptive statistical analysis of 638 sampling points indicates all four nutrients demonstrate moderate variability intensity (coefficient of variation 44.80–97.08%), with readily available potassium exhibiting the highest coefficient of variation (97.08%). This suggests it is significantly influenced by random factors such as artificial fertilisation and irrigation practices. Although organic matter and total nitrogen exhibited relatively lower coefficients of variation, they still demonstrated pronounced spatial heterogeneity. Nutrient concentrations displayed wide distribution ranges—organic matter fluctuated between 1.28 and 78.40 g kg-1 reflecting the uneven distribution of soil fertility across Wuwei’s cultivated lands.

(2) Nutrients exhibit significant synergistic variation patterns, with organic matter and total nitrogen showing the strongest correlation: Correlation analysis indicates all four nutrients are highly significantly positively correlated (p<0.01), with organic matter and total nitrogen exhibiting a correlation coefficient as high as 0.954. This demonstrates their high consistency in soil formation processes and fertility evolution. Available phosphorus and readily available potassium also exhibited moderate correlations (r=0.334–0.395) with other nutrients, indicating that soil nutrient management must holistically consider synergistic effects between elements to avoid the limitations of single-nutrient regulation.

(3) Semi-variance function modelling revealed that nutrient spatial patterns predominantly exhibited moderate autocorrelation. The GS+9.0 model yielded optimal semi-variance functions: organic matter followed an exponential model (R2 = 0.974), total nitrogen a spherical model (R² = 0.988), while available phosphorus and readily available potassium both conformed to Gaussian models (R2 = 0.992). Blockwise coefficient analysis indicated that organic matter (58.3%) and available phosphorus (55.9%) exhibited moderate spatial correlation, suggesting their variation is jointly influenced by structural factors (e.g., parent material, topography) and random factors (e.g., fertilisation, tillage). Conversely, total nitrogen (69.7%) and available potassium (73.4%) showed weaker spatial correlation, with random factors predominantly driving their variation. Spatial distribution patterns reveal organic matter follows a ‘low in the north, high in the south; high in the west, low in the east’ configuration, while available phosphorus and readily available potassium exhibit west-high-east-low and north-low-south-high trends respectively.

5. Conclusions

This study employed geostatistical methods to quantify the spatial variation patterns of nutrients in Wuwei’s cultivated land, providing precise target zones for differentiated , soil improvement, and enhanced land quality. Future work may integrate remote sensing and machine learning technologies to establish a dynamic monitoring system, further the spatio-temporal allocation strategies for nutrient management.

This study provides a spatially explicit and independently validated assessment of soil nutrient heterogeneity across cultivated land in arid Wuwei, northwestern China. Based on 638 sampling sites, all four key plough-layer nutrients (SOM, TN, AP, AK) exhibited substantial variability, with coefficients of variation ranging from 44.8% to 97.1%. Available potassium showed the greatest dispersion and weakest spatial dependence, whereas soil organic matter and available phosphorus exhibited clearer regional-scale spatial structure. The extremely strong correlation between SOM and TN (r = 0.954) indicates near-collinearity and suggests that these two variables largely reflect the same underlying fertility gradient driven by long-term organic matter accumulation and mineralization processes. In contrast, AP and especially AK displayed greater short-range varonsistiability, cent with stronger sensitivity to localized fertilizer management, irrigation redistribution, and soil-texture constraints. Semivariogram modelling revealed predominantly moderate spatial dependence (nugget-to-sill ratios 0.559–0.734), indicating that nutrient patterns in Wuwei are shaped by both structural controls (soil type, landscape gradients) and field-level management heterogeneity. The high nugget proportion for AK suggests that finer-scale sampling or integration of environmental covariates would be required to improve predictive precision for potassium management zones. Collectively, these findings demonstrate that nutrient-specific spatial structures differ substantially within the same agroecosystem, and that uniform fertilization strategies are unlikely to be agronomically efficient. Spatially differentiated nutrient management, particularly for potassium, is therefore justified in arid irrigated croplands of the Hexi Corridor. Future research should integrate auxiliary predictors (e.g., soil texture, salinity indices, irrigation intensity, remote sensing covariates) and explore regression-kriging or machine learning approaches to enhance predictive accuracy and support dynamic nutrient management frameworks.