Preferred hierarchical control strategy of non-point source pollution at regional scale

26 Non-point source (NPS) pollution has wide range of sources. Under rainfall conditions, NPS pollution occurs mainly by overland flow, resulting in difficult governance. In this study, based on the cooperative 28 analysis of critical periods (CPs) and critical source areas (CSAs), a preferred hierarchical control 29 strategy of NPS pollution, which was connected with management units, was proposed in the 30 Danjiangkou Reservoir Basin (DRB) to improve the pertinence of NPS pollution control. The practicality 31 of the grid-based CSA identification results was improved by point density analysis (PDA). CPs, sub- 32 CPs, and non-CPs were identified on the temporal scale; CSAs, sub-CSAs and non-CSAs were identified 33 on the spatial scale. The results showed that CPs (July, April, and September), sub-CPs (May, March, 34 and August), and non-CPs contributed 62.8%, 31.1%, and 6.1% of the annual TP loads, respectively. 35 Furthermore, we proposed a hierarchical NPS pollution control strategy: class Ⅰ (CSAs in CPs) → class 36 II (sub-CSAs in CPs, CSAs in sub-CPs) → class III (non-CPs, non-CSAs, sub- and non-CSAs in sub- 37 CPs). Class Ⅰ covered the periods and areas with the highest NPS pollution loads, contributing 26.2% of 38 the annual load within 14.5% of the area and 25.0% of the time. This study provides a reference for the 39 targeted control of NPS pollution at regional scale, especially in environmental protection with limited

shown that more than 50% of the load is generated from CPs, accounting for 25% of the time (Ruan et 54 al. 2020), and more than 50% of the load is generated form CSAs, accounting for less than 30% of the  Soranno et al. 1996). The dual-structure export empirical model (DSEEM) 67

Simulation of NPS pollution loads 114
In this study, we adopted the DSEEM to estimate the monthly TP loads from the perspectives of 115 particulate and dissolved pollutant loads based on the collected environmental data (Table 1) (Shen et al. 116 2010, Shi et al. 2002). 117

Particulate pollutant load 118
The particulate pollutant load is calculated as: 119 where L s is the particulate pollutant load (kg/ha/month), α is the transformation coefficient (kg/t), A is 121 the soil erosion (t/ha/month), η is the enrichment ratio of the soil pollutant (dimensionless), C s is the 122 concentration of particulate phosphorus in soil particles (%), and S d is the sediment transport ratio 123 (dimensionless). The values of each parameter in the formula were determined based on the background 124 data of the study area and previous studies , Zhuang et al. 2016. A is calculated using 125 MUSLE model: 126 where A is the monthly soil erosion amount (t/ha/month), R is the rainfall erosivity factor 128 (MJ·mm/ha/h/month), LS is the slope factor of slope length (dimensionless), K is the soil erosion factor 129 (t·h/ MJ/mm), C is the land cover and its management factor (dimensionless) and P is the soil and water 130 conservation measurement factor (dimensionless). 131 132

Dissolved pollutant load
where L d is the dissolved phosphorus load of ground types (kg/ha), β is the conversion coefficient, C d is 136 the concentration of pollutants (mg/L), and Q is runoff depth (mm). 137 Algorithms of all these factors are based on the background information of the study area (Shen et al. 138 2010). Q is determined using the Soil Conservation Service model: 139 where Q is the runoff depth (mm), P is the rainfall (mm), S is the maximum infiltration (mm), and CN 142 is the number of curves (dimensionless). Based on the simulation results of the monthly TP load, we adopted the load-time accumulation 147 curve and load-area accumulation curve to identify the CPs and the CSAs of the phosphorus runoff loss 148 in DRB during typical years (Ruan et al. 2020, Zhuang et al. 2016). CPs were given priority because the 149 temporal variability of NPS pollution was greater than the spatial variability in the study area. First, the 150 TP load of each month was arranged in descending order to obtain the data set (l 1 , … l i , … l m ), where l 1 151 is the maximum and l m is the minimum value. Then, the cumulative load L j and cumulative months T j 152 were calculated as: where L j is the cumulative TP load under load grade j in a month (t), l i is the TP load of the ith month (t), 156 T j is the cumulative number of months under load grade j, and t i is the number of months. 157 The cumulative percentage of TP load (Pl j ) and the corresponding cumulative percentage of time (Pt j ) 158 were calculated as: 159 where L total is the annual total TP load of the study area (t) and T total is the total number of months. 162 The cumulative load-time curve was fitted by taking Pl j as the ordinate and Pt j as the abscissa, with f as 163 the fitting function: 164 CSAs were similarly identified. The dataset (l 1 , … l i , … l n ) was generated by arranging the TP load of 166 each grid in a descending order, where l 1 is the maximum and l n is the minimum value. The relevant 167 parameters were calculated as: 168 where Pa j is the cumulative percentage of grids under load grade j, A j is the cumulative number of grids, 170 and A total is the total number of grids. 171 The cumulative load-area curve was fitted as: 172 where f is the fitting function. 174 175

Criterion selection of CPs and CSAs 176
relationship between the growth rate of the TP load and time. When k t > 1, the load grows faster than 178 time, implying that when the time increment is 1%, the load increment is greater than 1%. Here, the value 179 of k t is 1 as the index to divide the CPs and sub-CPs of phosphorus loss, and the value of k t is 0.5 as the 180 index to divide the sub-CPs and non-CPs. For the load-area curve, k a values represent the growth rate of 181 the TP load along the area. Combined with the concentration of the load in the study area, regions with 182 k a values greater than 2 were divided into CSAs, those with k a values greater than 2 and less than 1 were 183 divided into sub-CSAs, and those with k' values less than 1 were divided into non-CSAs. In addition, the 184 CSAs in the CPs, sub-CPs, and the non-CPs were identified to compare the load contribution of CSAs 185 during different periods. 186 k t and k a values were calculated as: 187

PDA of CSAs 191
The identification results of the CSAs based on the grid scale are often discretized, leading to a lack 192 of guidance in the implementation of NPS control measures. Considering the difficulty in management, 193 the environmental quality requirements of surface water and the self-purification ability of water bodies, 194 some discretely distributed CSAs do not need to adopt centralized control measures. Therefore, 195 identification results based on the cumulative load-area curve were processed further in this study. Here, 196 the PDA was employed to describe the distribution density of the original CSA grids, followed by the 197 extraction of high-density areas as CSAs that required critical control. Thus, CSAs too discrete to control 198 were screened out. Few small CSAs scattered around large CSAs with high TP loads (such as large 199 combining the distribution of slope and land-use types in the study area, few CSAs that are remote, 201 difficult to reach, and small in size were screened to improve the practicality of the identification results. (2) Load-time/area curve 216 The parameters of the fitting curves were presented, including the load-time curve, load-area curve 217 in CPs, load-area curve in sub-CPs, and load-area curve in non-CPs. The R-square values of the four 218 fitting curves were all greater than 0.94, indicating that a good fit (Table 2). 219 [Insert Table 2 here] 220 (3) CSAs after PDA 221 of the CSAs before and after PDA were compared. The results showed that the distribution of CSAs and 223 sub-CSAs became more concentrated after PDA (Fig. 5). The area proportion of CSAs and sub-CSAs 224 increased by 36.0% and 29.9% respectively, because non-CSA grids in regions with high CSA 225 aggregation were also considered when identifying CSAs. CSAs after PDA accounted for up to 26.2% 226 of the annual load with 14.5% area, and the load/area value was 1.81, which was important for the control 227 of NPS pollution in the entire region.

Spatiotemporal distribution characteristics of TP load 231
The TP load in the DRB showed significant spatiotemporal heterogeneities (Fig. 6). The variation 232 and distribution of the TP load in the reservoir area are consistent with rainfall and NDVI, respectively 233 ( Fig. 3 and Fig. 4 in the Appendix). The monthly total TP loads ranged from high to low in July, April, 234 September, May, March, August, June, October, February, November, December, and January. The total 235 TP load in June, October, February, November, December, and January were less than 0.02  10 4 t. In

236
July, the TP loads accounted for 39.2% of the entire year because of high rainfall and low vegetation 237 coverage, especially in the central part of the reservoir (Fig. 6g). Although it differed slightly from the 238 total TP loads in April, September, May, March, and August, the distribution characteristics and the main 239 factors affecting the distribution were different. In March and April with low rainfall and vegetation 240 coverage, the spatial distribution of TP loads was relatively homogeneous, especially in forest and 241 shrubland around the DRB. In May, due to the variation in NDVI, the TP load was concentrated in the 242 middle of the reservoir. During high rainfall (August and September), areas around the reservoir, which 243 concentrated in farmland and shrubland in the middle of the reservoir, and the distribution was consistent 245 with the slope (Fig. 2b in the Appendix). Therefore, rainfall and vegetation coverage are the most 246 important factors that impact the seasonal distribution of TP loads. 247 [Insert Fig. 6 here] 248 249

Fitted curve of CPs 251
CPs were identified quantitatively using the fitted cumulative load-periods curve (Fig. 9a). The 252 slope k of the fitting curve represents the increasing load rate with periods. The load changes 253 heterogeneously as the area grows. When k was 1, 27.9% of the time contributed 70.8% of the load. identified as sub-CPs, accounting for 31.0% TP load with 25.0% proportion of the period, and the 258 remaining months (June, October, February, November, December, and January) were identified as non-259 CPs, accounting for 6.1% of the TP load with 50.0% proportion of the period. 260 261

Fitted curve of CSAs in CPs 262
CSAs during the different periods (CP, sub-CP, and non-CP) were confirmed based on the TP loads 263 and the identified CPs ( Fig. 9b; 9c; 9d). The uneven increase in the k' value indicates the heterogeneity 264 of the spatial distribution of TP loads. Compared with the load distribution in sub-CPs and non-CPs, the 265 change is faster). 267 Considering the spatial distribution characteristics of the TP load in the study area, the slope of the 268 fitting curve of 2 was used as the boundary to divide the CSAs and sub-CSAs, and the slope of the fitting 269 curve of 1 was used as the boundary to divide sub-CSAs and non-CSAs. For CPs, when k' was 2, the 270 corresponding critical value of TP load was 2.72 kg/ha, and 10.7% of the areas contributed 31.3% of the 271 annual TP load. When k' was 1, the corresponding critical value of TP load was 1.35kg/ha, and 21.4% of 272 the area contributed 43.3% of the load. For sub-CSAs, when k' was 2, 11.2% of the area contributed 13.6% The results showed that farmland, shrubland, and forest were the main sources of phosphorus in the 290 region, contributing up to 91.7% of the annual TP load (Fig. 11b). Of the total TP load, 62.3% originates 291 from the CP, and 23.8% of the load originates from the CSA during the CP. Approximately 94.2% of 292 farmland was dry land, with 42.9% located in an area with a slope of 5-15° and 18.7% with a slope 293 greater than 15°. In the TP load contributed by farmland, 63.4% of the load was generated during CPs 294

Preferred hierarchical control strategies for NPS 303
Prior to causing NPS pollution, nutrients are impacted by processes such as vegetation interception, 304 river sedimentation, and water self-purification (Frankenberger et al. 2015). Therefore, all NPS sources 305 do not need to be controlled effectively. Periods and areas that generate disproportionately high pollutant 306 loads, were prioritized using the cooperative analysis of CPs and CSAs. For this, a preferred hierarchical 307 control strategy was proposed to achieve efficient short-term and small-scale prevention and control of 308 NPS pollution. Three levels were classified according to the concentration of TP loads: class Ⅰ (CSAs in 309 CSAs in sub-CPs) ( Table 3). 311 [Insert Table 3 Table 2 Fit and accuracy of cumulative load-time curve and cumulative load-area curve 529      sub-CSA 0.005 1.5 13.8 Ⅲ non-CSA 0.006 1.8 72.7 Ⅲ sum --0.019 6.1 100.0 -Class is the NPS pollution control level.