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

Investigating the phosphorus (P) sources, pathways, and final sinks are important to reduce P pollution and improve P management. In this study, substance flow analysis (SFA) was performed for P flow analysis from 1995 to 2016 in different crops of Dongying District, a core region of the alluvial delta at the estuary of the Yellow River. The results showed that P input steadily increased from 1.48 × 104 t in 1995 to 2.16 × 104 t in 2007, and then decreased from 1.90 × 104 t in 2010 to 1.78 × 104 t in 2016. Chemical fertilizers made the highest contribution to P input. The cotton with the highest P load was on the top of P load risk ranks. More importantly, this study applied the Partial Least Squares Path Modeling (PLS-PM) model for P flow analysis and established the numerical relationship between the variables (including fertilizers, straws return-to-field, harvested grains, discarded straw, and P erosion and runoff), P use efficiency (PUE) and P load. The analysis revealed that fertilizer and crop production are the key factors affecting the PUE. Therefore, optimizing the use of P-fertilizer whilst maintaining yields can be an effective strategy to improve the local region PUE.


INTRODUCTION
Phosphorus (P), one of the key crop nutrients, is important for crop growth and reproduction (Smil, 2000). Concerning material flow, P is of great significance for the stable functioning of the agricultural ecosystem (Jiang et al., 2019a;Villalba et al., 2008). However, human activities (e.g., fertilizer utilization) have intensively been affecting P cycling, causing serious environmental problems, such as increasing pressure on P resources. Notably, phosphate is a non-renewable and geographically restricted resource (Jiang et al., 2019a;Simons et al., 2014). Given that agricultural production is highly dependent on P-based fertilizers, phosphate rock has been excessively exploited. Meanwhile, excessive P utilization increases its soil accumulation, and a part of the surplus P leaches into the watery areas causing environmental issues, such as eutrophication (Van Drecht et al., 2009). Recently, China has become the largest producer and consumer of P fertilizer in the world, accounting for 37.5% of global production and 30% of global consumption (Jiang et al., 2019a;Zhang et al., 2008). From the perspective of P resource conservation and environmental protection, it is crucial to better manage P-based fertilizer inputs and P use efficiency (PUE). P use efficiency (PUE), an important indicator of P usage efficiency, denotes the ratio between the P content in the harvested grains and the total P input in the agricultural systems (Wu et al., 2016). Improving PUE can better manage P resources and guide crops transit toward a more efficient and sustainable agricultural model (Zheng et al., 2017). Meanwhile, studies based on PUE can guide policymakers and local farmers to develop eco-agriculture and green agriculture.
To estimate PUE in agricultural soil, it is essential to quantify and trace P cycling in agricultural systems (Wu et al., 2016). There are several quantitative methods for calculating the P amount based on Geographic Information System (GIS) (Eastman et al., 2010) and P loss formulas (Leone et al., 2008). However, those methods only consider the P amount but not the P flows. Therefore, these calculations fail to trace P flows and interpret their relationship with human activities. The mass-balance methods, which can analyze and quantify the P flows, are appropriate for PUE estimation. In particular, substance flow analysis (SFA) systematically connects the P sources, pathways, and final sinks, and therefore can efficiently quantify the P balance within the agricultural system . Substance flow analysis (SFA) has been widely used for PUE estimation and P flow analysis to help policymakers or local farmers to develop the proper P management system. For example, Yuan et al. (2011) used SFA to establish an anthropogenic phosphorus flows model within a socioeconomic system in Chaohu City over 2008 and found that fertilizers utilization in agricultural soil was the most important source of P load on local surface water. Eventually, they suggested limiting fertilizer use for ecological agriculture to reduce the eutrophication of water bodies. In New Zealand,  used SFA to estimate P recovery from pollution sources and suggested measures to improve PUE.
The Yellow River Delta (YRD), one of the fastest-growing deltas in the world, covers an area of 5,400 km 2 , of which 5,200 km 2 area is located in the Dongying District, which is a major agricultural base in Shandong Province, China (He et al., 2020b). The soil in the Dongying District has high salinity and low nutrients (Guo et al., 2020); therefore, excess P fertilizer is applied to achieve optimal crop yields, causing P pollution and low PUE. Additionally, the crop-planting pattern in the Dongying District has changed dramatically in the past 20 years from the traditional crop-planting pattern of "grain crops-economic crops" to "grain crops-economic crops-feed crops-energy crops" (Gu & Wang, 2008). In general, the cultivation of feed and energy crops demands more fertilizer than grain and economic crops (Huan et al., 2020). Thus, the current crop-planting pattern in the Dongying District forces excessive P fertilizer usage, increasing the risk of P pollution in agricultural soils. Therefore, it is urgent to analyze the P inputs, outputs, losses, and accumulations in different crops of the Dongying District and find the strategies to improve the current situation.
Accordingly, our study structured the P flow frame in the regional agricultural system by systematic evaluation of P flows, PUE, and P load risk rank. By SFA approaches (Wu et al., 2016), we quantified the P inputs, outputs, losses, and accumulations in seven crops (wheat, maize, rice, soybean, peanut, cotton, and fruit-vegetable) from 1995 to 2016 (over 22 years) in the Dongying District. Our results can help manage the sustainable agricultural development in the Dongying District.

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 2 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 Figure 2 The diagram of phosphorus (P) flows in crops of Dongying district. The figure showed the frame of P flows in crops of Dongying district. 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 cropping systems. In addition, harvested grains, harvested straw, erosion and runoff were recognized as P outputs. The difference between P inputs and outputs was considered as soil P accumulation.
(12) For Eq. (12) The PUE was described as the ratio between the P content from the harvested grain and the total P inputs in studied system. P grain i : P content of the harvested grains of wheat, maize, rice, soybean, peanut, cotton, and fruit-vegetable, respectively (t); P input i : P input of wheat, maize, rice, soybean, peanut, cotton, and fruit-vegetable, respectively (t).
(13) For Eq. (13) P load is used to assess the health of the soil system, and it refers to the difference between P inputs and P outputs per area. P input i : P input of wheat, maize, rice, soybean, peanut, cotton, and fruit-vegetable, respectively (t); P output i : P output of wheat, maize, rice, soybean, peanut, cotton, and fruit-vegetable, respectively (t). A area i : Sown areas of wheat, maize, rice, soybean, peanut, cotton, and fruit-vegetable, respectively (ha); (14) For Eq. (14) P chem i : P content of chemical fertilizers used in wheat, maize, rice, soybean, peanut, cotton, fruit-vegetable, respectively (t); P livestock i : P content of the livestock excrement applied to field of wheat, maize, rice, soybean, peanut, cotton, fruit-vegetable, respectively (t); P resid i : P content of the residents' excrement applied to field of wheat, maize, rice, soybean, peanut, cotton, fruit-vegetable, respectively (t); Fp is the appropriate amount of P fertilizers (70 kg/ha); Lp is the maximum P content in excrement used (35 kg/ha). 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 . 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.
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 4 t in 1995 to 2.16 × 10 4 t in 2007, while it decreased from 1.90 × 10 4 t in 2010 to 1.78 × 10 4 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 3 t in 1995 to 1.17 × 10 4 t in 2007 and maintained high levels until 2013. Furthermore, the P input of the soybean decreased sharply from 4.33 × 10 3 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

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 3 t in 1998 to 4.06 × 10 3 t in 2004 and then increased to 5.70 × 10 3 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 3 t in 1998 to 4.65 × 10 3 t in 2004, and then increased to 6.56 × 10 3 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 3 t in 1995 to 1.02 × 10 3 t in 2004, and then gradually increased from 1.55 × 10 3 t in 2007 to 2.81 × 10 3 t in 2016. The P output of the maize showed a similar trend first decreasing from 1.29 × 10 3 t in 1995 to 1.19 × 10 3 t in 2004 and then increasing from 1.70 × 10 3 t in 2007 to 2.16 × 10 3 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.
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 9 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 . 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 −1 ) 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 . 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 returnto-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.

Optimizing fertilizer application
(1) Soil testing and fertilizer recommendation Soil testing can determine the soil's fertility status to regulate fertilizer demand and improve fertilizer efficiency (Liu et al., 2017;Sun et al., 2013). Soil testing and fertilizer recommendation were implemented in the Dongying District from 2004. From then on, fertilizer inputs reduced by 25.83, 83.33, and 19.66 kg/ha for wheat, maize, and cotton fields in 2013, respectively. However, the soil testing and fertilizer recommendation can be significantly affected by several factors, such as local soil condition, climate situation, and cotton varieties (Jordan-Meille et al., 2012). Moreover, the structure of fertilizer can negatively impact the effect of soil testing and fertilizer recommendation. For example, the local farmers chose chemical fertilizers over organic fertilizers to save manpower.
The excessive use of chemical fertilizers decreases soil fertility and in turn PUE. Therefore, we suggest that the soil properties (pH, clay content), climate situation, cotton species, and fertilizer structure should be considered to attain the optimum results (Jordan-Meille et al., 2012).
(2) Modifying the P-fertilizer treatment In the Dongying District, P is applied by the broadcast fertilization method. However, the method is not effective for P uptake by crops (Schroder et al., 2011). Thus, to increase PUE, novel scientific fertilization modes should be encouraged in the Dongying District. A study showed that lowland rice soil of Sub-Saharan Africa needed twice the amount of P fertilizer if broadcasted instead of P-dipping (Rakotoarisoa, Tsujimoto & Oo, 2020). Here, we recommend the root-zone P fertilization instead of broadcast fertilization to improve P fertilizer efficiency in the Dongying District.
(3) Combination of inorganic fertilizer and organic fertilizer In the Dongying District, the amounts of organic fertilizers (for example, excrement) reduced gradually, while the amount of chemical fertilizer increased rapidly from 1995 to 2016 (Fig. 4A). The organic fertilizers, such as excrement, increase carbon and organic matter in soil that can improve soil P binding, retention, and availability to crops (Yang, Chen & Yang, 2019). However, long-term application of organic fertilizer would promote the soil accumulation of heavy metals, e.g., Zn, Cd, and Cr, inducing adverse effects on food safety (Ning et al., 2017). Thus, we suggest that the farmers of the Dongying District should have the reasonable application of organic fertilizers.
(4) Optimizing the planting structure We found that the cotton and soybean had low PUE and high P load. To improve PUE, the growing area of the cotton and soybean should be reduced appropriately, whilst the growing area of the crops with high PUE and low P load should be increased. However, the local agricultural structure is mainly dependent on the market economy and local planting habits, and therefore cannot be fully controlled.

Increasing the harvest
(1) Utilizing P solubilizing microorganisms (PSMs) Many studies have reported the solubilization of insoluble P by PSM (Richardson, 2001;Sharma et al., 2013), which can convert unavailable P into available P (Olsen P) improving P uptake by crop and higher yield. Generally, PSMs can be divided into two classes; (1) inorganic P-solubilizing microorganisms secreting organic acid to dissolve inorganic P compounds, and (2) organic P-mineralizing microorganisms secreting phosphatase to enzymatically mineralize organic P compounds (Alori, Glick & Babalola, 2017). Both field and pot experiments with or without P-fertilizers showed that PSMs can improve crop yield and P uptake (Bolo et al., 2021;Jiang et al., 2021). For example, Noor et al. (2017) used the co-application of P fertilizer-PSM (Pseudomonas putida) to improve corn dry matter yield and P uptake by 12% and 33%, respectively.
Additionally, PSMs decrease soil pH and are better suitable for natural and alkaline soils. Therefore, PSMs can greatly improve available P and PUE without influencing the soil health in the Yellow River Delta region. Various soil species of PSMs (including both fungi and bacteria) can solubilize phosphorus. However, it is difficult to predict effective PSMs in each field. Therefore, screening effective PSMs is necessary for a specific location.
(2) Improving the mechanization level The crop harvest is significantly affected by the mechanization level. Specifically, in cotton planting, the increase in labor cost and lower production restricts productivity which can be improved with mechanization. The farmers in the Shawan county of Xinjiang Province have realized the cotton mechanization from sowing to harvesting. The local government and farmers introduced a new packing cotton picker, and the process of cotton harvesting and packaging has become highly mechanized. This greatly improved the cotton yield and quality. We suggest that the farmers in the Dongying District should establish similar cotton production mechanization, and learn from the cotton planting experience of Xinjiang Province.
The abovementioned strategies can provide the possibilities to reduce P load and improve PUE, however, with farmer's requirements and acceptance. The ecological agriculture development in the Dongying District is largely dependent on the collaborative partnerships of policymakers, researchers, and farmers. Additionally, the factors of government policy, price, cropping mechanization, and farmer's motivation would indirectly affect PUE by influencing the cultivation area. For example, the government policy of encouraging local farmers to extend cotton plantation since 2004, the supply of maize exceeds demand in Dongying District. However, only a few farmers were willing to plant soybean, because of tedious harvesting, storage and post-harvest management, as well as low economic efficiency (Bern, Hanna & Wilcke, 2008). Furthermore, the low prices and cropping mechanization muted farmers' enthusiasm to grow wheat and corn from the year of 1995 to 2004, reducing the total P output of these crops. Therefore, the researchers should fully consider those factors, e.g. government policy, price, cropping mechanization, and farmer's motivation, in further study of how to improve PUE.

CONCLUSIONS
In this paper, we used SFA to establish the dynamic model of phosphorus flow in different crops of the Dongying District from 1995 to 2016. We found that P input increased steadily from 1995 to 2007, and then decreased from 2010 to 2016. Specifically, compared with other crops, cotton had the highest P input from 2004 to 2013. Among the contributing sources of P input, chemical fertilizers contributed the highest. Meanwhile, the P amount from harvested grains accounted for 80-90% of the total P output. Also, the cotton had the highest P load and topped in the P load risk ranking. Excessive P load in cotton caused problems with the plant uptake of elements and potentially increased eutrophication risk. Additionally, the PUE was significantly different among distinct crops. The PUE of cotton was the lowest. Based on the PLS-PM structure equation model, fertilizers and harvested grains were the two important factors affecting PUE. Therefore, reducing the use of fertilizer in agricultural soils whilst maintaining crop yields can effectively improve PUE.
Finally, we provide some recommendations for improving PUE and reducing P load in agricultural soil. However, the soil physicochemical properties and biases of cotton planting in the Dongying District should be considered before implementing recommendations. Additionally, farmers' requirements and acceptance should be respected. In-field demonstrations and formulating the fertilization scheme can significantly improve policy implementation.

ADDITIONAL INFORMATION AND DECLARATIONS Funding
This study was supported by the National Key Research and Development Program of China (2018YFD0800303), the National Natural Science Foundation of China (41977144 and 31970545), and the Postdoctoral Science Foundation of Qingdao City (61460079311133). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.