Simulation of Flow and Agricultural Non-Point Source Pollutant Transport in a Tibetan Plateau Irrigation District

Flow and transport processes in soil and rock play a critical role in agricultural non-point source pollution (ANPS) loads. In this study, we investigated the ANPS load discharged into rivers from an irrigation district in the Tibetan Plateau and simulated ANPS load using a distributed model. Experiments were conducted for two years to measure soil water content and nitrogen concentrations in soil and the quality and quantity of subsurface lateral flow in the rock and at the drainage canal outlet during the highland barley growing period. A distributed model, in which the subsurface lateral flow in the rock was described using a stepwise method, was developed to simulate flow and ammonium nitrogen (NH4-N) and nitrate nitrogen (NO3-N) transport processes. Sobol’s method was used to evaluate the sensitivity of simulated flow and transport processes to the model inputs. The results showed that with a 21.2% increase of rainfall and irrigation in the highland barley growing period, the average NH4-N and NO3-N concentrations in the soil layer decreased by 10.8% and 14.3%, respectively, due to increased deep seepage. Deep seepage of rainfall water accounted for 0–52.4% of total rainfall, whereas deep seepage of irrigation water accounted for 36.6–45.3% of total irrigation. NH4-N and NO3-N discharged into the drainage canal represented 19.9–30.4% and 19.4–26.7% of the deep seepage, respectively. The mean Nash–Sutcliffe coefficient value, which was close to 0.8, and the lowest values of root mean square errors, the fraction bias, and the fractional gross error indicated that the simulated flow rates and nitrogen concentrations using the proposed method were very accurate. The Sobol’s sensitivity analysis results demonstrated that subsurface lateral flow had the most important first-order and total-order effect on the simulated flow and NH4-N and NO3-N concentrations at the surface drainage outlet.


Introduction
The use of fertilizers in irrigation districts is considered as one of the most common sources of environmental pollution [1,2]. With low fertilizer recovery efficiency of nitrogen in the field, a major portion of applied fertilizers is lost due to various processes [3].
To simulate the hydrological processes and transport processes of agricultural non-point source (ANPS) pollutants, several models have been developed, such as the two-dimensional upland erosion model CASC2D (Colorado State University, Fort Collins, CO, USA) [4], dynamic watershed simulation model [5], and soil and water assessment tool (SWAT) model (United States Department of Agriculture, Washington, DC, USA) [6,7]. The modules in these models have been modified to reflect the impacts of irrigation practices and agricultural managements on migration processes of water and fertilizer in irrigation zones. For instance, the soil water and groundwater evaporation module was improved in the SWAT model to reflect the impact of crops on the hydrological cycle [8]. Additionally, the terrain elevation was corrected in the spatial grid unit to describe the irrigation drainage system [9,10].
However, hydrological processes in agricultural irrigation areas are largely influenced by manual irrigation and drainage processes. Therefore, the underlying surface conditions and runoff generation and confluence processes in irrigation zones differ significantly with those in natural watersheds [11][12][13][14]. Additionally, frequent water exchange that exists both inside and outside the irrigation region leads to the hydrological and transport processes being different from those simulated by the abovementioned models [15][16][17][18]. More importantly, because watershed hydrological and ANPS models must consider the macroscopic flow properties of the watershed, these models often describe the microscopic soil hydrological processes using relatively simple methods. In agricultural irrigation districts, however, the process of pollutant discharge into a river is mainly influenced by soil hydrological and pollutant transport and transformation processes [19][20][21]. As a result, the successful application of these models in complex irrigation and drainage zones remains limited.
An important issue when applying distributed models to simulate ANPS loads (pollution discharged into river) in Tibetan Plateau irrigation districts is the physical basis of the subsurface flow and transport and transformation processes. For example, in the Nyingchi prefecture of Tibet, irrigation districts are often located in valley terraces, where the soil layer is commonly as thin as 40-60 cm and underlain by rock comprised of loose rubble and sand gravels. Because the soil layer is relatively thin and highly permeable, it is difficult to attain direct surface runoff into rivers during the periods of rainfall and irrigation [22]. ANPS pollution is discharged from the soil into the river through various processes, including deep seepage from soil to rock, subsurface lateral flow in the rock, groundwater outflow, and surface drainage, as well as pollution conversion and transformation processes. At present, little is known about the ANPS pollution discharge process from soil into the surface drainage canal, especially the transport and transformation processes in the rock.
A number of parameters and physical variables have been used to describe ANPS pollution transport and transformation in different flow processes from the soil to the drainage canal outlet. Sobol's method is a global sensitivity analysis method that determines the impacts of each parameter and its interactions with other parameters on the model output [23][24][25]. Previous research has demonstrated the benefits of Sobol's method for identifying different flow and transport processes in models [24,26]. It is essential to investigating the crucial information on the use and meaning of the model parameters.
The objectives of this study are: (1) to investigate the effects of subsurface lateral flow and ANPS pollution transport and transformation processes on the ANPS pollution load from an irrigation district in the Tibetan Plateau, and (2) to simulate ANPS pollution load using a distributed model including a detailed description of flow and pollution transport and transformation processes in the rock. The study was conducted in the Danniang irrigation district for the major ANPS pollutants of ammonium nitrogen (NH 4 + -N) and nitrate nitrogen (NO 3 − -N).

Study Area
The experiments were conducted in the Danniang irrigation district of Miling (in Nyingchi prefecture, Tibet) from May to September in 2014 and 2015. The total and actual irrigated area were 12.8 and 4.26 km 2 , respectively. Spring highland barley was the main crop planted in this district. From 1953 to 2014, the average annual temperature was 8.7 • C and the extreme maximum and minimum temperatures were 32.4 • C and −15.3 • C, respectively. The annual average relative humidity was 64.2%. The average annual precipitation was 664.5 mm, which was predominantly concentrated during the period from May to September (81.4% of the annual precipitation). The average annual pan evaporation was 1734.1 mm. The annual average wind speed was 1.68 m/s, and the annual average hours of sunshine were 2064.6 h.
The irrigation and drainage canals in the irrigation district are shown in Figure 1. The main irrigation channel was arranged along the contour lines of the piedmont. Irrigation subdistricts 1, 2, 3, 7, and 8 were located in the long and narrow region in the piedmont and irrigated with water directly from the main irrigation channel. Subdistricts 4-6 were irrigated with water from the branch irrigation channels. The average elevation difference between the main irrigation channel and drainage canal was 3.6 m.

Monitoring Hydrological and ANPS Pollution Transport and Transformation Processes in the Soil and Rock
Soil samples were taken in different layers (each 10 cm thick) from 0-50 cm depth during the growing period of highland barley. In each sampling event, the samples were obtained at 12 locations to determine soil water content and concentrations of the major ANPS nitrogen pollutants (i.e., NH4 + -N and NO3 − -N). All soil samples were collected using a cutting ring with a volume of 100 cm 3 . Fourteen and fifteen sampling events were conducted during the growing period of highland barley in 2014 and 2015, respectively. The soil water content and NH4 + -N and The soil thickness was 50 cm and underlain with a loose rock layer comprised of coarse sand and rubble 2-40 cm in size. The soil was classified as dark brown soil according to the standards of the second national soil survey in China and sandy loam according to the standards of the US Department of Agriculture (UNSADO). The rock layer was more than 20 m thick. Samples were taken at 12 locations ( Figure 1) to determine soil physical parameters (Table 1). Fertilization was performed twice with a dose of 18,000.0 and 6000.0 kg urea/km 2 for the first and second irrigation, respectively.

Monitoring Hydrological and ANPS Pollution Transport and Transformation Processes in the Soil and Rock
Soil samples were taken in different layers (each 10 cm thick) from 0-50 cm depth during the growing period of highland barley. In each sampling event, the samples were obtained at 12 locations to determine soil water content and concentrations of the major ANPS nitrogen pollutants (i.e., NH 4 + -N and NO 3 − -N). All soil samples were collected using a cutting ring with a volume of 100 cm 3 . Fourteen and fifteen sampling events were conducted during the growing period of highland barley in 2014 and 2015, respectively. The soil water content and NH 4 + -N and NO 3 − -N concentrations in these samples were measured at the Farmland Water Conservancy Laboratory of the Tibet Agriculture and Animal Husbandry College. A subsurface lateral flow experimental site was constructed in the irrigation district ( Figure 1) to measure subsurface lateral flow and nitrogen transport and transformation processes in the rock ( Figure 2). Two 78-m long and 2-m wide experimental plots were set up. The plot length was identical to those of ridged fields in the irrigation district. The borders of the experiment plot were isolated with steel plates compressed into the ground to a depth of 20 cm. A 0.5-m wide protected area was set outside the steel plates, wherein steel plates were arranged in parallel to the experimental plots. Two subsurface lateral flow monitoring sections (A and B) were set at 0.5 m and 20.0 m, respectively, behind the ridge of plots. A trench (5.0 m long and 2.0 m wide) was excavated mechanically to a depth of 4.0 m. Subsurface lateral flow collection devices were set at a depth of 3 m below the soil. The collection device was compressed into the rock to a depth of 20 cm. A tipping-bucket automatic water recorder was introduced to measure subsurface lateral flow rate. The bucket had a volume of 50 mL and was automatically tipped after being filled. The time and frequency of tipping were automatically recorded, enabling automatic monitoring of the subsurface lateral flow rate. Meanwhile, water samples were collected to measure NH 4 + -N and NO 3 − -N concentrations. The area between sections A and B was covered with a rainproof film to prevent water infiltration into the soil. The amount of irrigation in the experimental site was identical to that in the field during the growing period of highland barley.
Water 2018, 10, x FOR PEER REVIEW 5 of 24 A subsurface lateral flow experimental site was constructed in the irrigation district ( Figure  1) to measure subsurface lateral flow and nitrogen transport and transformation processes in the rock ( Figure 2). Two 78-m long and 2-m wide experimental plots were set up. The plot length was identical to those of ridged fields in the irrigation district. The borders of the experiment plot were isolated with steel plates compressed into the ground to a depth of 20 cm. A 0.5-m wide protected area was set outside the steel plates, wherein steel plates were arranged in parallel to the experimental plots. Two subsurface lateral flow monitoring sections (A and B) were set at 0.5 m and 20.0 m, respectively, behind the ridge of plots. A trench (5.0 m long and 2.0 m wide) was excavated mechanically to a depth of 4.0 m. Subsurface lateral flow collection devices were set at a depth of 3 m below the soil. The collection device was compressed into the rock to a depth of 20 cm. A tipping-bucket automatic water recorder was introduced to measure subsurface lateral flow rate. The bucket had a volume of 50 mL and was automatically tipped after being filled. The time and frequency of tipping were automatically recorded, enabling automatic monitoring of the subsurface lateral flow rate. Meanwhile, water samples were collected to measure NH4 + -N and NO3 − -N concentrations. The area between sections A and B was covered with a rainproof film to prevent water infiltration into the soil. The amount of irrigation in the experimental site was identical to that in the field during the growing period of highland barley. The kinetic coefficients of pollutant transformation were determined based on the mass balance method [27,28]. Undisturbed rock samples were collected using a stainless-steel sleeve (10 cm long and 5 cm in diameter). Ion exchange resins were placed at the bottom and top of the steel sleeve to eliminate the effects of the in-and out-fluxes of pollutants at the border on the The kinetic coefficients of pollutant transformation were determined based on the mass balance method [27,28]. Undisturbed rock samples were collected using a stainless-steel sleeve (10 cm long and 5 cm in diameter). Ion exchange resins were placed at the bottom and top of the steel sleeve to eliminate the effects of the in-and out-fluxes of pollutants at the border on the measurement of the transformation process. The steel sleeve was returned to the sampling point with backfilling. The samples were retrieved 7 days later. The NH 4 + -N and NO 3 − -N concentrations were measured on the sampling day and 7 days later. The mass difference was taken as the amount of pollutant transformation. Samples were taken at three positions in the soil layer (10-40 cm depth) and loose rock layer (70-90 cm depth). The measurements were made twice, three times, and once, respectively, during the sowing-tilling, jointing-heading, and flowering-filling stages of the highland barley. A monitoring section, i.e., section O in Figure 1, was set at the outlet of the surface drainage canal in the irrigation district. The drainage process was monitored automatically using a Doppler flowmeter, and water samples were taken to measure NH 4 + -N and NO 3 − -N concentrations in the surface drainage. Meteorological parameters, including radiation, temperature, humidity, wind speed, and air pressure, were measured at a meteorological station within the study area.

Simulation of Hydrological and ANPS Pollution Transport Processes in the Plateau Irrigation District
Water flow and transport processes were simulated independently in the irrigation subdistricts with different parameters related to the irrigation events and size of the irrigation subdistrict. The flow and pollution fluxes at the outlet of the drainage canal represent the confluence and superposition processes of flow discharged from the rock, including from subsurface lateral flow and groundwater.
Subsurface lateral flow in the rock is illustrated in Figure 3. Monitoring data from 2014 and 2015 showed that no significant surface runoff occurred during irrigation. A stepwise method was used to describe subsurface lateral flow processes in the rock. In the interval from i to i + 1, the deep seepage flux L i that enters the rock from the soil can be expressed as where I is the amount of irrigation (or rainfall) and ∆W i is the variation of soil water content before and after irrigation (or rainfall).
Water 2018, 10, x FOR PEER REVIEW 6 of 24 measurement of the transformation process. The steel sleeve was returned to the sampling point with backfilling. The samples were retrieved 7 days later. The NH4 + -N and NO3 − -N concentrations were measured on the sampling day and 7 days later. The mass difference was taken as the amount of pollutant transformation. Samples were taken at three positions in the soil layer (10-40 cm depth) and loose rock layer (70-90 cm depth). The measurements were made twice, three times, and once, respectively, during the sowing-tilling, jointing-heading, and flowering-filling stages of the highland barley. A monitoring section, i.e., section O in Figure 1, was set at the outlet of the surface drainage canal in the irrigation district. The drainage process was monitored automatically using a Doppler flowmeter, and water samples were taken to measure NH4 + -N and NO3 − -N concentrations in the surface drainage.
Meteorological parameters, including radiation, temperature, humidity, wind speed, and air pressure, were measured at a meteorological station within the study area.

Simulation of Hydrological and ANPS Pollution Transport Processes in the Plateau Irrigation District
Water flow and transport processes were simulated independently in the irrigation subdistricts with different parameters related to the irrigation events and size of the irrigation subdistrict. The flow and pollution fluxes at the outlet of the drainage canal represent the confluence and superposition processes of flow discharged from the rock, including from subsurface lateral flow and groundwater.
Subsurface lateral flow in the rock is illustrated in Figure 3. Monitoring data from 2014 and 2015 showed that no significant surface runoff occurred during irrigation. A stepwise method was used to describe subsurface lateral flow processes in the rock. In the interval from i to i + 1, the deep seepage flux Li that enters the rock from the soil can be expressed as where I is the amount of irrigation (or rainfall) and i W Δ is the variation of soil water content before and after irrigation (or rainfall). The flux Li seeps into the rock from the soil and moves laterally toward to the surface drainage canal as well as vertically into the groundwater. The lateral inflow flux from soil seepage Li−1 at section i is denoted as Si (Figure 3b). In the interval from section i to section i + 1, only a partial flux of Si continuously moves to the surface drainage canal due to high arbitrariness of flow in the loose fractured rock.
A subsection-fitting method was used for the calculation of the subsurface lateral flow coefficient where t is the relative time (seepage time/irrigation duration) and α and β are parameters related to the rock structure. At section i + 1, the lateral flow flux Si decreases to Si,i+1. The flux at section i + 1 also includes lateral flow from soil seepage Li, denoted as flux Si+1. At section i, the lateral flux to the surface drainage canal is the superposition of all fluxes and can be expressed as In the rock, the impact of the diffusion process on pollutant migration can be ignored in the case of a long flow path, and the convection (i.e., the product of water flux and pollution concentration) process can be used to describe the pollutant transport flux. The mean values of the measured soil water content and NH4 + -N and NO3 − -N concentrations in the soil profile were used as inputs in the subsurface lateral flow and transport simulation.
The amount of seepage water from the rock to the aquifer in the section from i to i + 1 is The flux L i seeps into the rock from the soil and moves laterally toward to the surface drainage canal as well as vertically into the groundwater. The lateral inflow flux from soil seepage L i − 1 at section i is denoted as S i (Figure 3b). In the interval from section i to section i + 1, only a partial flux of S i continuously moves to the surface drainage canal due to high arbitrariness of flow in the loose fractured rock.
A subsection-fitting method was used for the calculation of the subsurface lateral flow coefficient µ(G) in the rock [29]: where t is the relative time (seepage time/irrigation duration) and α and β are parameters related to the rock structure. At section i + 1, the lateral flow flux S i decreases to S i,i+1. The flux at section i + 1 also includes lateral flow from soil seepage L i , denoted as flux S i + 1 . At section i, the lateral flux to the surface drainage canal is the superposition of all fluxes and can be expressed as In the rock, the impact of the diffusion process on pollutant migration can be ignored in the case of a long flow path, and the convection (i.e., the product of water flux and pollution concentration) process can be used to describe the pollutant transport flux. The mean values of the measured soil water content and NH 4 + -N and NO 3 − -N concentrations in the soil profile were used as inputs in the subsurface lateral flow and transport simulation. The amount of seepage water from the rock to the aquifer in the section from i to i + 1 is Assuming that variations in groundwater flow are linearly related to the rate of change in water table height, we obtain dQ gwi,i+1 dt where Q gwi,i + 1 is the groundwater flow into the surface drainage canal, µ is the specific yield of the shallow aquifer (cm 3 /cm 3 ), and L g is the distance from the ridge of the groundwater to the surface drainage canal. Water flow and pollutant transport processes in the drainage canal can be characterized by the following Equations: ∂A ∂t where A is the area of the drainage canal cross section; Q is the surface drainage flow rate; q is the inflow during the calculation period (as determined using Equations (4) and (6)); h is the depth of the drainage channel; n is the Manning roughness coefficient; R and i 0 are the hydraulic radius and longitudinal slope of the drainage channel, respectively; and i f is the friction slope. Due to the prevention of field ridges and the high soil permeability in the irrigation zone, irrigation and rainfall did not typically produce runoff. The transformation processes in the rock and drainage water were described using the first-order kinetics described by Equation [9]: where c t and c t+1 are the concentrations of pollutants (NH 4 + -N or NO 3 − -N) at times t and t + 1 (mg/L).
The k is a lump factor to quantify pollutant transformation caused by various physical, chemical, and biological processes in the rock or in the drainage water (day −1 ).

Calibration of Model Parameters
Monitored flow rates and NH 4 + -N and NO 3 − -N concentration at sections A and B in the rock and at the outlet of the drainage canal during the highland barley growing period in 2014 were used to calibrate model parameters. The objective function was defined as where m q is the number of monitoring variables (flow rate, NH 4 + -N concentration, and NO 3 − -N concentration), n qj is the monitoring number of variable j, g j * (x, t i ) is the monitored value of variable j at time t i in the monitoring position, and g j ( Additionally, v j is the weight of variable j, which is used to reduce the magnitude effect of different variables and different monitoring events on the calculated results.
where σ j is the standard deviation of the monitored value for variable j.
The simulation results of the NH 4 + -N and NO 3 --N concentrations were evaluated using four indicators: the Nash-Sutcliffe coefficient (NSE) [30], the relative root mean square error (rRMSE), the fraction bias (F B ), and the fractional gross error (F E ) [31]; NSE, F B , and F E are calculated as follows: where O t and P t are the monitored and calculated values at time t, respectively, O is the mean of the monitored value, and n is the number of observation points. The ideal value of NSE was 1. F B and F E were used to describe the systematic error and total deviation between the observed and simulated values, respectively.

Sobol's Sensitivity Analysis
The effects of input parameters on model performance are evaluated using Sobol's method. The model can be represented in the following functional form: where Y is the model output and X = (x 1 , x 2 , . . . x n ) is the parameter and physical variable set. For all of parameters and physical variables used in the model, two low-discrepancy random sequences are generated using the random sampling method [32].
where m is the sampling repetition and n is the number of parameters. A new matrix was formed by replacing the parameters in the ith column in matrix B (i.e., sample values of parameter x i ) with the ith column data in matrix A. Finally, matrices A and B were stacked into one "combined" matrix, C, which has n + 2 rows and n columns.
The first-order sensitivity index, S i , and the total-order sensitivity index, S Ti , were calculated as follows: where y

(k)
A , y (k) B , and y (k) C are the model outputs with parameters in matrices A, B, and C, respectively. The variance of the output with the parameters in matrix A is calculated as The first-order index, S i , is a measure of the variance contribution of the individual parameter x i to the total model variance. The difference between the first-order sensitivity index and the total-order sensitivity index can be regarded as a measure of the interaction between parameter x i and the other parameters.

Characterization of Flow and NH 4 + -N and NO 3 − -N Transport
The total rainfall and irrigation were 712 mm and 863 mm for the growing period of highland barley in 2014 and 2015, respectively. Figure 4a

Simulation of ANPS Pollutant Transport and Transformation Processes
The calibrated parameters are listed in Table 3. Figure 6a-c compares the simulated and measured flow rate and NH 4 + -N and NO 3 − -N concentrations at the outlet (Section O) of the surface drainage canal. The accuracy of calculated flow and transport processes from soil to the outlet of the drainage canal was quantified with the indexes of NSE, rRMSE, F E , and F B listed in Table 4. The NSE value exceeded 0.6. The F B remained less than 15% throughout the growing period of highland barley, demonstrating no systematic error between the simulated and measured results. F B and F E indexes of the simulated values were within the range of ±0.15/0.3 (F B /F E ), indicating high accuracy of the simulation [31]. All the results of the Nash-Sutcliffe coefficient, rRMSE, and F E demonstrate that the model can effectively describe the transport and transformation processes of pollutants in the rock and estimate the pollution load accurately in this plateau irrigation area. Table 5

Sensitivity Evaluation of Model Parameters
The first-and total-order sensitivity indices of 12 parameters are shown in Figure 7. The x-axis represents the parameter number listed in Table 3, and the y-axis represents the first-and total-order sensitivity indices (black and gray bars, respectively). It is noted that the sum of all first-order indices is less than 1, which means that the model is nonadditive, as expected. The first-and total-order sensitivity indices of the parameters related to subsurface lateral flow in the rock (α and β) were significantly higher than the other parameters. The first-and total-order sensitivity indices of the parameters related to groundwater discharge into surface drainage were the lowest of all parameters. These results were attributed to the lower discharge of water, NH 4 + -N, and NO 3 − -N from groundwater than from subsurface lateral flow.
The parameters related to subsurface lateral flow exhibited the maximum sensitivity. The value of S i for soil NH 4 + -N concentration was higher than the rates of deep seepage from the soil into the rock.
The difference in S i and S Ti values for the NH 4 + -N transformation coefficient in lateral transport and surface drainage transport was significantly higher than the differences between other parameters. The results also showed that 65.2% of the variations in the simulated NH 4 + -N concentrations at the surface drainage outlet were caused by variations of α and β in Equation (3)

Sensitivity Evaluation of Model Parameters
The first-and total-order sensitivity indices of 12 parameters are shown in Figure 7. The xaxis represents the parameter number listed in Table 3, and the y-axis represents the first-and total-order sensitivity indices (black and gray bars, respectively). It is noted that the sum of all first-order indices is less than 1, which means that the model is nonadditive, as expected. The first-and total-order sensitivity indices of the parameters related to subsurface lateral flow in the rock (α and β) were significantly higher than the other parameters. The first-and total-order sensitivity indices of the parameters related to groundwater discharge into surface drainage were the lowest of all parameters. These results were attributed to the lower discharge of water, NH4 + -N, and NO3 − -N from groundwater than from subsurface lateral flow.
The parameters related to subsurface lateral flow exhibited the maximum sensitivity. The value of Si for soil NH4 + -N concentration was higher than the rates of deep seepage from the soil into the rock. The difference in Si and STi values for the NH4 + -N transformation coefficient in lateral transport and surface drainage transport was significantly higher than the differences between other parameters. The results also showed that 65.2% of the variations in the simulated NH4 + -N concentrations at the surface drainage outlet were caused by variations of α and β in Equation (3)

Discussion
In the Danniang irrigation district, the experiment results indicated that NH4 + -N and NO3 − -N concentrations are directly affected by the subsurface lateral flow process. The NH4 + -N and NO3 − -N leaching per unit rainfall and irrigation were 4.25 and 18.45 kg km −2 mm −1 , and the ratio of the total N leaching to the fertilizer N application was 0.34. These results were substantially different with those of Ju and Zhang [33], who reported that the N leaching in hilly barley crop land was in the range of 1.05 × 10 −2 to 1.18 kg km −2 mm −1 , and the ratio of the N leaching to the fertilizer N application was in the range of 0.044 to 0.249. The thin soil layer, irregularities in rainfall distribution, and heavy irrigation attributed to the higher N leaching in the Danniang irrigation

Discussion
In the Danniang irrigation district, the experiment results indicated that NH 4 + -N and NO 3 − -N concentrations are directly affected by the subsurface lateral flow process. The NH 4 + -N and NO 3 − -N leaching per unit rainfall and irrigation were 4.25 and 18.45 kg km −2 mm −1 , and the ratio of the total N leaching to the fertilizer N application was 0.34. These results were substantially different with those of Ju and Zhang [33], who reported that the N leaching in hilly barley crop land was in the range of 1.05 × 10 −2 to 1.18 kg km −2 mm −1 , and the ratio of the N leaching to the fertilizer N application was in the range of 0.044 to 0.249. The thin soil layer, irregularities in rainfall distribution, and heavy irrigation attributed to the higher N leaching in the Danniang irrigation district, especially the NH 4 + -N leaching, was mainly attributed to heavy irrigation because fertilizer N was applied during the first and second irrigation and soil NH 4 + -N concentrations were maintained at high levels after fertilization. In total, 45.8 kg km −2 of NH 4 + -N, or 84.7% of total NH 4 + -N leaching, was caused by the irrigation. Given that the soil thickness is sufficient, the spatial and temporal distribution of NO 3 − -N in the soil profile may be significantly changed by fertilization, irrigation, and precipitation events [34,35]. NH 4 + -N can be converted to NO 3 − -N rapidly under favorable water and heat conditions. As a result, less NH 4 + -N leaching occurs [36]. In the studied irrigation district, 28.7 and 72.3% of the total NO 3 − -N leaching was caused by irrigation and rainfall, respectively, due to the thin and highly permeable soil. Although the mean value of the flow rates and NH 4 + -N and NO 3 − -N concentrations at the outlet of the drainage canal were accurately simulated, the coefficient of variation (the ratio of standard deviation to mean) of flow rates and NH 4 + -N and NO 3 − -N concentrations were overestimated. Due to a 21.2% increase in total rainfall and irrigation in 2015, the monitored coefficients of variation of the flow rate and NH 4 + -N and NO 3 − -N concentrations at the outlet of the drainage canal increased by 15.7, 18.9, and 8.2%, respectively, while the simulated increase of the coefficients of variation were 24.2, 28.4, and 10.8%, respectively. The deviations between the measured and the simulated flow rates and NH 4 + -N and NO 3 − -N concentrations might be attributed to using the lumping factor to characterize the various flow and transformation processes. The storage effect and the nonlinear flow processes in the rock also affected the accuracy of the calculated concentrations. As shown in Figure 6, the flow rates simulated with the calibrated parameters were lower than the values measured in July 2015. The mean values of the measured and simulated flow rates were 0.0338 and 0.0402 m 3 /s, respectively, at the outlet of the drainage canal. Soil seepage into the rock was significantly reduced because rainfall was only 14.5 mm in June 2015 (Figure 4), and part of the seepage was stored in fractures and pores in the rock instead of being discharged into the surface drainage canal. As a result, the model underestimated flow rates and overestimated NH 4 + -N and NO 3 − -N concentrations by 4.71 and 5.64%, respectively, for this period. The parameters of the dominant flow processes, i.e., the subsurface lateral flow process, have the most first-and total-order sensitivities. This result is similar to those found at natural basins, in which the runoff is the dominant flow process and related parameters have the most first-and total-order sensitivities [25,27,37,38]. The differences between the first-and total-order sensitivity indices of the parameters used to describe subsurface lateral flow process were significantly less than those for the parameters used in simulating the dominant flow process in the watershed, since the range of flow velocity in the rock was limited and fewer factors caused variation in the subsurface lateral flow. There were two parameters, i.e., α and β, in Equation (3) that were sensitive to the both subsurface lateral flow and nitrogen (NH 4 + -N and NO 3 − -N) transport and transformation processes. The Sobol's sensitivity indices of α were significantly greater than the indices of β for the flow process, whereas the indices of α decreased and the indices of β increased for the nitrogen transport and transformation processes. The result showed that the peak flow rate after rainfall or irrigation has more influence on the estimation of flow process, and the description of the nonlinear subsurface flow process was more influential on the simulated output of the NH 4 + -N and NO 3 − -N concentrations.

Conclusions
Field experiments were conducted to investigate the hydrology and ANPS transport and transformation processes in the Danniang irrigation region, Tibet, during the highland barley growing periods in 2014 and 2015. A distributed model was constructed to simulate flow and ANPS (NH 4 + -N and NO 3 − -N) transport and transformation from the soil to the outlet of the surface drainage canal.
Water flow and transport processes were simulated in irrigation subdistricts with different parameters related to the irrigation events and the size of the irrigation subdistricts. A stepwise method was used to describe subsurface lateral flow in the rock. The flow and pollution fluxes at the outlet of the drainage canal represented the confluence and superposition processes of flow discharged from the rock, including subsurface lateral flow and groundwater.
The simulated values of flow rate, NH 4 + -N concentration, and NO 3 − -N concentration exhibited systematic deviations of less than 15%. Across the growing period of highland barley, the mean values of the Nash-Sutcliffe coefficient, rRMSE, the fraction bias, and the fractional gross error between the simulated and observed flow rates and NH 4 + -N and NO 3 − -N concentrations at the outlet of the surface drainage canal were 0.712, 7.01%, 0.04, and 0.255, respectively. These indicate that the proposed method can effectively simulate the hydrological and ANPS pollutant migration processes in this plateau irrigation district. As a result of the 21.2% increase in the rainfall and irrigation amount during the highland barley growing period in 2015, the average NH 4 + -N and NO 3 − -N concentrations in the soil layer decreased by 10.8 and 14.3%, respectively, due to increased deep seepage. Deep seepage due to rainfall accounted for 0-52.4% of all rainfall, whereas deep seepage due to irrigation accounted for 36.6-45.3% of all irrigation. NH 4 + -N and NO 3 − -N discharge into the drainage channel represented 19.9-30.4% and 19.4-26.7% of the deep seepage, respectively. Parameters in the subsurface lateral flow simulation were shown to have the most important first-order and total-order effects on the simulated flow and NH 4 + -N and NO 3 − -N concentrations at the outlet of the surface drainage canal.