Annual nitrate load patterns in an agricultural watershed in consecutive dry years

Nitrate (NO 3 -N) load characteristics in consecutive dry years in the Huai River Basin (HRB), China, were examined using stream-ﬂ ow and NO 3 -N concentration data. The data set spanned 12 years including three consecutive dry years. Base ﬂ ow separation, load estimation, and nonparametric linear regression were applied to separate point source (PS), base ﬂ ow, and surface runoff NO 3 -N loads from the total load. The mean annual nonpoint source (NPS) load was 2.84 kg·ha (cid:1) 1 ·yr (cid:1) 1 , accounting for 90.8% of the total load. Base ﬂ ow contributed approximately one-fourth of the natural runoff and half of the NPS load. The base ﬂ ow nitrate index (i.e., the ratio of base ﬂ ow NO 3 -N load to total NPS NO 3 -N load) was 25.4% higher in consecutive dry years than in individual dry years. This study demonstrated that base ﬂ ow is the preferential hydrological pathway for NO 3 -N transport in the HRB and that base ﬂ ow delivers a higher NO 3 -N percentage to streams under long-term drought than under short-term drought. This study highlights the alarming evidence that continuous drought caused by climate change may lead to a higher rate of nitrogen loss in agricultural watersheds. nonparametric linear regression provides a convenient method to differentiate the nitrate load from different sources.


INTRODUCTION
The Huai River Basin (HRB), where the annual nitrogen fertilizer application rates amount to approximately 600 kg·N·ha À1 ·yr À1 , is one of the major grain-producing areas in China. Nitrogen fertilizers have been identified as the main contributor to the increases in the nitrate-nitrogen (NO 3 -N) concentration in groundwater in agricultural watersheds (Czekaj et al. 2016;Coskun et al. 2017;Paladino et al. 2020;Wu et al. 2021), and the HRB is no exception (Ju & Zhang 2017). The NO 3 -N concentration in groundwater below agricultural fields is generally much higher than in surface waters (Almasri & Kaluarachchi 2004). This NO 3 -N in groundwater may be primarily transferred to surface waters under baseflow conditions (He & Lu 2016;Richards et al. 2021). However, apart from Chen (2013) and Chen et al. (2017), who estimated the NO 3 -N load in groundwater discharge in a small subbasin, few researchers have reported the long-term baseflow contributions to the nonpoint source (NPS) NO 3 -N loads in streamflows in the HRB. Therefore, it is necessary to quantify the NPS NO 3 -N loads and baseflow contribution for NPS NO 3 -N assessment and management of the HRB.
Climate change is likely to make hydrological regimes more extreme, with wet seasons becoming wetter and dry seasons becoming drier. The runoff from the upper and middle HRB has decreased over the last 60 years because of regional climate change and human activities, implying that the risks and intensity of droughts have increased (An & Hao 2017). Consecutive dry years, such as those occurring from 2011 to 2013, will probably occur more frequently in the future. These changes may influence the temporal variations in NPS NO 3 -N loads. Many researchers have reported that the NO 3 -N loads transported from agricultural watersheds to surface STUDY AREA AND DATA

Watershed description
The study site is the watershed above the Bengbu Hydrometric Station (Figure 1). The area of the watershed is 121,330 km 2 . The study area mainly comprises hills and flatlands, which account for 71% of the watershed area. The quaternary layers have a thickness ranging from 50 to 200 m and mainly consist of fluvial and lacustrine loose sediments. The main aquifer is the quaternary unconsolidated rock, which is separated into shallow and deep aquifers by a thick layer of clayey soil. The thickness of shallow aquifer ranges from 50 to 60 m and the buried depth of shallow groundwater ranges from 2 to 6 m (Ge et al. 2006). The mean annual precipitation is 911 mm. The flood season lasts from June to September, and the dry season lasts from December to February. Rainfall accounts for 60% of the annual precipitation in the flood season (Figure 1(d)).
In general, most of the rivers are recharged by groundwater all year round (Pan 2014). Figure 2 shows the monthly river stages measured at four gauging stations (Figure 1(c)) from 2011 to 2012 and the nearby groundwater levels. It is evident that the groundwater levels are higher than surface water levels.
The agriculture, urban, and forest land account for 70.1, 14.2, and 6.8% of the watershed area, respectively. As a major grain-producing area in China, the study area is dominated by dry land crops and the application rate of chemical fertilizer (600 kg·ha À1 ·yr À1 ) ranks first in the country (Ju & Zhang 2017). Agricultural food processing, leather manufacturing, and papermaking are the main industries in this area.

Data
This study is based on daily stream discharge data collected from 2007 to 2018 at the Bengbu Station, which is operated by the Huai River Commission of the Water Resources Ministry of the P.R.C. The baseflow was extracted from the streamflow using the digital filter method. The wastewater discharge data originated from the Huai River Water Resources Bulletin (http://www.hrc.gov.cn). The precipitation data were downloaded from the National Meteorological Information Center (http://data.cma.cn/site).
A total of 128 surface water samples were collected at the Bengbu Station from 2007 to 2018, either monthly or bimonthly. Groundwater quality samples were collected from 122 phreatic wells (Figure 1   in the laboratory by ultraviolet spectrophotometry (detection limit ¼ 0.08 mg/L). The NO 3 -N concentration in the surface water samples ranged from 0.31 to 4.45 mg/L, with an average of 2.05 mg/L (Figure 3(a)). It was generally the highest when the streamflow was low. The result of a seasonal Kendall test (Helsel & Hirsch 1992) suggests that from 2007 to 2018, the streamflow NO 3 -N concentration exhibited no significant temporal trend (p . 0.05), indicating that the entire data set could be grouped together for analysis purposes. The streamflow NO 3 -N concentration on days when BFI ! 90% was assumed to be caused by the PS and baseflow NO 3 -N loads. This NO 3 -N concentration ranged from 1.52 to 4.45 mg/L with an average of 2.86 mg/L. The annual average groundwater NO 3 -N concentration over the 122 monitoring wells (Figure 1(c)) ranged from 4.12 mg/L in 2017 to 9.34 mg/L in 2008 (Figure 3(b)). The maximum groundwater NO 3 -N concentration reached 75.9 mg/L. This was much higher than the surface water NO 3 -N concentration.

METHODOLOGY Overall procedures
In this study, the total streamflow was assumed to be composed of discharge from sewage outfalls, baseflow, and surface runoff. The discharge from sewage outfalls was termed as anthropogenic PS discharge, and the total discharge subtracted by the anthropogenic PS discharge was termed as natural discharge (i.e., baseflow plus surface runoff). Anthropogenic PS discharge can obscure the contributions from surface runoff and baseflow to a streamflow hydrograph when using digital filter methods (Barlow et al. 2015). Therefore, the anthropogenic PS discharge was removed from the total flow before baseflow separation. Baseflow was separated from the natural streamflow using the recursive digital filter method proposed by Eckhardt (2005). The total NO 3 -N load consisted of PS and NPS contributions (i.e., anthropogenic PS discharge and the natural streamflow, respectively). First, the total load and the PS and baseflow loads were estimated using regression estimator LOADEST (Runkel et al. 2004). This regression estimator assumes that the relationship between the logarithms of the constituent concentration and discharge is approximately linear (Cohn et al. 1989). Second, the PS load was estimated using the method proposed by Albek (2003) with confining conditions. Finally, the surface runoff load and baseflow load were extracted from the total load and the PS and baseflow loads, respectively. For the purpose of comparison, the annual discharges were converted into millimeter, and the annual loads were calculated as the average daily load multiplied by the number of days per year and divided by the drainage area, including the PS discharge and PS load.

Baseflow separation
The natural flow was extracted from the total flow before baseflow separation: where Q n is the natural streamflow, k is the time step, Q t is the total streamflow, and Q p is the PS discharge, which remains constant over a given period. Assuming that the outflow from an aquifer is linear to its storage (Eckhardt 2008), the recursive digital filter method divides the natural flow into surface runoff and baseflow. The filter is expressed as follows: The filter is subject to Q b, k Q n, k , where Q b is the baseflow, Q s is the surface runoff, BFI max and α are filter parameters. The filter parameter BFI max was set to 0.8 in this study, which, according to Eckhardt (2005), is suitable for perennial streams and porous aquifers. Filter parameter α was set to 0.9982, which was determined by recession analysis according to Eckhardt (2008).

NO 3 -N load estimation
Two user-defined regression models were built in LOADEST to estimate the total NO 3 -N load (Model 1) and the PS and baseflow NO 3 -N loads (Model 2), as described in the following equations, respectively (Runkel et al. 2004): where a j and b j are model coefficients, L t is the estimated total NO 3 -N load, L pb is the estimated NO 3 -N load consisting of PS and baseflow contributions. g t and g pb are bias correction factors. g ln (Q t ), g ln (BFI), g ln (P L ), and g ln (Q p þ Q b ) are the centered values of the natural logarithms of the total streamflow, BFI, annual precipitation in the last year, and summation of the PS discharge and baseflow, respectively. g d time is the centered decimal time. The explanatory variables are centered to eliminate collinearity, and the center of a variable can be found in Runkel et al. (2004).
The explanatory variables in Model 1 (Equation (4)) are chosen to consider the influence of the total streamflow, BFI, antecedent precipitation, and seasonal factors, such as agricultural activities and temperature, which are major drivers of biological, chemical, and physical processes. Model 2 (Equation (5)) considers the influence of baseflow, antecedent precipitation, and seasonal factors. The model coefficients (a j and b j ) in Equations (4) and (5) were calculated using the maximum likelihood estimator (MLE) (Cohn et al. 1989). The MLE assumes that the model residuals follow a normal distribution with a constant variance. The calibration period ranged from 2007 to 2018. The coefficients in Model 1 were calibrated by the streamflow NO 3 -N concentration (n ¼ 128). The coefficients in Model 2 were calibrated by the streamflow NO 3 -N concentration on days when the BFI was between 90 and 100% (n ¼ 35).
The PS NO 3 -N load was estimated based on the mass balance (Albek 2003): where C t , C p , and C n are the NO 3 -N concentrations in the total streamflow, PS discharge, and natural runoff, respectively. If the total streamflow is much greater than the PS discharge (i.e., Q t ) Q p ), Q n can be replaced by Q t and Equation (6) can be written as follows: where L p, KTRL ¼ C p Q p , which is the initially estimated PS NO 3 -N load and is constant over a given period. This can be estimated by linear regression of C t against 1=Q t . For nonparametric linear regression, the KTR has been recommended to eliminate the effects of data outliers and to avoid any restrictions posed by normality assumptions (Helsel & Hirsch 1992). In this study, L p, KTRL was estimated in three time periods (2007-2010, 2011-2013, and 2014-2018) to account for the annual variations in the PS NO 3 -N load.
Assuming that the PS NO 3 -N load will remain constant over a given period and should be no higher than the NO 3 -N loads consisting of PS and baseflow contributions (L pb ), the upper limit of the PS NO 3 -N load for a given year can be written as follows: where L p, u is the upper limit of the PS NO 3 -N load, i is the year, and j is the jth day of year i. Finally, the PS NO 3 -N load in this study is estimated as follows: where L p is the final result of the PS NO 3 -N load estimation. The surface runoff, baseflow, and total NPS NO 3 -N loads are estimated as follows: where L s , L b , and L n are the estimated surface runoff, baseflow, and total NPS NO 3 -N loads, respectively. The baseflow nitrate index (BFNI, i.e., the ratio of baseflow NO 3 -N load to total NPS NO3-N load) is calculated as follows: where BFNI represents the baseflow contribution to the total NPS NO 3 -N load. The baseflow enrichment ratio (BER) is used to describe the preferential pathway of NO 3 -N transport (Schilling & Zhang 2004). The BER is calculated as follows: If NO 3 -N essentially followed the water flow, then the BFNI would be equal to the BFI and the BER would equal 1. A BER value greater than 1 implies that NO 3 -N is preferentially leached into groundwater and carried by the baseflow to streams.

Parameters calibration
The coefficients for Model 1 (Equation (4)) and Model 2 (Equation (5)) are presented in Table 1. The scatter plots of the observed load versus the estimated load are shown in Figure 4(a) and 4(b), and the correlation coefficients for Models 1 and 2 were all larger than 0.83. The scatter plots and normal probability plots for the model residuals (Figure 4(c)-4(f)) indicate that the residuals of the two models adhere to the normal distribution and constant variance assumptions of the MLE. The significance level (p . 0.05) of the independent-sample t-test for the model residuals in the consecutive dry years (2011)(2012)(2013) versus the other years indicates that climate conditions exhibited no significant differences between the consecutive dry years and the other years to the models. Therefore, it is appropriate to fit only one model for the entire study period when estimating the total NO 3 -N load or the PS and baseflow NO 3 -N loads. The bias percentage (Bp), Nash-Sutcliffe efficiency (NSE), and ratio of the root mean square error to the standard deviation of the measured data (RSR) were À2.53, 0.773, and 0.428 for Model 1, respectively, and 0.629, 0.755, and 0.443 for Model 2, respectively. The accuracy of the model simulations obtained in this study was satisfactory according to Moriasi et al. (2007), who recommended that model simulations could be considered satisfactory when NSE .0.50 and RSR 0.70.

Baseflow
The mean annual precipitation and runoff in the upper and middle reaches of the HRB were 932.7 and 237.2 mm, respectively, during our study period. The annual runoff series  were analyzed via frequency analysis based on the Pearson Type III distribution to identify the annual hydrological patterns in the study area. The results indicate that the period from 2011 to 2013 was continuously dry, 2009 was also dry, 2007 and 2017 were wet, and the other years were normal.
From 2007 to 2018, the annual PS discharge increased from 14.6 to 22.5 mm ( Table 2). The mean annual PS discharge was 17.5 mm, accounting for 9.0% of the total discharge (194.6 mm). The mean annual baseflow and surface runoff were 47.1 and 130.0 mm, accounting for 19.5 and 71.5% of the total discharge, respectively. Because of the reduced precipitation, the baseflow and surface runoff were lower in the consecutive dry years (2011)(2012)(2013) than in the individual dry year (2009) and the other normal or wet years ( Figure 5(a)). The coefficient of variance (CV) of the daily baseflow and surface runoff was 1.25 and 1.80, respectively. This result suggests that the baseflow was less variant than the surface runoff. Over the study period, the annual BFI ranged from 20.3 to 34.6% with an average of 26.6%. The BFI was the highest in the consecutive dry years, with an average of 33.0%, which was higher than the BFI in the individual dry year (24.6%) and the mean annual BFI (Figure 5(c)). In general, the BFI tended to increase in drier years because less precipitation was routed to streams as surface runoff, while the baseflow remained relatively steady.

PS NO 3 -N load
The PS NO 3 -N loads estimated by KTRL regression (L p , KTRL ) were 23.6, 10.5, and 4.1 t/d for the time periods of 2007-2010, 2011-2013, and 2014-2018, respectively (Figure 6). The upper limit of the PS NO 3 -N load is shown in Figure 7. It decreased from 17.5 t/d in 2007 to 11.1 t/d in 2018. The final PS load (L p ), estimated from Equation (9), ranged from 17.5 t/d (0.53 kg·ha À1 ·yr À1 ) to 4.1 t/d (0.12 kg·ha À1 ·yr À1 ), which also exhibited a decreasing trend (Table 2; Figure 5(b)). It decreased by 77% from 2007 to 2018. The mean annual PS load was 9.7 t/d (0.29 kg·ha À1 ·yr À1 ), which accounted for 9.2% of the total load, and the mean annual PS load in the consecutive dry years was 10.6 t/d (0.32 kg·ha À1 ·yr À1 ), which accounted for 22.2% of the total load. Because of the reduced natural runoff in the consecutive dry years, the NPS load decreased, but the proportion of the PS load increased.     The total NPS NO 3 -N loads ranged from 1.03 kg·ha À1 ·yr À1 in 2013 to 5.05 kg·ha À1 ·yr À1 in 2017 and averaged 2.84 kg·ha À1 ·yr À1 (Table 2). Over the study period, the NPS loads accounted for 90.8% of the total loads discharged from the study area. The temporal variations in the NPS loads were mainly controlled by the rainfallrunoff process, i.e., the NPS load often reached a maximum when the discharges were high (Figure 5(a) and 5(b)). The baseflow NO 3 -N load ranged from 0.67 to 2.14 kg·ha À1 ·yr À1 and averaged 1.37 kg·ha À1 ·yr À1 , while the surface runoff NO 3 -N load ranged from 0.25 to 3.08 kg·ha À1 ·yr À1 and averaged 1.47 kg·ha À1 ·yr À1 (Table 2). Baseflow and surface runoff contributed 43.8 and 47.0%, respectively, to the total load. The CV was 0.38 and 0.66 for the baseflow and surface runoff loads, respectively. This implies that the baseflow loads were less variable than the surface runoff loads. The BFNI represents the baseflow contribution to the NPS NO 3 -N load. It ranged from 37.5 to 79.2% and averaged 48.3% over the study period (Table 2; Figure 8(a)). The BFNI was the highest in the consecutive dry years with an average of 69.8%, which was higher than the BFNI in the individual dry year (44.2%) and the average over the study period (48.3%). The annual BER in the study area ranged from 1.42 to 2.29 and averaged 1.82 (Table 2; Figure 8(b)).

DISCUSSION PS and NPS NO 3 -N load patterns
Averaged over the 12-year study period, the PS, baseflow, and surface runoff NO 3 -N load contributions were 0.29, 1.37, and 1.47 kg·ha À 1 ·yr À1 , accounting for 9.2, 43.8, and 47.0% of the total load, respectively. The NO 3 -N exported from the upper and middle HRB was mainly attributed to NPSs.
Leather, paper, and chemical fertilizer industries are the main sources of PS NO 3 -N load in this area (Li et al. 2016), and the denitrification of ammonia nitrogen from domestic sewage is an important indirect source (Ding et al. 2018). With the improvement of PS pollution control measures, such as the increase of wastewater treatment plants, the improvement of sewage treatment process, and the upgrading of wastewater discharge standards (He et al. 2015), the PS NO 3 -N load declined by 77% from 2007 to 2018. This can also be interpreted from Figures 6 and 7. As shown in Figure 6, the slope of the KTRL and the Pearson correlation coefficient (r) between C t and 1/Q t decreased from 2007 to 2018, and the upper limit of the PS NO 3 -N load in Figure 7 also reveals a decreasing trend. Because the NO 3 -N concentration is not routinely monitored from the sewage outfalls, the estimated PS NO 3 -N load cannot be directly verified. According to the Huai River Water Resources Bulletin (http://www.hrc.gov.cn), both the chemical oxygen demand and ammonia nitrogen (NH 3 -N) emissions declined by approximately 70% from 2007 to 2018, which is close to the estimated decline of PS NO 3 -N load (77%). This provides indirect evidence for the correctness of the estimated PS NO 3 -N loads. In this study, the PS NO 3 -N load was estimated based on the data collected at the monitoring station. It represents the exported NO 3 -N load attributed to PSs, which was the result of PS NO 3 -N transportation and the chemical and biological transformations of PS nitrogen (e.g., nitrification and denitrification) between the source and monitoring station (Albek 2003).
The annual NPS NO 3 -N load ranged from 1.03 to 5.05 kg·ha À1 ·yr À1 and averaged 2.84 kg·ha À1 ·yr À1 . The mean annual baseflow NO 3 -N load was 1.37 kg·ha À1 ·yr À1 and accounted for 48.3% of the NPS load. This is similar to the result reported for a subbasin in the study area (0.95 kg·ha À1 ·yr À1 for the baseflow NO 3 -N load (Chen et al. 2017) and 40% for the BFNI (Chen 2013)). N fertilizers have been identified as the main source of NPS NO 3 -N pollution in this agricultural watershed (Chen 2013;Ma et al. 2019). According to Ju et al. (2010), soil residual N and nitrate leaching accounted for about 27 and 30% of the fertilizer N (approximately 600 kg·N·ha À1 ·yr À1 ), respectively. The large amount of NO 3 -N retained in the soil and that leached to groundwater will be transported to surface waters through surface runoff and baseflow. However, the NPS NO 3 -N load was much lower in this study area than that in other agricultural watersheds worldwide (Table 3). This may occur due to the high nitrogen uptake rate of crops (Ju & Zhang 2017) and the high denitrification rate in the subsoil zone (Li et al. 2020) in the study area. The denitrification in the groundwater (Kolbe et al. 2019) and hyporheic zones (Harvey et al. 2013) probably also imposes significant impacts and needs to be further investigated in this watershed.

Preferential hydrological pathway of NO 3 -N transport
Hydrological processes could be the dominant factors affecting N loss, such as baseflow (Schilling & Zhang 2004), tile drainage (Schilling & Wolter 2001;Deelstra et al. 2014), and wet deposition . Various reports have indicated that baseflow is one of the dominant hydrological pathways for NO 3 -N migration toward streams (He & Lu 2016). For example, the percentages of total NO 3 -N load contributed by baseflow were 66, 59, 52, and 62% for the Raccoon River watershed in the USA (Schilling & Zhang 2004), the Corberia   (Rasiah et al. 2013), and the Changle River watershed in China (He et al. 2020), respectively. In the upper and middle HRB, the baseflow only accounted for one-fourth of the total natural runoff, while it contributed half of the total NPS NO 3 -N load. This may result from the higher NO 3 -N concentrations in the groundwater (Figure 3). The excessive use of nitrogen fertilizers has caused extremely high NO3-N accumulation in the vadose zone of the study area (Ju 2017). Dominated by flatland, this area's groundwater buried depth is generally smaller than 6 m (Ge et al. 2006), which makes it easier for the retained NO3-N to leach to groundwater with the infiltration of rainfall. This results in much higher NO3-N concentrations in groundwater. The magnitude of BER (Table 2; Figure 8(b)) suggests that baseflow was the preferential hydrological pathway for NO 3 -N transport in the upper and middle HRB. The BER in this study area was relatively higher than that in other agricultural watersheds. For example, the BER in the Raccoon River Watershed, USA, ranged from 0.23 to 1.61 and averaged 1.23 (Schilling & Zhang 2004); the BER in the Daejeon region, South Korea, ranged from 0.6 to 1.5 (Kim et al. 2010). These results suggest that controlling the baseflow NO 3 -N concentration plays a key role in NPS NO 3 -N management of the upper and middle HRB.

Influences of the consecutive dry years on the baseflow NO 3 -N load
In the consecutive dry years, the reduced runoff resulted in a decrease in the NPS NO 3 -N load ( Figure 5(a) and 5(b)), while the proportions of the PS and baseflow NO 3 -N loads increased. The PS, baseflow, and surface runoff contributions were 22.2, 54.2, and 23.6%, respectively, to the total load. The baseflow was the major contributor to the total load in the consecutive dry years. The BFNI was 25.4% higher in the consecutive dry years than that in the individual dry year (2009), and the BER was 0.32 higher in the consecutive dry years than that in the individual dry year (Table 2). Furthermore, both the BFNI and BER are positively related to the BFI in the study area, as shown in Figure 9. These results suggest that a higher percentage of NO 3 -N will be delivered to streams through the baseflow in periods of long-term drought than in short-term drought. Climate change increases the probability of long-term drought (Paxian et al. 2019), during which more NO 3 -N retained in agricultural soils would leach to groundwater and then transfer to surface waters through baseflow. As a result, the risk of severe NO 3 -N pollution in this agricultural area may increase during the dry period. To minimize the NO3-N pollution in baseflow, appropriate measures should be applied, such as setting the threshold of N fertilizer application rate, limiting N fertilizer application on risk areas, and targeting fertilizer application to periods  , and vegetated buffer strips ( Janssen et al. 2018) are optional practical approaches to minimize the NO 3 -N pollution. Table 4 provides the mean annual runoff and NPS NO 3 -N loads for the periods before and after the consecutive dry years (Periods 1 and 2, respectively). The total natural runoff in Period 2 declined by 5.7% over that in Period 1, while the total NPS NO 3 -N load increased by 34.4%. This is mainly the result of the large amounts of NO 3 -N storage in the agricultural soils under long-term drought, which were later mobilized in periods of higher rainfall and runoff ( Ju & Zhang 2017). Compared to Period 1, the baseflow decreased by 12.2%, while the baseflow NO 3 -N load increased by 9.7% in Period 2. The surface runoff decreased by 3.4%, while the surface runoff NO 3 -N load increased by 61.1%. This indicates that more NO 3 -N was stored on the ground surface than in the groundwater under long-term drought. The surface runoff NO 3 -N loads should be controlled, especially after long-term droughts.

Impacts of the errors in PS NO 3 -N load estimation
In this study, the PS, baseflow, and surface runoff NO 3 -N loads were estimated using LOADEST and KTRL regression. Stenback et al. (2011) reported that LOADEST may overestimate the NO 3 -N loads for certain rivers, while it performs well at sites with large drainage areas (.20,000 km 2 ). The drainage area of the upper and middle HRB is much larger than 20,000 km 2 . Furthermore, Table 1 and Figure 4 reveal that LOADEST performed well in this area (i.e., NSE .0.75, Bp , +3, and RSR , 0.45).
Albek (2003) reported that for nonconservative substances, settling, chemical, and biological transformations between the source and monitoring station could induce uncertainty when estimating the PS NO 3 -N load based on the linear regression between the concentration and the reciprocal of the discharge. Because the total NO 3 -N load mainly stemmed from NPSs, we mainly analyzed the impacts of the errors in PS NO 3 -N load estimation on BFNI and BER. According to Equations (11)- (14), the upper and lower limits of BFNI and BER were calculated, when the PS NO 3 -N load ranged from 0 to its upper limit (Figure 7). A small influence interval (the difference between the upper and lower limits) indicates a nonsignificant impact of the PS NO 3 -N load errors on the BFNI and BER. It is evident from Figure 8 that the influence intervals for the BFNI (4.3-14.0%) and BER (0.17-0.52) are relatively small, and the characteristics of the interannual variations in the BFNI and BER have little changes. For example, the mean BFNI was higher in the consecutive dry years (64.9-76.3%) than that in the individual dry year (44.2-56.9%), while the BER, ranging from 1.31 to 2.45, was higher than 1 in all years. Therefore, it is reasonable to conclude that the errors in PS NO 3 -N load estimation imposed nonsignificant impacts on the BFNI and BER in this work.
Although a physical basis is lacking, the combined methods applied in this study could provide a rough estimation of the NO 3 -N loads consisting of PS, baseflow, and surface runoff contributions for regions where the PS load cannot be neglected and wastewater monitoring is not performed. Compared to physically based hydrological process models (He & Lu 2016), the above-combined methods require less data (mainly concentration and discharge data measured at monitoring stations). It is a good choice to use the combined methods to separate the PS, baseflow, and surface runoff NO 3 -N loads from the total load when few PS monitoring data are available.

CONCLUSIONS
Over a 12-year study period, the PS, baseflow, and surface runoff NO 3 -N loads were quantified in the upper and middle HRB. The PS, baseflow, and surface runoff contributions were 9.2, 43.8, and 47.0%, respectively, to the Uncorrected Proof mean annual total load (3.13 kg·ha À1 ·yr À1 ). The baseflow accounted for approximately one-fourth of the total natural runoff and contributed half of the total NPS NO 3 -N load, and the mean annual BER was 1.82. The BFNI and BER were higher in the consecutive dry years (69.8% and 2.12, respectively) than those in the individual dry years (44.2% and 1.80, respectively). We have demonstrated that NPSs are the major contributors to the total NO 3 -N load and that baseflow is the primary hydrological pathway for NO 3 -N transport in this area. Compared to short-term drought, a higher percentage of NO 3 -N will be delivered to streams through baseflow under long-term drought. This study highlights the alarming evidence that continuous drought caused by climate change may lead to a higher rate of nitrogen loss in agricultural watersheds. Therefore, controlling the baseflow NO 3 -N load (e.g., fertilizer application rate control, vegetated buffer strips, straw return, and permeable reactive barrier) plays a key role in NPS NO 3 -N management of the upper and middle HRB, especially in dry periods that span consecutive years.