Soil phosphorus composition, loss risk and contribution to the aquatic environment in a typical agricultural area

Phosphorus (P) is the primary limiting nutrient in most surface water, and a P concentration of over 20–35 μg/L will cause an algal bloom in a freshwater system (Ahmad et al. 2016). Non-point sources P in agricultural runoff contribute to a great portion of freshwater inputs, which accelerate eutrophication and arouse global environmental concern (Sharpley et al. 1994). It was believed that excess fertilization and manure production caused the accumulation of surplus P in soil, some of which is transported to aquatic ecosystems (Carpenter et al. 1998). Soil erosion from agricultural lands, which delivers large amounts of particulate P to surface water, has become the dominant source of P accumulation in the aquatic ecosystem (around 75%; Vilmin et al. 2018). Therefore, the identification of soil P levels, risk of soil P loss and the impact of soil P migration on water quality are required for making sustainable P management to ensure both environmental safety and crop production (Zhou et al. 2004; El-Nahhal et al. 2014a, 2014b).

Although accurate estimations of P leaching have been conducted for some soil types (Sharpley et al. 1994; Bai et al. 2013; Jalali & Jalali 2017; Khan et al. 2018), they are area-specific, time-consuming and can only provide limited information for soil P leaching potential. Therefore, soil P loss evaluation from soil properties at the basin or larger scale is very important for soil P management, utilization and environmental risk control (Amundson et al. 2015). Based on previously published research, soil P threshold levels for agronomic and environmental purposes have been established (Feagley & Lory 2005; Sharpley et al. 2008; NRCS 2011), which can be used for soil P loss risk evaluation in an agricultural area. In this case, a comprehensive soil P investigation at the basin scale was essential to provide sufficient basic information.

High-resolution mapping of soil P content is necessary to identify critical source areas where a large risk of loss coincides in agricultural landscapes. However, because of the high heterogeneity of soil type and property, it is difficult to precisely interpret soil P's content and distribution. Spatial variability analysis based on geostatistics has been widely applied on soil property survey, fertilization evaluation and farmland nutrient management (Ahmad et al. 2016; Denton et al. 2017; Vasu et al. 2017; Laekemariam et al. 2018), which facilitates reliable quantitative assessment of soil P spatial heterogeneity. This may allow a better understanding of soil nutrient evolution and its dynamics, which can further assist sustainable agricultural nutrients management, soil fertility maintenance and eutrophication prevention (Mousavifard et al. 2013; Vasu et al. 2017; de Oliveira et al. 2018; Wenhua et al. 2018).

Jialu River Basin is an agriculture intensive area, the largest tributaries of the Huaihe River Basin, China, with average P fertilizer input (2000–2010) of 135.5 kg/ha yr (China 2010). The annual runoff (13.41 and 26.09 m³/s at midstream and downstream, respectively) (Tianhai et al. 2016) contributes 2.93 t dissolved P and a considerable amount of particulate P to the Jailu River (unpublished data). As a result, P concentration in the river water exceeds by 1–2 times China's environmental quality standard for surface water (Ministry of Environmental Protection of P. R. China 2002), which has greatly impacted the water quality and ecological stability. On the other side, there are 11 main sluices along the Jialu River (The Huaihe River Commission of the Ministry of Water Resources, P. R. China 2010), resulting in low velocity and relatively retained hydrodynamic state, which aggregated the risk of eutrophication. Therefore, it is necessary and urgent to thoroughly investigate the soil P loss risk and its contribution to aquatic P accumulation in this representative area.

The objectives of this research were to evaluate the soil P content and loss potential in the Jialu River Basin, as well as its contribution to aquatic P. To achieve these goals, the soil P composition and spatial variability in the Jialu River Basin were depicted using geostatistics. At the same time, the correlation between soil P and aquatic P was analyzed to interpret the P transport pathway and accumulation pattern at the basin scale. This research could provide a reference for the adjustment of agricultural management practices such as fertilization recommendations and non-point source pollution control in similar agriculture intensive areas.

Jialu River, one of the largest tributaries of the Huaihe River, with a total length of 255.8 km, and a catchment area of 5,896 km². The geographical range of the Jialu River was 113°6 ′8″–114°38′43″ E and 33°38′9″–34°55′55″ N, with a decline of altitude toward the southeast (Figure 1). The average annual precipitation is 633 mm, most of which concentrates in July–September, accounting for about 60% of the annual amount. The average annual evaporation is about 1,700 mm. Jialu River Basin contains rich soil types, including KG soil, neutral calcareous soil, sand, tide soil, wet tide soil, salinization wet soil, most of which are KG soil and calcareous moist soil. Land coverage types in the basin are mainly farmland, forest and grassland, among which farmland occupies 75.1% of the total area. The main crops in the Jialu River Basin include wheat, maize, peanut and so on, with wheat–maize rotation as the main tilling method.

133 soil samples (1–133) and 10 sediment samples (c1–c10) within the basin were collected according to the principles of mesh points in October 2013, and the monthly river water TP data were collected from the Huaihe River Water Resources Commission (Figure 1 and Supplementary Material, Table S1). One kg of surface soil (0–20 cm depth) was collected using a ring knife. The bulk density of the sample at each site was measured immediately. Soil samples were removed from stones, weeds and other debris, stored in plastic bags and sent back to the laboratory. All samples were air-dried, sieved through 100 mesh for later use (Fu et al. 2010). Ten sediment samples (0–10 cm depth) along the Jialu River and its tributaries were collected, air-dried, passed through 100 mesh sieve and stored in plastic bags for later use.

Soil P speciation was done using the SMT (Standard Measurements and Testing) method (Hua 2010) developed under the framework of the European Commission (Figure 2). The amount of dissolved phosphate in the extract (SRP) was determined by the molybdenum antimony anti-spectro-photometric method undertaken on a UV spectrophotometer (UV-2100, Rayleigh, China). The measurement process was as follows. The concentration gradient of 0, 0.5, 1.0, 3.0, 5.0, 10.0 and 15.0 mL were taken for the standard curve. Take 0.04 g of potassium peroxydisulfate in the colorimetric tube, sealing and digesting under 120 °C for 30 min. Inject 1 mL of ascorbic acid solution in the colorimetric tube, mixing for 30 s. Inject 2 mL of molybdate solution and mix for 15 min. Test the absorbance at the wavelength of 650 nm, with a blank sample for contrast. Each sample was measured in triplicate. Soil total organic carbon (TOC) was measured by an elemental analyzer (CHN-O-Rapid, Heraeus, Germany).

The summary statistics of soil P, including minimum, maximum, mean, standard deviation and coefficient of variation were calculated using SPSS 21. Additionally, a geostatistical method was applied on spatial variability analysis of soil P by the geostatistical tool in ArcGIS 9.2 for windows (Guoan & Xin 2006).

Based on the fitted semivariogram model, spatial parameters such as nugget (C₀), sill (C+C₀), range and nugget coefficient (C₀/C+C₀) were calculated. Nugget represents variation caused by stochastic factors, such as an error in measurement. Sill represents the overall sample variability. The range represents the distance at which one variable becomes spatially independent of another (Reza et al. 2015). Nugget coefficient represents the degree of spatial dependence, whereby a value <0.25 indicates variables have a strong spatial autocorrelation, 0.25–0.5 indicates a significant spatial autocorrelation, 0.5–0.75 indicates a moderate spatial autocorrelation and >0.75 indicates a weak spatial autocorrelation, under which condition random variation takes the main role (Cambardella et al. 1994).

In the adopted speciation method in this study, the soil P species were classified into Al/Fe/Mn hydroxide bonded phosphorus (non-apatite phosphorus, NAIP), calcium phosphate (apatite phosphorus, AP), inorganic phosphorus (IP), organic phosphorus (OP) and total phosphorus (TP). The sum of IP and OP was theoretically equal to the concentration of TP, and the sum of AP and NAIP was equal to the concentration of IP. The difference between the analytical values and the calculated values for IP and TP was within the analytical error (3–4%). The soil TP concentration in the Jialu River Basin was 0.56 mg/g soil on average and ranged between 0.14 and 1.418 mg/g. The concentration of IP ranged from 0.12 to 0.96 mg/g soil and was the main component of soil TP. The average concentration of soil OP ranged from less than 0.01 mg/g soil to 1.32 mg/g soil with an average of 0.12 mg/g soil. The coefficient of variation of OP is relatively large (1.98), showing a significant variability. On the contrary, coefficients of variation of IP and TP were all less than 0.4, representing a low variability (Table 1). The average concentration of soil AP and NAIP was 0.26 mg/g soil and 0.23 mg/g soil, respectively, with an intermediate coefficient of variation of 0.96 and 0.82.

In this research, the soil TP was at the same level as the Wenyu River Basin, China (550 mg/kg; Jianling Liu & Zhang 2000), Brittany region, France (296–2,393 mg/kg; Matos-Moreira et al. 2017) and Sub-Saharan Africa (130–4,400 mg/kg; Magnonea et al. 2017), while it appeared lower than that in the sediments of Yellow River Delta (around 600 mg/kg; Liu et al. 2018) and higher than that of Danjiangkou reservoir, China (0.38–176 mg/kg; Li et al. 2020). IP was also the main component of soil TP in other areas (Liu & Zhang 2000).

In Jialu River, sediment TP (n = 10) ranged from 0.271 to 0.867 mg/g (Figure 3(a)). The proportion of IP was between 77.5% and 88.9% (average 82.13%), as the main form of sediment P. The proportion of OP was less than 10%, composing only a small part of the sediment P. The proportion of NAIP varied between 14.6% and 57.1% with an average of 28.7%, indicating an obvious heterogeneity of NAIP content in sediment along the Jialu River.

P content differed significantly in each soil type (Table 2), which was closely related to soil parent material (Mocek & Owczarzak 2011) and human activities (Amundson et al. 2015). NAIP was the main component of IP in meadow aeolian sandy soil with an average of 89.8% (0.336/0.374 mg/g), as well as calcareous cinnamon soil with an average of 88.1% (0.333/0.378 mg/g). Meadow aeolian sandy soil composed around 7.7% of the basin area, which was formed in the ancient riverbed, receiving mainly terrigenous detrital P that led to a higher NAIP proportion, while the overall P loss was severe after long-term weathering and leaching (Squires 2009). It was the other way around in fluvo-aquic cinnamon soil where AP took about 68.7% of the IP on average (0.325/0.473 mg/g). This soil took about 2.3% of the basin area, distributing mainly in the hawthorn complex alluvial plain when coupled with subsequent long-term agricultural activities, it led to high TP content and a high proportion of AP. The proportion of AP and NAIP in mixed soil and cinnamonized fluvo-aquic soil were roughly equal. Mixed soil (including middle mixed soil) was the main soil type in the study area (34.6%) presenting mainly in the slope area of plains. Fostering from the sandy silt of riverbed deposits and the upper humus with well-sorted particle size, mixed soil led to a stable P content.

Pearson Correlation Coefficients analysis (Supplementary Material, Table S2) indicated that soil IP and TP had a significant correlation (r = 0.825, p < 0.01, n = 132), as well as NAIP and TP (r = 0.530, p < 0.01, n = 132) in all soil samples. This result was caused by the high proportion of NAIP and IP in soil TP.

The spatial distribution of soil P showed that AP, IP and TP concentrations gradually increased downstream (Figure 4). At the same time, there was a significant symmetrical pattern cross river, among which IP was the most obvious one. OP showed no discernible pattern with the middle area concentration lower than that of upstream and downstream areas, resembling the NAIP. IP had the lowest concentration in the middle reach of the Jialu River between 0.10 and 1.00 mg/g. AP concentrations ranged from 0.03 to 1.00 mg/g, which was slightly higher than NAIP. Affected by the adsorption/desorption process, soil P content is related to particle size, pH, organic matter content and moisture (Kuo and Lotse 1974; Parfitt 1977). At the same time, NAIP can release bonded P under disturbed redox conditions, which may be transformed to low soluble and inactive AP (Wang et al. 2017). This process may explain the high content of AP in the soil.

In the Jialu River Basin, agriculture is the main land use type (75.1%). According to the 10-year fertilizers application in the Jialu River Basin (2000–2010), the input of agricultural P fertilizer in the high TP area was over 10,000 tons/yr. Compared with other less intensive agricultural exploitation areas (China 2010), the significantly higher application indicated that P fertilizers were the main source of soil TP.

The content of soil OP in the research area was generally low, and there was no obvious distribution pattern, which may be due to the mineralization of biogenic OP in the soil and lack of OP input. TOC and IP were negatively correlated (r = −0.554, p < 0.01, n = 70) (Supplementary Material, Table S3), which indicated that soil organic carbon could influence the soil P form significantly. It was reported that organic acid could increase the available phosphorous in soil and accelerate the accumulation of OP in specific soil type (Zhang et al. 2009), which may explain the negative association of the two parameters. This may also be explained by the simple dilution of IP in the soil by a poor phosphorous organic matter.

The loss of soil P in dissolved and particulate forms is influenced by topography, soil type, soil P content and soil hydrology (McDowell & Sharpley 2001), resulting in high uncertainty. Many studies have been devoted to estimating the P release amount and mechanism from agricultural soil to surface and subsurface runoff (McDowell & Sharpley 2001; Frossard et al. 2014), and it was found that P concentration in runoff was positively related to extractable soil P. Thus, the soil P surplus control was the most efficient way to release soil P loss through runoff. 50 mg/kg Mehlich-3 P for agronomic and 190 mg/kg for environmental threshold has been established in the USA (NRCS 2011). Soil available P above 41.2 mg/kg can pose leaching risk (Bai et al. 2013). Xi et al. reported that the agronomic and environmental P threshold were 15 and 30 mg/kg, respectively, in the semi-humid plain area of Yellow River and Huaihe River Basin (Xi et al. 2016). In this research, based on the average total dissolvable P (TDP) proportion of 69.1% in soil TP (Xi et al. 2016), TDP was from 94 to 622 mg/kg (393 mg/kg in average), which all exceeded (on average 13-fold) the environmental P threshold of 30 mg/kg. TDP in mixed soil (including middle mixed soil) and sandy soil (area ratio of 49.8%) were the highest in all soil types with an average of over 400 mg/kg, indicating a high P loss risk (Table 2). Therefore, stricter soil nutrient management strategies are recommended in the Jialu River Basin to control the loss of soil P.

The contrast between the semivariograms of different soil P types was carried out to analyze the spatial correlation and variability of soil P (Figure 5 and Table 3). The cross-validation results are shown in Table 4.

In this study, the nugget coefficient of IP was larger than 0.75, indicating a weak spatial correlation. The nugget coefficient of TP, AP and NAIP was between 0.5 and 0.75, indicating a moderate spatial correlation. The nugget coefficient of OP was between 0.25 and 0.5, indicating a significant spatial correlation (Table 3). The observed nugget coefficients of all soil P types were classified as moderate or weak, except for OP, indicating a slight spatial dependence and low heterogeneity in the Jialu River Basin. Particularly, soil OP showed a significant spatial correlation, which was probably influenced by agricultural activity. Thus, it can be concluded that the soil P spatial distribution of the study area was the result of a combination of structural and random factors, which is highly complex (Goenster-Jordan et al. 2018). Eljebri reported a moderate to weak spatial variance between soil P₂O₅ content in the irrigated plain of Doukkala, Morocco, which might be related to unbalanced P fertilizer application (Eljebri et al. 2018). The spatial dependence and heterogeneity of soil P in the Jialu River Basin were relatively low compared with the Danjiangkou basin (Li et al. 2020) where soil TP had moderate and high variability in soil with a CV of 40%. Similarly, in the study of spatial variability of the soil pedoindicators with different textures, the spatial dependence of almost all soil properties was classified as moderate to weak (de Oliveira et al. 2018; Goenster-Jordan et al. 2018), which is consistent with this study. It was reported that most soil physical properties, such as particle size distribution, water content at field capacity, permanent wilting point and available water content of alluvial floodplain soil, were also moderately spatial-dependent (Reza et al. 2015). The research in a hot and humid tropical region of India reported a varied spatial distribution pattern with moderate to strong spatial dependence for most of the soil properties, including pH, electrical conductivity, soil organic carbon, available P, K, S and B, exchangeable Ca and Mg, etc. (Behera et al. 2018). In another study, most variances in extractable soil P for the 0–15 cm depth were associated with differences among fields. However, the significant variance was associated with differences among sampling plots within each field (Wilson et al. 2016). It was shown that Particular P in rivers is primarily correlated to suspended solid concentrations, which in turn are correlated to average soil clay content and land use (Sandström et al. 2019). These results indicate that the spatial variability differed between soil parameters and was greatly influenced by the bedrock property, hydrologic background, weather condition and agricultural management.

The spatial variability of soil parameters can be influenced by many factors. It was reported that soil biogeochemical background corresponding to P inherited from natural soils at the conversion to agriculture and farming practices were the main drivers of the spatial variability in cropland soil P content (Ringeval et al. 2017). A structured spatial variability of soil parameters at different scales and magnitudes of strength was found in an Alfisol soil catena (Rosemary et al. 2017), and land use history showed a significant impact on the soil spatial variability. In this research, most of the samples were collected from farmland, so the influence of biogeochemical process (denitrification) and farming practices (wheat and corn rotation) on spatial variability of soil P might outweigh that of natural property. Besides, sampling depth also influenced the spatial variability as reported in other studies. Soil sampled deeper (from 3 to 15 cm) would reduce the variability of the measured P values, because soil P in 0–3 cm included a high level of statistical ‘noise’ (Kaul & Grafton 2017). It was reported that the chemical properties of soil have greater spatial dependence at 0–0.1 m, which are mostly influenced by an intrinsic factor, such as mineralogy and texture, whereas extrinsic variables such as tillage, fertilizer application, soil and water management and other management practices may control the variability of the moderately spatial-dependent properties (Ramzan & Wani 2018). In this research, a mixed soil sample from 0 to 20 cm, which was in the optimal range, was collected for P measurement. In conclusion, the spatial prediction of soil properties using the geostatistical approach is an alternative for the ordinary difference method, which will help in site-specific farming in the study area. In future research, a better sampling strategy should be developed to minimize the influence of soil property uncertainty.

The spatial range for different components of soil P indicated that in a specific area, a greater number of samples are necessary for parameter determination to acquire appropriate statistical precision. In this research, the range of each soil P component was around 1 m, which indicated that observed values of soil P are influenced by environmental parameters over a smaller distance. A much larger range for soil micronutrients (495–2,110 m and 2,200–7,364 m) was reported by Foroughifar et al. (2013) and Dharejo et al. (2011). Generally, soil P can be influenced by many factors such as characteristics of soil mineralogy and weathering history, and particularly by human activities especially intensive agricultural fertilization, while micronutrients are mostly determined by bedrock type, which could be rarely influenced by other factors. Therefore, a precise soil P profile needs a dense sampling network.

In this study, IP was the main P form in sediments (77.0–96.4%), most of which was composited with AP (70.4–95.6%), resulting in a significant correlation (r = 0.910, p < 0.01, n = 7) between them (Table 5). In the meantime, both IP and AP showed a significant correlation with TOC in the river sediments (r = 0.795/0.838, p < 0.05, n = 7). The possible sources of sediment IP include self-generated AP and other bio-particles from the soil. Although the content of NAIP content was low in the research area, it can still be converted into soluble form under reductive conditions and enter the water column, which is also an important endogenous source of P load (Xiaona 2006). This soil P leaching and NAIP dissolving process might explain the dominant P form shift from NAIP to AP and from terrestrial soil to the aquatic environment.

It was worth noticing that sediment AP had a significant positive correlation with soil AP (r = 0.842, p < 0.05, n = 7), while sediment OP was positively related to soil AP (r = 0.0.841, p < 0.05, n = 7) and negatively related to soil NAIP (r = −0.812, p < 0.05, n = 7) (Table 5). AP is mainly formed during the deposition and early diagenesis of calcareous bio-particles in the upper layer of water, as well as the mineral debris produced by the rock weathering in the basin. Therefore, it can be inferred that soil P is the main source of sediment P, mainly in the form of AP. According to the binding state, Fe/Al/Mn hydroxide bonded phosphorous (NAIP), Ca-bound authigenic apatite P and detrital apatite P (AP), and exchangeable phosphorous (Ex-P) is the main form in sediment P, among which Fe/Al/Mn bonded P mainly comes from terrestrial input, reflecting the input intensity of terrestrial P sources (Yu et al. 2011), while in sediment mainly comes from aquatic organisms, reflecting the amount of aquatic biomass in the river ecosystem. These distinct sources might explain the negative correlation between soil NAIP and sediment OP content.

Based on the TP accumulation pattern analysis, TP in upstream and sub-streams sediments were lower than that of the downstream and the mainstream (Figure 6), revealing an accumulation behavior of sediments P along the Jialu River. At the same time, a good consistency of TP in the sediment and water can be seen, which was also subject to soil TP content, and further proved the extensive P input from soil to river water. Ning et al. (2021) used the interactive evaluation index method to evaluate the risk of soil P loss, but this method could not identify the source of P in rivers. It was reported that P-rich soils lost substantial amounts of P stocks as high as 70% under the influence of cultivation (Alvarez et al. 2018), most of which might enter the aquatic ecosystem through runoff. Under this instruction, critical source areas for P loss could be identified. Many states in the USA considered the development of recommendations for P applications and watershed management based on the potential for P loss in agricultural runoff, ranging from 50 mg/kg in Delaware to 200 mg/kg in Texas (Sharpley 2001). The high risk of soil phosphorus loss in the Jialu River Basin indicates the urgency and necessity of applying the best management practices.

IP was the main form of soil P in the Jialu River Basin (82.13%). The soil P composition differed significantly in different soil types. For example, AP was the primary form of IP in fluvo-aquic cinnamon soil (the area ratio of 2.3%), while NAIP dominated in meadow aeolian sandy soil (the area ratio of 7.7%). In contrast, proportions of AP and NAIP in most of the soil types such as mixed soil were roughly equal (the area ratio of 34.6%). Calculated soil TDP varied from 94 to 622 mg/kg (393 mg/kg on average), which all exceeded the environmental threshold of 30 mg/kg. TDP in mixed soil and sandy soil (the area ratio of 49.8%) was the highest (average >400 mg/kg), indicating a high P loss risk. The spatial variability of soil P was moderate or weak, indicating a low heterogeneity in the Jialu River Basin.

In the river sediment, IP (77.0–96.4% of TP) had a significant correlation with TP (r = 0.766, p < 0.01, n = 10). Sediment IP, as well as AP (70.4–95.6% of IP), both showed a significant correlation with sediment TOC (r = 0.795/0.838, p < 0.05, n = 7). At the same time, sediment AP had a significant positive correlation with soil AP (r = 0.842, p < 0.05, n = 7), which confirmed that soil AP as the main source of sediment P was associated with organic matter from weathered soil. Furthermore, an accumulation behavior of sediment P along the Jialu River and a consistent trend with water TP along the flow path was revealed.

This study was funded by The Major Projects of Science and Technology for Water Pollution Control and Management (Grant Nos. 2017ZX07602-003 and 2012ZX07204-003).

All relevant data are included in the paper or its Supplementary Information.