The model simulation was conducted for Quzhou Experimental Station of China Agricultural University (36.9°N, 115.0°E) in Hebei Province, China, covering the period from 2007 to 2017. The experimental site, which located in the core agricultural region of the North China Plain, features calcareous fluvo-aquic soil, and the top 20-cm soil layer exhibited a loamy-silty texture comprising 14.7% clay, 74.0% silt, and 11.3% sand. The chemical properties of the soil at the start of APSIM simulation in 2006 were as follows: soil pH (1:2.5 w/v H₂O) measuring 8.7, extracted mineral N (N _min) 19.8 mg kg⁻¹, labile P 6.6 mg kg⁻¹, and organic matter content 10.7 g kg⁻¹, respectively. Daily climate data from 2007 to 2017 including daily maximum and minimum temperature, rainfall, and sunshine hours were obtained from the weather station of Quzhou County. Daily solar radiation was calculated with daily sunshine hours using the Angstrom equations (Angstrom, 1924).

The management practices implemented in the APSIM model precisely replicated those applied in the field experiments from 2016 to 2017. Both the simulation period (2007–2017 and field experiment (2016–2017) were single-season maize cultivation. The maize cultivar Zhengdan 958 (ZD 958) was sown on May 26 at a final planting density was 75,000 plants ha⁻¹ with 60 cm row spacing of. All treatments were received with 225 kg N ha⁻¹ as urea, 60 kg K₂O ha⁻¹, and P fertilizer as calcium superphosphate. P and potassium were broadcast and incorporated into 0–20 cm soil before sowing. Urea was applied in three applications during whole maize growth stage: 90 kg N ha⁻¹ at sowing, 60 kg N ha⁻¹ at the six-leaf stage, and 75 kg N ha⁻¹ at the 12-leaf stage. The automatic- irrigation switch was turned on to avoid water stress throughout the simulations.

The APSIM model version 7.9 was used to simulate above-ground biomass, grain yield, and P use efficiency of maize at the study site. In APSIM, phenological development of maize from emergency towards maturity is driven by the accumulation of thermal time. Above-ground biomass is simulated with the intercepted radiation.

The APSIM model has demonstrated capability in simulating maize aboveground biomass, crop yield, and nitrogen/phosphorus use efficiencies (Lai et al., 2025; Wang et al., 2014; Liu et al., 2012). Validation studies by Chen et al. (2010), Liu et al. (2012), Wang et al. (2012), Zhang et al. (2012), and Wang (2009) confirmed the model’s robust performance in simulating maize growth dynamics in the North China Plain and Northeast China, particularly regarding dry matter accumulation, yield formation, and nutrient uptake responses to nitrogen/phosphorus supply. Notably, Chen et al. (2010) identified systematic underestimation of aboveground biomass and yield in North China Plain simulations. This limitation was addressed by modifying the radiation use efficiency (RUE) parameter from its default value of 1.6 g/MJ to 1.8 g/MJ, based on methodologies established by Bastiaanssen and Ali (2003) and Tao et al. (2005), resulting in significantly improved biomass and yield predictions.

Within the APSIM framework, maize phenological progression from emergence to physiological maturity is governed by thermal time accumulation. Aboveground biomass production is determined by the interaction between canopy light interception and RUE across developmental stages, while P uptake dynamics are regulated by organ-specific P concentration thresholds during critical growth phases. The APSIM model, incorporating both the maize module and the soil P module (Figure 1), was calibrated using field experiment data (2016–2017) from Quzhou County, Hebei Province, China. Subsequently, the model was validated with data published in the literature, which was also sourced from the same site (Zhang et al., 2018). Comprehensive calibration of these parameters, including thermal time requirements, light interception algorithms, RUE adjustments, and P allocation thresholds of this research program was shown in Table 1, ensuring model fidelity to regional agroecological conditions.

The calibrated APSIM model was employed to evaluate long-term P fertilization management effects on maize growth and soil P pool dynamics in Quzhou County. Scenario analyses of field P management practices focused on crop-soil system responses under varying P application rates. Eight P fertilization rate levels (0–300 kg P₂O₅ ha⁻¹) were simulated: 0 (P0), 25 (P25), 50 (P50), 75 (P75), 100 (P100), 125 (P125), 150 (P150), and 300 (P300) kg P₂O₅ ha⁻¹. Single superphosphate (SSP) was selected as the P source, with full-dose basal application at a 5 cm soil depth. To eliminate water stress effects, automated irrigation was triggered when soil moisture dropped below 85% of field capacity. Model outputs encompassed interannual variations in crop biomass, grain yield, P uptake, and soil P pools (labile P and unavailable P fractions). This configuration enables systematic quantification of legacy P effects, fertilizer utilization efficiency, and environmental risks under continuous P fertilization regimes.

The influence of soil P supply intensity on aboveground biomass and yield of maize exhibited significant interannual variations (Figure 2). From 2007 to 2017, maize yield demonstrated a gradual increasing trend with elevated P supply intensity. However, during 2007–2009, no significant differences in maize biomass or yield were observed among P fertilizer treatments, indicating that soil P supply intensity did not substantially affect maize growth during the initial 3 years of the experiment. From 2010 to 2013, both biomass and yield under the P75 treatment were significantly higher than those under the P0 treatment. Notably, maize yield in 2014 under equivalent P fertilizer treatments significantly higher than that of other years with. Thus, the P75, P100, P125, P150, and P300 treatments resulted in significantly greater maize biomass and yield compared to P0, and P0 showed significant decrease both in shoot biomass and yield when compared to P25 and P50 treatments (Figure 2).

Aboveground biomass, grain yield, and P uptake in maize shoots increased significantly with increased P application rates (Figure 3). The linear-plateau model was used to describe the relationship between P application rates and maize shoot biomass (R² = 0.607), yield (R² = 0.627) and shoot P uptake (R² = 0.860) (Figure 3). The critical P fertilization rate was different among maize shoot biomass, yield and shoot P uptake. For shoot biomass and yield, the critical P fertilization rate was ranged from 64.5 to 71.0 kg ha⁻¹. In contrast, the critical P application rate for shoot P uptake was 193.5 kg ha⁻¹.

The incremental gains in maize biomass and yield across P application treatments exhibited diminishing trends with increasing P inputs (Figures 4A,B), accompanied by a gradually decline in P use efficiency (PUE) (Figure 4C). When P fertilizer application exceeded 75 kg P₂O₅/ha, the biomass and yield increments per kilogram of applied P₂O₅ significantly decreased from 92 kg/kg P₂O₅ and 52 kg/kg P₂O₅ to 25 kg/kg P₂O₅ and 15 kg/kg P₂O₅ at P300, respectively (Figures 4A,B). Meanwhile, PUE declined from 17% at 75 kg P₂O₅/ha to 10% under the 300 kg P₂O₅/ha treatment (Figure 4C).

Under long-term P fertilization, both soil available P (labile_P) and steady-state P concentrations (unavail_P) exhibit progressive accumulation. Figure 5 illustrates the divergent temporal trajectories of labile_P and unavail_P concentrations in calcareous soils under long-term P fertilization regimes during 2007–2017. Over this 11-year period, all P-fertilized treatments except P0 and P25 demonstrated significant increases in soil available P concentrations: P50 (18.1 → 22.2 mg/kg), P75 (20.0 → 28.1 mg/kg), P100 (22.2 → 35.2 mg/kg), P125 (24.6 → 42.8 mg/kg), P150 (26.9 → 50.8 mg/kg), and P300 (43.2 → 121 mg/kg). Conversely, P0 treatment showed a decline in topsoil available P concentration from 14.4 to 11.6 mg/kg. A parallel accumulation pattern was observed for unavailable P concentrations across all fertilized treatments except P0 and P25: P50 (205 → 250 mg/kg), P75 (212 → 316 mg/kg), P100 (218 → 389 mg/kg), P125 (225 → 467 mg/kg), P150 (231 → 549 mg/kg), and P300 (274 → 1,117 mg/kg). The P0 treatment exhibited a marked reduction in topsoil steady-state P concentration from 193 to 131 mg/kg. Notably, Bai et al. (2013) established environmental thresholds for soil available P concentrations across Chinese agricultural soils, ranging from 39.9 mg/kg (Yangling Lou soil) to 90.2 mg/kg (Qiyang Red soil). Given the calcareous soil characteristics of Quzhou, the regional environmental threshold for available P should not exceed 39.9 mg/kg. While the P75 treatment maintained available P concentrations below this critical threshold (32.0 mg/kg in 2017), persistent P accumulation poses ongoing environmental risks. Of particular concern is the P100 treatment, where available P concentrations surpassed the threshold (40.8 mg/kg in 2017). The 23% steady-state P increase observed between P75 and P100 treatments demonstrates a substantial reservoir for available P replenishment, significantly elevating the probability of threshold exceedance. This phenomenon highlights the critical need for optimized P management strategies to mitigate environmental risks while maintaining agricultural productivity in calcareous soil systems.

Figure 6 illustrates the dynamic changes in soil available P and steady-state P pools across various P-fertilization treatments following prolonged application (2007–2017) and a subsequent 22-year cessation period. During the 22-year post-application phase, all P-treated soils exhibited continuous declines in available P concentrations: P0 (11.2 →7.3 mg/kg), P25 (18.1→10.6 mg/kg), P50 (23.3→12.4 mg/kg), P75 (32.3→15.8 mg/kg), P100 (41.2→19.8 mg/kg), P125 (55.6→29.8 mg/kg), P150 (71.3→40.7 mg/kg), and P300 (134.6→102.1 mg/kg). A parallel decreasing trend was observed for unavailable P concentrations across all treatments: P0 (149→74.7 mg/kg), P25 (253 →115 mg/kg), P50 (327→142 mg/kg), P75 (447→195 mg/kg), P100 (567→254 mg/kg), P125 (703→383 mg/kg), P150 (869→516 mg/kg), and P300 (1,504→1,145 mg/kg). These systematic reductions in both available P and steady-state P pools highlight the gradual depletion of legacy P reserves in calcareous soils under extended fertilization discontinuation, emphasizing the critical hysteresis effect between historical P inputs and long-term soil P dynamics. The differential residual P persistence across treatment gradients (e.g., P300 maintaining SSP >1,000 mg/kg post-cessation) underscores the nonlinear relationship between initial P loading intensity and environmental legacy duration.

Figure 7 illustrates the dynamic changes in maize aboveground biomass, grain yield, and P uptake under different long-term P fertilization regimes following P withdrawal. After 22 years of continuous P deprivation, both biomass production and grain yield in P0, P25, and P50 treatments exhibited progressive decline (Figure 7A). In contrast, the P75 maintained biomass and yield levels comparable to those of P100, P125, P150, and P300 treatments, demonstrating sustained high productivity (Figures 6A,B). Notably, long-term fertilization (75 kg P₂O₅ ha⁻¹ yr⁻¹ for 11 years) builds up significant legacy P reserves, sustaining maize yields at 8–10 t ha⁻¹ for 12–13 years after stopping P applications without yield decline (Figure 7B). Shoot P content increase along with the increase of P fertilization rate, and treatments P0-P100 displayed continuous diminishing trends, while P125-P300 treatments maintained elevated P uptake capacities without observable downward trajectories over the experimental duration (Figure 7C). In the maize shoots, P content exhibited a consistent decreasing trend with treatments P0, P25, P50, P75, and P100. In contrast, treatments P125, P150, and P300 maintained elevated P uptake (approximately 40 kg ha⁻¹) without showing a significant downward trend (Figure 7C).

Figure 2 demonstrates that maize biomass and grain yield exhibited no significant declines during the initial three experimental years (2007–2010) under zero or low P supply regimes, indicating that initial soil P concentration sufficiently met crop demands during this period. This observation aligns with Jiao’s (2016) findings that residual soil P can sustain maize productivity for three consecutive years without P fertilization. The residual P pool in agricultural soils has been recognized as a critical potential P resource for sustainable crop production (Sattari et al., 2012). Empirical evidence from multiple long-term studies consistently confirms that legacy P in cultivated soils can effectively satisfy crop P requirements under optimized management conditions (Aulakh et al., 2007; Valkama et al., 2009, 2011). Based on the integrated assessment of aboveground biomass, grain yield, and P uptake, the optimal P fertilizer application rate was determined to be 71 kg P₂O₅/ha (Figure 3). At 71 kg P₂O₅ ha⁻¹ application rate, both maize biomass and grain yield were maintained at highest level with 17.5 t ha⁻¹ and 9.28 t ha⁻¹, respectively. Maize P uptake was only 23 kg ha⁻¹ (Figure 3C). Therefore, 71 kg P₂O₅ ha⁻¹ application rate can help to achieve the dual objectives of ensuring food security and conserving P rock resources.

Following 25 years of continuous P fertilization, both peanut and rapeseed systems demonstrated the capacity to sustain yields for three subsequent years without P application by utilizing residual soil P pools (Aulakh et al., 2007). In high P -fixing soils such as red soils, a single high-dose P application enabled crop production to remain stable for 7–9 years under subsequent P withdrawal regimes (Kamprath, 1967). This phenomenon extends beyond field-scale observations to regional patterns: Japan achieved maintained crop yields from 1985 to 2005 despite progressive reductions in both mineral and organic P fertilizer inputs (Sae and Kohyama, 2010). Similarly, European Union nations observed stabilized or even enhanced agricultural productivity post-1980s alongside declining total P fertilizer usage, attributable to efficient legacy P mobilization (Sattari et al., 2012). These findings collectively demonstrate crops’ ability to exploit residual soil P reserves to buffer yield declines during P fertilizer application reduction.

The APSIM model has incorporated a soil P module (soilP module) and coupled it with crop modules to simulate crop growth responses to soil P availability (Delve et al., 2009; Wang et al., 2014). This study represents the first application of this framework to simulate maize growth responses to P fertilization in fluvo-aquic soils of the North China Plain. Soil organic and inorganic P pools exhibit fundamentally distinct transformation pathways (Hansen et al., 2004; Turner and Leytem, 2004), accounting for 30%–65% and 35%–70% of total soil P respectively (Condron et al., 2005; Shen et al., 2011). The organic P pool predominantly exists in stabilized forms and contributes to available P through mineralization processes (Shen et al., 2011), with mineralization rates being regulated by soil moisture, temperature, and chemical properties. Inorganic P speciation varies significantly between soil types. In acidic soils, P primarily associates with Fe/Al oxides or forms complexes with clay minerals through adsorption processes, while desorption mechanisms can release P into soil solution. Conversely, in neutral/calcareous soils, phosphate ions tend to precipitate on calcium carbonate surfaces. The dissolution of these P precipitates becomes enhanced under decreasing soil pH conditions, thereby increasing labile P availability (Wang and Nancollas, 2008). These transformation pathways-encompassing adsorption-desorption equilibria, precipitation-dissolution reactions, and mineralization processes—constitute a complex dynamic system that requires comprehensive evaluation to advance our understanding of P cycling in agricultural ecosystems.

In the soil P module of the APSIM model, two principal P transformation processes are incorporated: 1) the mineralization/immobilization process between organic P pools and labile P pools, and 2) the sorption/desorption process between stable inorganic P pools and labile P pools (Delve et al., 2009; Wang et al., 2014). These core processes enable the model to accurately simulate P stress conditions in soils of eastern Kenya and the effects of P fertilizer types (chemical or organic) on maize biomass accumulation and grain yield (Kinyangi et al., 2004; Micheni et al., 2004; Probert, 2004). The module has also demonstrated robust performance in simulating crop rotation systems in Australia (Wang et al., 2014) and P dynamics across diverse soil types for maize and soybean cultivation (Delve et al., 2009). Our findings confirm the model’s capability to realistically simulate maize biomass growth and P uptake patterns (Zhang et al., 2024). However, it should be noted that the current APSIM P module calculates labile P content primarily through sorption/desorption mechanisms while neglecting the predominant precipitation/dissolution processes that govern P availability in neutral/calcareous soils. This limitation may constrain the model’s accuracy in simulating labile P dynamics in soils where precipitation-dissolution equilibria prevail. To enhance the model’s performance and broaden its applicability across pedologically diverse systems, further experimental data quantifying P fractionation and transformation kinetics across different soil types are critically required. We employed the APSIM model in conjunction with a scenario analysis approach to systematically evaluate maize productivity, P utilization efficiency (PUE), and the dynamics of soil P pools across different long-term P fertilizer application strategies. This integrated methodology establishes a robust quantitative framework to inform future long-term experimental research.