Rain Pattern Deeply Reshaped Total Phosphorus Load Pattern in Watershed: A Case Study from Northern China

Excessive phosphorus in aquatic systems poses a threat to ecosystem stability and human health. Precipitation has a profound influence on the phosphorus biogeochemical process; however, it has been inadequately considered at the watershed scale. In this study, the Bayesian latent variable regression model was utilized to quantify the impact of rainfall on the concentration of total phosphorus using daily monitoring data from 2019 to 2021. The results revealed a piecewise linear relationship between total phosphorus concentration and precipitation. It was further inferred that the breakpoint (The total rainfall during a single rainfall event equal to 39.4 ± 0.45 mm) represented the tipping point for the transformation of the primary river runoff generation mechanism. Subsequently, the excess phosphorus load caused by rainfall events was estimated in the Shahe basin by combining the regional nutrient management approach with the results of the Bayesian latent variable regression model. The results indicated that rainfall events were one of the most significant sources of TP loads from 2006 to 2017, accounting for 28.2% of the total. Non-artificial land, including farmland, forests, and grasslands, serves as the primary source of the excess phosphorus load resulting from rainfall events. This study provides insights into how to quantify the phosphorus load resulting from rainfall events at the basin scale, which is valuable for phosphorus management. Environmental managers should prioritize the regulation of phosphorus in non-artificial land moving forward. Implementing hierarchical management based on calibrated curve numbers and soil phosphorus content could prove to be an efficient approach for regulating phosphorus in the watershed.


Introduction
Phosphorus is a vital element in the biogeochemical cycle [1]. Almost no creature can survive without phosphorus. The distribution and total amount of phosphorus profoundly influence ecosystem dynamics [2]. Phosphorus is primarily stored in soil, water, and organisms. Insufficient phosphorus can significantly impact regional biomass and diversity [3,4], whereas excessive phosphorus can disrupt ecosystem balance [5]. Current research has advanced our understanding of the phosphorus biogeochemical process and aided in phosphorus regulation.
Simulating the phosphorus biogeochemical process in a basin remains an effective method for enhancing our understanding. Recently, several modeling methods have been developed, including the Spatially Referenced Regressions On Watershed attributes (SPAR-ROW) model [6], the Soil and Water Assessment Tool (SWAT) [7,8], and the Generalized Watershed Loading Function (GWLF) [9], among others. Additionally, simplified models such as the output coefficient method have also been widely employed [10]. These models are used to calculate the spatial and temporal distribution of phosphorus.
Based on previous research, human activities have been the primary contributor to phosphorus in rivers and lakes in recent decades. Agricultural phosphorus is consistently The Shahe basin covers an area of approximately 866.16 km 2 , with latitude and longitude ranging from 40 • 1 1 to 40 • 23 3 and from 117 • 35 10 to 118 • 5 14 ( Figure 1). The basin falls within the temperate, semi-humid monsoon climate zone. It has an average elevation of 196 m, an average annual rainfall of 754 mm, and an average annual temperature of 11.7 • C ( Figure 2). The Shahe basin exhibits a general trend of higher altitudes in the north and lower altitudes in the south, resulting in a north-to-south flow direction for the river. The northern region of the Shahe basin comprises mountainous terrain, while the central and southern regions are predominantly plains. Zunhua City is situated in the eastern part of the basin. The primary land use types in the Shahe basin include agricultural land, forest land, grassland, construction land, and water bodies ( Figure 3). Environmental managers have been grappling with the deterioration of river water quality caused by excessive nutrient discharge in recent years. Consequently, daily water quality monitoring has been implemented at the outlet of the Shahe basin since 2019 to facilitate effective water environment management. source of drinking water in Tianjin city. The Shahe River, being a crucial source for the Yuqiao Reservoir, has a direct impact on its water quality. Consequently, this poses a threat to the water security of Tianjin City downstream [20]. The Shahe basin covers an area of approximately 866.16 km 2 , with latitude and longitude ranging from 40°1'1" to 40°23'3" and from 117°35'10" to 118°5'14" (Figure 1). The basin falls within the temperate, semi-humid monsoon climate zone. It has an average elevation of 196 m, an average annual rainfall of 754 mm, and an average annual temperature of 11.7 °C (Figure 2). The Shahe basin exhibits a general trend of higher altitudes in the north and lower altitudes in the south, resulting in a north-to-south flow direction for the river. The northern region of the Shahe basin comprises mountainous terrain, while the central and southern regions are predominantly plains. Zunhua City is situated in the eastern part of the basin. The primary land use types in the Shahe basin include agricultural land, forest land, grassland, construction land, and water bodies ( Figure 3). Environmental managers have been grappling with the deterioration of river water quality caused by excessive nutrient discharge in recent years. Consequently, daily water quality monitoring has been implemented at the outlet of the Shahe basin since 2019 to facilitate effective water environment management.   source of drinking water in Tianjin city. The Shahe River, being a crucial source for the Yuqiao Reservoir, has a direct impact on its water quality. Consequently, this poses a threat to the water security of Tianjin City downstream [20]. The Shahe basin covers an area of approximately 866.16 km 2 , with latitude and longitude ranging from 40°1'1" to 40°23'3" and from 117°35'10" to 118°5'14" (Figure 1). The basin falls within the temperate, semi-humid monsoon climate zone. It has an average elevation of 196 m, an average annual rainfall of 754 mm, and an average annual temperature of 11.7 °C (Figure 2). The Shahe basin exhibits a general trend of higher altitudes in the north and lower altitudes in the south, resulting in a north-to-south flow direction for the river. The northern region of the Shahe basin comprises mountainous terrain, while the central and southern regions are predominantly plains. Zunhua City is situated in the eastern part of the basin. The primary land use types in the Shahe basin include agricultural land, forest land, grassland, construction land, and water bodies ( Figure 3). Environmental managers have been grappling with the deterioration of river water quality caused by excessive nutrient discharge in recent years. Consequently, daily water quality monitoring has been implemented at the outlet of the Shahe basin since 2019 to facilitate effective water environment management.

Data Source
This study utilized a diverse range of data, and Table 1 provides a summary of the information related to these datasets. TP concentrations exhibit two distinct temporal resolutions. Prior to 2018, TP concentrations were monitored on a monthly basis. However, since 2019, daily water quality monitoring has been conducted.

Bayesian Latent Variable Regression (BLVR)
BLVR has been widely utilized in various fields, including computer science [21], medical informatics [22], political and social research [23], and environmental sciences [24]. Bayesian regression is flexible in accommodating different types of data, allows for the reuse of experimental data, and mitigates the risk of overfitting. The experiences gained from previous studies are incorporated into the model as prior distributions. Moreover, it excels in assessing uncertainty. In accordance with the BLVR model, latent variables are utilized to describe hidden states that are challenging to detect, thereby enhancing accuracy and interpretability.
In this study, we used a piecewise linear model to describe the relationship between ∆c(TP) and PreciT. The BLVR is designed as Equations (1)- (4).

Data Source
This study utilized a diverse range of data, and Table 1 provides a summary of the information related to these datasets. TP concentrations exhibit two distinct temporal resolutions. Prior to 2018, TP concentrations were monitored on a monthly basis. However, since 2019, daily water quality monitoring has been conducted.

Bayesian Latent Variable Regression (BLVR)
BLVR has been widely utilized in various fields, including computer science [21], medical informatics [22], political and social research [23], and environmental sciences [24]. Bayesian regression is flexible in accommodating different types of data, allows for the reuse of experimental data, and mitigates the risk of overfitting. The experiences gained from previous studies are incorporated into the model as prior distributions. Moreover, it excels in assessing uncertainty. In accordance with the BLVR model, latent variables are utilized to describe hidden states that are challenging to detect, thereby enhancing accuracy and interpretability.
In this study, we used a piecewise linear model to describe the relationship between ∆c(TP) and PreciT. The BLVR is designed as Equations (1)- (4).
where ∆c(TP) indicates the increase in TP concentration during each rainfall event. PreciE indicates the total rainfall during a single rainfall event. PreciT indicates the sum of PreciE and the rainfall in the first 30 days. The intercepts (a 0 and a 1 ), slopes (b 0 and b 1 ) and stand error (σ 0 and σ 1 ) are different in every single state (Equations (1)- (3)). In the BLVR model, TP input states (k) are latent variables that are estimated by piecewise functions (Equation (4)). The ε represents the threshold for state transition.

ReNuMa (Regional Nutrient Management)
The ReNuMa model is a hydrologically driven and quasi-empirical model that builds upon two previous approaches: GWLF and NAPI (Net Anthropogenic Phosphorus Inputs). The GWLF is a lumped-parameter model that addresses hydrology, sediment, and nutrient transport at the watershed scale. The NAPI is an accounting methodology for quantifying phosphorus inputs to watersheds. It was developed for estimating nutrient fluxes in large watersheds ( Figure 4).
where ∆c(TP) indicates the increase in TP concentration during each rainfall event. PreciE indicates the total rainfall during a single rainfall event. PreciT indicates the sum of PreciE and the rainfall in the first 30 days. The intercepts (a 0 and a 1 ), slopes (b 0 and b 1 ) and stand error ( 0 and 1 ) are different in every single state (Equations (1)- (3)). In the BLVR model, TP input states (k) are latent variables that are estimated by piecewise functions (Equation (4)). The ε represents the threshold for state transition.

ReNuMa (Regional Nutrient Management)
The ReNuMa model is a hydrologically driven and quasi-empirical model that builds upon two previous approaches: GWLF and NAPI (Net Anthropogenic Phosphorus Inputs). The GWLF is a lumped-parameter model that addresses hydrology, sediment, and nutrient transport at the watershed scale. The NAPI is an accounting methodology for quantifying phosphorus inputs to watersheds. It was developed for estimating nutrient fluxes in large watersheds ( Figure 4). The GWLF model provides estimates for monthly loads of dissolved and total nitrogen and phosphorus in streamflow from complex watersheds, incorporating surface runoff, groundwater sources, and nutrient loads from point sources and on-site wastewater disposal (septic) systems.
According to the ReNuMa model, the impact of rain pattern on river TP concentration is negligible. Rain patterns are only considered by first-order wash-off functions when simulating nutrient loads in urban runoff [28]. For the other nutrient sources, monthly river TP concentrations are given by NAPI or directly defined by the user. Therefore, daily TP concentration in a single month is constant for those sources.
The ReNuMa model did not adequately consider the effect of rain patterns on nutrient loads. Thus, it can be combined with the results of the BLVR to estimate the effect of rain patterns on nutrient loads in the watershed ( Figure 5). According to ReNuMa, the basin TP load is calculated by Equation (5).
where Load TP indicates TP load in the watershed, which is calculated by flow depth (Q), TP concentration (c), and area (A). Those parameters are different for every land use and land cover (LULC). For further details, please refer to the ReNuMa model manual. Then, the rain pattern was taken into account to calculate the TP load (Equation (6)).
Load TP ' = ∆ event (Load TP )+Load TP (6) The GWLF model provides estimates for monthly loads of dissolved and total nitrogen and phosphorus in streamflow from complex watersheds, incorporating surface runoff, groundwater sources, and nutrient loads from point sources and on-site wastewater disposal (septic) systems.
According to the ReNuMa model, the impact of rain pattern on river TP concentration is negligible. Rain patterns are only considered by first-order wash-off functions when simulating nutrient loads in urban runoff [28]. For the other nutrient sources, monthly river TP concentrations are given by NAPI or directly defined by the user. Therefore, daily TP concentration in a single month is constant for those sources.
The ReNuMa model did not adequately consider the effect of rain patterns on nutrient loads. Thus, it can be combined with the results of the BLVR to estimate the effect of rain patterns on nutrient loads in the watershed ( Figure 5). According to ReNuMa, the basin TP load is calculated by Equation (5).
where Load TP indicates TP load in the watershed, which is calculated by flow depth (Q), TP concentration (c), and area (A). Those parameters are different for every land use and land cover (LULC). For further details, please refer to the ReNuMa model manual.
where Load TP ' indicates TP load after adjustment, which was the sum of Load TP and increased TP load from rainfall events (∆ event (Load TP )).
where d event indicates duration of rainfall influence. The A ' indicates land use area except water. Follow the assumption that the rain pattern is vital for TP load at the watershed scale. The influences of the rain pattern are quantified by comparisons between Load TP and Load TP ' . Nash-Sutcliffe efficiency coefficient (NSE) and R-squared (R 2 ) were employed to describe simulation accuracy (Equations (8) and (9)).

BLVR Modeling Results
Daily monitoring has been conducted at Shaheqiao station, which serves as the outlet of the Shahe basin, since 2019. TP concentration is one of the water quality indices monitored daily. The daily monitoring data for TP concentrations over a period of approximately three years has been recorded (Table 1 and Figure 6a). Subsequently, the influence of rainfall patterns on river TP concentrations and TP load was analyzed using the daily TP concentrations and precipitation data. Then, the rain pattern was taken into account to calculate the TP load (Equation (6)).
where Load TP indicates TP load after adjustment, which was the sum of Load TP and increased TP load from rainfall events (∆ event (Load TP )).
where d event indicates duration of rainfall influence. The A indicates land use area except water. Follow the assumption that the rain pattern is vital for TP load at the watershed scale. The influences of the rain pattern are quantified by comparisons between Load TP and Load TP . Nash-Sutcliffe efficiency coefficient (NSE) and R-squared (R 2 ) were employed to describe simulation accuracy (Equations (8) and (9)).

BLVR Modeling Results
Daily monitoring has been conducted at Shaheqiao station, which serves as the outlet of the Shahe basin, since 2019. TP concentration is one of the water quality indices monitored daily. The daily monitoring data for TP concentrations over a period of approximately three years has been recorded ( Table 1 and Figure 6a). Subsequently, the influence of rainfall patterns on river TP concentrations and TP load was analyzed using the daily TP concentrations and precipitation data. The seasonal distribution of TP concentrations at Shaheqiao station is shown in Figure 6a. The TP concentration was highest during the winter, subsequently decreasing to its lowest level in the spring. The fluctuation may be attributed to the death and growth of organisms. Furthermore, the river TP concentration increased during the summer due to the growth facilitated by biogeochemical and hydrological processes. Subsequently, it stabilized and declined during the autumn (Figure 6a).
Based on daily precipitation monitoring, approximately 85 precipitation events occurred in the Shahe basin from 1 January 2019, to 14 October 2022. The number of precipitation events decreased significantly with an increase in precipitation intensity. More than half of the precipitation events had an intensity below 10 mm (52.9%), while only 21.2% of the events had an intensity exceeding 40 mm (Figure 6b).
The relationship between ∆c(TP) and PreciT is shown in Figure 7. Both PreciT and PreciE contributed to an increase in river ∆c(TP). Therefore, river TP concentrations exhibited a nonlinear response to precipitation. Subsequently, BLVR was employed to further investigate the fluctuations of river ∆c(TP) in response to precipitation.  The seasonal distribution of TP concentrations at Shaheqiao station is shown in Figure 6a. The TP concentration was highest during the winter, subsequently decreasing to its lowest level in the spring. The fluctuation may be attributed to the death and growth of organisms. Furthermore, the river TP concentration increased during the summer due to the growth facilitated by biogeochemical and hydrological processes. Subsequently, it stabilized and declined during the autumn (Figure 6a).
Based on daily precipitation monitoring, approximately 85 precipitation events occurred in the Shahe basin from 1 January 2019, to 14 October 2022. The number of precipitation events decreased significantly with an increase in precipitation intensity. More than half of the precipitation events had an intensity below 10 mm (52.9%), while only 21.2% of the events had an intensity exceeding 40 mm (Figure 6b).
The relationship between ∆c(TP) and PreciT is shown in Figure 7. Both PreciT and PreciE contributed to an increase in river ∆c(TP). Therefore, river TP concentrations exhibited a nonlinear response to precipitation. Subsequently, BLVR was employed to further investigate the fluctuations of river ∆c(TP) in response to precipitation. The seasonal distribution of TP concentrations at Shaheqiao station is shown in Figure 6a. The TP concentration was highest during the winter, subsequently decreasing to its lowest level in the spring. The fluctuation may be attributed to the death and growth of organisms. Furthermore, the river TP concentration increased during the summer due to the growth facilitated by biogeochemical and hydrological processes. Subsequently, it stabilized and declined during the autumn (Figure 6a).
Based on daily precipitation monitoring, approximately 85 precipitation events occurred in the Shahe basin from 1 January 2019, to 14 October 2022. The number of precipitation events decreased significantly with an increase in precipitation intensity. More than half of the precipitation events had an intensity below 10 mm (52.9%), while only 21.2% of the events had an intensity exceeding 40 mm (Figure 6b).
The relationship between ∆c(TP) and PreciT is shown in Figure 7. Both PreciT and PreciE contributed to an increase in river ∆c(TP). Therefore, river TP concentrations exhibited a nonlinear response to precipitation. Subsequently, BLVR was employed to further investigate the fluctuations of river ∆c(TP) in response to precipitation.  distribution for all parameters, with no clear trend observed. TheR value serves as a quantitative convergence diagnostic method that compares the parameter estimates of each chain. If the chain has successfully converged and mixed, theR value should be close to 1. TheR values in Table 2, all below 1.01, indicate that the results have achieved convergence (Figure 9a). values of the parameters during the MCMC iteration process. These trajectory plots indicate that both models have reached convergence ( Figure 8). The five chains show a mixed distribution for all parameters, with no clear trend observed. The R value serves as a quantitative convergence diagnostic method that compares the parameter estimates of each chain. If the chain has successfully converged and mixed, the R value should be close to 1. The R values in Table 2, all below 1.01, indicate that the results have achieved convergence (Figure 9a).   The results of the parameter posterior distributions are summarized in Table 3. The mean of the posterior distribution for ε is 39.4 mm, based on a prior mean of 40 mm. It should be noted that the posterior distribution can be influenced by prior assumptions. To assess this, we conducted simulations with different prior means for ε and observed the resulting posterior distribution (Figure 9b). From the comparison with the posterior standard deviation of ε, the prior mean of 40 mm appears to be reasonable. Consequently,  values of the parameters during the MCMC iteration process. These trajectory plots indicate that both models have reached convergence (Figure 8). The five chains show a mixed distribution for all parameters, with no clear trend observed. The R � value serves as a quantitative convergence diagnostic method that compares the parameter estimates of each chain. If the chain has successfully converged and mixed, the R � value should be close to 1. The R � values in Table 2, all below 1.01, indicate that the results have achieved convergence (Figure 9a).   The results of the parameter posterior distributions are summarized in Table 3. The mean of the posterior distribution for ε is 39.4 mm, based on a prior mean of 40 mm. It should be noted that the posterior distribution can be influenced by prior assumptions. To assess this, we conducted simulations with different prior means for ε and observed the resulting posterior distribution (Figure 9b). From the comparison with the posterior standard deviation of ε, the prior mean of 40 mm appears to be reasonable. Consequently, The results of the parameter posterior distributions are summarized in Table 3. The mean of the posterior distribution for ε is 39.4 mm, based on a prior mean of 40 mm. It should be noted that the posterior distribution can be influenced by prior assumptions. To assess this, we conducted simulations with different prior means for ε and observed the resulting posterior distribution (Figure 9b). From the comparison with the posterior standard deviation of ε, the prior mean of 40 mm appears to be reasonable. Consequently, the threshold value for PreciE is determined to be 39.4 mm, which also serves as the breakpoint for the Bayesian regression analysis. It is important to note that the latent states depend on this threshold value, although their meanings are not clearly defined. The piecewise regression results are shown in Figure 10, suggesting that the effect of precipitation on river TP concentrations is not consistent. This phenomenon is evident in the nonlinear relationship between ∆c(TP) and PreciT as shown in Figure 8. When PreciE exceeds 39.4 mm, the slope of the regression line noticeably increases. Moreover, ∆c(TP) exhibited a moderate increase with PreciT when PreciE was below 39.4 mm. the threshold value for PreciE is determined to be 39.4 mm, which also serves as the breakpoint for the Bayesian regression analysis. It is important to note that the latent states depend on this threshold value, although their meanings are not clearly defined. The piecewise regression results are shown in Figure 10, suggesting that the effect of precipitation on river TP concentrations is not consistent. This phenomenon is evident in the nonlinear relationship between ∆c(TP) and PreciT as shown in Figure 8. When PreciE exceeds 39.4 mm, the slope of the regression line noticeably increases. Moreover, ∆c(TP) exhibited a moderate increase with PreciT when PreciE was below 39.4 mm.

Hydrological and Nutrient Modeling
The TP load was influenced by several factors in the basin. In order to assess the impact of rainfall patterns on TP loads, it was necessary to simulate the TP load production process. The ReNuMa was utilized in this study to trace nutrient sources in the basin and simulate the distribution of nutrients across different land areas. Monthly streamflow depth was simulated using the Soil Conservation Service-Curve Number (SCS-CN) method embedded within the GWLF. The Nash-Sutcliffe Efficiency (NSE) for the simulation of monthly streamflow depths was 0.84 during the training period and 0.74 during the prediction period (Figure 11a). The three reservoirs constructed upstream of the Shahe basin exert control over the streamflow. Therefore, artificial control may be the main factor driving the uncertainty in hydrological simulations. The simulation accuracy throughout the entire period was deemed acceptable (NSE = 0.78), indicating the reasonableness of the hydrological simulation. The subsequent monthly TP load simulation was conducted based on the hydrological simulation.

Hydrological and Nutrient Modeling
The TP load was influenced by several factors in the basin. In order to assess the impact of rainfall patterns on TP loads, it was necessary to simulate the TP load production process. The ReNuMa was utilized in this study to trace nutrient sources in the basin and simulate the distribution of nutrients across different land areas. Monthly streamflow depth was simulated using the Soil Conservation Service-Curve Number (SCS-CN) method embedded within the GWLF. The Nash-Sutcliffe Efficiency (NSE) for the simulation of monthly streamflow depths was 0.84 during the training period and 0.74 during the prediction period (Figure 11a). The three reservoirs constructed upstream of the Shahe basin exert control over the streamflow. Therefore, artificial control may be the main factor driving the uncertainty in hydrological simulations. The simulation accuracy throughout the entire period was deemed acceptable (NSE = 0.78), indicating the reasonableness of the hydrological simulation. The subsequent monthly TP load simulation was conducted based on the hydrological simulation. The TP load simulation exhibits low accuracy in the original ReNuMa model (Load TP ). The NSE values were only 0.48 during the training period and 0.18 during the prediction period (Figure 11b). While some uncertainty may arise from hydrological simulation, the main cause lies in the absence of specific simulations of biogeochemical processes in the Shahe basin, particularly regarding rainstorms. Therefore, more comprehensive TP load simulations were conducted during the same period, taking into account the additional TP input from rainstorms (Load TP ' ). The NSE of the simulations increased noticeably during different time periods to varying degrees ( Figure 11b and Table 4). These improvements were primarily observed during the rainy season (Figures 11 and 12). The modified TP simulations are relatively acceptable, despite their modest level of accuracy. The TP concentration observation in March 2016 was abnormal, resulting in an extremely high TP load and low precision. Furthermore, when excluding the abnormal data, the NSE of the modified TP simulations improved from 0.61 to 0.73 throughout the entire simulation period. Consequently, the modified TP load simulations are also considered reasonable. The TP load simulation exhibits low accuracy in the original ReNuMa model (Load TP ). The NSE values were only 0.48 during the training period and 0.18 during the prediction period (Figure 11b). While some uncertainty may arise from hydrological simulation, the main cause lies in the absence of specific simulations of biogeochemical processes in the Shahe basin, particularly regarding rainstorms. Therefore, more comprehensive TP load simulations were conducted during the same period, taking into account the additional TP input from rainstorms (Load TP ). The NSE of the simulations increased noticeably during different time periods to varying degrees (Figure 11b and Table 4). These improvements were primarily observed during the rainy season (Figures 11 and 12). The modified TP simulations are relatively acceptable, despite their modest level of accuracy. The TP concentration observation in March 2016 was abnormal, resulting in an extremely high TP load and low precision. Furthermore, when excluding the abnormal data, the NSE of the modified TP simulations improved from 0.61 to 0.73 throughout the entire simulation period. Consequently, the modified TP load simulations are also considered reasonable.   The ReNuMa model utilizes coupled simulations of watershed hydrological processes and nutrient processes to provide robust nutrient simulation results. In the Re-NuMa model, different LULCs are associated with distinct hydrological and biogeochemical processes. Therefore, the simulation results are capable of estimating nutrient loads from various LULCs.
Here, we calculated the TP load contributions from each type of landscape and conducted further pooled analysis ( Figure 13). The primary sources of TP load are farmland and artificial land, which have contributed 44.9% and 36.9%, respectively, over all years. The contribution from waterbodies remains relatively stable throughout the year. The contribution of artificial land to the TP load decreases over time, indicating progress in reducing urban pollution.
The annual TP load contribution pattern has been reshaped when considering TP loads associated with rainfall events. Rainfall events are among the most significant sources of TP load, accounting for 28.2% of the total. Furthermore, rainfall events have resulted in extreme TP loads during high flow years (2012-2013 and 2016-2017). The average contributions of farmland and artificial land have decreased to 31.2% and 33.3%, respectively. The ReNuMa model utilizes coupled simulations of watershed hydrological processes and nutrient processes to provide robust nutrient simulation results. In the ReNuMa model, different LULCs are associated with distinct hydrological and biogeochemical processes. Therefore, the simulation results are capable of estimating nutrient loads from various LULCs.
Here, we calculated the TP load contributions from each type of landscape and conducted further pooled analysis ( Figure 13). The primary sources of TP load are farmland and artificial land, which have contributed 44.9% and 36.9%, respectively, over all years. The contribution from waterbodies remains relatively stable throughout the year. The contribution of artificial land to the TP load decreases over time, indicating progress in reducing urban pollution.

Discussion
According to the model results and the observations, TP load is concentrated in the summer for all years ( Figure 11). Thus, the TP load in the summer has directly controlled the yearly load. Moreover, the proportions of additional TP input due to rainfall in the summer fluctuate with the precipitation. Over 40% of TP derives from rainfall events in the summer in partial years with high rainfall. However, the greatest proportion still floats around 40%, even though the TP load caused by rainfall increased sharply ( Figure 11).
There has been a noticeable reduction in urban pollution sources ( Figure 13). However, there has been a simultaneous increase in non-point source pollution caused by heavy rainfall events (Figure 14), which can potentially be attributed to the escalating frequency of extreme precipitation due to climate change. This phenomenon is significantly altering the distribution of total phosphorus (TP) sources within the basin. As a result, the The annual TP load contribution pattern has been reshaped when considering TP loads associated with rainfall events. Rainfall events are among the most significant sources of TP load, accounting for 28.2% of the total. Furthermore, rainfall events have resulted in extreme TP loads during high flow years (2012-2013 and 2016-2017). The average contributions of farmland and artificial land have decreased to 31.2% and 33.3%, respectively.

Discussion
According to the model results and the observations, TP load is concentrated in the summer for all years ( Figure 11). Thus, the TP load in the summer has directly controlled the yearly load. Moreover, the proportions of additional TP input due to rainfall in the summer fluctuate with the precipitation. Over 40% of TP derives from rainfall events in the summer in partial years with high rainfall. However, the greatest proportion still floats around 40%, even though the TP load caused by rainfall increased sharply (Figure 11).
There has been a noticeable reduction in urban pollution sources ( Figure 13). However, there has been a simultaneous increase in non-point source pollution caused by heavy rainfall events (Figure 14), which can potentially be attributed to the escalating frequency of extreme precipitation due to climate change. This phenomenon is significantly altering the distribution of total phosphorus (TP) sources within the basin. As a result, the non-linear influx of non-point source pollution is assuming a progressively crucial role. The additional load of total phosphorus (TP) originating from non-artificial land constitutes a significant portion of non-point source pollution within a watershed. This observation underscores the crucial role of meteorological factors in the management of the water environment within river basins, aligning with the findings of TP traceability analysis conducted by Zhang et al. (2020) in the Yuqiao Reservoir basin [29]. Consequently, it is imperative for the pertinent environmental management authorities to devote due attention to this emerging concern.

Discussion
According to the model results and the observations, TP load is concentrated in the summer for all years ( Figure 11). Thus, the TP load in the summer has directly controlled the yearly load. Moreover, the proportions of additional TP input due to rainfall in the summer fluctuate with the precipitation. Over 40% of TP derives from rainfall events in the summer in partial years with high rainfall. However, the greatest proportion still floats around 40%, even though the TP load caused by rainfall increased sharply (Figure 11).
There has been a noticeable reduction in urban pollution sources ( Figure 13). However, there has been a simultaneous increase in non-point source pollution caused by heavy rainfall events (Figure 14), which can potentially be attributed to the escalating frequency of extreme precipitation due to climate change. This phenomenon is significantly altering the distribution of total phosphorus (TP) sources within the basin. As a result, the non-linear influx of non-point source pollution is assuming a progressively crucial role. The additional load of total phosphorus (TP) originating from non-artificial land constitutes a significant portion of non-point source pollution within a watershed. This observation underscores the crucial role of meteorological factors in the management of the water environment within river basins, aligning with the findings of TP traceability analysis conducted by Zhang et al. (2020) in the Yuqiao Reservoir basin [29]. Consequently, it is imperative for the pertinent environmental management authorities to devote due attention to this emerging concern. The inclusion of urban surface erosion runoff loads in GWLF models allows for an assessment of the primary source of additional TP input resulting from rainfall, which is predominantly found to originate from non-artificial land. This characteristic makes the GWLF model an essential choice for conducting research in this area. Moreover, the construction of an urban drainage system in Zunhua City, situated within the Shahe Basin, has been enhanced since 2018. As a result, urban rainwater is predominantly discharged into the river through the drainage network. The calibrated threshold ε for the BLVR model is determined to be 39.4 mm, which may need to be reduced if urban runoff is found to primarily contribute to a non-linear increase in TP concentration. Thus, the non-artificial land is further verified as mainly sources to TP concentration nonlinear increase.
Some field experiments have concluded that rain patterns change nutrient concentrations in the output runoff by changing the runoff generation mechanism (RGM). In the study of river basin scale, the change in RGM is likely to be an important reason for the increase in TP in rivers during heavy rainfall. When PreciE is greater than 39.4 mm, more rainwater enters the river channel through surface runoff and also washes various substances from the rural land (e.g., farmland, woodland, and grassland) into the river. The rate and amount of phosphorus collected in this process may be higher than during periods of low runoff. The change in the TP concentration of runoff caused by the production flow mechanism is often nonlinear, which corresponds to the nonlinear relationship found in the paper.
The findings from field experiments to watershed simulations, including the present study, have consistently demonstrated the promoting effect of rain patterns on nutrient output [18,30]. It is important to note that rainfall does not exert a uniform impact on nutrient output, as it not only influences water quantity but also elevates nutrient concentration in runoff ( Figure 15). Regrettably, conventional models frequently overlook the potential influence of rainfall patterns on the concentration of total phosphorus. However, the objective of this study is to address this gap by undertaking a speculation-and simulation-based approach to unravel the intricate processes through which rainfall patterns indirectly influence the total phosphorus concentration in runoff. By comprehensively exploring the mechanisms governing nutrient production and transfer, our aim is to provide a comprehensive understanding of how these patterns ultimately shape the total phosphorus load in rivers. By solely considering a singular approach in the assessment of pollution sources, there is a risk of underestimating nutrient output from agricultural land, which indirectly impacts the formulation of effective nutrient management policies within the watershed.
trations in the output runoff by changing the runoff generation mechanism (RGM). In the study of river basin scale, the change in RGM is likely to be an important reason for the increase in TP in rivers during heavy rainfall. When PreciE is greater than 39.4 mm, more rainwater enters the river channel through surface runoff and also washes various substances from the rural land (e.g., farmland, woodland, and grassland) into the river. The rate and amount of phosphorus collected in this process may be higher than during periods of low runoff. The change in the TP concentration of runoff caused by the production flow mechanism is often nonlinear, which corresponds to the nonlinear relationship found in the paper.
The findings from field experiments to watershed simulations, including the present study, have consistently demonstrated the promoting effect of rain patterns on nutrient output [18,30]. It is important to note that rainfall does not exert a uniform impact on nutrient output, as it not only influences water quantity but also elevates nutrient concentration in runoff ( Figure 15). Regrettably, conventional models frequently overlook the potential influence of rainfall patterns on the concentration of total phosphorus. However, the objective of this study is to address this gap by undertaking a speculation-and simulation-based approach to unravel the intricate processes through which rainfall patterns indirectly influence the total phosphorus concentration in runoff. By comprehensively exploring the mechanisms governing nutrient production and transfer, our aim is to provide a comprehensive understanding of how these patterns ultimately shape the total phosphorus load in rivers. By solely considering a singular approach in the assessment of pollution sources, there is a risk of underestimating nutrient output from agricultural land, which indirectly impacts the formulation of effective nutrient management policies within the watershed. Moreover, previous studies have provided evidence for the influence of rainfall on total phosphorus (TP) loads at different basin scales. For instance, Yin et al. employed the SPARROW model to simulate TP concentrations in the expansive Poyang Lake basin. By incorporating the seasonal variation of rainfall through a Bayesian approach, they achieved significant improvements in simulation accuracy [31]. The Poyang Lake basin, encompassing an area of approximately 160,000 km 2 , signifies that the non-linear contribution of total soil phosphorus to rivers due to intense rainfall could potentially manifest on a larger spatial scale. Consequently, the impact of this mechanism on TP loads within Moreover, previous studies have provided evidence for the influence of rainfall on total phosphorus (TP) loads at different basin scales. For instance, Yin et al. employed the SPARROW model to simulate TP concentrations in the expansive Poyang Lake basin. By incorporating the seasonal variation of rainfall through a Bayesian approach, they achieved significant improvements in simulation accuracy [31]. The Poyang Lake basin, encompassing an area of approximately 160,000 km 2 , signifies that the non-linear contribution of total soil phosphorus to rivers due to intense rainfall could potentially manifest on a larger spatial scale. Consequently, the impact of this mechanism on TP loads within rivers across extensive river basins becomes non-trivial. Notably, even within smaller watersheds, the aforementioned non-linear relationship might be of greater magnitude. This is exemplified by the Lake Mendota watershed in southern Wisconsin, which spans a mere 604 km 2 . Carpenter et al. unveiled the correlation between TP loads and extreme precipitation events within this compact watershed [32]. Thus, the inference of rain patterns asserting a universal regulatory influence on TP loads can be extended across varying spatial scales.
However, it is important to note that rainfall does not always lead to an increase in TP concentrations in rivers. An example of this can be observed in Kasumigaura Lake, located in Ibaraki Prefecture, Japan, where the TP concentration in the rivers entering the lake does not show a clear relationship with rainfall in the preceding ten days. And in some cases, it is even negatively correlated [33]. This suggests that the flushing effect of rainfall events may be outweighed by their dilution effect on the TP content of the water body. In the Shahe Basin, there is a weak inverse relationship between the total monthly rainfall and the average monthly TP concentration, with the highest TP concentrations observed in winter. This indicates that the relationship between rainfall and river TP concentration may vary depending on the temporal scale, as also mentioned in existing literature [34][35][36][37][38]. These findings suggest that different driving forces may be responsible for changes in TP concentration at different time scales.
On a timescale of several days, the TP concentration in the Shahe River exhibits an increase on the day of rainfall or a few days after rainfall, with the magnitude of the increase corresponding to the intensity of the rainfall, as discussed earlier. This short-term increase in TP concentration during summer contradicts the trend observed on a monthly scale, suggesting that heavy rainfall events serve as the primary driver of TP concentration in the river within a few days, despite an overall downward trend. Similarly, on an annual scale, the patterns of variation and driving factors influencing TP concentration in the river may differ from those observed on monthly and daily scales. The regularity observed on an annual scale may be attributed to the long-term management of TP nutrient loads within the basin over the course of several years [34,39].
The daily concentration of TP in rivers is influenced by a combination of various drivers. Theoretically, any driver has the potential to be the primary driver, even if its impact may be short-lived. However, in practice, different areas exhibit unique characteristics, and not every factor within a specific watershed has the opportunity to become the main driver. This highlights the fact that in certain river basins, rainfall may not hold the potential to be a significant driver. Consequently, the concentration of TP in rivers does not always exhibit an increase in response to rainfall, even on a daily scale. This suggests that while short-term increases in nutrient load within river channels during the rainy season may be driven by rainfall, it is still crucial to establish multi-year targets for managing nutrient loads within watersheds.
The GWLF model used in this study takes into account the urban rainfall erosion effect, which implies that the contribution of TP from rainfall events primarily stems from non-point sources such as farmland, woodland, and grassland. Consequently, the soil phosphorus content in non-artificial land plays a crucial role in determining the TP load resulting from rainfall. Based on our findings, we propose that rainfall patterns indirectly influence the TP load from non-point sources and its spatiotemporal distribution by impacting the flow generation pattern. Furthermore, it can be inferred that the runoff generation capacity of different land types indirectly affects their contribution to the TP load. In light of this, managers can classify non-urban land based on soil phosphorus content and the runoff generation capacity of different land types for effective management strategies.
Human land use practices have a direct impact on runoff yield capacity (RYC) and soil TP content, as depicted by the black solid line in Figure 16. The red zone represents areas with high RYC and soil TP content, indicating a risk zone. Conversely, the green zone represents a safe zone with lower values for both RYC and soil TP content. The yellow zone corresponds to areas with a single factor exhibiting low values, designating it as a control zone. Generally, forest and grassland exhibit smaller TP content and RYC compared to farmland. In a risk zone characterized by forest and grassland, the management strategy should prioritize the reduction of RYC through increased vegetation cover, followed by implementing non-point source management measures to decrease soil phosphorus content, as illustrated by the cyan dotted line in Figure 16. In risk areas with farmland, the management approach should focus on reducing soil TP content through minimizing fertilizer application intensity, followed by regulating RYC using measures such as constructing water conservancy facilities and improving farming methods, as indicated by the blue dotted line in Figure 16. This management strategy aligns with Zhang et al., who employed the Universal Soil Loss Equation (USLE) to classify different phosphorus pollution control zones based on soil phosphorus content [40]. However, in this study, soil TP management zones are determined by calibrated CN values representing RYC and total soil TP content. The calibrated CN values in the GWLF incorporate information from the USLE and meteorological data, rendering the classification more specific compared to solely relying on the USLE. by the blue dotted line in Figure 16. This management strategy aligns with Zhang et al., who employed the Universal Soil Loss Equation (USLE) to classify different phosphorus pollution control zones based on soil phosphorus content [40]. However, in this study, soil TP management zones are determined by calibrated CN values representing RYC and total soil TP content. The calibrated CN values in the GWLF incorporate information from the USLE and meteorological data, rendering the classification more specific compared to solely relying on the USLE. Figure 16. TP management zones. (The risk zone is red; the control zone is yellow; the safe zone is green).

Conclusions
The nonlinear relationship between rain pattern and TP concentration is quantified using the BLVR model. And then we estimated the effect of this nonlinear relationship on the TP load in the Shahe basin. Furthermore, we obtained the following results:

Conclusions
The nonlinear relationship between rain pattern and TP concentration is quantified using the BLVR model. And then we estimated the effect of this nonlinear relationship on the TP load in the Shahe basin. Furthermore, we obtained the following results: