The experimental site was located at the National Monitoring Base for Purple Soil Fertility and Fertilizer Efficiency, which is in Beibei District (29°48′38″N, 106°24′40.6″E), Chongqing Municipality, China (Figure 1A). This field is in the typical subtropical monsoon climate zone with a mean annual temperature of 18.3 °C and a mean annual precipitation of 1105.4 mm. The rainy season is mainly from May to October, overlapping with the hot season. This region is rich in low hilly landforms and sloping farmlands. The soil is mainly purple, classified as Pup-Orthic Entisols in Chinese Soil Taxonomy or Regosols in Food and Agricultural Organization (FAO) Taxonomy or Entisols in United States Department of Agriculture (USDA) Taxonomy. The soils are shallow (approximately 30 cm–80 cm depth), and the parent bedrock has poor water permeability. The soil for this study is developed from the Jurassic Shaximiao Formation. The basic physicochemical properties of each layer’s soil are shown in Table 1. The common crop rotation in this region is summer maize rotated with winter mustard.

Each runoff field was constructed as a separate plot with dimensions of 9 m (length) × 3 m (width) × 0.8 m (depth). Each runoff plot’s slope and soil depth were 15° and 60 cm, respectively. Adjacent ones were hydrologically isolated with cement walls to avoid unexpected interplot cross-flow between plots. A reinforced concrete floor was made to imitate the parent bedrock with poor water permeability. The test soil in each runoff field was filled in different soil layers (L₁: 0 cm–20 cm, L₂: 20 cm–40 cm, and L₃: 40 cm–60 cm) according to the technical code for the newly cultivated land in Chongqing. The open trench collected surface (L₀) flow and subsurface flow from three different soil layers were collected via U-shaped PVC conflux troughs settled on the bottom of each soil layer and conveyed to buckets by using water discharge tubes (Figure 1B). Each conflux trough (rectangle without cover, 10 cm width, 5 cm depth) was of the same length, parallel to the experimental field, and then, filled with clean quartz sand and covered by two layers of 200-mesh gauze to enhance the trough’s hydraulic conductivity (Wang et al., 2021a). Runoff with sediment samples could be obtained in each corresponding bucket. The possible overflow was collected into a sub-bucket by using a pipe (nominal diameter 32 mm) to measure the overflow from the bucket to obtain the accurate runoff flow of each plot. Specifically, to improve the test accuracy, two identical collection troughs for each soil layer of a single runoff plot and three runoff plots with the same treatment (as three replicates) were conducted. Conventional along-slope cropping and fertilization were conducted according to farmers’ planting habits. The fertilization amounts and dates for different crops during the test period (3/1/2018-2/29/2020) are shown in Table 2.

The observation of P migration from different soil layers by runoff and sediment was from March 2018 to February 2020. Hourly meteorological data were obtained from a small automatic-collection weather station near the research plots. Soil samples were taken with auger and ring samples before planting and after harvesting. The basic physicochemical properties of the soil and sediment samples were determined by following Yang et al.’s (2008) methods.

After each rainfall-runoff event, the volumes of surface flow and subsurface flow from each soil layer in each runoff plot were monitored. Subsequently, water samples of surface and subsurface flow in the buckets were collected separately and then analyzed promptly in the laboratory. Water samples from surface and subsurface flows were used to measure sediment yield after oven drying at 105 °C for 24 h and the concentrations of P fractions. The total phosphorus (TP) concentration in the water samples was measured without filtration. Water samples used to analyze the phosphate–phosphorus (PO₄ ^3-–P) concentration and total dissolved phosphorus (TDP) were filtered through 0.45 μm filtration membranes. After digestion with potassium persulfate, TP and TDP concentrations were determined through the molybdenum blue colorimetry method (Huang et al., 2020). The concentration of PO₄ ^3-–P was measured using the same method used for TDP without digestion. The PP concentration was calculated by the following equation: PP = TP−TDP (Li et al., 2017). The TP content of each soil layer was determined by the molybdenum blue colorimetry method after melting with NaOH, while the available phosphorus (AP) was measured by the Olsen method (Wang, 2021). The sediment samples were dried in an oven at 105 °C to a constant weight. The conventional pipette method determined the sediment particle size distribution (PSD) (Zuo et al., 2020).

Each fluid sample’s runoff volume was the water collected in the corresponding bucket of each collection device. The runoff flux of surface runoff and subsurface runoff under the ith rainfall event was calculated by the following formulas, modified from the work of Peng et al. (2016): R s u r   i = r s u r   i A ∗ COS 15 ° , (1) R n i = r n i a ∗ COS 15 ° , (2) ∆ R n   i = R n − 1   i − R n   i , (3) where R s u r   i is the surface runoff flux under the ith rainfall event (mm); r s u r   i is the collected surface runoff volume (L); A is the area of a single test plot (m²); 15° is the slope of the test plot; R n i is the subsurface runoff flux at the nth soil layer (mm); r n i is the collected subsurface runoff volume at the nth soil layer (L); n is the number of soil layers (L₁: 0 cm–20 cm depth; L₂: 20 cm–40 cm depth; and L₃: 40 cm–60 cm depth); a is the area of a single collection device (m²); and ∆ R n − 1   i is the retention amount of subsurface runoff at the layer of n-1 (mm). When n = 1, R 0   i is the whole permeation flux of rainfall through the soil surface, which is approximately equal to the gap of precipitation and surface runoff; when n = 4, R 4   i is the final subsurface runoff flux lost from the whole soil layers.

The following formulas were used to calculate the P migration load of the surface and subsurface under the ith rainfall event: Q s u r   i j = r s u r   i * C s u r   i j A ∗ COS 15 ° ∗ 10 4 / 10 3 , (4) Q n i j = ( r n i * C n i j a * COS 15 ° * 10 4 / 10 3 , (5) ∆ Q n   i j = Q n − 1   i j − Q n i j , (6) where Q s u r   i j is the jth nutrient migration load in surface runoff under the ith rainfall event (g/ha); C s u r   i j is the concentration of the jth nutrient in surface runoff (mg/L); Q n i j is the jth nutrient migration load (g/ha); C n i j is the concentration of the jth nutrient at the nth soil layer (mg/L); and ∆ Q n   i j is the retention load of the jth nutrient at layer n-1 (g/ha). When ∆ Q n   i j is positive, it indicates that this layer intercepts N under the rainfall event; otherwise, it plays a releasing role in P migration. When n = 1, Q 0   i j represents the overall leaching amount of the jth nutrient through the soil surface, approximately equal to the jth nutrient gap in rainfall and surface runoff. When n = 4, ∆ R 3   i is the jth nutrient final loss load through subsurface runoff from the whole soil layer.

This study analyzed the migration loads of P under different planting stages (S₁: maize season, S₂: summer fallow period, S₃: mustard season, and S₄: winter fallow period) and different rainfall patterns (TR: torrential rain, HR: heavy rain, MR: moderate rain, and LR: light rain, classified by the Classification Standard of China Meteorological Administration), as shown in Figure 2.

IBM SPSS Statistics 26 was the data analysis tool utilized in this study, and GraphPad Prism 8 was the mapping software used in this paper. Pearson correlation was used for the correlation analysis. The least significant difference (LSD) method at p = 0.05 was performed to elucidate any significant differences between different levels.

There were differences in the annual cumulative rainfall and each soil layer’s annual cumulative runoff loss between different rainfall patterns in the test period. Figure 3A shows that the annual cumulative rainfall amounts of TR, HR, MR, and LR were 137.9, 246.1, 338.9, and 286.3 mm, accounting for 13.7%, 24.4%, 33.6%, and 28.4% of the annual cumulative rainfall, respectively. However, the differences in annual cumulative rainfall amounts between different rainfall patterns were nonsignificant. Figure 3B suggests that, except for the annual loss of surface runoff, the annual cumulative subsurface runoff loss amount measured from each soil layer decreased with increasing soil depth. The annual surface runoff losses under TR, HR, MR, and LR were 19.4, 26.0, 22.3, and 2.5 mm, respectively. The corresponding rainfall-surface runoff rates were 14.1%, 10.6%, 6.6%, and 0.9%, respectively, indicating that the lower the rainfall grade was, the more difficult it was to produce surface runoff. The annual cumulative amounts of intercepted subsurface flow were 63.2 mm (TR), 144.4 mm (HR), 220.6 mm (MR), and 270.1 mm (LR), accounting for 45.8%, 58.7%, 65.1%, and 94.3% of the annual cumulative rainfall, respectively. The lower the rainfall grade was, the higher the proportion of the intercepted flow of supplying soil water to the total rainfall. The subsurface runoff loss amounts from the L₁ layer under different rainfall patterns were 102.7 mm (TR), 166.7 mm (HR), 165.1 mm (MR), and 28.4 mm (LR), which meant that the intercepted amounts of the L₁ layer were 25.0%, 37.0%, 68.7%, and 94.6% of the total intercepted amounts under different rainfall patterns, respectively. The subsurface flow flux decreased significantly with increasing soil depth. In contrast, the intercepted subsurface runoff decreased, which meant that the interception ability of the soil layer for permeated subsurface runoff weakened with increasing soil depth. The eventual annual loss amounts of subsurface runoff under various rainfall patterns were only 46.7% (TR), 34.4% (HR), 30.3% (MR), and 4.8% (LR) of the annual total amounts permeated into the soil and only 40.1% (TR), 30.8% (HR), 28.4% (MR), and 4.8% (LR) of the annual cumulative rainfall.

The distributions of the rainfall and runoff loss amounts of each planting stage were uneven during the whole test period. Figure 3C shows that 56.0% and 29.1% of the annual rainfall amounts were from the maize season (S₁: 571.3 mm) and summer fallow period (S₂: 281.5 mm), respectively; only 12.9% and 2.0% were from the mustard season (S₃) and winter fallow period (S₄), respectively. Figure 3D shows that 67.7% of the annual surface runoff loss amount and 79.1% of the annual subsurface loss amount occurred in the maize season; 25.9% of the annual surface runoff loss amount and 12.0% of the annual subsurface loss amount occurred in the summer fallow period. No surface or subsurface runoff loss phenomenon was found due to the slightly scattered rainfall in the winter fallow period. The annual surface runoff loss and each soil layer’s annual subsurface runoff loss amount decreased from S₁ to S₄ seasons. The annual subsurface runoff-permeated loss amount in each planting season significantly decreased with the increase in soil depth (p < 0.05).

The difference between the average migration loss intensities of P per rainfall-runoff event under different rainfall patterns was significant. As shown in Figure 4, each soil layer’s average P migration loss intensity per rainfall-runoff event increased remarkably with increasing rainfall grades. Notably, the average surface migration loss intensities of TP and PP measured under the TR pattern were exceptionally higher than those measured under other rainfall patterns. The average leached loss intensities of each form’s P decreased sharply from soil depths of 20 cm–40 cm under each rainfall pattern. The average leached loss intensities of TP and PP at the 40 cm depth were significantly lower than those at the 20 cm depth, indicating that the 20 cm–40 cm soil layer intercepted most of the leached TP and PP from the 0 cm–20 cm soil layer under all rainfall events. The rainfall patterns significantly affected the average migration loss intensities of TP and PP in the surface runoff. The change characteristics of the average migration loss intensities of TDP and PO₄ ^3-–P per rainfall-runoff event under each rainfall pattern showed the same trends. In summary, the higher the rainfall intensity was, the greater the TP and PP migrating loads through runoff were.

Figure 5 shows that each form of P’s leached loss intensities from each soil layer presented positive linear relationships with rainfall. The leached loss intensities of each form of P from each soil layer increased significantly with rainfall. The regression coefficient between each soil layer’s leached loss intensities of various P and rainfall decreased with increasing soil depth, which meant that the effect of rain on each leached P loss intensity decreased with the increase in soil depth. Notably, the gaps of each form’s P leached loss intensities between the L₁ layer and L₂ layer were significantly higher than those between the L₂ layer and L₃ layer, indicating that the interception ability of the L₂ layer was significantly higher than that of the L₃ layer.

As shown in Table 3, the average annual migration loss loads of each form’s P from different soil layers under different rainfall patterns were significantly different. The TR, HR, MR, and LR patterns caused 1041.09 g*ha⁻¹*yr⁻¹, 1884.35 g*ha⁻¹*yr⁻¹, 1551.64 g*ha⁻¹*yr⁻¹, and 153.99 g*ha⁻¹*yr⁻¹ of TP yearly migration loss from the L₁ layer, 29.0% (302.27 mm), 25.2% (474.67 mm), 18.4% (285.11 mm), and 12.3% (18.95 mm) of which was lost through surface runoff, respectively. It was concluded that the higher the rainfall grade was, the higher the proportion of the annual surface loss to the L₁ layer’s total annual migration loss was. Additionally, 71.9% of the TP’s annual surface loss occurred under the TR and HR patterns. The annual leached loss loads of TP from L₁ and L₂ under TR, HR, MR, and LR were 738.82 g*ha⁻¹*yr⁻¹, 1409.68 g*ha⁻¹*yr⁻¹, 1266.53 g*ha⁻¹*yr⁻¹, and 135.04 g*ha⁻¹*yr⁻¹ and 259.05 g*ha⁻¹*yr⁻¹, 295.52 g*ha⁻¹*yr⁻¹, 332.70 g*ha⁻¹*yr⁻¹, and 30.15 g*ha⁻¹*yr⁻¹, respectively, indicating that the L₂ layer intercepted 64.9%, 79.0%, 73.7%, and 77.7% of the L₁ layer’s leached TP. The annual leached loss loads of TP from the L₃ layer under the TR, HR, MR, and LR grades were 163.96 g*ha⁻¹*yr⁻¹, 195.53 g*ha⁻¹*yr⁻¹, 165.07 g*ha⁻¹*yr⁻¹, and 17.36 g*ha⁻¹*yr⁻¹, respectively, indicating that the L₃ layer intercepted 12.9%, 7.1%, 13.2%, and 9.5% of the L₁ layer’s leached TP, respectively.

The annual surface migration loss loads of PP under different rainfall patterns were 91.8% (TR), 94.6% (HR), 88.7% (MR), and 83.6% (LR) of the surface migration loss loads of TP. The annual leached loss loads of PP underground under different rainfall patterns were 70.4% (TR), 83.4% (HR), 74.1% (MR), and 58.2% (LR) of the annual leached loss loads of TP underground. According to one-way ANOVA, except for TDP, the annual migration loss loads of each form’s P under different rainfall grades showed significant differences. The annual surface migration loss loads of each form’s P and the annual leached loss loads of each form’s P at each soil layer under the TR, HR, and MR grades were significantly higher than those under LR. Notably, TR and HR were the major rainfall patterns that caused P surface migration losses; however, HR and MR were the primary rainfall grades that induced underground leached P losses.

The annual migration loss loads of different P forms from the same soil layer in different planting stages were significantly different. The annual migration loss loads of each form’s P in different planting stages mostly occurred in the S₁ and S₂ stages. There was no observed migration loss of each form’s P from every soil layer in the S₄ stage. Table 4 shows that the annual migration loss loads of surface TP in different planting stages were 850.03 g*ha⁻¹*yr⁻¹ (S₁), 193.39 g*ha⁻¹*yr⁻¹ (S₂), 37.58 g*ha⁻¹*yr⁻¹ (S₃), and 0.00 g*ha⁻¹*yr⁻¹ (S₄), accounting for 78.6%, 17.9%, 3.5%, and 0.0% of the annual total surface migration loss load of TP, respectively. There was no observed migration loss of each form’s P from every soil layer in the S₄ stage. The annual leached loss loads of TP from different soil layers were 3550.07 g*ha⁻¹*yr⁻¹ (L₁), 917.43 g*ha⁻¹*yr⁻¹ (L₂), and 541.91 g*ha⁻¹*yr⁻¹ (L₃). The sequence of contribution rates of each soil layer’s TP leached loss loads in different planting stages to the corresponding soil layer’s annual total leached loss loads were 60.7% (S₁), 31.9% (S₂), 7.4% (S₃), and 0.0% (S₄) at the L₁ layer, 62.8% (S₁), 30.1% (S₂), 7.1% (S₃), and 0.0% (S₄) at the L₂ layer, and 73.0% (S₁), 14.7% (S₂), 12.3% (S₃), and 0.0% (S₄) at the L₃ layer. Except for the S₄ stage, the percentage of PP surface loss loads to the TP surface loss load in each planting stage gradually decreased from 93.0% (S₁) to 83.3% (S₃), while the percentages of the leached PP loss load to the leached TP loss load in each planting stage were 75.0% (S₁), 79.4% (S₂), and 76.0% (S₃). The percentages of the leached PO₄ ^3--P loss load to the leached TP loss load in each planting stage gradually increased from 5.4% (S₁) to 14.7% (S₃).

Figure 6 shows that soil layers significantly differed in each form of P’s annual migration loss loads. The annual surface loss loads of different forms of P were 1.08 kg*ha⁻¹*yr⁻¹ (TP), 0.99 kg*ha⁻¹*yr⁻¹ (PP), 0.09 kg*ha⁻¹*yr⁻¹ (TDP), and 0.07 kg*ha⁻¹*yr⁻¹ (PO₄ ^3-–P). The eventual annual leached loss load of each form’s P was 0.54 kg*ha⁻¹*yr⁻¹ (TP), 0.41 kg*ha⁻¹*yr⁻¹ (PP), 0.13 kg*ha⁻¹*yr⁻¹ (TDP), and 0.10 kg*ha⁻¹*yr⁻¹ (PO₄ ^3-–P), accounting for 50.1%, 41.3%, 152.3%, and 151.5% of their annual total surface migration loss loads, respectively. TP’s annual leached loss load from the L₁ layer was up to 3.55 kg*ha⁻¹*yr⁻¹, among which 65.2% was PP. TP’s annual leached loss loads from the L₂ and L₃ layers were 0.92 kg*ha⁻¹*yr⁻¹ and 0.54 kg*ha⁻¹*yr⁻¹, among which 68.2% and 75.8% were PP, respectively. The surface loss loads of PP and PO₄ ^3-–P accounted for 92.0% and 6.3% of the annual total surface TP loss load, respectively, while the underground eventual leached loss loads were 75.8% and 19.0% of the annual total underground leached loss load of TP, respectively. According to one-way ANOVA and the LSD test, the annual leached loss loads of TP and PP weakened significantly with the increase in soil depth (p < 0.05).