Prediction of regional-scale groundwater recharge and nitrate storage in the vadose zone: A comparison between a global model and a regional model

Extensive nitrogen loads at the soil surface exceed plant uptake and soil biochemical capacity, and therefore lead to nitrogen accumulation in the deep vadose zone. Studies have shown that stored nitrogen in the vadose zone can eventually reach the water table and affect the quality of groundwater resources. Recently, global scale models have been implemented to quantify nitrate storage and nitrate travel time in the vadose zone. These global models are simplistic and relatively easy to implement and therefore facilitate analysis of the considered transport processes at a regional scale with no further requirements. However, the suitability of applying these models at a regional scale has not been tested. Here, we evaluate, for the first time, the performance and utility of global scale models at the regional scale. Applied to the Loess Plateau of China, we compare estimates of groundwater recharge and nitrate storage derived from global scale models with results from a regional scale approach utilizing the Richards and advection-dispersion equations. The estimated nitrate storage was compared to nitrate observations collected in the deep vadose zone (>50 m) at five sites across the Loess Plateau. Although both models predict similar spatial patterns of nitrate storage, the recharge fluxes were three times smaller and the nitrate storage was two times higher compared with the regional model. The results suggest that global scale models are a potentially useful screening tool, but require refinement for local scale applications.

. Uncertainties in predictions derived from global scale hydrological models are attributed to the coarse-resolution of the data that underpins these models (López et al., 2016). Sources of uncertainties in nutrient transport models are associated with lack of data, gaps in databases, or inconsistency of the data (Johnes & Butterfield, 2002;Leip et al., 2011;Oenema et al., 2003). Consequently, previous research has identified that comparing outputs from models of different complexities and at various scales is a major research need (Haddeland et al., 2011;Johnes & Butterfield, 2002;Koch et al., 2016;López et al., 2016).
The vadose zone has an essential role in partitioning between the precipitation, infiltration, runoff, evapotranspiration and groundwater recharge (e.g., Rempe & Dietrich, 2018). Despite this, biogeochemical processes in the vadose zone remain poorly understood (Rempe & Dietrich, 2018), and there is a dearth of quantitative information of how this compartment affects legacy nutrient dynamics (Chen et al., 2018;Marçais et al., 2018). In the past, global and regional water and nutrient transport models considered the vadose zone as a dimensionless system and represented it through lumped parameter models (Harter & Hopmans, 2004;Ronen & Sorek, 2005). Soil physicists, on the other hand, have applied non-linear models, such as the Richards equation, to model flow in the vadose zone at small scales (e.g., lab or field scales). Recently, global and regional hydrological models have been developed to quantify water cycle components employing distributed models of varying complexity, some of which include the Richards equation and the advection-dispersion equation (ADE; Ascott et al., 2017;Beusen et al., 2015;Castillo, Castelli, & Entekhabi, 2015;Keese, Scanlon, & Reedy, 2005;Koch et al., 2016;López et al., 2016;Turkeltaub, Jia, Zhu, Shao, & Binley, 2018). A commonly used global approach is the PCR-GLOBWB model (Van Beek et al., 2011) that calculates, for each time step and map cell, the water storage and water exchange in the soil in a simplistic manner, while considering the variations of elevation, land cover, vegetation, and climate. Ascott et al. (2017) utilized the PCR-GLOBWB model for estimation of recharge to derive patterns in global storage of nitrate in the vadose zone.
The accessibility of global models has resulted in extensive application in hydrological studies at different scales (e.g., Emerton et al., 2016;Hoch et al., 2017;Lehner, Döll, Alcamo, Henrichs, & Kaspar, 2006;Straatsma et al., 2020). Recently, a special issue elaborated the challenges and opportunities in the development of global scale models (Hofstra, Kroeze, Flörke, & van Vliet, 2019). In the same issue, a study presented an analysis regarding the missing linkages between global and basin/local-scale water quality models (Tang et al., 2019) Despite these recent analyses, no studies have evaluated global scale models of water flow and nitrate transport in the unsaturated zone. Evaluation of global scale models to simulate unsaturated processes is not trivial, as there are no observations at global scales.
Regional-scale models of the unsaturated zone are more detailed and based on intensive local data, which yield high credibility to their outputs (Assefa & Woodbury, 2013;Keese et al., 2005;Turkeltaub et al., 2018). In addition, the regional models cover large areas that overlap with the global model scales. Therefore, a regional model could be used to examine the outputs of global models of the unsaturated zone. The objective of this study is to evaluate groundwater recharge fluxes and nitrate storage predicted by global and regional models. For this analysis, we compare regional recharge and nitrate storage in the Loess Plateau of China (LPC) calculated by a regional model based on well-established approaches utilizing the Richards and ADE (Turkeltaub et al., 2018) and a global modelling approach (Ascott et al., 2017). Simulations from the models are compared to local soil sampling investigations. Following the inferences of this comparison, we assess the suitability of applying this global model for estimating processes at a regional scale.
2 | METHODOLOGY 2.1 | The loess deposit and the LPC Loess is an aeolian deposit that evolved mainly during the Quaternary and covers 10% of the Earth's surface (Pye, 1995;Smalley, Markovi c, & Svirčev, 2011, Figure 1). The loess sediments are dominated by silt grain size, often resulting in a limited spatial variability of the soil properties (Smalley & Markovi c, 2014). The abundance in silt particles facilitates exceptional conditions for agriculture cultivation, which place the loess sediments among the most fertile and productive soils (Catt, 2001). Unconfined groundwater systems can be found under the loess sediments or the loess sediment could be divided to saturated (groundwater) and unsaturated zones (El Etreiby & Laudelout, 1988). Therefore, the loess is also an important recharge source for most of these groundwater resources. Intensive agricultural cultivation of the loess sediments has raised concern about degradation of groundwater quality due to nitrate and other agrochemical sources in different parts of the world (Baran, Richert, & Mouvet, 2007;El Etreiby & Laudelout, 1988;Huang, Pang, & Yuan, 2013;Isla, Londoño, & Cortizo, 2018;Keller, Butcher, Smith, & Allen-King, 2008;Wagner & Roberts, 1998).
The LPC is the thickest and largest loess deposit in the world (Kukla, 1987). Groundwater supplies 22% of the total water supply in the LPC, which accommodates more than 100 million people (Li & Qian, 2018;Zhao, Mu, Wen, Wang, & Gao, 2013). The LPC is comprised of arid, semiarid and semi-humid regions in the north of China.
Most precipitation occurs in the rainy season from June to September (55-78%) in the form of high intensity rainstorms, ranging, per annum, from 226 mm in the northwest to 683 mm in the southeast (Xin, Yu, Li, & Lu, 2011). The annual estimated evaporation is 650-1,200 mm and the mean annual temperature ranges from 3.6 C in the northwest to 14.3 C in the southeast. An unconfined aquifer is embedded within the loess sediments and the water table is located on average at 52 m depth, but can vary between 0 and 233 m according to the model suggested by Fan, Li, and Miguez-Macho (2013). This groundwater resource has been overexploited, and the regional water table is in rapid decline (Huang & Pang, 2011;Li et al., 2014). Additionally, for soil stabilization reasons, many soil conservation measures were applied across the LPC, which included significant land use changes (Jia, Shao, Yu, Zhang, & Binley, 2019;Zhang, Zhang, Zhao, Rustomji, & Hairsine, 2008). The land use/cover distribution over the outcrops of the LPC aquifer are shown in Figure 1. Irrigated agriculture is mainly in the south outcrops of the aquifer, while rainfed agriculture covers the western parts and the northern parts of the outcrops. Forest and grassland coverage dominates the central and east of the aquifer's outcrops. The LPC is a unique environment given the relatively insignificant soil texture variability. Large databases of soil properties and climate variables exist, making the LPC an ideal focus for comparing modelling approaches.
In order to characterize the potential accumulation of nitrate in the deep loess vadose zone across the LPC, loess samples were collected at five sites from land surface to bedrock (Figure 1, Jia et al., 2018). The observed nitrate storage profiles are used to evaluate the performances of the global and regional models' predictions. These profiles were not used to calibrate either model, and thus serve as independent data.

| Regional and global model approaches
Two approaches, a global approach and a regional approach were implemented to calculate groundwater recharge and nitrate storage in the vadose zone across the LPC. The local approach is based on a large database of local soil properties, climate variables, vadose zone thickness and land use/land cover (Turkeltaub et al., 2018). These data are prepared in gridded format (raster maps), and the daily climate data are interpolated with the inverse distance weighting method from local meteorological stations to the specific cell on the gridded map. Note that for the current study, the maps are reconstructed to a 48 km × 48 km pixel size; such significant coarsening (upscaling) of the grid was implemented in order to use a comparable grid scale to that in the global model: a much 16× finer grid (3 km × 3 km pixel size) was implemented in Turkeltaub et al. (2018). The region being investigated is discretized into multiple 1D columns, with no exchange between columns. Groundwater recharge fluxes in the vadose zone of a column are calculated using the Richards equation, which is coupled with the ADE to simulate ammonium and nitrate transport throughout the vertical profile of the column. These equations are solved with the Hydrus 1D code (Šimů nek, Šejna, & van Genuchten, 2005). The nitrogen processes that are included in the multicolumn model are as follows: ammonium volatilization, adsorption and nitrification, nitrate denitrification, and nitrate passive root uptake. Nitrogen inputs are wet nitrogen deposition, which was assumed to follow mean values reported by Zhu et al. (2015) and anthropogenic nitrogen fertilizer.
Earlier studies indicated that the anthropogenic nitrogen input (fertilizers) did not exceed plant uptake until the 1980s, and since then, nitrogen fertilizer consumption in China has increased substantially (Huang, Pang, & Yuan, 2013;Zhang, Tian, Zhang, & Li, 1996). Therefore, as nitrogen applications before the 1980s were insignificant, fertilization application was added from the 1980s to the model simulations. The average pore water velocities were calculated by dividing the yearly recharge with the yearly mean water content for each cell.
F I G U R E 1 Global loess distribution (a) derived from the Global Unconsolidated Sediments Map database (GUM) as was reported by Börker, Hartmann, Amann, and Romero-Mujalli (2018) source: https://doi.org/10.1594/PANGAEA.884822). USGS Land Use/Land Cover (b) for the Loess Plateau of China (https://lta.cr.usgs.gov/glcc/eadoc2_0), the approximate location of the unconfined groundwater system modified from Huang, Pang, and Edmunds (2013). Pixel size 1,000 m × 1,000 m. The black circles represent the locations of five study sites where the nitrate storage of deep vadose zone profiles was assessed and investigated Subsequently, the vadose zone thickness was divided by the average pore water velocities to derive the yearly average nitrate travel times in the vadose zone.
For the global modelling approach used by Ascott et al. (2017), a number of existing models are integrated to estimate nitrate stored in the vadose zone. The PCR-GLOBWB model, which is a "leaky bucket" type of model, is used to estimate the regional groundwater recharge distribution across the LPC aquifer's outcrops for each grid cell (0.5 × 0.5 discretization) (Wada et al., 2010). PCR-GLOBWB calculates the water storage in two vertically stacked soil layers, as well as the water exchange between the layers and between the top layer and the atmosphere (precipitation, evaporation, and snow melt). Additionally, short vegetation extracts water from the upper layer only, while tall vegetation extracts water from both soil layers. Using a piston-flow assumption, recharge estimates are then combined with depth to water table data estimated by Fan et al. (2013) and nearsurface porosity estimates by Gleeson, Moosdorf, Hartmann, and van Beek (2014) to derive an estimate of the travel time for nitrate in the vadose zone. Nitrate leaching from the base of the soil zone was derived from the IMAGE model (Beusen et al., 2015) on an annual basis on a 0.5 grid. The nitrogen input of the IMAGE model is based on nutrient data covering the period 1900-2000 presented by Bouwman et al. (2013). This study illustrated, by subdividing the 20th century to two periods, that between 1900 and 1950, soil N surplus almost doubled compared to the period before 1900, and between 1950 and 2000, soil N surplus was nearly eight times more than before 1900.
For each grid cell, nitrate leaching estimates were combined with the derived travel times to calculate nitrate stored in the vadose zone, considering a simulation period of 1958-2000. In Tables 1 and 2, there are additional details regarding the different parameters and components of the regional and global models. For further description of the regional and the global approaches, the reader is referred to the publications by Turkeltaub et al. (2018) and Ascott et al. (2017). The models' performances are evaluated by comparing between models' predictions and local observations ( Figure 1b). We recognize that a significant contrast in spatial scale between observation and model state, for both models. However, there is an information scarcity regarding vadose zone nitrate storage beyond 4 m depth in the LPC . In addition, an earlier regional (km scale) study indicated that soils in the LPC exhibit limited textural variation horizontally and vertically, which allows such observations to be effective data for comparison between nitrate vadose zone storage predictions that were estimated by the two models (Zhao, Shao, Jia, & Zhang, 2016). An additional comparison was conducted to evaluate the differences between the model inputs, local data of climate and soil parameters were compared with the PCR-GLOBWB model inputs (see Supplementary Information). The meteorological data includes the mean monthly temperature, monthly Penman-Monteith potential evapotranspiration (PET), which is implemented in both model approaches, and monthly precipitation. The soil data include the saturated water content and the saturated hydraulic conductivity at different depths. Figure 2 shows the long-term average annual groundwater recharge for both the global (PCR-GLOBWB) and regional approaches. There are significant differences in the simulated recharge spatial variability and magnitude between the two methods ( Figure 2). According to the PCR-GLOBWB model predictions, the perennial average recharge flux is about 12 mm/year (Figure 2c). Additionally, the intensive recharge rates occur mainly in the southern part of the LPC outcrops (concentrated with agriculture activity) and very low recharge fluxes elsewhere ( Figure 2a). The perennial average recharge flux calculated by the regional approach is about three times larger (38 mm/year) than that from the PCR-GLOBWB model (Figure 2c). Moreover, the predicted recharge fluxes from the regional model exhibit high fluxes in the central-north of the LPC outcrops, which according to the land use map (Figure 1), is covered mainly with grass ( Figure 2b). Wu, Si, He, and Wu (2019) reported average annual groundwater recharge rates of 39.9 ± 26.5 mm/year and 48.3 ± 12.5 mm/year according to local investigations and satellite information, respectively. Both the satellite data and local methods indicate similar recharge rates to the recharge rates predicted by the regional approach here. It appears, therefore, that the PCR-GLOBWB model, when applied at a regional scale, underestimates the recharge rates in the LPC.

| Groundwater recharge
Groundwater recharge is controlled by climate, soil properties and vegetation. These variables effect the groundwater spatial and temporal distribution. Therefore, it is challenging to determine a dominant factor that causes the differences between the regional model estimations and the PCR-GLOBWB model estimations. To elucidate the similarities and differences between the models, the meteorological and soil parameters inputs of the PCR-GLOBWB model were compared to local observations, which are the inputs for the regional approach (see Supplementary Information). The analysis indicates that the monthly mean temperature of the PCR-GLOBWB model inputs are very similar or almost identical to the observed monthly mean temperature values (r = .99). Relatively high correlations (r = .96) were calculated between the PET inputs of the PCR-GLOBWB model and local observations. However, the PET inputs of the PCR-GLOBWB model are generally 20% lower than the measured values. Moderate correlation (r = .75) was calculated between local precipitation measurements and the precipitations inputs of the PCR-GLOBWB model.
In addition, the precipitation inputs are 10% higher than the local observations. This is an unexpected result considering that higher recharge rates should be calculated under conditions of smaller PET and higher precipitation. Nevertheless, the comparison between the soil parameters of the two models show very poor correlations, where the PCR-GLOBWB model inputs do not capture the LPC regional soil variability (see Supplementary Information). Generally, the saturated hydraulic conductivity (Ksat) values of the PCR-GLOBWB model inputs are lower than those observed. It is possible that the low Ksat values in the PCR-GLOBWB model might encourage higher runoff and evaporation rates compared with the regional approach. Note T A B L E 1 Summary of the water flow parameters as integrated to the global and regional models that other parameters and factors in both models are incomparable due to the differences in assumptions, structure and equations that construct the models.

| Nitrate storage in the vadose zone
To illustrate the consequences of the regional and global model They concluded that the higher masses of STN occur under croplands and in regions with higher precipitation and temperatures. These conditions are mainly located in the south central parts of the aquifer's outcrops ( Figure 1). Hence, the predictions by the global approach and the regional approach agree with these regional investigations.
The discrepancies between the two approaches are illustrated by the wider spatial distribution of nitrate storage and larger magnitude computed from the global approach, in comparison with the regional model output (Figures 3 and 4). According to the global model output, an intensive nitrate accumulation started in the 1950s and showed a rapid increase in the mid-1960s (Figure 3). In contrast, the predictions T A B L E 2 Summary of the nitrogen fate parameters as integrated to the global model and regional model F I G U R E 2 The groundwater recharge maps of the Loess Plateau of China predicted with the (a) global approach (PCR-GLOBWB model), (b) regional approach with daily climate inputs (pixel size 48 km × 48 km), and (c) boxplot of the groundwater recharge. The horizontal line shows the median groundwater recharge, the box shows the second and third quartile range and the whiskers show the first and fourth quartiles calculated with the regional approach indicate intensive nitrate accumulation started at the beginning of the 1980s (Figure 3). Further comparison between the simulated nitrate storage of the global and regional approaches and local observations indicates that the global approach overestimates the nitrate storage, while the regional model simulations are comparable to the nitrate storage field-based observations ( Figure 4). This is an indication that the intensive nitrogen input in the global approach starts too early, well before it actually started in the LPC region, and in China (Huang, Pang, & Yuan, 2013;Zhang et al., 1996). Only very specific locations show short travel times and are not part of the general trend ( Figure 5). The travel times estimated by the global approach are four times larger than the estimated travel times of the regional approach. This is to be expected since the recharge fluxes calculated by the global model are three times smaller than the estimated recharge fluxes by the regional approach.

| Nitrate travel time in the vadose zone
The magnitude of travel times shown in Figure 5 reveal that the nitrate storage maps in Figure 4, on the whole, represent an accumulation of nitrate over the complete simulation period, since very little will have reached the regional groundwater body. Therefore, the storage maps in Figure 4 should be reasonably similar. However, given the contrast between travel times estimations obtained by the global model and the regional model ( Figure 5), simulating for a longer time period (centuries) would lead to a greater contrast in nitrate storage

| DISCUSSION
An investigation of groundwater recharge and nitrate storage in the vadose zone at a regional scale is necessary for appropriate management of unconfined aquifers. This study provides a first test of a global scale model at the regional scale for investigation of vadose zone flow and transport processes. Various studies that were conducted over loess sediments under different climate conditions across the world have reported similar recharge and nitrate leaching fluxes magnitude (e.g., Baran et al., 2007;Green et al., 2008;O'Geen, McDaniel, & Boll, 2002;Sophocleous, 2005;Wu et al., 2019). Moreover, O'Geen, McDaniel, Boll, and Keller (2005) illustrated that the level of heterogeneity of loess will control the recharge rates. Higher recharge rates were estimated for locations with homogeneous loess. This might suggest that the water and nitrate fluxes in loess sediments are largely controlled by the loess physical properties. Furthermore, as is indicated in this study, the discrepancies of the groundwater recharge predictions of the PCR-GLOBWB model is an outcome of a combination of factors that might be not well represented because of the coarse resolution, especially with regard to the soil parameters. Therefore, global models should be adjusted according to the findings of this study (i.e., input of detailed soil properties and climate), before being applied to other loess regions in the world. Notably, similar methods to the one presented in the current study, are implemented to improve global soil maps. A "bottom-up" approach that involves the collection of soil profile data is combined with a "top-down" approach to produce gridded maps by using global modelling (Arrouays et al., 2017). Here as well, the "bottom-up" approach (local modelling) provides better results than global modelling, due to generalization of different relations between covariates and soil properties. Nevertheless, there is no doubt regarding the benefit of using top-down products: they provide early proof of concept (a screening tool), can be updateable according to local data, and combined with lower scale methods using ensemble approaches (Arrouays et al., 2017).
van  suggested the use of consistent spatial/ temporal resolutions of the inputs and output when comparing between models for uncertainty analysis. However, in the current study, the harmonization of the input and output datasets of the models was limited by the differences in the structure of the models. These distinctions between the models challenged the search for a dominant process that leads to the uncertainties in the outputs of the global model. Although mechanistic models might improve our understanding of the processes occur at the basin/local scales (Tang et al., 2019), it is still unclear how to integrate these inferences to the global approaches.
Despite the discrepancies between the models' inputs, due to relatively low recharge rates in the loess vadose zone, the impact of the variability of vegetation and soil properties on nitrate transport might be less significant in the LPC. In addition, the combination of a thick unsaturated zone and low rates of downward movement in the LPC results in long nitrate travel times. Nevertheless, previous studies indicated intensive water fluxes in locations with coarser soil types or vegetation with shallow root systems, which facilitated nitrate leaching (e.g., Green et al., 2008;Turkeltaub et al., 2015). For environments with soil type other than loess, and with larger soil and vegetation variability, the nitrate leaching predictions presented in this paper cannot be directly replicated. An additional challenge is the implementation of preferential water flow in vadose zone models. Our current understanding, observations and number of studies regarding preferential flow are limited.
Therefore, both the regional and global approaches cannot account for the contribution of preferential flow.
Clearly, the global model used here could be improved by local calibrations, although this could be argued as undermining a key value of such approaches. As a first step to improve global nitrogen F I G U R E 5 N-NO3 travel time maps in the Loess Plateau of China as predicted by the (a) global approach and (b) regional approach with daily climate inputs (pixel size 48 km × 48 km) prediction, finer temporal (of decades) subdivisions of the anthropogenic nitrogen application should be implemented, instead of the coarse subdivision of the 20th century for two periods. Moreover, previous studies that presented different calibration procedures for nutrient and hydrological global models relied on observations that were obtained at catchment, watershed and regional scales (e.g., Beusen et al., 2015;López et al., 2016). Mostly, these observations are obtained from surface water resources, for example, river discharge, temperature and nutrient concentrations or databases obtained from topsoil. Currently, the deep vadose zone database, especially with regard to nutrient inventories, is limited and fragmented. The establishment of a global dataset of local groundwater recharge fluxes and borehole information could contribute to the improvement of the vadose transport simulations by global models.

| CONCLUSION
This study is the first to evaluate the performance and utility of a global scale model of vadose zone nitrate storage and transport at the regional scale. Relatively large differences in predicted recharge rates between the approaches are related to over/under estimation of the meteorological conditions, mainly PET and rainfall, and inadequate representation of the soil parameters (Ksat and saturated water content). This is probably a consequence of the coarse resolution in the global scale model. In our application to the LPC, the total nitrate storage in the vadose zone is overpredicted by the global approach. However, the nitrate storage maps produced by the global and the regional approaches show similar spatial patterns: large nitrate inventories in the south central parts of the aquifer's outcrops and a decreasing trend in nitrate storage from south to north.
The results obtained in this study could be implemented for other loess environments worldwide. In regions with larger soil and vegetation variability, detailed information regarding these variables is required and as well as investigations that include sensitivity analysis of the possible impact of the soil and vegetation variability on the predicted fluxes.
Other issues such as the contribution of preferential flow to recharge and nitrate fluxes at the regional scale should get more attention. Further work should include the implementation of local scale data in the global model for better representation of the vadose zone processes.
Ultimately, this study benefitted from an extensive database of observations. In an absence of these types of data, global models could be used only as a primary step, in decision of recognizing locations where investigations of plot, field, and regional scales are required. https://www.bgs.ac.uk/services/ngdc/citedData/catalogue/800dd09b-5848-4803-a70c-cd15d5590f16.html. The regional model is available for public download through the following link: 10.17635/lancaster/ researchdata/316.
The regional model is available for public download through the following link: 10.17635/lancaster/researchdata/316.