A Comparative Study of Conceptual Model Complexity to Describe Water Flow and 1 Nitrate Transport in Deep Unsaturated Loess 2 3

28 Understanding nitrate migration through the deep vadose zone is essential for aquifer 29 vulnerability assessments. The effect of variability of physical properties of the deep vadose zone 30 on nitrate transport has been scarcely explored. Recently, deep nitrate storage profiles were 31 determined in the vadose zone of the Loess Plateau of China. Using these observations along 32 with measured soil properties, this study investigates the effect of loess vertical heterogeneity on 33 water movement and nitrate transport through the deep vadose zone. Models of different 34 complexity were established and calibrated. First, a simple piston flow and nitrate mass balance 35 approach was calibrated to the observed nitrate storage. The results indicate that the total nitrate 36 storage is estimated well, while the estimation of the distribution of nitrate is relatively poor. 37 Subsequently, Richards’ equation and the Advection-Dispersion equation were evaluated. Three 38 different conceptualizations of the numerical models were calibrated against deep vadose zone 39 nitrate and water content observations: (1) one-layer model assuming homogenous loess vadose 40 zone; (2) a model that considers a hydraulic conductivity (Ks) decay function and (3) a model 41 where the Miller-Miller scaling factors are prescribed to account for changes of the hydraulic 42 functions with depth. Accounting for the vertical Ks decay in the numerical models improved 43 water flow performances. The study reveals the adequacy of implementing water flow and nitrate 44 transport numerical models together with a simple representation of the vertical loess variability, 45 for simulating nitrate migration in loess deep vadose zone environments. 46

. The green circles represent the four study 172 sites where deep boreholes were drilled (Jia el al., 2018).  175 A full description of the soil sampling and soil physical analysis relating to data used here can be 176 found in Jia et al. (2018), and, therefore, only a brief explanation of the database is given here.

177
The deep vadose zone profiles were drilled (15 cm diameter borehole) from the land surface to 178 bedrock between May and June 2016, using the under-reamer method (Overburden Drilling Exploration;Izbicki et al., 2000). As this study focuses only on the unsaturated zone, data collected below the water table is not considered. The water tables are located at 81, 96, 141 and 181 54 m depth in Yangling, Changwu, An'sai and Shenmu, respectively. For the Changwu, An'sai 182 and Shenmu sites, soil cores were sampled at 1m intervals, while at the Yangling site soil cores 183 were taken every 0.5 m to a depth of 10 m, and at 1 m intervals for depths beyond 10 m. The 184 samples were analyzed for gravimetric water content, particle size distribution, bulk density, pH, 185  O in nitrate. For further details, the reader is referred to Jia et 186 al. (2018). Note that the volumetric water contents were computed from the gravimetric water 187 content and bulk density. For clarity, the data are elaborated below briefly.
188 Figure 2 shows the vertical variability of the particle size distribution (PSD) at the four study 189 sites, organized from south to north of the LPC (Figure 1). The PSD data clearly confirm the 190 predominance of the silt fraction for all sites. Nevertheless, there is a distinguished increase in 191 sand fraction and decrease in silt and clay fractions with depth from Yangling (south) to Shenmu 192 (north), which is accordance with previous regional studies that investigated distribution of loess 193 particle size extensively (e.g. Derbyshire, 2001). 194 The bulk density measurements show a more complex behavior compared to the textural 195 information ( Figure 3). An increase in bulk density from soil surface to about 10 m depth can be 196 seen at all sites. Note that the depth where the bulk density stabilizes is different for each site. 197 Moreover, some transitions in bulk density can be observed at depth (e.g., 40 m at Yangling). 198 This phenomenon was not reported previously and its effect on the water flow in the vadose zone 199 is unknown. As is illustrated by Derbyshire (2001)  Soil retention curves and unsaturated hydraulic curves are commonly described according to the 223 van Genuchten-Mualem (VGM) model (Mualem, 1976;van Genuchten, 1980):

228
Hydraulic conductivity is often described by:

230
where Ks [L T -1 ] is the saturated hydraulic conductivity and l is the pore connectivity parameter 231 prescribed as 0.5.

232
Three undisturbed soil cores were collected from the upper 10 cm of the soil at each of the four 233 deep vadose zone study sites ( Figure 1, green circles point of previous regional studies (Wang et al., 2013b;Jia et al., 2015).  zone, an objective function (Φ(p)) was establish to minimize the differences between the 280 reference hydraulic curves and the hydraulic curves estimated from Rosetta3 for each depth i.

281
Note that the fitted hydraulic parameters in Table 1 were assumed to represent the reference 282 hydraulic curves. The objective function is as follows (Clausnitzer et al., 1992;Nasta & Romano, 283 2016): where p is the parameter vectors that includes all scaling factors, N denotes the number of values 287 of the degree of saturation, Se, ranging from 0 to 1. Note that only the soil hydraulic functions 288 parameters were available. Therefore, the hi,η and Ki,η were estimated by implementing a range of  Here, the water balance was calculated at a daily resolution.

307
The transpiration rate is controlled by two mechanisms: the atmospheric demand and the supply 308 of water in the soil. Assuming that there is no effect of salts, the uptake function can be described where Tp represents the maximum transpiration rate as dictated by the atmospheric demand, S* is 312 a threshold value and Sw is the saturation at which the uptake is zero and the plant wilts.

313
Evaporation from the root zone is described similar to the root water uptake approach: where Sh is the hygroscopic saturation, at which evaporation ceases. Note that the average 316 saturation at which E reaches its maximum, S*, is the same as that for transpiration (equation 9).

317
As described by Laio et al. (2001), drainage of the root zone (potential recharge) is set equal to 318 the unsaturated hydraulic conductivity, described with an exponential form: where Ks is the saturated hydraulic conductivity, β is a parameter of the soil and Sfc is the field 321 capacity.

322
The nitrate leaching rate is calculated simultaneously to the water balance approach according to 323 the following mass balance approach: The nitrate mass balance approach is based on the assumption that three main processes control 326 the nitrate dynamics in the root zone; passive nitrate uptake (NO3 Uptake), and nitrate leaching 327 (NO3 Leaching). Note that the NO3 Soil represent the residual nitrate in the root zone and NO3 Input is 328 the fertilizer amount.

329
The passive nitrate uptake is assumed to be proportional to the transpiration rate, T(s) in equation where kl is the leaching coefficient and equals to 0.02 in the current study.

337
To calculate vertical nitrate displacement in the vadose zone using the piston flow approach, 338 pore water velocities are estimated using the recharge fluxes (R). Water velocities are linearly 339 related to the calculated water fluxes by a coefficient that represents the specific volume through 340 which the water and solutes are transported. Assuming only vertical advective transport of nitrate 341 in the loess unsaturated zone, the displacement ΔZ of the nitrate is given by: where R(t) is the (time-dependent) recharge, and Δt is the time step. The np parameter in this 344 study was set to the average measured water contents in the vadose zone of each of the study 345 sites ( Figure 1).

346
The piston flow and nitrate mass balance model were calibrated against the total nitrate storage 347 calculated from the nitrate profiles in Figure 4 and the observed water contents (see Supporting 348 Information). Only the β parameter in equation 11 was modified during the calibration process.  Simulation of the root water uptake rate (the sink term) was conducted according to the model suggested by Feddes et al. (1978). The parameters used for the different plant types were 362 obtained from the Hydrus 1D database.

363
The following set of equations was used to model the 1D vertical transport of NO3: to the latter suggestion for simulations with no model calibration (see model calibration section).

376
The spatial distribution of root density with depth was assumed to follow the exponential model 377 presented by Vrugt et al. (2001): The three conceptual models were calibrated against water content profiles (Supporting 394 Information) and nitrate profiles ( Figure 4) obtained from the loess vadose zones. An inverse 395 problem was formulated to find an optimum combination of parameters that minimizes the 396 following objective function:   parameter was considered to be a fitting parameter in the current study.

477
In Yangling, most recharge and nitrate leaching occurred during the cultivation periods (Table   478 2). These intensive fluxes during the crop period are probably due to irrigation that exceeds plant  Yangling and Changwu are substantially over-fertilized.   therefore, seems to be site specific ( Figure 6). In addition, the cessation of the decay in Ks might 526 indicate the reduction of the compaction effect; the reason why this is site specific is not clear. The vertical variability of soil properties can also be described by Miller-Miller scaling factors.

535
Firstly, the VGM parameters for soil samples from all depths were estimated using Rosetta3 536 based on the texture and bulk density data for each sampling point along the four vertical profiles 537 (Figures 2 and 3). A single set of scaling factors for each site were calculated by applying 538 equations 4 to 7. The unsaturated hydraulic functions in Table 1 were used as the reference 539 curves.

540
The scaling factors mostly ranged between 0.2 and 0.8, which indicates that there is no strong 541 contrast of the soil physical properties with depth ( Figure 7). In general, the scaling factors 542 display higher values close to the loess surface and a decrease with depth, which is similar to the observed trend of the bulk density measurements along the profiles. The scaling factors of the 544 Yangling site show slightly higher values and variability compared to the other sites (Figure 7a).

545
Note that the Yangling site is differ from the other sites by the climate conditions (higher 546 temperature, ET and rain) and this site is also intensively cultivated (double cropping system).

547
The agriculture cultivation might affect the variability of the near surface loess (e.g., HET) are examined and compared. All three models were further calibrated against the water 561 content and nitrate observations (Figures 9 and 10). The calibrated parameters are documented in 562 Supporting Information (Table S1). Sensitivity tests showed that the parameters n and θs in 563 equations 1 and 2, are the most sensitive for all sites. However, to improve model performances, the Ks parameters were also included in the calibration process. In Changwu and Shenmu sites, 565 the longitudinal dispersivity (λ) parameter was also included as a calibration parameter.

566
Statistical evaluation of model performance is summarized in Tables 3 and 4. To explore the 567 contribution of accounting for the loess vertical variability in the modeling, the statistics of 568 model performances prior to calibration are also shown (Table 3 and 4).

569
The observed water content profiles from the four sites are plotted together with the simulated 570 values from the three model configurations in Figure 8. In general, the simulated water content is 571 within the same range as the observed volumetric water contents. Nevertheless, according to the 572 evaluation results in Table 3, the EXP model that includes the Ks decay function produces a 573 relatively better performance compared with the HOM and HET models for all sites. In fact, the 574 HET model shows the poorest performance (Table 3). A comparison between evaluation results 575 before and after calibration illustrates that calibration is necessary for most cases to improve 576 model performance.

577
The results of the nitrate simulations are less conclusive (Figure 9, Table 4). Essentially, all 578 models show vertical nitrate distributions that are similar to the observed nitrate profiles ( Figure   579 9). Here, as before, the HET model displays the poorest performance. The HOM shows similar 580 or slightly better performance relatively to the EXP model for all sites (Table 3, Figure 9). When 581 comparing between the evaluation results of calibrated and uncalibrated models, for the 582 Changwu, An'sai and Shenmu sites, the calibration procedure improves model simulations.

583
However, for Yangling site no changes in the simulation results are apparent.  Table 3. Evaluation results of the water content profiles as were simulated by the three 590 conceptual models (Figure 8)   implies that the nitrate arrives as one concentrated cluster, while the numerical models account 658 for nitrate arrival that is spread over time ( Figure 10)