Nitrate Runoff Contributing from the Agriculturally Intensive San Joaquin River Watershed to Bay-Delta in California

Abstract

Nitrogen loading from agricultural landscapes can trigger a cascade of detrimental effects on aquatic ecosystems. Recently, the spread of aquatic weed infestations (Eichhornia crassipes, Egeria densa, Ludwigia spp., and Onagraceae) in the Sacramento-San Joaquin Delta of northern California has raised concerns, and nitrogen loading from California’s intensive farming regions is considered as one of the major contributors. In this study, we employed the Soil and Water Assessment Tool (SWAT) to simulate nitrogen exports from the agriculturally intensive San Joaquin River watershed to the Delta. The alternate tile drainage routine in SWAT was tested against monitoring data in the tile-drained area of the watershed to examine the suitability of the new routine for a tile nitrate simulation. We found that the physically based Hooghoudt and Kirkham tile drain routine improved model performance in representing tile nitrate runoff, which contributed to 40% of the nitrate loading to the San Joaquin River. Calibration results show that the simulated riverine nitrate loads matched the observed data fairly well. According to model simulation, the San Joaquin River plays a critical role in exporting nitrogen to the Delta by exporting 3135 tons of nitrate-nitrogen annually, which has a strong ecological implication in supporting the growth of aquatic weeds, which has impeded water flow, impairs commercial navigation and recreational activities, and degrades water quality in Bay-Delta waterways. Since nitrate loadings contributed by upstream runoff are an important nutrient to facilitate weed development, our study results should be seen as a prerequisite to evaluate the potential growth impact of aquatic weeds and scientific evidence for area-wide weed control decisions.

1. Introduction

Agriculture plays a vital role in feeding the growing world population of 7.6 billion people, and nitrogen fertilizer is essential for promoting and enhancing crop growth and development. However, agricultural activities continue to alter the global nitrogen cycle by increasing the input of bioavailable N to the environment [1]. There are growing concerns about the negative impacts of anthropogenic nitrogen inputs on aquatic ecosystems and human health. The United States Environmental Protection Agency (USEPA) reported nutrient enrichment as a primary reason for the pollution in streams, lakes, estuaries, and other water bodies in the U.S. [2]. Excess nitrogen can lead to water quality degradation via surface runoff into water bodies and leaching through the soil profile to groundwater. High nutrient loading can contribute to eutrophication that triggers a cascade of detrimental effects on aquatic ecosystems, including hypoxic waters, altered food webs, and harmful algal blooms [3,4,5,6]. For example, recently, the spread of aquatic weed infestations in the Sacramento-San Joaquin Delta of northern California has raised concerns. Aquatic weeds and associated pests threaten Delta’s water supply to 25 million people and the statewide irrigated agriculture of $27 billion, a $300 million recreational boating industry, and habitats for more than 750 species [7]. The Delta is also impacted by pollution from agricultural areas upstream. Given this urgent need, the United States Department of Agriculture (USDA) Agricultural Research Service has initiated the Delta Region Areawide Aquatic Weed Project (DRAAWP) to develop and implement integrated aquatic weed control programs (https://ucanr.edu/sites/DRAAWP/About_us/) [7]. To reach this goal, it is critical to understand the nitrogen loading from upstream, which is intensively cultivated in the San Joaquin River watershed and is the prerequisite to link nitrogen exposure with aquatic weed growth in Bay-Delta waterways.

Ecohydrological models are increasingly used to investigate the off-site movement of agrochemical pollutants [8,9,10,11,12,13], since continuous water quality monitoring is labor expensive and cannot cover the entire space for any time range. Statistical models, for example SPARROW (SPAtially Referenced Regressions On Watershed attributes), could also be used to simulate nitrogen loads for California streams. However, this model provides simulation results relevant to an average condition and does not capture the dynamics of nitrogen loading [14]. Another advantage of ecohydrological models is their capability to simulate the effects of various management practices in large and complicated basins on pollutant loading to receiving water bodies under various scenarios, thereby assisting with area-wide management and planning [15,16,17].

The Soil and Water Assessment Tool (SWAT) developed by USDA in the early 1990s is one of the most universally utilized water quality models [18]. It is a continuous, semi-distributed, and mechanism-based watershed-scale hydrologic model designed to investigate the effects of nonpoint-source pollution and alternative management scenarios in large river basins. Many studies reported that SWAT performed well in simulating nitrogen loads in surface water on an annual or monthly basis [19,20,21] or on a daily basis [22]. Other studies suggested that SWAT overestimated nitrate export during wet periods [23].

Although there is scientific need to understand the riverine nitrogen dynamics exported from the San Joaquin Valley, which is one of California’s most intensive farming regions, nitrogen runoff modeling is usually not the main focus for previously published studies. Most modeling studies in this region focused on other agrochemicals, such as pesticides or salinity [24,25,26]. A mechanical based nitrate simulation for the entire San Joaquin river basin is either briefly mentioned in previous studies [8,20], or usually conducted at subbasin level, for example, in the Merced river basin and Mustang river basin [27,28], which are located in the eastern mountain side of the San Joaquin river basin. Nitrogen runoff loads from the western flat region of the San Joaquin valley, for example, Orestimba Creek, Mud and Salt Sloughs, are rarely reported. However, monitoring studies indicated that almost half the nitrate loading in the San Joaquin River is contributed by west-side agricultural regions, where agricultural tile drainage is installed, contributing increasing nitrate loadings since 1970 [29].

A subsurface tile drainage system (tile drain) is an agricultural practice used to maintain acceptable root zone aeration and salinity levels for agricultural production through control of shallow water tables [30,31]. Ultimately, tile drains convey drainage water to surface water outlets. Tile drains are installed in considerable portions of the western San Joaquin River watershed, making the receiving water bodies vulnerable to nitrate pollution [14]. Therefore, when using the modeling method to investigate nitrate fate and transport in areas with tile drainage installations, it is important to first configure the model’s ability to capture the nitrate yield contributed by the tile-installed western San Joaquin River watershed, and then provide predictive modeling and/or scenario analysis for the entire basin. SWAT currently supports two methods to simulate tile drainage. The original tile drain module is contributed by three empirical parameters. The alternate tile drainage routine results from physical-based Hooghoudt and Kirkham equations [32], but has only been tested with a few datasets. Moriasi et al. [33] and Boles et al. [34] evaluated the alternate tile drainage equations based on measured nitrate-nitrogen in tile flow in the Midwestern U.S. at three 13.5 × 15 m continuous corn plots and in a small (47 km²) agricultural watershed, respectively. Evaluation against additional datasets is necessary in order to examine the suitability of the alternate tile drainage routine for tile nitrate simulation. The monitoring data in the tile-drained area of the San Joaquin River watershed in California is a good dataset to test the new routine.

This paper aims to model nitrogen loading delivered from the San Joaquin River watershed into the Delta using SWAT to advance scientific understanding of nitrogen transport processes in a partially tile-drained agricultural basin. Our primary objectives are to (1) examine the performance of the physically based tile drainage routine in nitrate simulation in the agriculturally intensive San Joaquin River watershed, (2) quantify nitrate prediction uncertainty at the tile drainage installed area and entire watershed outlet, and (3) investigate the nitrate export from the upstream San Joaquin River watershed to the Bay-Delta estuary and discuss the implications for aquatic weed growth.

2. Materials and Methods

2.1. Site Description

The San Joaquin River watershed is one of the most productive agricultural regions in California and worldwide. It is located in the middle of California’s Central Valley (Figure 1). About two-thirds of the watershed is agricultural land [35]. The gross value of the agricultural production is over $25 billion annually in the study region [36]. Its agricultural abundance includes more than 250 commodities. The major crop types include almonds, vineyards, alfalfa hay, oats, corn, cotton, and tomatoes [37]. Pasture is also a major land use type, as three-quarters of the dairy herd in California is located in the San Joaquin River watershed [35].

The study watershed is drained by the San Joaquin River that originates in the Sierra Nevada Mountains, flows northwest across the central valley, and eventually drains to the San Francisco Bay-Delta estuary. The three major tributaries to the San Joaquin River are the Merced, the Tuolumne, and the Stanislaus Rivers, all located on the east side of the mainstem San Joaquin River (Figure 1). Tile drainage is mostly installed at the west side flat agricultural land, with a covered area of 27,721 hectares. Regional observation reported a 25 mg/L N concentration in the tile flow to the Mud Slough [38]. A long-term United States Geological Survey (USGS) monitoring site (#11303500) on the San Joaquin River near Vernalis was chosen as the entire watershed outlet with the four inlets defined at the USGS monitoring sites on the major eastern rivers (Figure 1). Another USGS site (#11261500) is located in the middle of the watershed, which collects outflows mostly from the western flat region. This region has a typical Mediterranean climate, with a wet and cool winter, but a dry and hot summer. Previous studies reported 10–12 inches of mean annual precipitation in the San Joaquin river basin [39].

2.2. Data Sources and Preprocessing

2.2.1. Input Data

Required inputs to the SWAT model include topography, stream network (in flat terrain), weather, land use, soil type, and management data. Topographic data were obtained from the 3D Elevation Program [40]. Retrieved data are the 1/3 arc-second (10 m) digital elevation models, with a vertical accuracy of 2 m (expressed as the root mean square error). Stream network data were collected from the 1:100,000 scale National Hydrography Dataset [41]. We used the stream network data to match the SWAT subbasin reaches with known stream hydrography [42]. This is most useful when dealing with flat terrain, where delineation based on topography data could be difficult. Daily climate data (precipitation, temperature, solar radiation, humidity, and wind speed) were obtained from the California Irrigation Management Information System (CIMIS) and the Climate Forecast System Reanalysis (CFSR) global meteorological dataset [43,44]. Land-use data were retrieved from the 2014 Cropland Data Layer [37]. Soil properties were obtained from the Soil Survey Geographic (SSURGO) database [45].

Nitrogen fertilizer information (application rates at each crop) was estimated based on the information from literature review and use rate recommendations from farm advisors [46,47,48,49,50,51,52,53,54,55]. Corn and oats were defined to have 80% of the N fertilizer as manure [53,56]. For winter wheat, the value was designated as 40% [53]. Continuous fertilizer operation was enabled for pastures where fresh manure was applied daily. For the remaining crop types, only inorganic fertilizer was applied.

2.2.2. Monitoring Data

The monitoring data for nitrate concentrations in surface water were obtained from two sources, the National Water Information System (NWIS) through the Water Quality Portal [57], and the California Environmental Data Exchange Network (CEDEN), created by the State Water Resources Control Board [58]. Nitrate concentrations were available for the San Joaquin River at Fremont Ford Bridge (USGS #11261500) and for the watershed outlet, the San Joaquin River near Vernalis (USGS #11303500). The rloadest package, the USGS water science R functions for Load Estimator (LOADEST), was utilized to estimate monthly nitrate loads [59].

2.3. Model Setup

ArcSWAT (version 2012.10_2.18) was employed to organize and prepare spatial input for the SWAT model. We partitioned the watershed into 27 subbasins and 647 hydrologic response units (HRUs) where the HRUs were defined based on an overlap of topography, soil, and land use maps using a 5% coverage threshold. The study period was from 1/1/2001 to 7/31/2014, as the most recent data in the CFSR database were available up to 7/31/2014. The years of 2001 and 2002 were used as the model initialization period. We used the next 6 year period of data for model calibration and data for the period 1/1/2009 to 7/31/2014 for model validation.

Both the default and alternate tile drainage routines of SWAT were evaluated in this study. For the original tile drainage routine, tile flow is computed by water table depth, tile drain depth, and the time draining soils to field capacity. The alternate tile drainage routine employs the Hooghoudt and Kirkham tile drain equations to calculate tile flow, which considers water table depth, tile drain depth, size, and spacing [33]. It also introduces a drainage coefficient that confine the flow rate similar to that discovered in actual tile drainage systems [34]. Details of the alternate tile drainage routine can be found in Moriasi et al. [32].

We assumed an equivalent network of parallel tile drains based on the tile-drained areas in the western San Joaquin River watershed derived from the National Resources Inventory [60], following Saleh and Domagalski [14]. 28 tile-drained HRUs in Los Banos Creek subbasin (Figure 1) were identified by overlaying the tile drainage map with the HRU layer. Crop production in the tile-drained areas were corn, cotton, tomatoes, alfalfa hay, and almond orchards, consistent with a previous report [61]. Most of the soils in the tile-drained areas are hydrologic soil groups D and C (slow infiltration rate), except for two almond HRUs where the soils were classified as group B. The tile drain parameters were set at the starting point of calibration (

Table 1. Parameters for the original and alternate tile drainage routines of the Soil and Water Assessment Tool.

Parameter	Description	Value
Parameters for the original tile drainage routine
ITDRN.bsn	Tile drainage equations flag	0 = original, 1 = DRAINMOD
DDRAIN.mgt	Depth to the subsurface drain (mm)	1510
TDRAIN.mgt	Time to drain soil to field capacity (h)	24
GDRAIN.mgt	Drain tile lag time (h)	96
DEP_IMP.hru	Depth to impervious layer in soil profile (mm)	Approximated by depth to the bottom of the soil profile
Additional parameters for the alternate tile drainage routine
DRAIN_CO.sdr	Daily drainage coefficient (mm/day)	35
LATKSATF.sdr	Multiplication factor to determine lateral ksat (conk(j1,j)) from SWAT ksat input value (sol_k(j1,j)) for HRU	1
RE.sdr	Effective radius of drains (mm)	20
SDRAIN.sdr	Distance between two drain tubes or tiles (mm)	30,000

) to compare the performance of the two tile drainage routines. The default values for tile drain parameters were determined based on the SWAT manual, literature review, and reports on the tile-drained areas in the San Joaquin River watershed [32,61,62]. Later, we used only the alternate tile drainage routine in subsequent analyses due to its better performance in simulating nitrate loads, which is considered in more detail in the Results and Discussion sections.

2.4. Calibration and Uncertainty Analysis

Chen et al. [63] determined the SWAT hydrologic parameters for the San Joaquin River watershed previously; therefore, our calibration here focused on the nitrogen parameters only while keeping the hydrologic parameters unchanged. The Sequential Uncertainty Fitting version 2 (SUFI-2) algorithm embedded in the SWAT-CUP program was used for auto-calibration and uncertainty analysis. Aggregated uncertainties came from model input, monitoring data, model structure, and parameterization [64,65]. The parameter uncertainty was propagated through SWAT to generate a 95% prediction uncertainty band (95 PPU) in the output by SUFI-2. In order to evaluate the prediction uncertainty band, two indices were developed to evaluate observation versus simulation, the P-factor and the R-factor. The P-factor is the percentage of the observed data bracketed by the 95 PPU. A large P-factor indicates that a large portion of the model uncertainty is being accounted for [66]. The R-factor is the ratio between the average thickness of the 95 PPU uncertainty band and the standard deviation of measurement. The R-factor implies the width of the uncertainty interval. Therefore, small values (close to 1) are desirable. A trade-off should be sought between the two factors, since the larger P-factor can be obtained when the R-factor is allowed to increase [67].

In this study, initial one-at-a-time sensitivity analysis indicated that nutrient transport from the landscape was the most sensitive process. SWAT was first manually calibrated using nitrate loads measured at the two monitoring stations beginning with the upstream station and subsequently moving to the watershed outlet. Initial ranges of the relevant parameters in the upstream and downstream subbasins were determined based on the results of manual calibration. Automatic calibration was performed thereafter, which included both monitoring stations in the objective function $\mathsf{\varphi }$:

$$\mathsf{\varphi }=\{\begin{array}{cc}\left|\mathrm{b}\right|{\mathrm{R}}^2\hfill & \mathrm{if} \left|\mathrm{b}\right|\le 1\hfill \\ {\left|\mathrm{b}\right|}^{-1}{\mathrm{R}}^2\hfill & \mathrm{if} \left|\mathrm{b}\right|>1 \hfill \end{array}$$

where $\mathsf{\varphi }$ is the weighted R², b is the gradient of the regression line, and R² is the coefficient of determination. Considering both b and R², $\mathsf{\varphi }$ quantifies systematic over- or under-prediction together with dispersion [68], with ranges only between 0 and 1, so model performance is not governed by the worst events [66,67]. After each iteration, SWAT-CUP calculated the 95 PPU, the P-factor, and the R-factor, and new parameter ranges for the next iteration. Iterations were not finished until satisfactory P-factor and the R-factor values were obtained. Here, we performed four iterations with 500 SWAT runs each to achieve satisfactory modeling results.

2.5. Performance Measures

Deterministic simulation results (SWAT with the best parameter estimates) are also presented so that the established performance measures (PMs) and performance evaluation criteria (PEC) can be applied to evaluate simulation results. The following PMs (

Table 2. Performance measures and associated performance evaluation criteria.

Performance Ratings	Nitrogen
	${\mathbf{R}}^2=(\frac{{{\displaystyle \sum }}_{\mathbf{i}=1}^{\mathbf{n}}\left({\mathbf{y}}_{\mathbf{i}}-\overline{\mathbf{y}}\right)\left(\widehat{{\mathbf{y}}_{\mathbf{i}}}-\overline{\widehat{\mathbf{y}}}\right)}{\sqrt{{{\displaystyle \sum }}_{\mathbf{i}=1}^{\mathbf{n}}{\left({\mathbf{y}}_{\mathbf{i}}-\overline{\mathbf{y}}\right)}^2}\sqrt{{{\displaystyle \sum }}_{\mathbf{i}=1}^{\mathbf{n}}{\left(\widehat{{\mathbf{y}}_{\mathbf{i}}}-\overline{\widehat{\mathbf{y}}}\right)}^2}}{)}^2$	$\mathbf{NSE}=1-\frac{{{\displaystyle \sum }}_{\mathbf{i}=1}^{\mathbf{n}}{\left({\mathbf{y}}_{\mathbf{i}}-\widehat{{\mathbf{y}}_{\mathbf{i}}}\right)}^2}{{{\displaystyle \sum }}_{\mathbf{i}=1}^{\mathbf{n}}{\left({\mathbf{y}}_{\mathbf{i}}-\overline{\mathbf{y}}\right)}^2}$	$\mathbf{PPBIAS}=\frac{{{\displaystyle \sum }}_{\mathbf{i}=1}^{\mathbf{n}}\left(\widehat{{\mathbf{y}}_{\mathbf{i}}}-{\mathbf{y}}_{\mathbf{i}}\right)}{{{\displaystyle \sum }}_{\mathbf{i}=1}^{\mathbf{n}}{\mathbf{y}}_{\mathbf{i}}}\times 100 (\%)$
Very good	${\mathrm{R}}^2>0.70$	$\mathrm{NSE}>0.65$	$\left|\mathrm{PBIAS}\right|<15$
Good	$0.60<{\mathrm{R}}^2\le 0.70$	$0.50<\mathrm{NSE}\le 0.65$	$15\le \left|\mathrm{PBIAS}\right|<20$
Satisfactory	$0.30<{\mathrm{R}}^2\le 0.60$	$0.35<\mathrm{NSE}\le 0.50$	$20\le \left|\mathrm{PBIAS}\right|<30$
Unsatisfactory	${\mathrm{R}}^2\le 0.30$	$\mathrm{NSE}\le 0.35$	$\left|\mathrm{PBIAS}\right|\ge 30$

Note: ${\mathrm{y}}_{\mathrm{i}}$ and $\widehat{{\mathrm{y}}_{\mathrm{i}}}$ are the ith observed and modeled values, $\mathrm{y}¯$ and $\overline{\widehat{\mathrm{y}}}$ is the average of the observed and modeled values, and $\mathrm{n}$ is the sample size.

) were chosen in this study: (1) coefficient of determination (R²), (2) Nash-Sutcliffe efficiency (NSE), and (3) percent bias (PBIAS). The corresponding PEC for these PMs were adopted from a recent meta-analysis on PMs and PEC for widely used watershed-scale models including SWAT [69]. The statistical analysis and data visualization were performed in R [70].

3. Results

3.1. Tile Drainage Simulation

Virtually all the tile drainage from the western San Joaquin River watershed has been discharged directly to the San Joaquin River since 1985 [38]. The Fremont Ford Bridge monitoring station (USGS #11261500) along the San Joaquin River is located immediately downstream from the discharge point, making it the perfect location to use in evaluating the capability of the alternate tile drainage routine in SWAT. Under default settings, the two tile drainage routines produced almost identical streamflow at the Fremont station, while peak nitrate loads were reduced by half when the alternate tile drainage routine was used (Figure 2). By incorporating the more physically based drain equations (Hooghoudt and Kirkham tile drainage algorithm) with a drainage coefficient, the alternate tile drainage routine led to better match to measured nitrate loads. According to the model simulation, tile nitrate accounted for 40% of nitrate yields in the San Joaquin River watershed from 2003 to 2015, which is consistent with previous published studies [29,71].

3.2. Sensitive SWAT Parameters

Table 3. Sensitive parameters, calibrated ranges, and optimal values.

Parameter	Description	Lower Limit	Upper Limit	Optimal Value
CMN.bsn	Rate factor for humus mineralization of active organic nutrients (N and P)	0.000145	0.000249	0.000169
CDN.bsn	Denitrification exponential rate coefficient	1.1	3	2.357197
SDNCO.bsn	Denitrification threshold water content	0.57	1.1	1.001510
NPERCO.bsn	Nitrate percolation coefficient	0	0.2	0.176200
ANION_EXCL.sol	Fraction of porosity from which anions are excluded	0.01	0.79	0.446
HLIFE_NGW.gw	Half-life of nitrate in the shallow aquifer (days)	33	200	90.870819
DEP_IMP.hru	Depth to impervious layer in soil profile (m)	0	6	1.23500
BC1.swq	Rate constant for biological oxidation of NH4 to NO2 (day−1)	0.21	1	0.909150
BC2.swq	Rate constant for biological oxidation of NO2 to NO3 (day−1)	0.37	2	1.789952
BC3.swq	Rate constant for hydrolysis of organic N to NH4 (day−1)	0.24	0.4	0.267680

lists the set of sensitive parameters for nitrogen, their calibrated ranges and values of the optimal performance. Simulated nutrient fate and transport in SWAT is affected by the land phase controlling nutrient loading from HRUs and subbasins, and the channel phase regulating the movement of nutrients via the main channel and in-stream nutrient cycling [72].

We found that the parameters affecting nitrogen transport from the landscape were the most sensitive in SWAT predictions. These included the rate factor for humus mineralization of active organic nutrients (CMN.bsn), the denitrification exponential rate coefficient (CDN.bsn), the denitrification threshold water content (SDNCO.bsn), the nitrate percolation coefficient (NPERCO.bsn), the fraction of porosity from which anions are excluded (ANION_EXCL.sol), and the half-life of nitrate in the shallow aquifer (HLIFE_NGW.gw). Three channel parameters were found to be sensitive for nitrate load simulation at the watershed outlet. These parameters were the rate constants for biological oxidation of NH4 to NO2 (BC1.swq), NO2 to NO3 (BC2.swq), and N to NH4 (BC3.swq). Since tile drainage is only installed in the Los Banos Creek subbasin, the DEM_IMP parameter is only assigned to HRUs in this subbasin.

3.3. Simulation of Riverine Nitrate Loads

According to the model simulation, an average of 3135 tons of nitrate were exported annually through the San Joaquin River near the Vernalis outlet over the study period (Figure 3). Annual nitrate loads at the Vernalis outlet ranged from 1657 tons/year in 2013 to 6283 tons/year in 2006. The 95 PPU band (shaded area in Figure 3) was calculated based on the final parameter ranges presented in Table 3. For the San Joaquin River at Fremont Ford Bridge, the 95 PPU band bracketed 56% and 67% of the measured data during calibration and validation, respectively (

Table 4. Calibration and uncertainty analysis for monthly nitrate loads in the San Joaquin River watershed. NSE, Nash-Sutcliffe efficiency; PBIAS, percent bias.

Station	P-Factor	R-Factor	R2	R2 Rating	NSE	NSE Rating	PBIAS (%)	PBIAS Rating
Calibration of nitrate simulation (2003–2008)
Vernalis	0.53	0.63	0.68	Good	0.45	Satisfactory	−22	Satisfactory
Fremont Ford Bridge	0.56	1.72	0.67	Good	−0.11	Unsatisfactory	16	Good
Validation of nitrate simulation (2009–2014)
Vernalis	0.51	0.86	0.71	Very good	0.25	Unsatisfactory	24	Satisfactory
Fremont Ford Bridge	0.67	1.82	0.52	Satisfactory	−0.73	Unsatisfactory	29	Satisfactory

). For the San Joaquin River near Vernalis, more than 50% of the measurements were captured by the 95 PPU. The average values for the R-factor were 0.74 and 1.77 for the Vernalis and Fremont stations, respectively. For deterministic simulations, nitrate simulation results were rated as at least “satisfactory” based on R2 and PBIAS during the calibration and validation periods at both stations. In particular, the performance rating was at least “satisfactory” according to all three metrics during the calibration period at the watershed outlet. NSE ratings were consistently lower as compared to the other two metrics. Illustrations of the observed and simulated nitrate loads indicate that some peak nitrate loads were overestimated, especially at the Fremont station (Figure 3). The poor performance in NSE is caused by the mismatch of peaks in 2005 and 2011. Since NSE computes the differences between the observation and simulation as square values (Table 2), it is sensitive to peak values [73]. However, SWAT still tracked the monthly nitrate loadings in the San Joaquin River reasonably well for other months.

4. Discussion

4.1. Tile Drainage Simulation

The default tile drainage routine calculates drainage flow timing based on only one parameter, TDRAIN (time to drain soils to field capacity). This method does not differentiate TDRAIN between large storms and small events [34]. The alternate tile drainage routine in SWAT incorporates the more physically based Hooghoudt and Kirkham tile drain equations, instead of a simple parameter TDRAIN. It also incorporates a drainage coefficient that limits the possible drainage volume in 24 h, therefore reducing the allowed maximal nitrate loading rates. The alternate tile drainage routine changes the shape of the nitrate runoff hydrograph. In that way, the pollutographs with abrupt peaks are altered to ones with smooth, more widely spread “bumps” (Figure 2). As a result, model predictions compared fairly well with the observed data. Our results show that the physically based Hooghoudt and Kirkham tile drain equations together with a drainage coefficient are able to reliably predict nitrate losses delivered by tile drains. No significant differences in the predicted streamflow compared to the measured values were found between the original and alternate tile drainage routines. This result occurs because the flow volume at the Fremont Ford Bridge is heavily dependent on upstream reservoir releases [24]. As a result, tile drainage has little impact on the overall streamflow trends.

We found that about 40% of the simulated nitrate losses were from tile nitrate in the San Joaquin River watershed over the study period. This finding is consistent with our knowledge that tile drains are critical for nitrate movement in surface water within the San Joaquin River watershed [14]. Results from a previous study suggested that tile drainage contributed 22% of the total nitrogen loads in the San Joaquin River at Vernalis for the years 1985 to 2004 [38]. This percentage was derived from measured nitrogen loads and a methodology for estimating tile drain discharge and concentrations [39]. Since total nitrogen is the sum of nitrate-nitrogen, nitrite-nitrogen, ammonia-nitrogen, and organic nitrogen, we consider our estimates comparable to the previous study.

4.2. Simulation of Riverine Nitrate Loads

In general, monthly nitrate loads in the San Joaquin River were well captured by SWAT, except for the overestimation of peak loads at the Fremont station (Figure 3). As a result, NSE showed less correspondence between observed and simulated values (Table 4), despite R² and PBIAS ratings of at least “satisfactory” that indicated both strong linear relationships between the observed and simulated loads and low average model simulation bias (over- vs. under-estimation). One of the major disadvantages of NSE is its sensitivity to large values in a time series [68]. Therefore, overestimation of peak loads is likely to lead to poor NSE performance ratings.

In hydrologic modeling, uncertainties in predicted output may arise from model structure, model parameters, and spatial and temporal variability in the input data [64,74]. SUFI-2 maps all parameter uncertainty sources and aims to achieve the smallest range of parameter uncertainty while bracketing most of the measurements in the 95 PPU band. These two aspects were quantified by the P-factor, which measures to what extent model uncertainties are being considered, and the R-factor, which indicates the strength of calibration [66]. Results from uncertainty analysis indicate that there was a balance between uncertainty measures and the strength of calibration (Table 4). Although the R-factor was slightly greater than the desired value of 1, we did not continue to refine the uncertainty analysis because a further decrease in the parameter ranges would lead to a P-factor of less than 0.5.

In this study, the major sources of uncertainty in nitrogen simulation include the lack of information about the rate and date of fertilizer and manure application at the HRU level. Fertilization application rates and timing vary among regions and growers, depending on weather patterns, water availability, crop prices, yield goals, etc. (M. Lundy, personal communication, March 2017). Uncertainty also exists in manure application rates and areas. Manure generated by confined animal feeding operations is typically applied to adjacent forage cropland. Application rates are dictated by a variety of factors, including the capacity of the manure storage pond and the type of crop [54]. In preparing the nitrogen input data, we used average nitrogen application rates and a single schedule. We also assumed that a fixed portion of the nitrogen input was derived from manure for forage crops. While such simplifications do not strictly represent the actual distribution of nitrogen inputs across the San Joaquin River watershed, the general lack of explicit fertilizer and manure application data requires the use of a simplified, standard nitrogen management routine. Such approximations are frequently adopted in nutrient modeling [19,75,76,77,78]. Previous studies [14,79,80] derived estimates of nitrogen inputs in the San Joaquin River watershed from reports compiled by USGS in support of its National Water-Quality Program (NWQP) [2,81,82]. These data are too coarse for our modeling purposes, since they are annual county-level estimates. There is clearly a need for better estimates of fertilizer and manure application data with improved spatial and temporal resolution in intensive nutrient input areas, such as the San Joaquin River watershed.

4.3. Riverine Nitrate Exports, Aquatic Weed Infestation, and Future Climate Change Impact

The San Joaquin River plays a critical role in exporting nitrogen to the Delta, as 3135 tons of nitrate were exported annually through the Vernalis outlet (Figure 3). Loading from the San Joaquin River is particularly important during the summer months when the San Joaquin River watershed contributed 50% of the total nitrogen loading to the Delta [38]. Results in this study demonstrate the impacts of tile drainage on nitrogen loading in the San Joaquin River watershed. Therefore, conservation practices for reducing tile nitrate loss could be helpful in mitigating aquatic weed infestations in the Delta. For example, previous research shows that a controlled drainage system, which contains moveable flash boards to control the subsurface water table, is reliable for reducing nitrate loss through tile drains [83].

River nutrient inputs are one of the major drivers of the summertime growth of invasive aquatic vegetation, especially floating aquatic vegetation like water hyacinth, which obtains nutrients primarily from the water column. In general, increasing nitrogen levels lead to increases in biomass and density of aquatic macrophytes [3,84]. However, the relative contribution of nutrients in facilitating the growth and spread of invasive aquatic vegetation is not well understood, especially in the context of the complex Delta ecosystem [3]. Ongoing efforts to understand the response mechanisms of invasive aquatic weeds to nutrient inputs include the development of aquatic weed growth models [7,85]. Responses to nitrogen, as well as other limiting factors, such as temperature and light, are explicitly defined based on data collected in lab and field experiments. Simulated nitrogen loadings and concentrations from the outlet of the San Joaquin river basin could support aquatic weed growth models with important water quality parameters, such as nitrate loading to the Delta.

Climate change in the study region indicates a drier tendency with higher temperature and lower precipitation [86], which reduces snowpack in the eastern mountain region, and hence reduces the water storage in the state and Federal reservoirs [87]. Brekke, et al. [88] evaluated the potential climate change impacts on regional water resources in the San Joaquin river basin. Streamflow shifts in both the amount and timing are predicted by the CalSim model [89]. A generalized tendency of decreased reservoir inflows with earlier peaks in the eastern San Joaquin region is expected, indicating substantial impacts on the operation and planning of water resources, especially irrigated agriculture with heavy reliance on released reservoir water. This tendency will produce adverse effects on water quality, especially water salinity and nitrate concentration. For example, reduced water deliveries to the west agricultural regions from reservoirs encourage water reuse and groundwater exploitation, resulting in tailwater and tile drainage with higher concentrations of salts and nitrates conveyed to the San Joaquin River, eventually draining to the Bay-Delta Estuary.

Historical observation records already indicate upward trends of nitrate concentrations in the San Joaquin river basin. For example, nitrate concentration has increased from 1 mg/L to around 1.4 mg/L from 1975 to 2004 near Vernalis. For a specific monitoring site (Mud Slough) at western agricultural region, the increasing trend is as high as 493% [38]. Drier and warmer climates in the future has potentials to accelerate this increasing trend in San Joaquin water ways, facilitating the rapid growth of water hyacinth, especially under warmer conditions [90]. Previous studies indicated that higher density and rapid growth rates (10 kg/m² in 50 days) are predicted for water hyacinth in a eutrophic environment (dissolved nitrogen more than 1 mg/L) under warm conditions (30 °C) [91].

Additional modeling efforts are still required to characterize the fate and transport of nitrogen within the Delta [7,92]. In the future, coupling watershed and Delta models with aquatic weed growth models would allow us to explore the impact of agricultural nitrogen loading on the downstream aquatic environment. The coupled model could be used to assess alternative agricultural management practices and climate scenarios (e.g., drought) [6,93,94,95,96], leading to proactive management strategies for aquatic weed control.

5. Conclusions

In this study, we used the SWAT model to investigate nitrogen exports from the San Joaquin River watershed to Bay-Delta, California. We evaluated the tile drain routines of SWAT and found that the more physically based Hooghoudt and Kirkham tile drain equations, and the inclusion of a drainage coefficient, better represent the measured nitrate loads. Simulations suggested that 40% of nitrate loading to the San Joaquin River came from tile nitrate, consistent with previous studies, which showed that tile drains are important for nitrate delivery in the San Joaquin River watershed. SWAT tracked the monthly riverine nitrate dynamics reasonably well, given the uncertainty in the nutrient management data. More accurate and consistent data on nitrogen use at finer temporal and spatial resolutions are clearly needed. The San Joaquin River plays an important role in supplying nitrogen to the Bay-Delta by exporting 3135 tons of nitrate-nitrogen annually, which is an important nutrient source for aquatic weed growth. Future climate change indicated an increasing trend of nitrate concentration exported from San Joaquin river basin, which facilitates the rapid growth of water hyacinth under warmer conditions, leading to a more severe infestation situation in the Bay-Delta waterways. Our modeling results could be further used to drive downstream aquatic weed growth models to assess the relative importance of nutrient loading on aquatic weed growth, potentially yielding insights into sustainable management practices for improved weed control.