Phosphorus flow analysis of different crops in Dongying District, Shandong Province, China, 1995–2016

Study area

Dongying District, located in the northeast of Shandong Province, is a core region of the alluvial delta at the estuary of the Yellow River. The region has a warm temperate continental monsoon climate with an average temperature of 14.3 °C and annual rainfall of 684.7 mm (Boss et al., 2011). Presently, the total area of the city is ~8,243.26 km² including 51.72% as agricultural land (China Statistics Press (2018)). The agricultural soils mainly contain tidal and saline-alkali soils (Ottinger et al., 2013). The soil of the region has high salinity and low nutrients (Guo et al., 2020). Additionally, the planting structure in the Dongying District had significantly changed greatly between 1995 and 2016 (Fig. 1). The traditional crops of the region include the grain (e.g., wheat and maize) and economic (e.g., cotton) crops. Since the year 2000, several newly formed lands have been developed and utilized as agricultural soils and the farmers are encouraged to cultivate economic crops (e.g., cotton) (Gu & Wang, 2008). The crops analyzed in our study, including wheat, maize, rice, soybean, peanut, cotton and fruit-vegetable, are under conventional farming, which mainly relies on chemical intervention to provide crop nutrition.

System definition

We used SFA following a previous work of Wu et al. (2016) with minor modifications to analyze the P flows in seven crops of Dongying District, including wheat, maize, rice, soybean, peanut, cotton, and fruit-vegetable. The frame of P flows is shown in Fig. 2. The studied period covered from 1995 to 2016. For the simplification of statistical results, the years 1995, 1998, 2001, 2004, 2007, 2010, 2013, 2016 were chosen to represent change every three years.

In our investigation, we assumed that the agricultural ecosystem in the Dongying District was in a steady state. Concerning the available data, atmospheric deposition, seeds, chemical fertilizers, pesticides, straw return-to-field, irrigation, resident and livestock excrement were considered as P inputs to the crops. In addition, harvested grains, discarded straw, erosion and runoff were recognized as P outputs. The difference between P inputs and outputs was considered as soil P accumulation.

Quantification

Based on the mass balance principle, the equation “P accumulation = P input − P output” was used to describe P flows in the above crops. Additionally, a dynamic model was developed to quantify annual P flows and distinguish the changes during 22 years. The detailed model equations of the P inputs and outputs of the system are described in Table 1.

In this study, we quantified the PUE, P nutrient load, and the risk index of P load in the seven crops of the Dongying District from 1995 to 2016. Specifically, PUE described as the ratio between the P content from the harvested grain and the total P inputs (Senthilkumar et al., 2012; Wu et al., 2016), was used to assess the P use intensity and efficiency in the agricultural soil systems. The P nutrient load and the risk index of P load can evaluate the healthy status of soil in the agricultural systems (Table 1). Here, P load refers to the difference between P inputs and outputs per area (Dambeniece-Migliniece, Veinbergs & Lagzdins, 2018). The risk index of P load was estimated according to the ratio between P content of fertilizers (chemical fertilizer and excrement) used and the appropriate amount of P fertilizers used in local farmland (Table 1, Eq. (14)) (Oenema et al., 2004).

Structural equation model

Partial least squares path modeling (PLS-PM) is a variance-based structural equation modeling technique that relies on an alternating least squares algorithm (Henseler, 2018). It is widely used in various disciplines as an effective tool to analyze the relationships and influence of different aspects on an observed phenomenon. Recently, the improved PLS-PM method has been used for the analysis of cause-effect relationships in confirmatory and explanatory research (Benitez et al., 2020). In this study, factors with higher P input proportions, i.e., chemical fertilizers, resident and livestock excrement, and straw return-to-field, were selected as block variables. Furthermore, the factors of P output including harvested grains, discarded straw and P erosion and runoff were considered as block variables. Thus, the block variables of the PLS-PM in our study included “fertilizers (residents’ excrement, livestock excrement, and chemical fertilizers)”, “straws return-to-field”, “harvested grains”, “discarded straw” and “P erosion and runoff”.

In the aspect of accuracy of PLS-PM, the loadings of those block variables were selected with the threshold 0.7, because the indicators with loadings >0.7 were regarded as adequate indicators for the corresponding block variables (Wang et al., 2016). In our study, the loadings of the five block variables were >0.7. Moreover, the PLS-PM model was performed in the plspm package (version 0.4.9) in R (version 4.0.3) and was validated by 1,000 bootstraps. The relationships between the five block variables, PUE and P loads were described by path coefficients, which was considered significant when P < 0.05.

Data collection

To quantify the P flows in agricultural soil systems of the Dongying District, the data was collected from the official statistical database, interviews, questionnaires, and published data. The statistical Yearbooks of Dongying District (China Statistics Press, 1995, 1998, 2001, 2004, 2007, 2010, 2013, 2016) were used to obtain the precipitation, crop areas, harvested grains, chemical fertilizers, pesticides used in the field, and the amount and area of straw returning to the field. These data are listed in Table 1 and the File S1. The Table S1 provides the details about data acquisition.

In addition, parameters such as the average wind erosion intensity $(w^{wind}$; 0.05 Kg/Ha; Smil, 2000), P content in the rail fall ( $r^{rain}$; 0.01 mg/L; Smil, 2000), P-containing rate of the compound fertilizer $(r^{com}$; 11.71%; Wu et al., 2016), P-containing rate of the phosphate fertilizer ( $r^{phos}$; 43.66%; Wu et al., 2016), P-containing rate of the pesticide $(r^{pest}$; 3.30%; Yan et al., 2008), P content of residents and livestock excrement ( $w^{resid}$ = $w^{livestock}$; 0.73 kg; Wu et al., 2016), the proportion of the excrement applied to the field ( $r^{resid}$ = $r^{livestock}$; 95%; Wu et al., 2016), and P content of irrigation per sown area of crop ( $w^{irrig};$ 0.45 kg ha⁻¹; Chen, Chen & Sun, 2008) were determined. The P-containing rates of wheat, maize, rice, soybean, peanut, cotton, and fruit-vegetable seeds or grains ( $r_i^{seed}=r_i^{grain})$ were assumed to be 0.5%, 0.4%, 0.4%, 0.6%, 0.5%, 0.78%, and 0.13% (Wu et al., 2016), respectively. The P-containing rates of straw $(r_i^{straw})$ of wheat, maize, rice, soybean, peanut, cotton, fruit-vegetable were 0.08%, 0.15%, 0.13%, 0.20%, 0.16%, 0.15%, and 0%, respectively (Smil, 2000; Wu et al., 2016). The grain to straw ratios ( $r_i^{grain-straw})$ of wheat, maize, rice, soybean, peanut, cotton, and fruit-vegetable were 1.10, 3.00, 1.00, 1.60, 1.50, 3.00, and 0.00, respectively (Smil, 2000; Wu et al., 2016). From 1995 to 1998, straw return-to-field ratios ( $r^{returedstraw})$ of wheat, maize, rice, peanut and cotton were 0.5, which later increased to 0.85 from 2001 onwards. Comparatively, the straw return-to-field ratio ( $r^{returedstraw})$ of soybean was 0.16, 0.16, 0.56, 0.56, 0.56, 0.56, 0.56, and 0.56 in years of 1995, 1998, 2001, 2004, 2007, 2010, 2013, and 2016, respectively (File S1). These data and respective sources are listed in Table S1.

Furthermore, the P loss data was variable from 1995 to 2016. Notably, the P loss is associated with P erosion and runoff and is tough to determine (Jiang et al., 2019b; Schroder et al., 2011). Therefore, this study used the P loss rates based on the published literature and inferred that it changed with the applied fertilizers (Chen, Chen & Sun, 2008; Fan et al., 2009; Senthilkumar et al., 2012; Wu et al., 2016; Yan et al., 1999). The values related with P loss were shown in File S1.

Statistical analysis

In this study, the parameters of P input, P output, PUE, P nutrient load, and the risk index of P load in the crops of the Dongying District were estimated. The model equations of these parameters mentioned above are described in Table 1 and File S1. In addition, the model PLS-PM was performed in plspm package in R (version 4.0.3). The block variables from P inputs were selected based on their P input proportions. Thus, the factors with higher P input proportions, i.e., chemical fertilizers, resident and livestock excrement, and straw return-to-field, were selected as block variables. Furthermore, the correlation between P input from fertilization and PUE of different crops was studied by “cor” function (method = “pearson”) in stats package in R (version 4.0.3).

P inputs

Atmospheric deposition, seeds, chemical fertilizers, pesticides, straw return-to-filed, resident and livestock excrement, and irrigation were considered as P inputs to study the different crops. Of the Dongying District, P input of agricultural soil systems increased steadily from 1.48 × 10⁴ t in 1995 to 2.16 × 10⁴ t in 2007, while it decreased from 1.90 × 10⁴ t in 2010 to 1.78 × 10⁴ t in 2016 (Fig. 3). Specifically, compared with other crops, cotton had the highest P input from 2004 to 2013 (Fig. 3). It rose by 5.36 folds from 2.19 × 10³ t in 1995 to 1.17 × 10⁴ t in 2007 and maintained high levels until 2013. Furthermore, the P input of the soybean decreased sharply from 4.33 × 10³ t in 1998 to 65 t in 2016, which was caused by decreasing soybean cultivation. Only a few farmers were willing to plant soybean in the Dongying District, due to tedious harvesting, storage and post-harvest management of and low economic efficiency (Bern, Hanna & Wilcke, 2008). Since then, the gap between supply and demand of soybean has been expanding, while the supply of maize exceeded demand in the Dongying District. We recognized that the crop structure of the Dongying District is unbalanced making it unfavorable for sustained agricultural development.

Concerning the source of P inputs, the chemical fertilizers contributed the highest from 1995 to 2016, followed by the excrement and straw return-to-field (Fig. 4A). Moreover, the proportion of P input from chemical fertilizers increased from 55.4% in 1995 to 78.1% in 2016, while the excrement proportion decreased from 39.1% to 12.2% in the same period (Fig. 4A).

P outputs

Harvested grains, discarded straw and erosion and runoff were considered as P outputs. Phosphorus (P) amount from harvested grains accounted for 80–90% of the total P output (Fig. 4B). It decreased from 5.59 × 10³ t in 1998 to 4.06 × 10³ t in 2004 and then increased to 5.70 × 10³ t in 2016 (Fig. 5). The trend of change in total P output is consistent with that of harvested grains, which decreased from 6.81 × 10³ t in 1998 to 4.65 × 10³ t in 2004, and then increased to 6.56 × 10³ t in 2016 (Table 2).

Compared with other crops, wheat and maize had the highest P output during the 22 years (Fig. 5). Specifically, the total P output of the wheat had the highest proportion 21–43%, followed by the maize 20–33%. The P output of the wheat first dropped from 2.47 × 10³ t in 1995 to 1.02 × 10³ t in 2004, and then gradually increased from 1.55 × 10³ t in 2007 to 2.81 × 10³ t in 2016. The P output of the maize showed a similar trend first decreasing from 1.29 × 10³ t in 1995 to 1.19 × 10³ t in 2004 and then increasing from 1.70 × 10³ t in 2007 to 2.16 × 10³ t in 2016.

P load and risk index

In this study, we used the P load index to indicate the soil health status. We found that the cotton had the highest P load value ranging from 63.42 kg/ha (in 1998) to 104.62 kg/ha (in 2004), followed by the soybean with a P load value of 49.47 kg/ha (in 1998)–87.23 kg/ha (in 2004) (Fig. 6). Next, we estimated the P load risk index of these crops in the study areas (Table 3), and found that it followed the trend of cotton > soybean > wheat > Fruit-vegetable > maize, rice, and peanut. The cotton was on the top of the P load risk ranks, in the II and III classes during the 22 years study period.

PUE

Next, we examined the PUE of the studied crops. The PUEs of fruits and vegetables, and wheat and maize were higher than that of the other crops (Fig. 4C), while the PUE of cotton was the lowest (~6.9–10.2%). From 1998 to 2004, the PUEs of studied crops showed a slow downward trend 36.9% to 19.6% (Fig. 4C). However, after 2004, it steadily increased to 24.1%, 28.4%, 25.2%, and 32.1% in 2007, 2010, 2013, and 2016, respectively (Fig. 7).

Here, we performed PLS-PM structure equation models to explore the effects of P inputs and outputs on PUE and P load (Fig. 8). In the final PLS-PM model after removing the variables with loading values <0.7, the fertilizer block included the factors of residents’ excrement, livestock excrement, and chemical fertilizers. The harvested grains was the most important factor affecting the PUE of crops (path coefficient = 0.934, P < 0.05), followed by fertilizers (path coefficient = −0.597, P < 0.05). Moreover, the fertilizers positively affected the P load, harvested grains, and P erosion and runoff (P < 0.05). Overall, fertilizer block is the key element in the PLS-PM structure equation model as it not only affects the PUE and P load but also affects the factors such as harvested grains and P erosion and runoff. Therefore, to improve the PUE, the comparatively effective method is to optimize fertilizer usage in agricultural soils whilst maintaining crop yields.

The application of PLS-PM in substance flow analysis

Phosphorus is not only a crop growth-limiting nutrient in estuaries but also an essential element for eutrophication in estuarine delta ecosystems (Koh, 2019; Wu et al., 2019). The P enrichment in sediments from eroded upstream water and upland sources can significantly affect the crop yields and eutrophication risks of coastal areas (Qu et al., 2021).

Notably, the Yellow River carries 1.08 × 10⁹ t of fluvial sediment annually to the Bohai Sea (Milliman & Meade, 1983). The sediment deposition has formed the Yellow River Delta which is recognized as the most active coastal zone with land-ocean interaction (Qu et al., 2021). The Dongying District, the core region of the Yellow River Delta, is located at the intersection of the Shandong Peninsula and the Bohai Sea and has been developed into a major agricultural production base (He et al., 2020a). In order to increase crop yield, P fertilizers were excessively used which increased the P accumulation in soils, seawater, and ultimately led to eutrophication of the Bohai sea (Xu et al., 2020). The typical soil type in the Dongying District is classified as saline alluvial soil (Meng et al., 2020). Notably, excessive application of P fertilizers aggravated soil salinization and ultimately reduced PUE in the Dongying District. The Dongying government firmly implements the guidelines of green agriculture development which encourage the farmers to reduce chemical fertilizers and improve PUE (Xu et al., 2020). To improve PUE, it is necessary to understand its causal relationship with P inputs and outputs. Here, we used the PLS-PM model to quantify the contributions of P inputs and outputs to PUE in the agricultural system of Dongying District.

PLS-PM is a variance-based structural equation modeling technique that relies on an alternating least squares algorithm (Henseler, 2018). It is widely used in various disciplines as an effective tool to analyze the relationships and influence of different aspects on an observed phenomenon. Recently, the PLS-PM method has been further improved for the analysis of cause-effect relationships in confirmatory and explanatory research (Benitez et al., 2020). In our study, the fertilizers, straws return-to-field, harvested grains, discarded straw, and P erosion and runoff were defined as block variables to estimate the relationship with PUE and P load using the improved PLS-PM method. The PLS-PM model was run based on the bootstrap procedure to validate the path coefficients and the coefficients of determination.

Our study is the first to use the PLS-PM model for phosphorus flow analysis. Moreover, we established the numerical relationship between the abovementioned five block variables, PUE, and P load, which is the most significant improvement over other studies. In agreement with previous studies (Wu et al., 2016), we found that fertilizer and crop production are the key factors affecting the PUE of the Yellow River Delta region. We suggest that optimizing fertilizer application rate and increasing harvest can effectively increase the PUE. This is a feasible approach for soil P management. We summarize some recommendations with discussion in “Data uncertainties”.

Comparison of PUE between Yellow River Delta and other deltas

The PUE of crops in the Yellow River Delta was compared with the reported PUEs of other river deltas (Table 4). The crop species, research methods, boundaries and research periods in those studies are listed in Table 4. We found that the PUEs of crops in the Yangtze and Pearl River Deltas were generally lower than that of the Yellow River Delta. This phenomenon can be attributed to P sorption and P forms in soils. Soil inorganic P could account for up to 50–90% of soil total P (Feng et al., 2016). The soil inorganic P includes the compounds Al-P and Fe-P in acidic soils and the compounds Ca-P in calcareous soils (Meena et al., 2018). According to the P transforming dynamics, the soil inorganic P could be fractionated into the available P (Olsen P), moderate-cycling P (variscite (Al-P), strengite (Fe-P), dicalcium phosphate (Ca2-P) and octacalcium phosphate (Ca8-P)), and recalcitrant P fractions (hydroxyapatite (Ca10-P) and occluded P (O-P)). Olsen P, which can be directly absorbed by crops, is closely related to PUE. The Olsen-P content showed great differences among deltas (Li et al., 2015). It was higher in the Yellow River Delta than in the Yangtze and Pearl rivers deltas. Specifically, the Olsen-P content in >90% of agricultural soils in Pearl River Delta is lower than the recommended value (39 mg kg⁻¹) for crop production (Li et al., 2015). In Pearl and Yangtze rivers deltas, the P-deficient land accounts for >50%, compared with 46% in the Yellow River Delta. That explains the reason for lower PUE in Yangtze and Pearl rivers deltas than in the Yellow River Delta.

The distribution and morphology of soil P are the key factors affecting PUE. Although the mineral forms of phosphate are not readily available to crops, Al-P, Fe-P, Ca2-P and Ca8-P can be transformed into free P forms as important buffering pools for Olsen P. Al-P and Fe-P are the main P forms in acidic soil of Pearl River Delta, while Ca2-P and Ca8-P are prominent in calcareous soils of Yellow River Delta (Xu et al., 2019). Notably, dissolution of Ca-P compounds is the main approach for P release under acidic conditions, which is related to ion exchange between OH- and the Fe-P and Al-P compounds (Huang et al., 2021). Furthermore, the Fe-P/Ca-P ratio significantly affects the PUE. Soils with higher Fe-P/Ca-P ratios release more P under alkaline conditions, and in turn higher PUE (Huang et al., 2005). The soil system is a complex dynamic ecosystem (Chen et al., 2002). Thus, to improve PUE, a better understanding of the distribution and morphology of P is important to frame scientific and sensible guidelines.

Data uncertainties

The data and parameters used in the SFA model on P flows were mainly obtained from the statistical yearbooks of the Dongying District, face-to-face interviews, questionnaires, and research articles. The data from the official yearbooks are publicly available and are frequently used in the research of material flow analysis (Han et al., 2021). The local statistician may have changed the statistical methods during the period of 22 years. Consequently, the uncertainty level is greater for the data from the older yearbooks.

The parameters such as the P-containing rates of crop seeds (r^seed ), grain to straw ratio of crops (r^grain-straw ), P-containing rate of crop straws (r^straw ), P content of residents’ excrement (w^resid), P content of livestock excrement (w^livestock), P-containing rate of grains (r^grain), were obtained from previous research articles (Table S1), with the limitation of space and time, which could inevitably introduce errors into our study (Wu et al., 2016).

Moreover, the parameters including average wind erosion intensity per sown area (w^wind ), P content in the rainfall (r^rain), P loss content per unit of crop area (r^loss) are mainly based on the local natural environment and climate. These data are difficult to get and relatively uncertain. Furthermore, the data obtained from the interview and questionnaire, for example, amounts of crop seeds per sown area (w^seed), amount of compound fertilizer used (B^com), amount of phosphate fertilizer used (B^phos), straw return-to-field ratio (r^{retured straw}), proportion of residents’ excrement applied to the field (r^resid), and proportion of livestock excrement applied to the field (r^livestock) might only have reflected the farmland situation of interviewee rather than the whole field information of Dongying District. Although we consider these uncertainties because of data limitations, the uncertainties could not be fully eliminated but were largely minimized by comparing the related studies and broadening the investigation. Additionally, the local monitoring and field experimentation could have improved the parameters and data accuracy (Wu et al., 2016).

Suggestions for improving PUE

We found that the PUEs of studied crops were low (~20–40%) and differed among different crops. The PUE of cotton was the lowest (~6.9–10.2%). Improving the PUE of cotton can improve P management in the agricultural systems of the Yellow River Delta. Additionally, the relationship between P input from fertilization and PUE was studied (Fig. S1). The PUE of wheat, rice and soybean positively correlated with P contents of fertilizers, the Pearson correlation coefficients (r) were 0.047, 0.16, and 0.25 for wheat, rice, and soybean respectively. However, the PUE of maize, peanut, cotton, and fruit-vegetables was negatively related with the P content of fertilizers (r = −0.62, −0.58, −0.29, −0.44 for maize, peanut, cotton, fruit-vegetables, respectively). The farmers of maize, peanut, cotton and fruit-vegetable should make the best use of fertilizers. Furthermore, our research suggests that the low PUE in the Yellow River Delta region is attributed to heavy application of fertilizers and low harvest. Here, we suggest the following measures to improve PUE.