Recharge and Nitrate Transport Through the Deep Vadose Zone of the Loess Plateau: A Regional‐Scale Model Investigation

The Loess Plateau of China (LPC) is suffering from the impacts of overexploitation of groundwater, leading to a decline in the regional water table. This is compounded by relatively recent large‐scale transformation of land use, which may impact groundwater recharge and threaten future water quality from intensive nitrogen fertilizer application. Understanding the regional unsaturated water flow and nitrate travel time in the deep vadose zone is crucial for the sustainable management of the LPC groundwater system. We develop here a regional‐scale model that exploits an extensive database of soil properties derived from recent intensive soil profile sampling over the LPC, together with climate, vadose zone thickness, and land use data. The model solves the Richards equation and a conservative solute transport model using a multiple 1D column approach. A comparison between model simulations and observations at five sites reveals good model performance. Application over the LPC indicates that areas with dense irrigated agricultural land use contribute significantly to groundwater recharge. By using land use cover information from 1975 to 2008, simulations reveal that the transition from agriculture to forest and grassland accounts for a small reduction of the total annual recharge (6.1%). Modeled nitrate travel times through the deep unsaturated zone range between decades and centuries; however, recent loading of anthropogenic nitrate is expected to reach the aquifer in the near future under irrigated areas or where the vadose zone is relatively thin. The model simulations provide valuable information for assessing future vulnerability of groundwater resources over a regional scale.


Introduction
Water resources management is a great challenge in China due to its large population and rapid economic growth (Currell et al., 2012). Moreover, Qiu (2010) indicated that China faces a groundwater crisis due to overexploitation and water quality degradation. In the north of China, where surface water availability is limited, the groundwater systems are subject to even higher pressure to meet water demand (OECD, 2007). A specific case is the groundwater system under the Loess Plateau of China (LPC), located in arid and semiarid regions of northern China. Water resources in the LPC have been overexploited, and the regional water table declines rapidly (Huang & Pang, 2011;C. Li et al., 2014). Furthermore, as a result of high soil erosion rates in this region, soil conservation measures were applied (Zhang et al., 2008). Land-use changes have been linked with the reduction of runoff, depletion of soil water, and decrease of groundwater recharge fluxes (e.g., Gates et al., 2011;Huang & Pang, 2011;Huang, Pang, & Edmunds, 2013;Jia et al., 2015;Wang et al., 2012Wang et al., , 2013Zhang et al., 2008). Another factor that is likely to impact the water quality of the already vulnerable groundwater system of the LPC is the excessive application of fertilizers (Huang, Pang, & Yuan, 2013). Zhang et al. (1996) reported that since the 1980s, nitrogen fertilizer inputs in China have exceeded plant use, thus establishing the conditions for enhanced nitrate leaching.
Groundwater recharge and nitrate transport and travel time in the unsaturated zone are controlled by many factors (e.g., soil, climate, land use, vadose zone thickness, vegetation, and water input) and therefore are spatially and temporally distributed across a landscape. An accurate representation of these processes is necessary for water resources management (Narula & Gosain, 2013). Moreover, the groundwater response to an input of water or solute at the ground surface is critically dependent on the travel time through the unsaturated zone (Jackson et al., 2008;Rossman et al., 2014;Sousa et al., 2013;Vero et al., 2014). Travel time (or lag-time) estimations are required to get a better assessment of the timescale for expected changes in groundwater quality and quantity due to, for example, changes in climate agricultural management practices (Rossman et al., 2014;Vero et al., 2014). Earlier studies suggested that estimates of water and nitrate fluxes from deep unsaturated zone observations could better recognize the sites that are "hot spots" of high recharge rates and vulnerable to nitrogen (N) leaching (e.g., Green et al., 2008). It is becoming more common to utilize so-called "physically based" models coupled with ArcGIS (Environmental Systems Research Institute, 2016) to describe regional groundwater recharge, nitrate and unsaturated vertical travel time distributions across multiple agricultural settings, soil types, and climate conditions (Assefa & Woodbury, 2013;Narula & Gosain, 2013;Rossman et al., 2014).
Application of physically based unsaturated flow and transport models requires soil hydraulic functions, such as the retention curve and unsaturated hydraulic conductivity properties; however, direct measurements of these parameters are often rare or obtainable locally. Consequently, there has been much use made of soil hydraulic functions predicted by pedotransfer functions (PTFs), which convert relatively easily obtainable soil properties (e.g., particle size distribution [PSD]) to soil hydraulic functions. There are inevitable uncertainties associated with soil hydraulic functions predictions by PTFs, which therefore impact on model predictions of states (e.g., water content; Faust et al., 2006;Loosvelt et al., 2011;Sinowski et al., 1997;Wang et al., 2009). However, various studies have demonstrated the reliability of the Rosetta for estimation of soil parameters estimation after evaluating model performance with independent measurements (Assefa & Woodbury, 2013;Jana & Mohanty, 2012). In this study we exploit the PTF approach in a study of unsaturated zone flow and transport since it allows us to examine processes in the LPC in a regional context. Our aims are to develop a regional-scale model of flow and solute transport in the vadose zone of the LPC and use this to examine impacts of land use change and assess likely solute travel times to groundwater.
Recently, an intensive soil profile sampling was conducted across the LPC Zhao, Shao, Jia, Nasir, et al., 2016). The analysis of the spatial variability of soil PSD suggests that clay, silt, and sand contents vary moderately in space and that soil texture changes with depth are insignificant. The availability of a detailed continuous map of soil properties has enabled the establishment of a regional-scale hydrological model at a 3-km resolution. The Richards equation and the advection-dispersion equation are implemented in a multiple 1D column model using interpolated climate at a daily resolution, gridded soil parameters, gridded vadose zone thickness, and land use data distributed over the LPC. Thus, the regional model combines all the factors controlling recharge and nitrate fluxes to groundwater. The regional model simulations are compared with local-scale measurements. The model is then used to investigate locations and dominant factors controlling the spatial variability of groundwater recharge and the travel time of nitrate in the unsaturated zone over the LPC.

Study Area
The LPC region covers a total area of 0.64 × 10 6 km 2 and is subject to a semiarid to subhumid climate ( Figure 1). According to Jia et al. (2015), the subregion of the LPC with continuous loess soil has an area of 0.43 × 10 6 km 2 (Figure 1). The principal unconfined aquifer in the region is embedded within the bounds of continuous loess soil, in Shaanxi, southern Ningxia, western Shanxi, and south-eastern Gansu provinces (Jakeman et al., 2016) (Figure 1). Loess is a Quaternary windborne deposit characterized by predominance of silt-sized particles (Lin & Wei, 2006). In the LPC, various developed geomorphologic units are present with similar grain-size distributions. The major particle fraction of loess is 0.05-0.01 mm, which accounts for about 45-60%; the remainder includes a <0.005-mm fraction, accounting for about 15-25% (Lin & Wei, 2006). The soil typically lacks grains greater than 0.25 mm, consistent with aeolian deposition (Lin & Wei, 2006). The Quaternary loess is composed of lower Pleistocene deposits (Wucheng Loess, Q1), middle Pleistocene deposits (Lishi Loess, Q2), and upper Pleistocene deposits (Malan Loess, Q3) (Huang & Pang, 2011;C. Li et al., 2014). The hard and compacted Wucheng Loess has a low permeability and is usually considered as an aquitard. The Lishi Loess has a thickness of 120-150 m, is unconsolidated and has a relatively high porosity, and is considered to be a good aquifer. The Malan Loess has a thickness of 10-15 m and is distributed above the Lishi Loess.

Soil Properties, Rosetta, and Parameter Regionalization
We utilized soil data collected from 243 locations over an area of 43 × 10 5 km 2 in the LPC (Figure 2). For further information concerning the soil sampling strategy, data collection, and the laboratory 10.1029/2017WR022190 Water Resources Research measurements, the reader is referred to  and Zhao, Shao, Jia, Nasir, et al. (2016). In summary, undisturbed soil cores were excavated at depths of 0-10, 10-20, and 20-40 cm to measure bulk density (BD), saturated soil water content, field capacity (FC), and saturated hydraulic conductivity (K s ; Figure 2). In addition, disturbed soil samples were collected at the same corresponding soil depths and from 40-60, 60-80, 80-100, 100-150, 150-200, 200-300, 300-400, and 400-500 cm for soil particle distribution analyses ( Figure 2). Note that in this study, the soil properties sampled at one of the depths (10-20 cm) were used as a representative deep unsaturated profile at the sampling location. It is debatable whether it is necessary to exploit all data (e.g., to average the sampled values) or select only one of the depths to obtain an optimal characterization of deep soil profile properties. However, the narrow PSD of the loess soil and the minor vertical variability indicate that no significant layering exists across a regional scale. Therefore, the use of uniform representative profiles in a regional-scale model is justified.
Soil retention curves and unsaturated hydraulic conductivity were estimated according to the van Genuchten-Mualem (VGM) model (Mualem, 1976;van Genuchten, 1980). Kriging interpolation was applied to regionalize the soil parameters over the study area. The saturated hydraulic conductivity (K s ) and saturated soil water content (θ s ) were interpolated directly from measurements. Since the shape parameters (α and n) and the residual soil water content were to be determined by the Rosetta PTF (Schaap et al., 2001), regionalization of the parameters was performed differently. Previous studies have suggested interpolating first the input soil properties and then applying the PTF computed values (Heuvelink & Pebesma, 1999;Sinowski et al., 1997). The different soil properties input to the PTF display a different correlation structure, which might be lost or corrupted when applying PTFs first. Therefore, the measured soil properties including BD, FC (water content at pressure head value of 330 cm), sand percentage, silt percentage, and clay percentage were interpolated, and subsequently, the shape and residual water content parameters were predicted by Rosetta.

Climate Data
Daily climate data for the time period 1955-2016 (22,646 days) were acquired from 56 stations within and near the Loess Plateau region (Figures 1 and 3; State Bureau of Meteorology, 2018; http://cdc.cma.gov.cn/ home.do). This climate database contains daily measurements of mean air temperature (°C), maximum temperature (°C), minimum temperature (°C), precipitation (mm), relative humidity (%), sunshine hours (hr), and wind speed (m/s). Using this database, daily reference evapotranspiration (ET ref ) was calculated according to the Penman-Monteith equation (Allen et al., 1998) for each meteorological station. The precipitation and ET ref mean annual and standard deviation values display temporal variability with no noticeable trend, while the mean air temperature shows minor inclination ( Figure 2). A long-term average of the annual precipitation and reference evaporation over the outcrops of the LPC aquifer is estimated at 437 ± 73 and 1,060 ± 48 mm, respectively. The error bars (standard deviations) indicate a relatively large spatial variability of the precipitation but smaller relative variability of ET ref .

Land Use Data
The dominant land uses over the LPC aquifer outcrops are agriculture, followed by grassland and forest. Because of the high rates of agricultural soil erosion in the LPC, a wide range of soil conservation measures were implemented (Zhang et al., 2008). Hence, these measures were

10.1029/2017WR022190
Water Resources Research expressed as land use changes at a regional scale ( Figure 4). Close examination of land use cover from the years of 1975 and 2008 indicates the reduction in agricultural land and increase in areas of grassland and forestry ( Figure 4). In 1975, 41% of the aquifer outcrop was covered with agricultural land, 39% was covered by grassland, and 15% was covered by forests. By 2008, the agricultural land uses decreased to 34% of the area, and simultaneously, there was an increase in grassland and forest cover to 43% and 18%, respectively. Moreover, the irrigated area, which was 24% of the agricultural land use in 1975, increased to 28% of the agricultural land use by 2008 (Zhang et al., 2008, for example, reported similar trends). Thus, the recharge flux responses to such changes in land use were examined in the regional model.

Unsaturated Zone Thickness
The unsaturated zone thickness was derived from a global hydrologic model suggested by Fan et al. (2013). This model provides the global distribution of water table depth constrained to observations. The derived unsaturated zone thickness ranges between 0 and 233 m with an average of 52 m. Close examination of the unsaturated zone thickness map indicates that thicker unsaturated zones are located at the center of the aquifer outcrops ( Figure 5), where the land cover is mainly forest (Figure 4). The south and west-south areas, which are characterized by high density of agricultural activity, also overlie a relatively thin to moderately thick unsaturated zone ( Figure 5). The unsaturated zone thickness values were incorporated into the simulations as described below and to allow the estimation of nitrate travel time to groundwater (see section 5.3).

Model Setup
Previous studies have reported on successful approaches for obtaining the regional spatial distribution of recharge fluxes and unsaturated nitrate transport by using GIS, Richards' equation, and the advectiondispersion equation (e.g., Assefa & Woodbury, 2013). Here unsaturated water flow and nitrate transport models were implemented as multiple 1D columns over a large area using interpolated climate with the inverse distance weighting method and gridded land use and soil parameter data. The representative grid cell spacing (3,200 m × 3,200 m) was estimated according to the suggested approach by Hengl (2006).

Water-Flow and Solute Transport Models
The unsaturated flow and nitrate transport were simulated with the Hydrus 1D code (Šimůnek et al., 2005). This code implements the Richards equation to account for water flow and the advectiondispersion equation for solute transport in the unsaturated zone. The 1D vertical Richards equation with a root water uptake sink can be written as follows: the unsaturated hydraulic conductivity function, is a function of the matric potential head, and S is a root water uptake sink term [T À1 ]. Simulation of the root water uptake rate (the sink term) was conducted according to the model suggested by Feddes et al. (1978); parameters used for the different crop (winter wheat, maize, alfalfa, and apple trees) and plant types (natural grass and conifer trees) were obtained from the Hydrus 1D database and Lv et al. (2014).
The unsaturated hydraulic functions of the different cells in the gridded data were described by the VGM hydraulic functions (see section 3.1). Moreover, a scaling factor that describes the hydraulic conductivity decay with depth was implemented according to Wang et al. (2017).
The advection-dispersion equation for conservative transport was implemented to describe chloride and nitrate transport. Since nitrogen fertilizer in China is primarily (95%) ammonium bicarbonate (Bouwman & Boumans, 2002), the following set of equations was used to model the 1D vertical transport of NH 4 and NO 3 : where C NH4 and C NO3 [M/L 3 ] are concentrations of the nitrogen species in the pore-water solution, A representative value of the longitudinal dispersivity parameter was estimated with a "trial and error" approach, comparing against field data (see section 5.1.2). There are no explicit daily nitrogen uptake rate values reported in literature for the simulated plants in the current study; therefore, the uptake values were prescribed according to general guidelines (Hanson et al., 2006;Kurtzman et al., 2013;Ramos et al., 2011). To verify that the prescribed values are within an acceptable range, simulated and observed nitrate profiles were compared (see section 5.1.2). The nitrification rates and ammonium partition coefficient values were taken from Hanson et al. (2006). The choice of a representative NH 4 volatilization rate value was based on a number of studies (Bouwman & Boumans, 2002;Kurtzman et al., 2013;Y. Li et al., 2015;Rochette et al., 2014). Most of these studies indicated that the volatilization process occurs at the top of the soil profile and ceases at 10 cm depth (e.g., Rochette et al., 2014). Therefore, for each simulation the soil profile was divided into two layers,

10.1029/2017WR022190
Water Resources Research the upper layer being 10 cm thick. Identical properties were assigned to each layer, except for volatilization rate, which was prescribed to the upper layer only.

Crop Root and Growth Models
The spatial distribution of root density with depth was assumed to follow the exponential model presented by Vrugt et al. (2001): where β(z) denotes the dimensionless spatial root distribution with depth, z m is the maximum rooting depth [L], and p z [À] and z* [L] are empirical parameters. Earlier studies in the Loess Plateau gave a range of maximum root depth values for different plant types (Canadell et al., 1996;Fan et al., 2016;Huang et al., 2004;Huang & Gallichand, 2006;Wang et al., 2013). Hence, representative values of the z m parameter could be directly estimated. The model performance using initial estimates of the fitting parameters (p z and z*), based on literature sources, was examined as discussed in section 5.1.2.
The increase in leaf area index (LAI) during the growing season for winter wheat, maize, alfalfa, grassland, and an apple orchard was estimated with the model of Leenhardt et al. (1998), where the increase in LAI is assumed a function of temperature according to the following: where LAI max is the maximum LAI of the crop, T i (°C) is the sum of temperature at the inflexion point of the curve, and b is a curvature parameter. The LAI max and the b parameters were estimated using the temperature database and reported LAI curves (Kang et al., 2003

Initial and Boundary Conditions
For each simulation, atmospheric boundary conditions with surface runoff were prescribed at the upper boundary (land surface) as water input (precipitation and irrigation), LAI, potential ET 0 , chloride and ammonium concentrations, and the minimum allowed pressure head at the soil surface (hCritA) (Šimůnek et al., 2005) at a daily temporal resolution. To estimate the potential ET 0 values, the ET ref values were multiplied with the single crop coefficients (K c ). K c values for the different crop or plant types were retrieved from the literature (Allen et al., 1998;Kang et al., 2003;Liu et al., 2017;Zheng et al., 2012). Applied irrigation rates, following the work of Kang et al. (2017) and Huang et al. (2004), were assigned to the model. Lower boundary conditions were prescribed as free drainage, and the length of the simulated soil column was defined according to the gridded unsaturated zone thickness data.
Various figures were reported for the concentration of chloride in precipitation over the LPC: 1.7 mg/L by Huang, Pang, and Edmunds (2013) Li et al. (2017), there are rain events with much higher chloride concentration than the range above. Therefore, a chloride concentration for each year was randomly sampled from a possible range of uniformly distributed values (0.25-2.5 mg/L), rather than prescribing one deterministic value.
The N input to the regional model was defined according to the recommended N fertilization inputs: 145 kg/ha/year for wheat and 180 kg/ha/year for corn, in the study region (Fan et al., 2005). Zhang et al. (2012) estimated that one third of the farmers overfertilize, while one third of the farmers underfertilize. Thus, for the local-scale simulations (see section 5.1.2), ammonium bicarbonate concentration was sampled randomly from a predefined range of values for both crops: 120-300 kg/ha/year for wheat and 160-380 kg/ha/year for corn. Another source of nitrogen that was implemented in the model is wet nitrogen deposition, which was assumed to follow mean values reported by Zhu et al. (2015). Earlier studies indicated that the anthropogenic nitrogen input (fertilizers) did not exceed plant uptake until the 1980s, and since 10.1029/2017WR022190

Water Resources Research
then, nitrogen fertilizer consumption in China has increased substantially (Huang, Pang, & Yuan, 2013;Zhang et al., 1996). Therefore, as nitrogen applications before the 1980s were insignificant, fertilization application was added from the 1980s to the model simulations.
For the crop rotation simulations (wheat-corn, wheat-corn-apple, wheat-corn-grassland, and wheat-corn-forest), the initial conditions of the first simulation were: water content at FC, and zero concentration for chloride, ammonium and nitrate. The simulated water content, chloride, ammonium, and nitrate concentration values at the end of each crop case were implemented as initial conditions for the following crop cover conditions. Each simulation had different crop parameters as described above. To minimize the influence of the initial conditions on results, the atmospheric boundary conditions were replicated three times. Therefore, each simulation ran for 67,938 days.

Examination of the Regional Model
Implementation of complex models to describe unsaturated flow and nitrate transport requires the input of many variables (e.g., soil "effective" parameters, and climate). Since employment of many parameters might lead to bias in calculations, it is necessary to evaluate the capacity of the regional model to simulate these processes. The model's performance and sources of errors in inputs were examined comparing against observations obtained from local and regional-scale studies. An extensive literature survey was conducted to search for local-scale (plot-scale) studies on groundwater recharge and nitrate fluxes across the LPC. Five sites, which are distributed over the LPC aquifer outcrops, document one measured retention curve and seven measured chloride and nitrate profiles under rain-fed agriculture cultivation (Figures 6, 7, and 8). Also, an average retention curve obtained in a large-scale study on the VGM parameters distribution over the top 5 cm of the LPC was included.

Soil Hydraulic Parameters
The gridded maps of the VGM parameters were sampled from the soil parameter maps (see section 3.1) according to the coordinates that were specified in Gates et al. (2011). The calculated retention curve using the sampled parameters as an input for Rosetta was plotted against the observed retention curve data and the two estimated retention curves reported by Gates et al. (2011; Figure 7). One retention curve was estimated using only PSD data, which was used as an input to Rosetta to predict soil hydraulic parameters. The second curve was fitted to an observed retention curve. For further comparison, an average retention curve of the entire LPC suggested by Wang et al. (2015) is plotted in Figure 7.
The estimated retention curve using only the PSD data shows poor agreement with the observed retention curve relative to the other estimated curves (Figure 7). This discrepancy is mainly related to the α parameter in the VGM. Using only the PSD data as an input to the ROSETTA, the α parameter is underestimated by an order of magnitude. When including the BD and FC (θ 330 ) data, the estimated α parameter is comparable to the estimated α values reported by Gates et al. (2011) and Wang et al. (2015; Figure 7). Note that the other parameters display comparable values. A low value of the α parameter indicates high airentry value (higher matric head) and the retention curve shifts upward (Figure 7). By using the FC value as an input in Rosetta, the retention curve is forced to pass through the FC point, which reduces the air entry value (i.e., increases the α parameter). Moreover, Schaap et al. (2001) showed that the accuracy of the estimation of the α parameter  increases considerably when one or two retention points are added to the Rosetta. Hence, following the described results, the soil parameter estimation procedure for the regional model was restricted to one of the three layers, which contained more data alongside the PSD (Figure 2). The suitability of the soil hydraulic parameters from a shallow depth layer as a representative of the deep vadose zone of the LPC was tested and compared with deep vadose zone observations (see next section).

Local Deep Vadose Zone Simulations
Previously reported chloride and nitrate profiles were used to identify the capacity of the regional unsaturated water flow and transport models to simulate processes at a local scale and for the deep vadose zone ( Figure 8). Moreover, this stage enabled a search for representative values for the following parameters: the longitudinal dispersivity, the nitrogen root uptake rates, and the root distribution model (see sections 4.1 and 4.2). Note that these values were initially prescribed according to the literature and were manually adapted in a "trial and error" approach. Five of the profiles were sampled under a wheat-corn rotation system, one profile was sampled under a field where the wheat-corn rotation system was replaced with an apple orchard, and one sampled under wheat cultivation (Figure 8). For the local-scale models, the climate was interpolated, and soil parameters were extracted from the gridded data of the regional scale by using the reported coordinates of the local-scale study sites ( Figure 6). The initial conditions were implemented as described above (section 4.3). For each location, the model was run 25 times to include the variability in chloride and nitrogen inputs. For example, the wheat-corn rotation ran for 67,938 days with a set of random values of chloride and ammonium bicarbonate concentration (1980, from day number 54,423) that were sampled from the predefined range. Subsequently, another simulation of rotation ran with different set of randomly sampled values. The simulation results are presented as ranges of the nitrate and chloride concentrations.
All profiles revealed higher chloride and nitrate concentrations at the shallower depths (root zone) than at greater depths of the vadose zone. The root zone is subjected to various factors that change seasonally (month to year): precipitation, fertilizer application, root uptake, and deep drainage. Therefore, the high observed concentrations in the shallower depths might be related to salt accumulation over a long period, or a snapshot of the soil profile just after a dry period/crop harvest or specific fertilization input, or both.

Water Resources Research
The processes in the deep vadose result in less temporal variability, compared to the root zone. The chloride concentrations in the deeper parts are mostly lower than in the root zone and display a narrower range of values. This might indicate relatively small variability of deep drainage fluxes that leave the root zone. Rainy years might lead to lower chloride concentrations and higher deep drainage rates, whereas dry periods may be characterized by higher chloride concentrations and low deep drainage rates. In contrast to the chloride profiles, the observed nitrate concentrations decrease to zero at shallower depths than the chloride. According to the analysis of Huang, Pang, and Yuan (2013), the nitrate accumulation at shallower depths of the unsaturated zone and its absence at deeper depths relate to the late overuse of nitrogen fertilizer (from the 1980s) in China and in the LPC specifically and low deep drainage fluxes.
The predicted ranges of the nitrate and chloride concentrations were compared with available observed profiles (Figure 8). Most of the simulated chloride and nitrate concentrations showed similar trends and values as the observed. However, discrepancies between observed and simulated exist in each of the chloride profiles. For example, the deep observed chloride concentrations at the Guyuan site falls within the simulated ranges, while the simulated chloride concentrations in the upper profile show lower values compared with the observed concentrations ( Figure 8a). In contrast, the simulated nitrate concentration at the same site shows a good fit with the observed nitrate concentrations (Figure 8e). Although it is difficult to determine the cause of the chloride misfit in the upper profile without full knowledge of the site's history, the discrepancies in chloride concentration might be related to a change of crop type to one with higher water demand, or a change of the fertilizer type, or both. The simulated chloride concentrations in the Zhifangou site adequately fit the chloride concentrations in the upper profile, while the simulated chloride concentrations at depth are slightly higher than the observed chloride concentrations (Figure 8d). Here the discrepancies are a result of vertical heterogeneity of the soil properties, which is not accounted for in the regional model.
The observed chloride concentrations were used to estimate recharge fluxes according the chloride mass balance method. The simulated recharge fluxes for the Zhifangou, Xifeng, and Guyuan sites are 30, 25, and 60 mm/year, respectively, and are within the range of recharge fluxes estimated with the chloride mass balance as reported in the literature (Gates et al., 2011;Huang & Pang, 2011). Note that the regional model was modified for the Guyuan site to account only for winter wheat crop with no rotation for validation purposes. The land use information for the specific cultivation in each location was not available. However, most studies indicated the dominance of the wheat-corn system across the LPC. Therefore, the winter wheat-corn rotation system was implemented in the regional model.
The aforementioned simulations provided a representative longitudinal dispersivity parameter of 7.5 cm. To obtain better fitting results, the root parameters were slightly modified for the different sites. For example, at the Zhifangou and Luochuan sites, the z m and z* parameters (equation (5)) are 2.5 and 1 m respectively, while for the Heihe (before apple plantation) and Xifeng sites, the z m and z* parameters are 2 and 0.45 m, respectively. Nevertheless, for simplicity purposes, 2 and 0.45 m were prescribed as representative values for the wheat-corn rotation system in the regional model. While the choice of the value of the maximum root depth was based on literature sources, the value of the z* parameter was chosen arbitrarily with the assumption that most of the root water uptake occurs in the upper soil. For the other plant types, the maximum root depths were prescribed according to the literature (see section 4.2).
Ultimately, the results presented above have shown the ability of the regional model to capture processes that occur at local scales and the land use conversion from a wheat-corn rotation system to tree plantation. The assumptions made for the regional model appear reasonable and indicate only minor deviations at the local scale.

Groundwater Recharge and Land-Use Change at the Regional Scale
The regional model was implemented for simulating the magnitude and spatial distribution of the mean annual groundwater recharge fluxes across the LPC (Figure 9). Two land use raster maps for 1975 and 2008 were used for analyzing the impact of land use change on recharge. The results indicate significant spatial variation across the region. Relatively high simulated recharge rates (>120 mm/year) are located in the south area of the LPC aquifer outcrops and some other locations. These areas are characterized by concentrated irrigated agricultural fields (Figure 4). More moderate recharge rates are distributed across the north parts of the outcrops, while the center of the region is characterized by low recharge rates (Figure 9).

10.1029/2017WR022190
Water Resources Research The regional model simulates rotation between periods of corn and winter wheat, and only one crop is cultivated per year. Winter wheat is cultivated between October and June when rain is low and irrigation is required to achieve better yields. Applications of irrigation keep the soil under relatively wet conditions or with no major soil moisture deficit. Therefore, during summer (June-September), when most of the rain falls (>50%), the soil conditions are preferable for deep percolation. Furthermore, the moderate recharge rates in the north part of the aquifer outcrop are related to the assumption in this study that crop development (LAI growth, equation (6)) is temperature dependent. Since the north part is characterized by low temperatures, crop development is limited. This limitation decreases the root water uptake potential and generates conditions for soil evaporation and deep drainage. For locations with higher temperatures and no irrigation, soil water replenishment by rain is not enough to ensure deep percolation.
For nonirrigated areas, the average simulated recharge flux was 64 ± 22 mm/year, which is in the range of previously reported recharge fluxes in the LPC (Gates et al., 2011;Huang & Gallichand, 2006;Huang & Pang, 2011;Huang, Pang, & Edmunds, 2013;Huang, Pang, & Yuan, 2013). Interestingly, recharge fluxes under similar climate conditions were reported in studies conducted over loess soils in other places in the world (Baran et al., 2007;Green et al., 2008;Sophocleous, 2005). The median of the simulated recharge fluxes under the irrigated areas is 150 ± 38 mm/year. Unfortunately, there are no reported studies of plot-scale recharge estimations under irrigated agricultural fields in the LPC, and the results cannot be verified with local-scale studies. However, Sophocleous (2005) and Green et al. (2008), for example, reported comparably high recharge values once irrigation was introduced to these types of soils. Previous studies on recharge in the LPC suggest different and contradicting sources for recharge, focused or diffused sources, based upon stable isotopes of unsaturated pore water and springs (Gates et al., 2011;Huang & Pang, 2011;Huang, Pang, & Edmunds, 2013;Huang, Pang, & Yuan, 2013;Z. Li et al., 2017). According to the estimates in this study, although the irrigated areas only cover 30% of all the agricultural land, their contribution to groundwater recharge is of the same magnitude as the total recharge that arrives from the rain-fed agricultural fields. Moreover, irrigated land uses are rather concentrated at specific locations (e.g., the southern parts) and contribute significantly to groundwater recharge in these areas (Figure 9).
To analyze the effect of land use change on the spatial distribution of groundwater recharge, both recharge flux maps of 1975 and 2008 are presented (Figure 9). Areas that were affected by the land use changes and displayed reduction in recharge rates are marked with dashed circles (Figure 9b). In the 1975 flux map, many locations at the top dashed circle display recharge rates between 60 and 90 mm/year (yellow, Figure 9a). By 2008, following the land use changes, the recharge fluxes in many areas declined to rates ranging between 0 and 60 mm/year (blue and green, Figure 9b). The total groundwater recharge estimates are 13.1 × 10 9 and 12.3 × 10 9 m 3 /year for 1975 and 2008, respectively, a reduction of 6.1% in groundwater recharge. Many studies reported the substantial effect of cropland conversion to pasture and forests on the soil water in the LPC

Water Resources Research
(e.g., Jia et al., 2015). However, at a regional groundwater perspective, these land use changes only slightly impact recharge fluxes.
Alongside the impact of land use change on the spatial distribution of recharge, there is also a temporal aspect. The change in recharge flux does not always occur simultaneously with the land use conversion; that is, a delay occurs. This phenomenon is influenced by many factors, such as soil type, vadose zone thickness, vegetation type, water input frequency and quantity, and climate. For further elaboration, the daily simulated recharge fluxes under two locations that have similar climate conditions are plotted in Figure 10. Both locations were cultivated as a wheat-corn rotation system before conversion to a forest and an alfalfa field. The difference between the recharge flux magnitudes at the two sites is related to the soil parameters ( Figure 10). Moreover, the vadose zone thickness under the alfalfa field is larger than that under the forest. Close examination of the recharge fluxes indicates that under the forest land use with a thinner vadose zone, the recharge decreases at a faster rate than under the alfalfa field. As shown in Figure 10, recharge fluctuation continues for a period of 5 years after the alfalfa plantation. Moreover, it appears that the recharge fluxes under both locations continue to decrease. Therefore, it is challenging to assess the full impact of change in crop cover on reduction in recharge. Further elucidation of the recharge rates requires several runs of future scenarios.
A decrease of 11% (cf. 6.1% estimated here) in groundwater recharge was previously reported by Li et al. (2014). The authors concluded that the reduction in groundwater recharge is related to the increase in ET following land use change and reduction in precipitation. Their model was calibrated according to observation well data at a specific location on the Loess Plateau assuming that the groundwater levels were affected by the climate change and land use change. This assumption might be justified in their specific location. However, in many other sites across the LPC, the vadose zone is relatively thick, which indicates long travel times of water and solutes. Any change in climate variables that might result in reduction or increase of the fluxes to groundwater would be noticeable only after a substantial period of time. To examine the potential effect of climate change on groundwater fluxes, a long-term climate database is required.

Nitrate Storage and Travel Time
The 1D vertical transport of NH 4 and NO 3 was coupled with the water flow simulations for grassland and agricultural land covers. Wet deposition of NH 4 and NO 3 were introduced to the time variable boundary conditions of both land covers, and fertilizer applications were introduced to the agricultural land uses. Subsequently, the simulation results enabled mapping of the spatial distribution of N-NO 3 leaching beyond the root zone (<2 m) of the LPC (Figure 11a). N-NO 3 leaching occurs across the entire outcrops of the LPC (Figure 11a). Relatively high leaching fluxes are located in the east side of the aquifer outcrop. As expected, Figure 10. The simulated recharge fluxes at two locations that were converted, in 1975, to alfalfa field and forest. The shaded areas designate the period before (1955)(1956)(1957)(1958)(1959)(1960)(1961)(1962)(1963)(1964)(1965)(1966)(1967)(1968)(1969)(1970)(1971)(1972)(1973)(1974) and after the conversions .

Water Resources Research
leaching rates under agricultural land uses are more intensive than under grassland and ranged between 0 and 75 kg/ha/year, with a mean N-NO 3 leaching rate of 10 kg/ha/year (Figure 11a). There is not a significant difference in N-NO 3 leaching rates between irrigated and nonirrigated fields. As discussed above regarding recharge fluxes, comparable leaching rates were reported in earlier studies conducted over loamy soils with similar crops and fertilization applications (Baran et al., 2007;Green et al., 2008;Knappe et al., 2002;Kucke & Kleeberg, 1997). Knappe et al. (2002) showed that loess soils can temporarily buffer high N fertilizer inputs and the nitrate leaching occurs in a sporadic manner. This might indicate the dominance of this soil type over other factors in dictating water and solute fluxes.
In addition to the nitrate flux quantification, elucidation of the time lag or residence time of nitrate in the unsaturated zone is necessary for monitoring, planning, and, if appropriate, for developing strategies for remediation of the groundwater system, especially for long residence times (Sousa et al., 2013;Vero et al., 2014). Earlier studies highlighted the long time lag for groundwater quality change due to the relatively thick unsaturated zone of the LPC (Huang, Pang, & Edmunds, 2013;Huang, Pang, & Yuan, 2013). However, these studies focused on rain-fed land uses. As mentioned before, there has been limited investigation of irrigated agricultural land uses in the LPC, where the correspondence between simulated high recharge rates and large N-NO 3 storage in the deep unsaturated zone might shorten the nitrate residence time in the unsaturated zone. Moreover, the map of unsaturated zone thickness indicates areas with relatively shallow groundwater ( Figure 5). For each cell, N-NO 3 velocity and travel time were estimated in two stages. First, the velocity and travel time between each pair of nodes in a 1D column were calculated, and subsequently, the calculated velocities and travel times were summed to representative velocity and travel time values (Figures 11b and  11c). The estimated average vertical velocity is 0.59 ± 0.48 m/year. The large standard deviation value indicates significant variability across the LPC (Figure 11b). Under rain-fed agricultural fields, the average velocity is 0.54 ± 0.28 m/year. Only the lower values in this range are similar to reported velocities in loess soils (Baran et al., 2007;Huang, Pang, & Yuan, 2013). However, this gap is likely due to the differences in velocity calculation methods. The shortest estimated travel times (between 5 and 100 years) are distributed throughout the outcrops of the LPC's aquifer, mainly under irrigated land uses or where the vadose zone is relatively thin (Figure 11c). A comparison between the velocity map and the travel time map reveals that travel time is not only dependent on velocity. Lastly, according to the travel time map, since the overapplication of nitrogen fertilizer started in the early 1980s, in the near future groundwater quality at some locations would be affected by anthropogenic N-NO 3 (Figure 11c).

Discussion and Future Research
There are a number of drawbacks or challenges in implementing the multiple column approach used here. A major challenge is the reliance on a wide range of data for the establishment of the model. As stated earlier, knowledge of soil properties is essential for reasonable soil parameter estimation. Moreover, there are other parameters, such as crop parameters and solute reaction parameters, which are required for utilizing this approach. However, these limitations could be addressed by extended literature surveys and field experiments. Another potential limitation with the multiple column approach is the assumption of vertical downward flow, neglecting horizontal flow and preferential flow paths (Sousa et al., 2013). Interestingly, Jana and Mohanty (2012) simulated recharge fluxes with the Hydrus 3D code and indicated the potential to include nonvertical flow paths in regional model applications. However, application of more complicated methods should be verified with deep vadose zone field data. A more fundamental challenge in using the multiple column approach over a regional scale is the debatable physical meaning of the models' calculations. The assumption that, for example, a 10-km 2 area is represented by one soil column is questionable when one examines the local-scale variability of a natural soil. However, such simplifications are essential as we move to such large scales. In this particular study we utilized an extraordinary dense array of soil sampling data.
Due to the high computational demand of the multiple columns modeling approach, the possible use of simpler nitrate quantifications should be considered. Previous studies have suggested different types of nitrate budget quantification at a regional scale, which include the transient nature of nitrate storage in deep vadose zones (Roelsma & Hendriks, 2014). Green et al. (2018) introduced a mechanics-based approach to account for reactive nitrate transport in the unsaturated zone at a regional scale. Subsequently, this method was implemented within a statistical framework based on GIS variables to reduce model complexity .

Water Resources Research
Recently, Ascott et al. (2017) adopted a simplistic nitrate budgeting model that permitted the computation of vadose zone nitrate storage on a global scale. The results calculated from such budget models could be tested against the multiple column approach calculations. Thus, two independent approaches could elaborate the necessity in highly detailed regional-scale data for understanding water flow and nitrate travel time processes in deep vadose zone.

Summary and Conclusions
Many factors, such as climate, soil properties, land use, vegetation, and vadose thickness control recharge and N-NO 3 fluxes and travel time in the unsaturated zone. Each of these factors might affect water flow and nitrate transport differently, and therefore, the quantifications of these processes are complicated. The ability to apply the Richards equation and the advection-dispersion equation helps address the challenge of combining all these processes.
A unique dense soil properties sampling campaign over the LPC aquifer's outcrops elucidated the spatial distribution of soil properties and formed the basis of a regional-scale model developed here. The data enabled us to distribute, at some level of reliance, the parameters required for applying models of unsaturated flow and nitrate transport. Raster maps of the land uses and soil parameters were produced, and climate data were interpolated. Subsequently, by assuming that each cell on the gridded landscape can be represented by a 1D soil column, groundwater recharge, nitrate storage, and travel time calculations were made. Model simulations were verified against plot-scale measurements.
Mapping the simulated groundwater recharge distribution revealed intensive recharge fluxes mainly in the southern parts of the LPC aquifer's outcrops. This indicates the significance of these parts in the total replenishment of the LPC aquifer. Furthermore, according to the model results, the land use change, transition from rain-fed crops to grassland and forest appears to have had a minor effect on total recharge. However, if the land use change trend continues toward the southern or the northern parts of the aquifer's outcrops, the effect may be significant. Another potential challenge is highlighted from nitrate transport model simulations, specifically the relatively large nitrate inventories and short travel times in locations with intensive recharge rates. Future management decisions in the southern and the northern parts of the LPC aquifer's outcrops may need to account for the high vulnerability of the groundwater system in this area.