E ﬀ ects of Land Cover and Atmospheric Input on Nutrient Budget in Subtropical Mountainous Rivers, Northeastern Taiwan

: The nutrient budget, the di ﬀ erence between the nutrient output via stream and input via precipitation, can provide insights into how environmental processes a ﬀ ect forested ecosystem biogeochemistry. In this study, ﬁeld measurements of the nutrient budgets—including Na + , Cl − , K + , Mg 2 + , Ca 2 + , NO 3 − , and SO 42 − —of 19 sites were conducted in Feitsui Reservoir Watershed (FRW) of northeastern Taiwan. A series of power-law regressions were developed to establish the relationship of the nutrient budget to the discharge, nutrient input, agricultural land cover, and slope. The result show that the weekly nutrient budget is signiﬁcantly a ﬀ ected by agricultural land and input via precipitation (R 2 of regression models ≥ 0.90), yet the relationship varies among di ﬀ erent nutrient elements. The agricultural land cover is the major factor, while the input via precipitation plays a relatively minor role in the budget of Cl − , Mg 2 + , Ca 2 + , and SO 42 − . These nutrients could be provisioned abundantly from the system, and thus the input via precipitation is not the predominant controlling factor. By contrast, the Na + and K + inputs via precipitation are indispensable for accurately estimating the riverine exports. Because weathering is a limited source of K + , the roles of agricultural activities and input via precipitation are likely decisive for transport. Besides, the NO 3 − budget reveals a strong interplay between the atmospheric input and agricultural land, as expected. Because the nutrient budget model of NO 3 − is strongly improved, the R 2 changes from 0.34 to 0.99 when a larger coe ﬃ cient in exponent term (10.2) for agricultural land cover (showing that NO 3 − export is strongly hydrologically controlled) and precipitation input are included. Our analysis is based on one year of data, so extrapolating the result to a long-term period should be done with caution, as there could be substantial inter-annual variation. The nutrient budget approach provides a preliminary assessment to evaluate the impacts of agriculture and atmospheric deposition on nutrient export, which can provide a precursory reference for watershed management for improving water quality and mitigating eutrophication.


Introduction
Nutrient input from atmospheric deposition and output through stream water are the most fundamental components of the ecosystem biogeochemical cycle, and the variation often reveals changes in ecosystem processes, especially for forested ecosystems [1]. The current urban expansion Water 2020, 12, 2800 3 of 18 upstream tributaries of the reservoir watershed. The annual mean temperature of the region is 18.9 • C, with monthly temperature ranging from 12.5 (January) to 26.5 • C (July) (Figure 1). The high mean annual precipitation of 3770 mm, with 68% taking place between May and September [30,31], is also highly variable spatially, increasing from 3500 mm in the southwestern part of the FRW to 5000 mm in the northeastern of the FRW during 2001-2010 due to the mountainous landscape. Typhoons, which mainly occur in summer and autumn, usually bring large amounts of precipitation in a few days. In the winter, the northeast monsoons could also bring substantial but more gentle rainfall. During the study period (2014), two typhoons made landfall to Taiwan: typhoon Matmo on 22-23 July and typhoon Fong-Wong on 20-22 September. The weekly precipitation (streamflow) associated with the two typhoon-affected weeks was 250 (200) mm and 180 (130) mm, respectively, based on the record of watershed 19 ( Figure 1). One strong thunderstorm event from 22 to 26 June also poured down 270 (140) mm of precipitation (streamflow) (Figure 1). During the winter monsoon period (December to February), there were 10 weeks with a weekly precipitation greater than 50 mm, showing the abundant winter precipitation in the region. Long-term observation that takes inter-annual climate variation into consideration is desirable, but it is difficult to carry out for multiple sites without institutional long-standing support, especially in mountainous regions. The number of typhoons, precipitation, and temperature were not atypical compared to the records from the past two decades ( Figure S1).

Study Sites
This study was conducted in Feitsui Reservoir Watershed (FRW) in northeastern Taiwan, which is the main domestic water supply for the mega city, Taipei ( Figure 1). Totally, 7 precipitation sites and 19 stream water sites were deployed along the Pei-Shi Creek and the Di-Yu Creek, two main upstream tributaries of the reservoir watershed. The annual mean temperature of the region is 18.9 °C, with monthly temperature ranging from 12.5 (January) to 26.5 °C (July) (Figure 1). The high mean annual precipitation of 3770 mm, with 68% taking place between May and September [30,31], is also highly variable spatially, increasing from 3500 mm in the southwestern part of the FRW to 5000 mm in the northeastern of the FRW during 2001−2010 due to the mountainous landscape. Typhoons, which mainly occur in summer and autumn, usually bring large amounts of precipitation in a few days. In the winter, the northeast monsoons could also bring substantial but more gentle rainfall. During the study period (2014), two typhoons made landfall to Taiwan: typhoon Matmo on 22-23 July and typhoon Fong-Wong on 20-22 September. The weekly precipitation (streamflow) associated with the two typhoon-affected weeks was 250 (200) mm and 180 (130) mm, respectively, based on the record of watershed 19 ( Figure 1). One strong thunderstorm event from 22 to 26 June also poured down 270 (140) mm of precipitation (streamflow) (Figure 1). During the winter monsoon period (December to February), there were 10 weeks with a weekly precipitation greater than 50 mm, showing the abundant winter precipitation in the region. Long-term observation that takes interannual climate variation into consideration is desirable, but it is difficult to carry out for multiple sites without institutional long-standing support, especially in mountainous regions. The number of typhoons, precipitation, and temperature were not atypical compared to the records from the past two decades ( Figure S1).  The elevation of FRW ranges from 180 to 1127 m, with a mean slope steepness of 42%. The dominant soils are Entisols and Inceptisols, with high silt contents developed from argillite and slate with sandstone interbeds [32]. The land use types within the FRW mainly include forests (84%), agricultural lands (8.2%), and other built-up areas (7.7%) [33]. The natural secondary forests are dominated by species of Fagaceae and Lauraceae, including Persea thunbergii (Sieb and Zucc.) kosterm, Castanopsis uraiana (Hayata) Kanehira and Hatusima, Cryptocarya chinensis (Hance) Hemsl., and Araliaceae such as Schefflera octophylla (Lour.) Harms [34,35]. The agricultural activities are restricted to pre-existing farmlands because the region was designated for water resource protection following the construction of the Feitsui Reservoir in 1987. Tea plantations (approx. 1200 ha) dominate the agriculture land, and the fertilizer applications can reach 786 kg-N ha −1 yr −1 [27]. The drainage area of the 19 watersheds varying from 1.4 (site 18) to 195.4 km 2 (site 16) presents a wide range of proportions of different land use types, with agricultural land ranging from 0.2% (site 13) to 22.1% (site 4) ( Figure 1 and Table 1).

Water Collections and Chemical Analysis
From January 2014 to December 2014, bulk precipitation and stream water were collected weekly ( Figure 1 and Table 1). The bulk precipitation samples were collected with a 20 cm-diameter polyethylene (PE) bucket (approx. 10 L), while the stream water samples were taken by plunging a 1-L PE bucket into the stream. For both precipitation and stream water, a subsample (approx. 600 mL) of water was taken to measure the pH and conductivity in situ; the rest was filtered with a 0.45 µm filter. A 100 mL filtrate was reserved within a PE bottle and transported to the laboratory in the National Taiwan University in Taipei. The filtered samples were kept at 4 • C in a refrigerator without any chemical preservatives before chemical analysis. were below the detection limits for more than 75% of the precipitation and stream samples, they were excluded from subsequent analyses.

Precipitation Estimation and Stream Flow Simulation
In FRW, there are in total 10 meteorological stations; five are maintained by the Water Resource Agency (WRA) and the other five are maintained by the Central Weather Bureau (CWB) of Taiwan ( Figure 1). The daily evaporation data were provided by the Taipei Feitsui Reservoir Administration (TFRA) and the station is located at the dam of the reservoir. The areal rainfall pattern over the study area was interpolated by the Thiessen polygon approach to present the spatial variation, from which the precipitation for individual watersheds can be derived from these weekly spatial rainfall distributions. The records of discharge were provided by the TFRA, and the inflow discharge is summed in the records from two main creeks nearby sites 16 and 19 ( Figure 1). Thus, the streamflow of the ungauged watersheds was simulated via a conceptual topographic-based hydrological model, topmodel [37][38][39], using the R package dynatopmodel Ver. 1.2.1 [39]. The topmodel was used to simulate the daily streamflow considering the daily rainfall, daily evaporation, and terrain data obtained from ALOS World 3D-30m (Ver. 2.1) [40]. The observed daily streamflow was utilized to train the best parameter set to fit low, normal, and extreme values of simulated flow particularly, and the criterion of best parameters was selected based on the Kling-Gupta efficiency (KGE), where the closer a value is to unity the better [41,42]. The simulation performed well (KGE = 0.88) and the parameter set (Table S1) was applied to all watersheds, but employing their own climatic inputs and terrain information to simulate their daily streamflow. The annual runoff ratio (estimated streamflow/precipitation) of all the watersheds ranged from 0.70 to 0.80, which was consistent with our previous estimation [5]. The paired weekly ion concentrations and water quantity of precipitation and streamflow were used to calculate the precipitation input, stream water output, and nutrient input-output budgets (stream water output−precipitation input). Table 1. Basic information on the sampling sites and watersheds (referring to the locations of the sites in Figure 1).

Sites
Watershed Area (km 3

Statistical Nutrient Budget Model
In this study, the widely used rating curve method (power-law regression between nutrient export and stream discharge) was used to estimate the nutrient export (Equations (1) to (8)). The constant in the rating curve is relevant to the nutrient source and represents the basal export, as the stream discharge is zero. Although it cannot perfectly explain the nutrient export at zero stream discharge, it would not be the case for perennial streams. The exponent coefficient is also relevant to the source, but it responds to the change in the stream discharge. Conceptually, a coefficient greater than one indicates that the discharge accelerates the nutrient export, which is often referred to as "enhancement", with sediment transport as a typical example. A coefficient that is close to one presents that the nutrient export linearly increases with stream discharge. In this case, it often implies that the nutrient source is abundant. A coefficient smaller than one represents that the nutrient has a limited source, so that the rate of nutrient export decreases as the stream discharge increases, which is also termed "dilution" [19]. We established regression models of the input-output nutrient budgets in which the agricultural land cover, discharge, and atmospheric inputs are regarded as source-associated factors and are then incorporated to determine their contribution to stream nutrient export [19,[43][44][45]. Five independent variables were included in the regression models-i.e., the fraction of agricultural land cover (Agr Frac ) of each watershed (from 0.2 to 22.1 here, Table 1), the weekly discharge (Q, mm w −1 ), the weekly precipitation input for analyzed ion i (Input i , kg ha −1 w −1 ), and topographical slope and elevation. Different multiple power-law regressions with various combinations of nutrient input, runoff, Agr Frac , slope, and elevation were examined to evaluate their predictability in nutrient budgets [7,[46][47][48]. Eight settings of multiple regressions are used in this study (Equations (1)-(8)): Model 3 : Model 4 : Model 5 : Model 8 : where y i is the weekly modeled or estimated nutrient budget (kg ha −1 w −1 ) for ion species i, and Q, Agr Frac , and Input i stand for weekly discharge (mm w −1 ), the fraction of agricultural land cover, and the precipitation input for species i (kg ha −1 w −1 ), respectively. The SLOPE and Elev are average values of slope (%) and elevation (m) for watersheds. The coefficients (a, b, c, and d) were parameters determined through performance measures. The weekly nutrient data derived from 12 of 19 watersheds across a wide range of proportions of agricultural land cover (Agr Frac ranging from 0.2 to 22.1) were used to establish the regression models, while the data of the other 7 watersheds (Agr Frac ranging from 1.0 to 18.9) were utilized to validate the models ( Table 1). The Akaike Information Criterion (AIC) was applied to select the best model from several alternative models, in which the model with the lowest AIC is the most cost-effective [49,50]. The calculation of AIC can be expressed as Equation (9): where m is the maximum likelihood and n is the number of parameters in the model. The model performance between the predicted and observed values was evaluated by the coefficient of determination (R 2 ) and the Nash-Sutcliffe efficiency (NSE), in which the closer a value is to unity, the better the model will be. The residual (i.e., the difference between the predicted and observed values) was evaluated by the root mean square error (RMSE) [51,52].

Nutrient Concentrations and Fluxes
The mean annual volume-weighted mean (VWM) pH (calculated from the annual VWM H + concentration) of the seven bulk precipitation sites ranged from 4.64 to 5.05, which was significantly lower than the mean annual VWM pH of the stream water of the 19 sites, ranging from 6.69 to 7.27 (p < 0.01, Figure 2 and Table S2). The VWM ion concentrations of stream water were considerably higher than those of precipitation for all ions except K + , which were not significantly different between precipitation and stream water ( Figure 2). The VWM concentrations of Na + and Cl − , the most abundant ions in sea salts, were 180-266 and 100-198 µeq L −1 , respectively, in stream water, and were more than two times higher than those in precipitation, 73-102 and 74-108 µeq L −1 (Figure 2 and Table S2). The annual VWM concentrations of Mg 2+ and Ca 2+ were more than five times higher in stream water than in precipitation. The SO 4 2− concentration in stream water was approximately 2-4 times of that in precipitation (Table S2) Table S2). In addition, there was a trend of decrease from northeastern sites to southwestern sites in the concentrations of Na + , K + , Cl − , and Mg 2+ in precipitation, whereas the Ca 2+ , NO 3 − , and SO 4 2− did not have such a spatial pattern ( Figure 1 and Table S2). The patterns of annual ion fluxes in precipitation and stream water were similar to the patterns of ion concentrations described above (Tables S2 and S3).

The Nutrient Budget Simulation
At an annual scale, the relationships between agricultural land cover (%) and annual nutrient budgets were significant (p < 0.01) for all ions with Na + (R 2 = 0. 58 (R 2 = 0.79, p < 0.01), and NO 3 − (R 2 = 0.24, p = 0.035), but not for K + , Mg 2+ , Ca 2+ , and SO 4 2− (data not shown). The average elevation of the watershed was only significantly related to the Na + (R 2 = 0.49, mboxemphp < 0.01) and Cl − (R 2 = 0.42, p < 0.01), whereas the average slope of the watershed was not significantly related to any of the annual budgets (data not shown).

The Nutrient Budget Simulation
At an annual scale, the relationships between agricultural land cover (%) and annual nutrient budgets were significant (p < 0.01) for all ions with Na + (R 2 = 0.58), K + (R 2 = 0.75), Mg 2+ (R 2 = 0.68), Ca 2+ (R 2 = 0.59), Cl − (R 2 = 0.38), NO3 − (R 2 = 0.86), and SO4 2− (R 2 = 0.46) (Figure 3). The annual discharge was also highly correlated with the annual nutrient budgets for Na + (R 2 = 0.73, p < 0.01), Cl − (R 2 = 0.79, p < 0.01), and NO3 − (R 2 = 0.24, p = 0.035), but not for K + , Mg 2+ , Ca 2+ , and SO4 2− (data not shown). The average elevation of the watershed was only significantly related to the Na + (R 2 = 0.49, p < 0.01) and Cl − (R 2 = 0.42, p < 0.01), whereas the average slope of the watershed was not significantly related to any of the annual budgets (data not shown).  However, the relationships between nutrient budget with discharge and agricultural land cover for ions were diverse at weekly scale. Based upon the various settings of multiple power-low equations and AICs, we tried to identify the best fit model using the data of 12 watersheds ( Table 1). The results showed that "Model 4", which simultaneously took precipitation input, discharge, and agricultural land cover into account, had the most satisfactory performances for Na + , K + , Cl − , and NO 3 − , as indicated by the lowest AICs (  Table 2). All the regressed lines paralleled with 1:1 lines (Figure 4) indicated that the simulations from the models were unbiased. For the validation using data of the other seven watersheds (  (Figure 5). The interception of fitting models ranged from −0.789 for Cl − to 0.119 for Mg 2+ . The coefficients c in the fitting models, between discharge (Q) and the fraction of agriculture (Agr Frac ), indicating the strength of the interactive effects between Q and Agr Frac on nutrient export, ranged from 0.991 for Na + to 4.491 for K + and 10.245 for NO 3 − ( Table 3). The coefficients d of precipitation input, standing for the contribution of precipitation to the nutrient budget estimation, varied from −0.961 for Ca 2+ to −1.016 for SO 4 2− (Table 3). Furthermore, the exports of Cl − , Mg 2+ , Ca 2+ , and SO 4 2− could be estimated promisingly by Model 1 (considering the discharge and fraction of agriculture, Table 3), and the performances increase slightly from Model 1 to Model 4. The export estimations of Na + and K + based on Models 2-4 (R 2 > 0.90) are much better than those of Model 1 without considering the precipitation input (R 2 < 0.30, Table 2).   Table 1). The modeled budgets were developed using a multiple powerlaw regression considering the precipitation input, discharge, and fraction of agricultural land cover. Please also refer to the regression models detailed in Table 3. The red lines are regressed lines. The R 2 , NSE, and RMSE stand for the coefficient of determination, Nash-Sutcliffe efficiency, and root mean square error, respectively.  Table 3 and the observed budgets for the other seven watersheds (validated data in Table 1). The red lines are regressed lines. The R 2 , NSE, and RMSE stand for the coefficient of determination, Nash-Sutcliffe efficiency, and root mean square error, respectively.  Table 1). The modeled budgets were developed using a multiple power-law regression considering the precipitation input, discharge, and fraction of agricultural land cover. Please also refer to the regression models detailed in Table 3. The red lines are regressed lines. The R 2 , NSE, and RMSE stand for the coefficient of determination, Nash-Sutcliffe efficiency, and root mean square error, respectively.
Water 2020, 12, x FOR PEER REVIEW 9 of 17 Figure 4. The weekly nutrient budget simulation from Model 4 (Na + , K + , Cl − , and NO3 − ) and Model 6 (Mg 2+ , Ca 2+ , and SO4 2− ) between the observed budgets and the modeled budgets based on 12 watersheds (trained data in Table 1). The modeled budgets were developed using a multiple powerlaw regression considering the precipitation input, discharge, and fraction of agricultural land cover. Please also refer to the regression models detailed in Table 3. The red lines are regressed lines. The R 2 , NSE, and RMSE stand for the coefficient of determination, Nash-Sutcliffe efficiency, and root mean square error, respectively.  Table 3 and the observed budgets for the other seven watersheds (validated data in Table 1). The red lines are regressed lines. The R 2 , NSE, and RMSE stand for the coefficient of determination, Nash-Sutcliffe efficiency, and root mean square error, respectively.  Table 3 and the observed budgets for the other seven watersheds (validated data in Table 1). The red lines are regressed lines. The R 2 , NSE, and RMSE stand for the coefficient of determination, Nash-Sutcliffe efficiency, and root mean square error, respectively.  Figure 4.

Characteristics of Precipitation and Stream Water Chemistry
The annual VWM pH between precipitation and stream water showed a striking difference across all sites, in which the precipitation annual VWM pH was consistently less than 5.0 (the criterion for acid rain), while the minimal annual VWM pH of stream water was 6.69, with 13 of the 19 sites > 7.0 ( Figure 2 and Table S2). This result was consistent with that at the Fushan Experimental Forest (FEF), a forested watershed 20 km southwest of the FRW. The long-term monitoring (1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013) at the FEF revealed that about 90% of the annual VWM pH of precipitation was lower than 5.0, but there was no sign of acidification in stream water with a 20-year annual VWM pH of 6.95 [4]. We had no measurements of groundwater across various watersheds, so we did not include it in the estimation of the nutrient budgets. The input from groundwater with a higher concentration of base cations, such as Na + , Ca 2+ , and Mg 2+ , was believed to play a critical role in neutralizing the acidity in the FEF [4,53], and likely also in our study sites (Figure 2). In many temperate forests, the acidification of stream water was usually followed by the acidification of precipitation [1,54]. However, the divergence of the precipitation and stream water acidity in this subtropical region was strikingly different from those in temperate forests. The high rainfall accompanies a higher concentration of ions such as SO 4 2− and NO 3 − , leading to high acidic depositions, which are comparable as findings in FEF and are higher than in most forests worldwide [55]. In spite of there continuously being a high deposition of sulfates and nitrates in Asia and Taiwan [26], there is no sign of harmful effects on the ecosystem structure and function except in China [56]. However, the Ca 2+ deposition in FEF showed a significant declining trend over the past 20 years, indicating that the soils could be further acidified and that the effects on the ecosystem remain a critical concern [55]. The ratio of Cl − /Na + was 1.05 in precipitation, which was close to 1.1 in sea water, indicating that the composition of Na + and Cl − did not change significantly during the aerosol formation and transportation before the collection of the deposition [57,58]. However, the ratio of Cl − /Na + in stream water decreased to lower than 0.80 because the enhancement of the Na + concentration in stream water was approximately 50% higher than that of Cl − (Table S2), which was attributed to the deep chemical weathering [59]. The differences in nutrient export between disturbed and forested watersheds were most evident for K + and NO 3 − , ions associated with fertilizers, which were approximately three and eight times higher, respectively, in watershed 4 than in watershed 13 (Table S3). The relationships between the annual nutrient budgets and agricultural land cover also revealed that when the agricultural land cover was higher than 5%, the net export of K + and NO 3 − was be remarkably increased ( Figure 3). A study that examined the effects of land cover on the sediment regime in the Upper Little Tennessee River a with gentle slope (5%) in North Carolina showed that the disturbed sites with >10% of non-forested land cover would have 5-to 9-fold more suspended sediment and bedload transport than less disturbed sites (<3% of non-forested land cover) [60]. A study in the coastal plain watersheds surrounding Chesappeake Bay also suggested that the nitrate concentrations of the stream water significantly increased when the non-forested land cover was higher than 10% [44]. In our study region (FRW) with a slope steepness of 42%, the hydrochemical responses are likely more sensitive to subtle conversion from forested to non-forested land cover, because a larger slope steepness will enhance erosion [26]. Many studies have indicated that water quantity (precipitation or streamflow) exerts a strong control on the export and budget of major ions-i.e., the precipitation control or hydrological control-and not just nitrogen in undisturbed forest watersheds [4,61,62]. Our study showed that, in watersheds disturbed by anthropogenic activities, the agricultural land cover is probably an even more important factor that contributes to the elevation of the concentration and flux of most ions (Figure 3).

Indications from the Nutrient Budgets Estimations
Land cover change, especially the conversion of natural forest to cropland and urban land use, is a critical issue for soil-water resource stabilization worldwide [25,63]. The issue is more urgent in humid tropical/subtropical mountainous regions, such as Taiwan, where the plentiful rainfall, at an annual average of 2500 mm, together with the rough topography control the soil erosion and nutrient budgets [4,64]. Previous studies conducted in steep-slope watersheds in Taiwan and South Korea have suggested that the sediments and nutrient exports, such as total nitrogen and total phosphorous, could be fitted well with discharge (R 2 > 0.8) using a logarithm-or exponent-type regression model [65,66]. We further considered the interaction between discharge and the agricultural land cover as a proportion of the watershed area and atmospheric input across 19 watersheds, which significantly improved the accuracy of the nutrient budget estimation (R 2 > 0.9).
Our power-law regressions presented the different levels of incorporation of three parametersprecipitation input, discharge, and fraction of agricultural land (Agr Frac )-in nutrient budgets. The estimations of Cl − , Mg 2+ , Ca 2+ , and SO 4 2− budgets can perform well using Model 1 (considering discharge and Agr Frac ), and the performances increase slightly from Model 1 to Model 4 (considering precipitation input, discharge, and Agr Frac ). Because the Cl − is the conservative trace and there is no significant evaporite in this region [67], the satisfactory estimations of Cl − from Model 1 to Model 4 indicate that the transport of Cl − from precipitation to streamflow is relatively unreactive in ecosystems under a high Cl − deposition, as in our study region [68]. The coefficient b of the discharge term, which is lower than unity, showed a significant dilution effect during chloride transport (because the Agr Frac must be less than one). In contrast to chloride, the other three ions (Mg 2+ , Ca 2+ , and SO 4 2− ) showed some improvement from Model 1 to Model 4 (including precipitation input into estimation), indicating that the Agr Frac within watersheds is a first control on the ion export, but the atmospheric input acted as a secondary control. For Na + and K + , the estimations of Models 2-4 considering precipitation input, discharge, and Agr Frac with different settings (R 2 > 0.90) are much better than those of Model 1 (R 2 < 0.30, Table 2), demonstrating that the two ion exports require taking atmospheric input into account. Even so, there are still some differences between the two export behaviors. The Na + export, in fact, is quite similar to that of Cl − in terms of the exponent coefficient, likely because sea salt is also the main Na + source. Another main source of sodium is silicate weathering and fertilizer, and therefore its export is more complicated than that of chloride. Comparing to Na + , the K + export is more significantly regulated by Agr Frac . The higher coefficient c of exponent between discharge and Agr Frac presents the magnitude of hydrologic control (Table 3). For NO 3 − export, it is the only ion which necessitates taking the enhancement effect into account (Table 3) [5,45,69]. The higher agricultural land cover could enhance the hydrologic control significantly under the same exponent coefficient, according to Model 4 ( Table 3). The higher coefficient, such as 4.491 in K + and 10.245 in NO 3 − , indicated that the potassium and nitrate would be transported more dramatically than the other ions ( Table 3). The highest exponent coefficient connected discharge with Agr Frac in nitrate, and also implied that fertilizer application could have dominant effects on the net export of nitrate from watersheds. A study conducted on watersheds along Lanyang-Shi in northeastern Taiwan indicated that the dissolved inorganic nitrogen (DIN) deposition likely led to a high background yield, and agricultural activity within the watershed could account for 49% of the DIN export [70]. Our result is similar to that when both the precipitation input and agricultural land cover are considered for comprehending the input-output budget of nitrogen. The effects will be likely exaggerated when the amount of fertilizer applied exceeds the watershed retention capacity. If the > 700 kg-N ha −1 yr −1 application of N-fertilizers reported previously is common in the FRW [27], it would be 3.5 times that of the maximum global nitrogen application rate, 200 kg-N ha −1 yr −1 [71,72], and is likely way more than is required for plant growth. It can be expected that the excessively added nitrogen will be rapidly flushed out with runoff, especially in regions with abundant rainfall, such as our study area. Other geographical factors, such as elevation and slope, could be important factors in nutrient export if the study is conducted on a global or continental scale. From our results, the performance of regression models shows some improvement in estimating the nutrient budgets of Mg 2+ , Ca 2+ , and SO 4 2− , but not for Na + , K + , Cl − , and NO 3 − when the slope was also included to the model in addition to the discharge, fraction of agricultural land cover, and precipitation input (i.e., Model 6 in Tables 2 and 3). Adding elevation to the models (Models 7 to 8 in Table 2) did not improve the models' performance in estimating the nutrient budgets. The limited effects of slope and no effect of elevation in estimating the nutrient budgets are possibly because the slope (31-45%) and elevation (484-647 m) do not vary much among the 19 watersheds in a mountainous region spanning over a few kilometers, which minimizes their effects on modeling nutrient exports ( Table 1). The higher annual net budgets (net export) of base cations of Na + , Ca 2+ , and Mg 2+ likely indicated that the rapid chemical weathering could provide a great quantity of these cations [59,67]. Chemical weathering that dissolves rock into ions is crucial for soil production and acidity mitigation and is strongly accelerated by runoff and physical erosion, and thus it prevails in wet young mountain environments [67,73]. The long-term observation in the FEF near to FRW revealed that both the higher pH and base cation concentrations in stream water were attributed to the high contribution of groundwater to the streamflow and the continuous exchange of cations between the groundwater and stream water [4]. The considerable net export of base cations during heavy rainfall events can be also found in the Leinhuachi and Fushan Experimental forest in central and northern Taiwan [4,74], and in Liwu and Beishi river in eastern and northern Taiwan [59,67]. These results suggested that these extreme events such as typhoons and thunderstorms had a great contribution to annual exports of base cations. Agricultural practices exposed fresh soil, and the application of nitrogenous fertilizers promoted acidity generation, which both accelerated weathering and increased the base cation export into stream water. The elevated erosion rates will deplete the soil productivity and may also generate an extra carbon source and consequent eutrophication downstream [75,76]. Our analysis of 19 watersheds from one complete year reveals several important insights into the estimation of various nutrient budgets, comprising atmospheric deposition, discharge, and agricultural land cover, yet it does not take inter-annual variations into consideration. Thus, extrapolating the result to the long-term situation needs to be done with caution. Although it is challenging to take long-term measurements for many years, efforts are strongly encouraged whenever possible, as such efforts, even rare ones, are highly valuable. Despite the unclear internal processes within watersheds, our nutrient budget approach based upon a series of watersheds provides the preliminary assessment to evaluate the impacts of agriculture and atmospheric deposition on water quality, which will be helpful for water resource management.

Conclusions
The dynamics of nutrient budget-the difference between output via stream water and input via precipitation-reveal the overall outcomes of interactive effects between internal ecosystem process and external environmental forces. It is critical to confirm the key control affecting water quality under a rapidly shifting climate and land cover change. Our simple series of regression models revealed that the weekly nutrient budget can be well estimated considering discharge, agricultural land cover, atmospheric inputs, and slope simultaneously (R 2 ≥ 0.95), yet the effects of these variables on the nutrient budgets varied among different ions. The role of atmospheric input of Cl − , Mg 2+ , Ca 2+ , and SO 4 2− in the nutrient budget is minor. The slope of the watershed has some improvement on simulations of nutrient budgets of Mg 2+ , Ca 2+ , and SO 4 2− , but not for those of Na + , K + , Cl − , and NO 3 − . The role of the atmospheric input of Cl − , Mg 2+ , Ca 2+ , and SO 4 2− in nutrient budget is minor.
Furthermore, the incorporation of atmospheric Na + and K + inputs is indispensable for accurately estimating the corresponding riverine exports, since the sources of Na + and K + are in part from sea salt, except for provisioning abundantly from the watershed. The NO 3 − budget presents a strong interplay with the atmospheric input and agricultural land. Moreover, the large coefficient (10.2) in the exponent term between discharge and the fraction of agricultural land cover for the NO 3 − export indicates that NO 3 − is the strongly hydrologically controlled. Due to the difficulty of sampling 19 watersheds for multiple years, our result is based on one complete year of data, without including inter-annual variation in the analysis, and so extrapolating the result to infer long-term patterns needs to be done with caution. Although the internal processes within watersheds are still not clear, our nutrient budget approach based upon a series of watersheds provides the preliminary calculation to assess the impacts of agriculture and atmospheric deposition on water quality, suggesting that strategies of both fertilizer reduction and emission control are effective ways to improve water quality and mitigate eutrophication. If projections of future interactions between climate change and land uses are available, they would provide an opportunity to realize the biogeochemical processes of the terrestrial ecosystem under different warming scenarios [77,78].

Supplementary Materials:
The following are available online at http://www.mdpi.com/2073-4441/12/10/2800/s1: Table S1: Input parameters in Hydrological Response Unit (HRU); Table S2: The estimated annual water quantity (mm), observed volume-weighted mean (VWM) pH and concentrations (µeq L −1 ) of major elements in precipitation and streamwater measured in 2014 across 19 watersheds; Table S3: The annual total fluxes of precipitation input and stream water output of major elements in 2014 across 19 watersheds. Figure S1. The annual typhoon occurrence (a), annual mean temperature and annual precipitation (b) during the period 2001-2019 based on data from the climate station of C0A530 maintained by Central Weather Bureau (between sites number 8 and number 9) in Figure 1. The typhoon occurrence was defined as those which were < 100 km in distance from the study region to their paths at the closest (Tu et al., 2009). Figure S2: The daily hydro-climate data of 18 watersheds in 2014 except for number 19 in Figure 1. The AT, AP, and AS indicate the annual mean temperature, annual precipitation, and annual streamflow, respectively.