## Abstract

Nitrate (NO_{3}-N) load characteristics in consecutive dry years in the Huai River Basin (HRB), China, were examined using streamflow and NO_{3}-N concentration data. The data set spanned 12 years including three consecutive dry years. Baseflow separation, load estimation, and nonparametric linear regression were applied to separate point source (PS), baseflow, and surface runoff NO_{3}-N loads from the total load. The mean annual nonpoint source (NPS) load was 2.84 kg·ha^{−1}·yr^{−1}, accounting for 90.8% of the total load. Baseflow contributed approximately one-fourth of the natural runoff and half of the NPS load. The baseflow nitrate index (i.e., the ratio of baseflow 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 baseflow is the preferential hydrological pathway for NO_{3}-N transport in the HRB and that baseflow 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.

## HIGHLIGHTS

Baseflow is the preferential hydrological pathway for nitrate transport in the Huai River Basin.

Baseflow delivers a higher percentage of nitrate to streams in consecutive dry years than in individual dry years.

The combination of baseflow separation, load estimation, and 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 sub-basin, 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 waters are low in low-flow periods, owing to the reduced runoff, but will increase after droughts end (Schilling & Zhang 2004). However, few studies focus on the impact of long-term droughts persisting over many years on baseflow NO_{3}-N loads. If carrying much NO_{3}-N, baseflow as the main source of environmental flow in dry periods (de Graaf *et al.* 2019) may threaten the health of riverine ecosystems under long-term drought. Within the context of climate change, to support the development of effective policies to manage the water quality, such as total maximum daily load programs and best management practices, the impact of consecutive dry years on baseflow NO_{3}-N loads needs to be evaluated.

Traditionally, constituent loads are calculated as the product of the discharge volume and constituent concentration in water. However, no direct approach exists at present to continuously measure baseflow and the constituent concentration in baseflow. One of the common methods is to estimate the baseflow NO_{3}-N load by assuming the streamflow NO_{3}-N concentration, with a baseflow index (BFI, i.e., the ratio of baseflow to total streamflow) ≥90 or 80% as the real baseflow concentration combined with baseflow separating methods and load estimation programs (e.g., LOADEST; Runkel *et al.* 2004; He & Lu 2016). However, this method may be problematic when applied to regions where the point source (PS) NO_{3}-N load cannot be neglected, because the assumed baseflow NO_{3}-N concentration (i.e., streamflow NO_{3}-N concentration at BFI ≥90 or 80%) consists of both PS and baseflow contributions. Therefore, the PS NO_{3}-N load needs to be quantified when applying this method to estimate the baseflow NO_{3}-N load. However, it is difficult to directly quantify the PS contributions for regions like the HRB, where the NO_{3}-N concentration is not a routine monitoring item of sewage outfalls, in addition to the influence of the chemical and biological transformations between the source and monitoring station (Albek 2003). Albek (2003) provided an indirect method – the Kendall–Theil robust line (KTRL) regression method (Helsel & Hirsch 1992) – to estimate the PS constituent load based on discharge and concentration data measured at monitoring stations. This approach assumes that the PS constituent load remains constant over a given period and is estimated as the slope of the KTRL. In this study, the NO_{3}-N loads consisting of PS, baseflow, and surface runoff contributions at the outlet of the upper and middle HRB were estimated via a combination of baseflow separation, load estimation, and KTRL regression, as described above.

The main objectives of this study are to (1) quantify the total NPS NO_{3}-N loads exported from the upper and middle HRB and quantify the baseflow contribution, (2) identify the preferential hydrological pathway of NO_{3}-N transport, and (3) evaluate the impact of consecutive dry years on the baseflow NO_{3}-N loads.

## 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(c)) in both the flood and non-flood seasons between 2007 and 2018, excluding 2013. The NO_{3}-N concentrations were determined 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 millimeters, 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

*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 , 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

_{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*and

_{j}*b*are model coefficients, is the estimated total NO

_{j}_{3}-N load, is the estimated NO

_{3}-N load consisting of PS and baseflow contributions. and are bias correction factors. , , , and 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. 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*) in Equations (4) and (5) were calculated using the maximum likelihood estimator (MLE) (Cohn

_{j}*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).

_{3}-N load was estimated based on the mass balance (Albek 2003):where , , and 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., ), can be replaced by and Equation (6) can be written as follows:where , 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 against . 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, 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.

_{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 (), the upper limit of the PS NO

_{3}-N load for a given year can be written as follows:where is the upper limit of the PS NO

_{3}-N load,

*i*is the year, and

*j*is the

*j*th day of year

*i*. Finally, the PS NO

_{3}-N load in this study is estimated as follows:where is the final result of the PS NO

_{3}-N load estimation.

_{3}-N loads are estimated as follows:where , , and 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.

## RESULTS

### 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–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.

Items . | Model 1 . | Model 2 . |
---|---|---|

Model coefficients | a_{0} = 9.651; a_{1} = 0.908;a_{2} = 0.006; a_{3} = 0.243;a_{4} = −0.256; a_{5} = 0.409;a_{6} = 0.017 | b_{0} = 13.453; b_{1} = 0.995;b_{2} = −0.434; b_{3} = 0.055;b_{4} = −0.014; b_{5} = 0.191 |

Regression statistics | ||

R ^{2} | 0.834 | 0.839 |

S ^{2} | 0.166 | 0.049 |

PPCC | 0.991 | 0.993 |

p^{a} | 0.967 | 0.907 |

p^{b} | 0.421 | 0.126 |

Bias diagnostics | ||

Bp | −2.53 | 0.629 |

NSE | 0.773 | 0.755 |

RSR | 0.428 | 0.443 |

Items . | Model 1 . | Model 2 . |
---|---|---|

Model coefficients | a_{0} = 9.651; a_{1} = 0.908;a_{2} = 0.006; a_{3} = 0.243;a_{4} = −0.256; a_{5} = 0.409;a_{6} = 0.017 | b_{0} = 13.453; b_{1} = 0.995;b_{2} = −0.434; b_{3} = 0.055;b_{4} = −0.014; b_{5} = 0.191 |

Regression statistics | ||

R ^{2} | 0.834 | 0.839 |

S ^{2} | 0.166 | 0.049 |

PPCC | 0.991 | 0.993 |

p^{a} | 0.967 | 0.907 |

p^{b} | 0.421 | 0.126 |

Bias diagnostics | ||

Bp | −2.53 | 0.629 |

NSE | 0.773 | 0.755 |

RSR | 0.428 | 0.443 |

^{a}The PPCC test is at a significance level of 0.05.

^{b}The independent-sample *t*-test for the model residuals in the consecutive dry years versus the other years at a significance level of 0.05.

*R*^{2}, coefficient of determination; *S*^{2}, residual variance; PPCC, probability plot correlation coefficient; Bp, bias percentage; NSE, Nash–Sutcliffe efficiency index; RSR, ratio of the root mean square error to the standard deviation of the measured data.

### 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 (1956–2018) 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–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.

Year . | Annual precipitation (mm) . | PS discharge (mm) . | Natural runoff (mm) . | PS nitrate load (kg·ha ^{−1}·yr^{−1})
. | NPS NO_{3}-N load (kg·ha^{−1}·yr^{−1}). | BER . | ||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|

Baseflow . | Surface runoff . | Total natural runoff . | BFI(%) . | Baseflow . | Surface runoff . | Total Nonpoint source . | BFNI (%) . | |||||

2007 | 989.9 | 14.6 | 87.4 | 219.2 | 306.6 | 28.5 | 0.53 | 2.14 | 1.35 | 3.49 | 61.2 | 2.15 |

2008 | 947.6 | 14.6 | 63.9 | 152.4 | 216.2 | 29.5 | 0.51 | 1.50 | 1.31 | 2.81 | 53.3 | 1.81 |

2009 | 887.4 | 16.0 | 27.9 | 85.5 | 113.4 | 24.6 | 0.47 | 0.70 | 0.88 | 1.58 | 44.2 | 1.80 |

2010 | 974.7 | 16.5 | 59.5 | 190.4 | 250.0 | 23.8 | 0.45 | 1.63 | 1.97 | 3.60 | 45.2 | 1.90 |

2011 | .8005 | .166 | .220 | .434 | .653 | .336 | .032 | .071 | .040 | .112 | .638 | .190 |

2012 | .7696 | .170 | .260 | .492 | .751 | .346 | .032 | .096 | .025 | .121 | .792 | .229 |

2013 | .7481 | .168 | .163 | .381 | .544 | .299 | .032 | .067 | .036 | .103 | .653 | .218 |

2014 | 939.6 | 17.5 | 34.0 | 99.4 | 133.3 | 25.5 | 0.12 | 1.26 | 0.93 | 2.19 | 57.4 | 2.25 |

2015 | 934.4 | 17.7 | 52.4 | 146.2 | 198.6 | 26.4 | 0.12 | 1.70 | 1.80 | 3.50 | 48.5 | 1.84 |

2016 | 1,114.7 | 18.4 | 42.4 | 166.5 | 209.0 | 20.3 | 0.12 | 1.39 | 2.32 | 3.71 | 37.5 | 1.85 |

2017 | 997.4 | 21.6 | 76.1 | 202.0 | 278.1 | 27.3 | 0.12 | 1.97 | 3.08 | 5.05 | 39.0 | 1.42 |

2018 | 1,053.0 | 22.5 | 57.1 | 168.1 | 225.2 | 25.4 | 0.12 | 1.86 | 2.98 | 4.84 | 38.4 | 1.51 |

Average over the consecutive dry years | .7728 | .168 | .214 | .435 | .649 | .330 | .032 | .078 | .034 | .112 | .698 | .212 |

Mean annual average | 929.7 | 17.5 | 47.1 | 130.0 | 177.1 | 26.6 | 0.29 | 1.37 | 1.47 | 2.84 | 48.3 | 1.82 |

Year . | Annual precipitation (mm) . | PS discharge (mm) . | Natural runoff (mm) . | PS nitrate load (kg·ha ^{−1}·yr^{−1})
. | NPS NO_{3}-N load (kg·ha^{−1}·yr^{−1}). | BER . | ||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|

Baseflow . | Surface runoff . | Total natural runoff . | BFI(%) . | Baseflow . | Surface runoff . | Total Nonpoint source . | BFNI (%) . | |||||

2007 | 989.9 | 14.6 | 87.4 | 219.2 | 306.6 | 28.5 | 0.53 | 2.14 | 1.35 | 3.49 | 61.2 | 2.15 |

2008 | 947.6 | 14.6 | 63.9 | 152.4 | 216.2 | 29.5 | 0.51 | 1.50 | 1.31 | 2.81 | 53.3 | 1.81 |

2009 | 887.4 | 16.0 | 27.9 | 85.5 | 113.4 | 24.6 | 0.47 | 0.70 | 0.88 | 1.58 | 44.2 | 1.80 |

2010 | 974.7 | 16.5 | 59.5 | 190.4 | 250.0 | 23.8 | 0.45 | 1.63 | 1.97 | 3.60 | 45.2 | 1.90 |

2011 | .8005 | .166 | .220 | .434 | .653 | .336 | .032 | .071 | .040 | .112 | .638 | .190 |

2012 | .7696 | .170 | .260 | .492 | .751 | .346 | .032 | .096 | .025 | .121 | .792 | .229 |

2013 | .7481 | .168 | .163 | .381 | .544 | .299 | .032 | .067 | .036 | .103 | .653 | .218 |

2014 | 939.6 | 17.5 | 34.0 | 99.4 | 133.3 | 25.5 | 0.12 | 1.26 | 0.93 | 2.19 | 57.4 | 2.25 |

2015 | 934.4 | 17.7 | 52.4 | 146.2 | 198.6 | 26.4 | 0.12 | 1.70 | 1.80 | 3.50 | 48.5 | 1.84 |

2016 | 1,114.7 | 18.4 | 42.4 | 166.5 | 209.0 | 20.3 | 0.12 | 1.39 | 2.32 | 3.71 | 37.5 | 1.85 |

2017 | 997.4 | 21.6 | 76.1 | 202.0 | 278.1 | 27.3 | 0.12 | 1.97 | 3.08 | 5.05 | 39.0 | 1.42 |

2018 | 1,053.0 | 22.5 | 57.1 | 168.1 | 225.2 | 25.4 | 0.12 | 1.86 | 2.98 | 4.84 | 38.4 | 1.51 |

Average over the consecutive dry years | .7728 | .168 | .214 | .435 | .649 | .330 | .032 | .078 | .034 | .112 | .698 | .212 |

Mean annual average | 929.7 | 17.5 | 47.1 | 130.0 | 177.1 | 26.6 | 0.29 | 1.37 | 1.47 | 2.84 | 48.3 | 1.82 |

*Note*: The bold and italic rows denote the consecutive dry years.

### 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.

### NPS NO_{3}-N load

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 rainfall–runoff 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.

Region . | Agricultural land (%) . | Annual runoff (mm) . | NPS NO_{3}-N load(kg·ha ^{−1}·yr^{−1})
. | Citation . |
---|---|---|---|---|

Walnut Creek, USA | 86 | 69–865 | 4–66 | Jaynes et al. (1999) |

Raccoon River, USA | 76.2 | 223 | 26 | Schilling & Zhang (2004) |

Upper Sangamon River, USA | 87 | 300 | 26.4 | Guo et al. (2002) |

Central Germany | 66.2 | 117 | 18.5–41.2 | Rode et al. (2009) |

Noridc and Baltic countries | 35–99 | 130–1,246 | 1.7–102.8 | Deelstra et al. (2014) |

Huai River, China | 70.1 | 237.2 | 2.84 | This study |

Region . | Agricultural land (%) . | Annual runoff (mm) . | NPS NO_{3}-N load(kg·ha ^{−1}·yr^{−1})
. | Citation . |
---|---|---|---|---|

Walnut Creek, USA | 86 | 69–865 | 4–66 | Jaynes et al. (1999) |

Raccoon River, USA | 76.2 | 223 | 26 | Schilling & Zhang (2004) |

Upper Sangamon River, USA | 87 | 300 | 26.4 | Guo et al. (2002) |

Central Germany | 66.2 | 117 | 18.5–41.2 | Rode et al. (2009) |

Noridc and Baltic countries | 35–99 | 130–1,246 | 1.7–102.8 | Deelstra et al. (2014) |

Huai River, China | 70.1 | 237.2 | 2.84 | This study |

### 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 (Povilaitis *et al.* 2014). 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 watershed in Spain (Rodríguez-Blanco *et al.* 2013), the Mulgrave River watershed in Australia (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 & Zhang 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 NO_{3}-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 NO_{3}-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 when crops require N (Andersen *et al.* 2014). In addition, straw return (Li *et al.* 2017), permeable reactive barrier (Zhang *et al.* 2018), 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.

Time span . | Mean annual precipitation (mm) . | Baseflow . | Surface runoff . | Total natural runoff . | |||
---|---|---|---|---|---|---|---|

Discharge (mm) . | NO_{3}-N load(kg·ha ^{−1}·yr^{−1})
. | Discharge (mm) . | NO_{3}-N load(kg·ha ^{−1}·yr^{−1})
. | Discharge (mm) . | NO_{3}-N load(kg·ha ^{−1}·yr^{−1})
. | ||

Period 1 (2007–2010) | 949.9 | 59.7 | 1.5 | 161.9 | 1.38 | 221.5 | 2.87 |

Period 2 (2014–2018) | 1,007.8 | 52.4 | 1.6 | 156.4 | 2.22 | 208.8 | 3.86 |

Change | 6.1% | − 12.2% | 9.7% | − 3.4% | 61.1% | − 5.7% | 34.4% |

Time span . | Mean annual precipitation (mm) . | Baseflow . | Surface runoff . | Total natural runoff . | |||
---|---|---|---|---|---|---|---|

Discharge (mm) . | NO_{3}-N load(kg·ha ^{−1}·yr^{−1})
. | Discharge (mm) . | NO_{3}-N load(kg·ha ^{−1}·yr^{−1})
. | Discharge (mm) . | NO_{3}-N load(kg·ha ^{−1}·yr^{−1})
. | ||

Period 1 (2007–2010) | 949.9 | 59.7 | 1.5 | 161.9 | 1.38 | 221.5 | 2.87 |

Period 2 (2014–2018) | 1,007.8 | 52.4 | 1.6 | 156.4 | 2.22 | 208.8 | 3.86 |

Change | 6.1% | − 12.2% | 9.7% | − 3.4% | 61.1% | − 5.7% | 34.4% |

### 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 < , 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 few 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 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.

## ACKNOWLEDGEMENTS

This work was supported by the Major Projects of Science and Technology for Water Pollution Control and Management (Grant/Award Nos 2012ZX07204-003 and 2017ZX07602-003). The authors thank Xushui Cheng, Kailai Qian, and Tiansong Qi for their assistance with data collection. We also thank Chengpeng Lu for his valuable suggestions.

## DATA AVAILABILITY STATEMENT

All relevant data are available from an online repository or repositories (https://doi.org/10.6084/m9.figshare.14167403.v2).

## REFERENCES

*Modeling of the Water Flow and Nitrate Transport in the Shallow Aquifer of the Shaying River Basin and Its Contribution to River Pollution*

*M.S. Thesis*

*Study on the Influence of Baseflow on the Environmental Flow in Huai River Basin, China*

*PhD Thesis*