Effects of upscaling temporal resolution of groundwater flow and transport boundary conditions on the performance of nitrate-transport models at the regional management scale

Managed aquifer recharge and changes in crop type or nutrient management on agricultural lands are promising approaches to address groundwater quality degradation by nitrate. Tools to assess nonpoint-source contaminant transport are needed to better understand the interaction between agricultural management practices and long-term nitrate dynamics in groundwater basins. This study investigates the impact of time-resolution upscaling of groundwater flow stresses (i.e., recharge, pumping, and evapotranspiration rates) on the long-term prediction of nitrate transport at the regional scale. A three-dimensional, monthly transient flow and nitrate-transport model using MODFLOW and MT3D is applied as the reference simulation. The reference model results are compared to temporally upscaled models with (1) upscaled annual-averaged flow and transport stresses and (2) steady-state flow stresses, across different management scenarios. Models with annual-averaged flow and nitrate-loading stresses were found to be the best alternative to the reference model. However, employing a steady-state flow field to parameterize transient transport models, using a time series of spatially variable annual total contaminant loading, provides a useful alternative to predict the trend and variability of nitrate-concentration breakthrough curves at wells across the regional scale and to differentiate the effects of various agricultural management scenarios, if the history of the source contaminant mass is known. The difference between concentrations resulting from steady-state-flow versus transient-flow models is less than 2 mgN/L for nearly 75% of shallow groundwater cells in the model. However, the steady-state-flow-model-based transport simulation does not capture short-term oscillations of nitrate concentrations in pumping wells at the local scale.


Introduction
Failure to preserve clean groundwater is of growing concern due to the important role of groundwater as a source of drinking water worldwide. Nitrate is among the most common pollutants, especially in regions with significant agricultural activities; these are regions that are also among the most likely to rely on domestic wells for drinking water (Johnson and Belitz 2017;Ransom et al. 2017). High levels of nitrate in groundwater degrades the quality of drinking water in municipal-supply and domestic wells and creates health issues. Different sources of nitrate pollution are distributed continuously across landscapes (nonpoint-source groundwater contamination) and nitrate in rural regions originates mainly from agricultural activities (Dalgaard et al. 2014;Harter et al. 2017;King et al. 2012;Van Grieken et al. 2019).
The nitrate contamination problem can be addressed in the long run by managing pollutant sources, while different strategies have been introduced in the short term to remove or attenuate nitrate in drinking water pumped from aquifers . Regulators, policy-and decision-makers, the public, and the regulated agricultural community want to understand the feasibility of long-term improvements in groundwater quality that can be achieved with a range of Electronic supplementary material The online version of this article (https://doi.org/10.1007/s10040-020-02133-x) contains supplementary material, which is available to authorized users. source management practices with respect to an ensemble of domestic and public water supply wells in a particular aquifer system of concern. Properly designed simulation models are needed to assess the impact of promising irrigation-, nutrient-, and recharge-management practices on long-term distributed groundwater quality (Fogg and LaBolle 2006). Nonpoint source pollutants originate from across a regional landscape, with spatially and temporally highly variable leaching to groundwater on a continuous basis (Kourakos et al. 2012;Hashemi et al. 2016;Henri and Harter 2019;Padilla et al. 2018). For point source contamination, e.g., by volatile organic compounds (VOCs), regulated maximum contaminant levels (MCLs) are commonly in the low parts per billion (ppb) range and concentrations are often several orders of magnitude above those MCLs. In contrast, nitrate source concentrations are typically less than one order of magnitude above the nitrate MCL (Burow et al. 2010;Hamilton and Helsel 1995;Nolan et al. 1997).
Hence, policy concerns about nonpoint sources differ from those for point source pollution. The latter typically encompass both source control and aquifer remediation, focusing on preventing or reversing the impact on specific wells at a specific site, while nonpoint source management often focuses exclusively on source control only, and does not entail active aquifer remediation (Dowd et al. 2008;Harter 2015). Nonpoint source policy necessarily is about ensemble behavior of the pollutant in hundreds to thousands of potentially impacted wells. Therefore, predicting exact local nonpointsource contaminant behavior is less important, while the more critical decision support will be to correctly predict regional distribution, variability, and trends of water quality across the wells of concern, and to explore how these may change over time without and with land management changes.
One of the most important challenges in numerical simulation of groundwater contaminant transport is to employ proper spatial and temporal scales for both flow and transport simulations (Cushman et al. 2002;Dentz et al. 2011;Scheibe et al. 2015). The scale at which fundamental processes are understood, parameters are measured, and data are collected (e.g. microscopic to meters, and seconds to days) is different from the scale at which model applications are discretized (e.g. meters to kilometers, and months to decades). Measurements taken at the meter scale suggest the need for high spatial and temporal resolution of flow and transport processes. However, given the size of regional aquifer systems, computational costs typically dictate much coarser resolution. Instead, upscaled or effective groundwater flow and transport models for practical scenario analysis and predictions at large scales are needed. Upscaled flow models can be developed, for example, through a process of replacing a heterogeneous domain by an equivalent homogeneous one that produces the same response to the upscaled boundary conditions imposed (Gelhar 1986;Rubin 2003). It is a simplification of a given problem that is strategically designed to allow for an appropriate yet efficient mathematical solution of the resulting model, while eliminating detailed information that does not affect the required accuracy of that solution. However, developing the spatial-upscaled-transport models that actually work in practice at basin-scale is a grand challenge among the hydrogeology community (Fogg and LaBolle 2006;Fogg and Zhang 2016;Soltanian et al. 2015a, b).
Besides spatial discretization, adequate representation of temporal variability in model stresses is a key challenge in groundwater flow and solute transport simulations (Blöschl and Sivapalan 1995;Jiang et al. 2019;Wallendar and Grismer 2002). Two main concerns regarding transient simulations are: (1) lack of time series associated with input flow-and-transportmodel boundary conditions, and (2) computational cost of the model. Groundwater-flow changes and contaminant-transport processes occur at time scales from seconds (e.g. geochemical interactions) to hours (e.g. precipitation, infiltration, changes in groundwater velocity, and transpiration) to days (e.g. well discharge). Transient boundary conditions used to define these changes and processes in a model include, but are not limited to, concentration rates of source/sink species, recharge rates, well pumping rates, and evapotranspiration rates. However, availability of information on these boundary condition data is usually restricted to monthly or annual measurements (or estimates) because monitoring data are often not collected more frequently. Furthermore, implementation of a highly detailed transient model is computationally expensive, rendering it less adaptable to support decision-making at regional scales (Turkeltaub et al. 2018).
For nonpoint-source management needs, upscaling issues arise that are not addressed in the existing upscaling literature, which has mostly considered point source pollution: at what temporal and spatial scale are continually discharging sources best resolved? The model domain of nonpoint source pollutants is much larger than that of point source contaminants (i.e., hundreds to thousands of square kilometers/miles instead of less than one to a few square kilometers/miles, and hundreds of impacted wells instead of a handful of impacted wells). Little work has been done to understand the degree to which nonpoint source pollution of groundwater can be upscaled in space or time.
The objective of this paper is twofold-first, to investigate the accuracy of temporal upscaling schemes for groundwater flow and transport boundary conditions in a numerical nonpoint-source transport model; and second, to test the temporal upscaling schemes specifically within the context of a decision-making tool for nonpoint source management. This is the first known work that attempts to extend available studies of managerial groundwater quality simulation methods by addressing these practical objectives for nonpoint source contaminants in agricultural areas at the scale of thousands of square kilometers. Here, the study evaluates the long-term impact of two nonpointsource management practices: nitrate load reduction (land use improvement) and agricultural-managed aquifer recharge, employed across a large regional agricultural landscape. Bastani and Harter (2019) showed the potential benefits of these practices when targeted at the source area of a specific production well. At the regional scale, recharge opportunities are limited in volume, but employing these nonpoint-source management practices may afford more wide-ranging water quality improvements in thousands of domestic wells.

Study area description
The study site is in the eastern San Joaquin Valley, in the Central Valley of California, USA ( Fig. 1), covering a region of about 2,700 km 2 surrounding two main cities (i.e., Modesto and Turlock), three main streams (i.e., Stanislaus, Tuolumne, and Merced), and bounded by the Sierra Nevada Mountains to the northeast and the San Joaquin River to the southwest. The land surface has a gentle slope (<0.5%) from east-northeast to westsouthwest, where groundwater discharges portions of the hydrological water to the San Joaquin River (Phillips et al. 2007).
The basin is composed of mostly unconsolidated fluvial deposits of medium-to-coarse grains with interlayered clay units. The aquifer system is described as unconfined in shallow deposits, becoming semiconfined with depth. The Corcoran Clay, a semi-permeable barrier, is a key hydrogeological feature that limits vertical groundwater flow below the base of this shallower portion of the aquifer system on which this paper focuses. In the study area, the top of the Corcoran Clay is between about 26 and 79 m below land surface (thickness of the unconfined to semi-confined aquifer overlying the Corcoran Clay), and the unit has a maximum thickness of about 58 m. A semi-confined aquifer exists below the Corcoran Clay layer in large parts of the San Joaquin Valley (Phillips et al. 2015).
Agriculture is the primary land use, covering about 65% of the Central Valley and about 80% of the study area. Agricultural land use is highly diverse, with hundreds of different crops grown, each with its own array of potential irrigation and nutrient management tools, adjusted to local soil conditions and agronomic considerations of individual landowners. The primary crops in the regional study area are almonds, peaches, grapes, grain, corn, pasture, and alfalfa. The cities of Modesto, Turlock, and a few smaller urban areas cover about 10% of the regional study area. The remaining area is predominantly natural vegetation near the foothills and in riparian areas- Fig. S1 of the electronic supplementary material (ESM).
Land use has varied over the decades throughout the study area ( Fig. 2) with almonds and corn making major gains in acreage. Increase in corn acreage (predominantly sileage corn) is due to an even larger increase in the dairy herd, milk production, and commensurate manure production, recycled mostly on corn fields (doublecropped with winter cereal). The increase in corn and almond acreage reflects significant intensification commensurate with increased land application of manure and commercial nitrogen fertilizers, respectively. The movement of water and solutes in the eastern San Joaquin Valley aquifer system is complex due to the heterogeneous distribution of alluvial fan sediments and the modification of the natural hydrologic system. Under natural conditions, recharge occurred at the mountain front and groundwater moved laterally until it discharged at streams near the center of the basin (southwestern boundary of the study area). However, diversion of surface water from streams and development of groundwater supplies significantly altered the natural flow system. Development of the groundwater basin began in the late 1800s. Groundwater pumping increased slowly until the 1940s and 1950s, when pumping increased significantly. Irrigation water (supplied from both surface water and groundwater) became the primary form of groundwater recharge, and irrigation pumpage became the primary form of groundwater discharge, causing water and solutes to also move vertically downward through the system .

Conceptual and numerical model overview
Groundwater flow setup Phillips et al. (2015) developed a three-dimensional (3D) calibrated transient (years 1960-2004) groundwater flow model for the regional study area using MODFLOW-OWHM (Hanson et al. 2014). The temporal dynamics in the boundary conditions are resolved at a monthly scale (i.e. 540 monthly stress periods). The model has 153 rows and 137 columns of square finite-difference cells with horizontal dimensions of 400 × 400 m. The aquifer is discretized into 16 layers vertically, with layer thickness increasing gradually with depth from one to 60 m. Layer eight represents the confining sub-Corcoran portion of the system. The model specifies general head boundary condition at all lateral boundaries. The exception is the northeastern boundary, where a no-flow boundary represents the contact between the alluvial aquifer system and the Sierra Nevada foothill bedrock complex. For the purpose of this study, the Phillips et al. (2015) model was used and converted to MODFLOW-2005 to be compatible with the MT3D platform for transport modeling (Bastani and Harter 2019). The same flow and transport model was used by Bastani and Harter to investigate the feasibility and longterm impacts of applying agricultural best management practices at the nitrate-contaminated source area of a pumping supply well to increase the quality of the drinking water at the well and the ambient groundwater. The study here uses these management concepts across larger scales while comparing the sensitivity of the solution to changes made in the time-resolution of the boundary conditions in the numerical model.
MODFLOW solves the 3D transient groundwater flow equation in a saturated heterogeneous porous medium (Harbaugh 2005): where K x , K y , and K z are aquifer hydraulic conductivities in the x, y, and z directions [L/T], h is hydraulic head [L], and S s is aquifer storage coefficient [L −1 ]. For steady state, the righthand side of Eq. (1) is zero. The heterogeneity of the study system was not represented in detail in this regional model. The sedimentary material hydrofacies was defined from the work of Phillips et al. (2007) based on the general approach of Laudon and Belitz (1991). Each hydrofacies represents a percentage of coarse grains in the sediment texture upscaled to   the size of the model cells. The sediment texture in each model cell was categorized into four materials based on the fraction of the coarse-grained materials available in that cell (Table 3). In this study, the temporal resolution of the boundary conditions to the transient flow model was upscaled to different levels. The highest frequency of available flow boundary conditions was a time series of monthly specified flows. The "reference" simulation takes advantage of the data resolution and uses monthly varying boundary conditions. Temporal upscaling ("scheme") simulations were defined by integrating monthly flow data over longer time periods (Fig. 3). Arithmetic mean averaging of monthly flows yields mass conservation at the larger time scale (e.g., averaging 12 months in a year or 540 months over the entire simulation period).
The highest level of temporal homogenization (upscaling), for conditions where long-term groundwater inflows and outflows are balanced, is a steady-state flow model. Two different approaches are used for developing steady-state boundary conditions. One alternative is to use a hydrologic year that represents balanced inflows and outflows over the study area. For example, Phillips et al. (2007) simulated a steadystate version of the study region using data from year 2000. Here, the year 2000 data from the transient model are used for upscaling flow boundary conditions to the steady state (scheme 6). For steady-state conditions, groundwater storage change over the study area must be negligible (inflowoutflow = storage change). In the transient model, actual total change of storage in year 2000 was -18 Mm 3 (Fig. 4), about 1% of the total pumped water during that year. For the steadystate flow simulation, outflows were proportionally scaled to match total inflows.
A second alternative approach to developing the steadystate flow model was to take the average of flow boundary conditions over the simulation period (i.e. 45 years) and import them to the model while solving the steady-state equation (scheme 5). The discrepancy between total inflow and outflow of the groundwater budget in this scheme after 45 years (Table 1) is less than 0.18% of the total pumping well extraction through the system. Again, a completely balanced inflow and outflow was achieved by proportionally scaling the pumping rate for all wells in the steady-state flow model. A total of six different upscaling schemes were developed and tested, four transient upscaling schemes, and two steady-state schemes ( Table 2).

Nonpoint-source-contaminant-transport simulation
The general form of the non-reactive groundwater contaminant fate and transport equation is (Zheng and Bennett 2002): where θ is porosity [unitless], C is the groundwater nitrate concentration [M/L 3 ], v is the actual pore velocity vector (equal to Darcy velocity divided by porosity) [L/T], D is the hydrodynamic dispersion tensor [L 2 /T], C s is the source or sink nitrate concentration [M/L 3 ], and q s is the source (e.g., recharge) or sink (e.g., discharging well) volumetric flow rate per unit volume of aquifer [T −1 ]. The right-hand side of the equation includes the dispersion term, advection term, and source and sink term, respectively. ∇ represents the gradient operation and ∇· represents the divergence. The MT3D (Zheng 1990) program was used to solve Eq. (2) for the 3D transient transport of nitrate in the region. A nonreactive transport model was developed because the denitrification reaction rate is assumed negligible in the aquifer due to its oxic conditions, although previous studies have shown denitrification occurs locally in some areas of this region where nonoxic conditions exist (Green et al. 2010(Green et al. , 2016McMahon et al. 2008).
Firstly, the porosity values were collected from literature review and distributed in the simulation cells considering the percentage of coarse-grained texture in each model cell categorized (Table 3). In the original flow model, each cell was associated with a number that represents the percentage of coarse sediment texture in that cell. Here, sediment size is assigned to specific fractions of coarse material (Table 3)-for instance, if the fraction of coarse sediments in a grid cell is more than 75%, it is classified as gravel and an effective porosity of 25% is assigned to that cell. This allows for appropriate representation of aquifer heterogeneity in both the flow model and the transport model. Secondly, the value for longitudinal dispersivity was estimated to be 11 m, using the empirical equation presented by Schulze-Makuch (2005) as where α is longitudinal dispersivity, c is characteristic parameter of the geological medium, L is the flow distance, which is the grid cell size in the model, and m is the scaling exponent.
The parameter values (m = 0.81, c = 0.085) were selected according to the type of the medium, which is unconsolidated sediments (Schulze-Makuch 2005). The ratio of transverse horizontal and transverse vertical dispersivity to longitudinal dispersivity is 0.1 and 0.01, respectively. Thirdly, initial nitrate concentration in the aquifer was estimated based on available distributed sampling-well data in the study area. Since there were not enough measured concentration data distributed in the study region in year 1960 to use as the initial condition, especially for depths above the Corcoran (layers 1-7), a combination of data distributed across years 1957, 1958, 1959, and 1960 was used for this purpose. The measured concentration data were categorized based on the perforated depth of the well screens (intervals of 30 m) and were assigned to the model layers. Ordinary kriging (Journel and Huijbregts 1978) was used to distribute the initial concentrations horizontally in each layer using available well location coordinates and the measured concentrations. Moreover, nitrate concentrations in streams were assumed to be negligible compared to other ground-surface spots in the model. The main source of nitrate in the area is nitrogen fertilizer. The nitrogen loads were equally applied over a period of 12 months on top of the first layer cells. In this study, it is assumed that the nitrate load travels through the vadose zone and reaches the water table in the same year when nitrate is lost from the root zone. Nitrate mass is not attenuated or delayed before percolating to the water table.
The next step is to define the nitrate mass load spatiotemporal distribution in the system. Harter et al. (2017) calculated the nitrate mass load for different crop types in the Central Valley using a mass balance method considering the amount of nitrogen fertilizer applied, harvested crop nitrogen, and other fluxes at 15-year intervals from 1945 onward (Table 4). For the purpose of these simulations, it is assumed nitrate load enters the groundwater system through deep percolation. Nitrate concentration (e.g. in g/m 3 ) of the percolated water is then calculated based on the volumetric recharge flow (e.g. in m 3 /day) and the rate of nitrate mass (e.g. in g/day) lost into that volume. The highest frequency of available nitrate mass data for individual crops in this study is every 15 years (Table 4) from which a linearly varying loading of nitrate was calculated to obtain varying annual nitrate loading.
The nitrate concentration in recharge at each field is determined by: (1) the agricultural field type, which may rotate at any given field over time, (2) a time-varying nitrate load associated with specific crops (Table 4), and (3) recharge fluxes that vary at monthly or temporally upscaled intervals. The specific combination of these three factors at each field creates a unique nitrate load and concentration history (Fig. 5).

Regional management scenarios
The six time-upscaling "schemes" are employed to several management "scenarios" designed to potentially improve groundwater nitrate in the study area. Four different scenarios were defined following the work of Bastani and Harter (Bastani and Harter 2019): (1) business as usual, (2) low-intensity crop, (3) winter recharge, and (4) low-intensity crop with winter recharge. In practice, the implementation of these practices is locally limited by: (1) availability of appropriate soil for agricultural groundwater recharge, (2) availability of excess stormwater in the area, and (3) suitable crops or management practices available for farmland (Fig. 6). The simulation period of the alternative management scenarios is the same as the reference scenario, from January 1960 to December 2004. The historic period for the alternative scenarios was being focused in this study to allow for an evaluation of the upscaling schemes against a well-defined reference case under various types of stresses. "Business as usual" (BAU) is the baseline scenario. Results of this scenario indicate how nitrate concentrations evolved in the region during the historic 45-year simulation period, given available data on historic agricultural management practices and land use distribution.
The "low-intensity crop" (LIC) scenario assumes a change of crop type or in management practices that significantly reduces nitrogen loading. The risk for nitrogen loading typically increases with the amount of nitrogen fertilizers applied, if management practices are not adjusted. Harter et al. (2017) estimated the nitrate leaching intensity through a mass balance approach for 58 different crops in the Central Valley. The benchmark of 35 kgN/ha/year for nitrogen loss to groundwater was used to distinguish between high-intensity and lowintensity crops. For instance, almond is a high-nitrogendemand crop (Table 4), but it is also capable of being recharged with stormwater. Alfalfa and vineyards are representatives of low-nitrogen-demand crops that are also suitable for recharge (Dahlke et al. 2018). In this scenario, all almonds were replaced by alfalfa (equivalent to treating almonds with highly efficient fertilization practices) from the beginning of the simulation or whenever land was historically converted to almonds with corresponding changes to the recharge and the nitrate concentration of those cells.
The third scenario is "winter recharge" (WR) which provides additional volume of water to farmlands during winter months of wet years. This scenario is also called agricultural managed aquifer recharge (Ghasemizade et al. 2019;Kourakos et al. 2019;Niswonger et al. 2017). It increases landscape recharge without increasing agricultural nitrate Additional recharge on almonds and alfalfa in the WR scenario is further restricted by four criteria: (1) floodwater that is made available only during winter months (i.e., November to April), (2) the volume of the floodwater available, (3) the maximum flood recharge amount tolerated by each crop, and (4) soil suitability for recharge. Availability of "high magnitude flows (HMF)" for on-farm recharge potential of groundwater banking within the Central Valley, California, has been assessed by Kocis and Dahlke (2017). The monthly total volume of available floodwater data for Stanislaus River, Tuolumne River, and Merced River were used as the source of water for "agricultural managed aquifer recharge" in the study area (Fig. 7).
A field's soil suitability for recharge was considered using the California Soil Agricultural Groundwater Banking Index (SAGBI; O'Geen et al. 2015). Soils rated "moderately good" to "excellent" were selected as suitable for winter recharge (Fig. 8). Due to historic landuse changes, the number of model cells with almond land-use on suitable soils gradually increases from 2,283 to 3,977 cells and decreases for eligible alfalfa lands from 1,037 to 455 cells during the simulation time.
The fourth scenario is "low intensity crop with winter recharge" (LICWR) which is the combination of the latter two scenarios. The LICWR scenario was developed based on the idea that the combination of the two practices further improves outcomes.

Comparative analysis
The performance of each scheme is evaluated by comparing the schemes' concentration results for each scenario with the reference simulation results at (1) shallow aquifer modeling cells, (2) domestic wells, and (3) agricultural and publicsupply wells in the study domain. A total of 28 simulations (1 reference plus 6 schemes, for 4 scenarios) were developed for the time-upscaling schemes and agricultural management scenarios in this study.

Shallow aquifer analysis
All active model cells in this study are classified as either "shallow zone" or "deep zone". Shallow zone is defined as the zone above the top elevation of the Corcoran Clay layer within the model which includes all cells of the first seven layers in the model. Nitrate concentrations in all the shallowzone cells at the last time-step of simulations are used as the sample population for statistical calculations to describe the degree of similarity between the reference and scheme outcomes in this study (Montgomery and Runger 2007).
Summary statistics include (1) the 99th percentile of nitrate in the sample population, (2) the arithmetic mean nitrate of the sample population, (3) the nitrate standard deviation of the sample population, and (4) the percentage of the sample population with concentrations below the nitrate MCL (i.e., 10 mg/L). Only cells with concentrations at and below the 99th percentile are used to calculate the mean and standard deviation to avoid outlier effects. To evaluate the accuracy of schemes in shallow groundwater, the root mean square error (RMSE) values to the reference simulation are calculated for all scheme simulations: The RMSE values are normalized (NRMSE) by the range of sample population data to further facilitate the comparison between different scenarios: where C s (x i ) is the nitrate concentration from the (upscaled) scheme simulation at the last time-step at location x i , C r (x i ) is the nitrate concentration from the reference simulation at the last time-step at location x i , P is the number of data points used for calculation, C max is the reference maximum nitrate concentration in the range of the sample, and C min is the minimum nitrate concentration in the same range of reference data. Residuals for each scheme simulation are defined as the difference between reference and scheme results (i.e., error = schemereference). The scheme simulation with the least RMSE is called "best fit" for the reference model. For the purpose of this study, this depth is being considered as the average depth of domestic wells in the study area. Layer seven of the model covers this depth; therefore, all active cells at layer 7 are each assumed to be occupied by one or several domestic well screens. Effectively this is a subset of the population investigated for the shallow aquifer analysis, representing the most likely location of domestic well screens. Breakthrough curves of nitrate concentrations at the depth of the domestic wells in the study region are compared between the reference and scenario simulations for the various upscaling schemes. It is assumed that each model cell represented all overlying households that extract their drinking water from a typical domestic well. For each time-step, exceedance frequency (%) is calculated as the ratio of the number of domestic well cells with nitrate concentrations above the MCL to the total number of domestic well cells in the study area (or subarea).

Agricultural and public-supply wells analysis
Agricultural-supply and public-supply wells are operating across the study area at various depths with discharging rates much higher than domestic wells. Using their historical information of pumping rates and their location coordinates, in the groundwater flow and contaminant transport model these wells are explicitly simulated with time-varying extraction rates. Nitrate breakthrough curves of all agricultural and public-supply wells at all layers of the scheme simulations are compared with the reference results at the scale of the regional study domain and the subregional water districts scale. Different percentiles of the ensemble of breakthrough curves are plotted as the statistical measurement for the purpose of the comparison.

Results and discussion
Historic nitrate contamination of groundwater in the BAU scenario Temporal dynamics of shallow groundwater nitrate in the reference (BAU) case is mainly controlled by the time-variant heterogeneity of the land use in the region. Changes in nitrate loading over time occur from crop intensification in general, but also from higher acreages of some of the crops having highest nitrogen losses. For example, almond acreage nearly doubles, rising from 13% to nearly one quarter of the total agricultural acreage in the study area, and corn, double-cropped with winter grain and both used as dairy forage, triples in acreage, from 3 to 9% (Fig.  2). Most crops see increased estimated nitrate leaching from 1960 through 1990, but corn and almond have the highest increases, more than doubling their estimated losses between 1975 and 1990, before improving nitrogen use efficiency after 1990 (Table 4). The combined impact of increased corn and almond acreage (Fig. 2) and increased corn and almond N losses between 1960 and 1990 (Table 4) accounts for nearly all of the total nitrate loading increase over that time period in the study region. The fraction of nitrate loading from corn and almonds increases from about half of all nitrate loading in 1960 to about 80% in 1990 (see Appendix). In contrast, estimates indicate no increase in N loading from alfalfa and a decrease in N loading from vineyards between the 1960s and the 1990s. In addition, their combined acreage decreases by nearly one half, from 12% in 1960 to only 7% in 2005 (Fig. 2). As a result, from 1960 to 2005, their contribution to N loading in the region decreases by over 70%.
In response to the nitrate loading dynamics, simulated mean shallow groundwater nitrate concentration under corn and alfalfa increases by 50%, between 1960 and 2005, from 6 mg N/L to 9 mg N/L. Even beneath the alfalfa or vineyard fields, nitrate concentration increases throughout the simulation period, but by less than 25%. The increase in groundwater nitrate underneath alfalfa fields and vineyards is partly due to their planting in areas with prior increases in N loading from other crops, and partly due to advection (and some dispersion) of nitrate from adjacent fields and orchards with higher nitrate losses. Elevated contamination at any sampling cell in the aquifer is the result of many years, even decades of increased loading and accretion of nitrate from the water table moving downward which is captured along well screens crossing multiple layers.
In domestic wells, deeper than the shallow-most aquifer region, the increase in nitrate is delayed. Simulated mean nitrate increases only slightly, from about 5.5 to 6.5 mg N/L, and no simulated increases occurs after the 1990s (Fig. 9). Deep production wells show a small decrease in mean concentration from the 1960s to the 1980s, from 4.8 to 4.2 mg N/L, then some increase into the 2000s (to 5.2 mg N/L). Across all groups of simulated measurement points (shallow groundwater, domestic wells, large production wells), there is significant variability among locations. By the 1990s/2000s, highest nitrate concentrations in all three groups exceed the nitrate MCL (10 mg N/L).
Simulated concentrations and their distribution are in reasonable agreement with measured nitrate data in the region, where available from public databases (Fig. 9). Simulated data are compared only at the specific locations (and depths) for which measured data are available. Overall trends and magnitude of measured nitrate concentrations for shallow groundwater, domestic, and large production wells are in good agreement with simulated values, showing the increasing trend in shallow and domestic well nitrate over the simulation period and the more stagnant nitrate concentrations in deeper production wells with only recent increases. The two most prominent differences between simulated and measured nitrate distributions are that measured data vary more widely than simulated data; and that both mean and upper quantile nitrate concentrations in the second half of the simulation period are measured to be significantly higher than simulated. In shallow and domestic wells, initial conditions appear to be higher than measured values which is mainly because of using a combination of available measured data from years 1957 to 1960 as the initial condition of the system. However, the trend of change in nitrate concentrations is the same for both measured and simulated values during the simulation period.
The presence of wider and more skewed distributions in the measured data would be consistent with the occurrence of non-Fickian transport processes (i.e., preferential transport ahead of the mean front and late-time tailing where nitrate comes out of the low-K layers by diffusion) and the ubiquitous effects of heterogeneity in the system such as the preferential flow through well-connected high-K zones not captured in the model grid scale, which produce a tail of high nitrate values. Also, geochemical heterogeneity, not represented in the model, may contribute to differences between measured and Moreover, the larger variability in measured data is thought to reflect between-farm and within-farm variability in loading due to soil-, irrigation-, and management-practice heterogeneity. This is not captured by the uniform loading assumed for each crop type (Table 4). Finally, the higher concentrations, especially in the 1990s and 2000s, stem from neglecting onfarm manure discharge, which is not accounted for in the cropcentric loading approach used here and which would have resulted in twice as much loading or more, on average across the study region, by 2005 (see Appendix). The comparison to measured data indicates that future simulation work must include that nitrogen source. For the purpose of this paper, the comparison to measured data indicates sufficient agreement of model results with actual data (Fig. 9), confirming that the overall conceptual and numerical approach taken is an appropriate basis for both, comparing to alternative management scenarios, and for evaluating the time-upscaling approach.

Effect of alternative scenarios on nitrate loading in the historic context
The alternative scenarios lead to various degrees in hypothetical groundwater nitrate improvement over the simulation period (compare the "Reference" cases among scenarios in Table 5).
The BAU scenario leads to shallow groundwater nitrate concentrations averaging just below the MCL, at 8.7 mg N/ L by 2004. During the 45 years of implementing the LIC scenario, continued degradation of shallow groundwater quality was halted and reversed. The average concentration above the Corcoran layer decreases by 1.1 mg N/L when all almond land uses are changed to alfalfa (or equivalent nutrient load reductions are achieved in almonds) and the number of exceedance cells in the model decreases by more than 10% at the end of the simulation. On the other hand, the winter recharge (WR) scenario has almost the same average nitrate value as the BAU. The frequency of available floodwater for the eligible lands is seasonally disrupted and also not available in all years. Implemented at the regional scale, the WR scenario is therefore not an effective tool to attenuate nitrate pollution, given that "agricultural managed aquifer recharge" in this WR scenario adds water to the system with irregular recharge frequency, in insufficient amounts to achieve significant dilution.
Similarly, the smallest fraction of wells in compliance with the MCL is observed in the BAU scenario (62%). Only slightly better results are observed for the WR scenario (64%). Importantly, the recharge does not lead to groundwater nitrate that is already stored in shallow groundwater to be flushed out. This is a concern that is sometimes raised by water purveyors that depend on groundwater. The simulations indicate that LIC does not lead to a faster increase in the nitrate trend observed for the BAU scenario. However, when WR is used in combination with LIC, recharge has a disproportionally larger effect: 71% MCL compliance for LIC, and 81% for LICWR (Table 5). The standard deviation in nitrate concentration, across wells, is similar between all scenarios, with a somewhat lower value for the LICWR scenario.
The main difference between nitrate dynamics in the shallow aquifer (as a whole) compared to domestic wells (assumed to be located in the lowest portion of the shallow aquifer) and public supply and irrigation wells is a delay in the nitrate impact observed, leading to the highest MCL exceedance frequency in 2004 for the shallow aquifer system. Because public supply and irrigation wells have longer and therefore deeper screens, exceedance frequencies are even BAU business as usual, LIC low-intensity crop, WR winter recharge, LICWR low-intensity crop with winter recharge lower in those than in domestic wells, consistent with observations (Fig. 9).

Performance of time-upscaling schemes at the shallow aquifer modeling cells
All upscaling schemes for all management scenarios have mean values close to their respective reference model which indicates that the selected upscaling schemes produce results that are statistically nearly identical to the reference schemes, for the shallow aquifer zone (Table 6). There is, however, a measurable but small increase in error with increased upscaling, the highest being for the steady-state upscaling scheme. This increased error occurs regardless of the type of management scenario. From a practical perspective, the somewhat higher simulated error appears of little consequence. For example, predicting the percentage of wells below the MCL, the maximum difference to the reference model is less than 7%, which occurs in the LICWR scenario with scheme 6 (steady state with year 2000 data). Even the extreme values of the 99th percentile nitrate value for all six upscaling schemes of scenario BAU are in close agreement with the corresponding value in the BAU reference model. They vary from 3 to 9 mg/L for the alternative management scenarios. It is always highest for scheme 6. In all management scenarios, the NRMSE increases from scheme 1 to scheme 6, i.e., with increasing temporal aggregation. The NRMSE values are larger for the BAU and WR scenarios when compared to the LIC and LICWR scenarios (at the last time step of the simulations). However, even in the steady-state simulation, the LIC and LICWR scenarios produce results closer to the reference case than in the BAU. In all upscaling schemes the overall error remains small, given greater uncertainty associated with the reference model inputs (see Fig. 9).

Water quality at the depth of domestic wells in the upscaled schemes
Capturing appropriate exceedance levels in domestic wells is critical in decision support tools. Upscaling schemes 1, 2, 3, and sometimes even scheme 4 were found to have nearly identical exceedance frequencies to the reference model, across all management scenarios. The model exceedance values in the steady-state schemes are slightly higher than in the reference model, but with less than 5% difference to the reference case, at the end of the simulation period (WR scenario). For the LIC and LICWR scenarios in particular, the steady-state flow models are in close agreement with the transient scenarios and the reference models (Fig. 10).
The temporal dynamics of the least accurate model (scheme 6) is compared in more detail with the reference cases. The temporal development of the median (50th percentile) as well as the 90th percentile nitrate concentration across the ensemble of domestic wells vary significantly over time, increasing for the BAU and WR scenarios, but mostly decreasing (after an early peak) for the LIC and LICWR scenarios. Even under steady-state flow boundary conditions (scheme 6), these temporal dynamics are well captured and the differences between management scenarios well articulated (Fig. 11).
Although not identical, the steady-state scheme closely follows the reference case dynamics (within about 10% of the reference value at any given time). Some short-term (less than 5 year) temporal variations expressed in the transient reference case are smoothed out by scheme 6 upscaling. At the end of the simulation period, the results of scheme 6 for BAU and WR overestimate the transient model concentration values with a discrepancy by 5%. Given the uncertainty associated with the reference model inputs, such errors appear rather acceptable from a decision-making perspective for domestic wells, indicating that temporal upscaling is appropriate for these types of nonpoint-source transport analyses. Water quality at the agricultural and public-supply wells in the upscaled schemes At even larger depth and after longer transport distances, nitrate reaches agricultural and public supply wells. Again, median and 90th percentile nitrate concentration is compared across irrigation and public-supply wells in the study region.
In the deeper irrigation and public supply wells with large pumping rates and long screens, median nitrate slightly increases over the first 10 years from its initial value of 4.6 mg/L, but then begins to decrease during the second decade, before leveling off and again increasing at 30 years through the end of the 45-year simulation period. The (reference) WR scenario leads to only slightly lower medians and follows a similar dynamic. In contrast, both LIC and LICWR scenarios continue the decreasing trend of the second decade, albeit at a decelerating trend, with LICWR medians about 0.4 mg/L lower than LIC, but almost 50% less than BAU at 45 years (Fig.  12). Similarly, the 90th percentile nitrate concentration increases from its initial value and continues a decreasing trend for the early period. Only after 20 years is the 90th percentile beginning to increase, quite quickly relative to the increase in median concentration. For the 90th percentile nitrate, differences of the WR, LIC, and LICWR scenarios to the BAU scenario are similar to those for the median. For the 90th percentile, the discrepancy between estimated nitrate concentration curves of scheme 6 and reference scenarios increases the most during the last decade (Fig. 12). The difference may be related to the dynamics of boundary conditions, especially recharge rate. More divergence of results is seen during highly transient conditions in the system. Change in recharge not only influences the flow direction, but also changes the nitrate concentration imported to the system (See Fig. 5). In this study, nitrate concentration has a generally declining trend for the first 15 years of the simulation following with an increasing and another decreasing pattern for the second and the third 15 years, respectively. However, the input nitrate concentration illustrates a sequence of sudden falling and rising shifts during the last decade of the simulation. The 90th percentile graph of the reference model manifests the change in the most recent 10 years (Fig. 12), but the steadystate flow model adds nitrate concentration to the system based on only one steady recharge value which causes a continuous increase of concentration at discharge locations.
The simulated breakthrough curve of the scheme 6 scenario ( Fig. 12) is consistent with previous nitrate concentration estimations in the study area. Green et al. (2016) applied a steady-state flow system (as estimated with agedistributions) combined with time-varying inputs in a datadriven approach to estimate historical trends of nitrate in this study domain considering regional oxygen reduction and denitrification rates. The late-time slope of the median nitrate concentration breakthrough curve in deep wells in their study ( Fig. 9 in Green et al. 2016) is about 1.3 mg/L per decade, which is similar to about 1 mg/L per decade in the current study (Fig. 12, 50th percentile, scheme 6 BAU).
Relative to these reference simulations, the upscaling schemes, even at this depth, and for wells with long screens, produce comparable results, although not quite as accurate as in the domestic well case. Again, scheme 6 produces the least accurate results. It overestimates reference concentrations for much of the simulation period, in all management scenarios. The overestimation of reference medians and 90th percentiles is larger for the BAU and WR scenarios than for the LIC and LICWR scenarios and increases over time. Median concentrations are overestimated by nearly 10 and 5%, respectively, whereas the 90th percentiles are overestimated by 15 and 10%, respectively. However, the steady-state model accurately captures the longer-term dynamic behavior of nitrate across these wells and appropriately differentiates the impact of management scenarios.
Distribution of the irrigation and public-supply wells in the study region is different for the eight subregions delineated by water district boundaries (Fig. 13). The 90th percentile of breakthrough curves of discharging wells are next evaluated at the subregional scale (Fig. 14). Comparing results of transient reference models at the two different aggregation scales (Figs. 11 and 14), the regional scale 90th percentile Fig. 12 Time series of water quality change in large production wells: a 50th percentile nitrate and b 90th percentile nitrate concentration timeline is much smoother than at the subregional scale, where the smaller sample population generates many sudden jumps and falls in this high percentile curve.
The steady-state model does not capture these oscillations of the 90th percentile breakthrough curves. Differences to the reference case are as much as 40-50%. It is not completely accurate to state that the steady model overestimates concentrations at the subregional scale, since the 90th percentile concentrations can sometimes be much lower than the reference value. Still, scheme 6 captures the general trend of this extreme value concentration change through time and recognizes the impact of management scenarios in a long-term analysis. The results indicate that the temporal upscaling schemes, while accurate for measures of the mode, mean, or median of the nitrate distribution across a group of wells, are not as reliable, but perhaps still useful, in predicting extreme values across a group of wells, including a smaller, subregional group of wells. However, even the reference model is not expected to predict actual extreme values well (see Fig. 9).
The RMSE between the reference model and the steadystate model captures the accuracy aggregated across all simulated wells (and cells in the case of the shallow aquifer), for all eight subregions (water districts), using the BAU scenario. Over time, the RMSE increases, although there is a decrease at early times, perhaps owing to errors in the initial concentrations. (Fig. 15). Each subregion has its own minimum and maximum concentrations through the simulation time. The minimum value of concentration is approximately zero in all Insights on the steady-state-flow-based-nonpointsource-transport modeling Few studies have discussed the impact of transientgroundwater-flow-model-boundary conditions, in particular the typical seasonal changes of pumping rates, on variation of nonpoint-source-contaminant fate-and-transport modeling results for domestic or larger production wells. The focus of previous works (Bexfield and Jurgens 2014;Chen 2010;Jeffrey et al. 2014;Yager and Heywood 2014) has mainly been on contaminants that originated from a point source. Libera et al. (2017) investigated the impact of temporally variable groundwater extraction rates on the solute concentration history at wells in a heterogeneous field. They showed that uncertainty associated with detected breakthrough curves at the well increases with the increase of pumping pattern oscillations. The findings concur with their conclusion when the steady-state flow model is used at a small subregional scale, even if the source is a nonpoint source. It was observed that ignoring detailed time-variation of flow boundary conditions in a steady-state model produces breakthrough curves that cannot track the immediate oscillation of a transient concentration change at the subregional scale (Fig. 14).
As the scale of the study area increases to many hundreds of square kilometers, the number of wells sampled encompasses hundreds rather than tens; hence, the percentile breakthrough curves resulting from transient and steady-state-flowbased transport models are more similar (Figs. 11 and 12). In fact, the steady-state-flow-based transport model provides a plausible prediction for the long-term variations of the concentrations across three sets of discharging points, even at the subregional scale (shallow aquifer, domestic wells, large production wells).
Moreover, change of the velocity field under transient flow conditions makes the solute dispersion development irregular and complicated, moving solutes into different directions during the simulation time. The transverse dispersion, specifically, controls the mixing and natural attenuation of the contaminant plume (Cirpka 2005;Rolle et al. 2009;Schirmer et al. 2000); however, the spreading of the plume under steady-state flow conditions in a heterogeneous media, in point-source contaminant transport, has larger dispersion in the longitudinal direction than the transverse (Cirpka and Attinger 2003). On the other hand, studies have shown that considering the detailed heterogeneous structure of the porous media has a more dominant impact on transverse mixing than temporally changing flow fields (Elfeki et al. 2012;Rolle et al. 2009). In this study, the same value of transverse dispersion was used in all transient and steady-state-flow-based transport models and explicitly simulated aquifer heterogeneity, so that the lateral dispersion of nitrate in the system has minor impact on the arrival concentration levels at the well. Finally, a significant fraction of variability in local concentration is lost when nonpoint-source pollutants are extracted from production wells, due to the large mixing in the well itself, while the pollutants occur at a relatively narrow concentration range compared to many point source pollutants (Henri and Harter 2019).
Transient flow changes the spreading direction and magnitude of the plume in the aquifer. The concentration level observed at targeted wells reflects the intensity of the stress changes in the groundwater flow system at each stress period of the model. Utilizing the steady-state flow does not capture sudden fluctuations in the concentration levels. Instead, this model often overestimates the simulated concentrations obtained as compared to the transient flow and transport model. Because of changes of flow conditions over time for the latter (i.e., recharge rates, pumping rates, evapotranspiration rates, and storage change) the transient model enhances contaminant dispersion (more so in a heterogeneous media), causing concentrations at different wells to be lower than those obtained from the steady-state flow model. The over-or underestimation tendency of concentrations at targeted wells in the steady-state-flow-based transport model is a site-specific behavior which depends on the variability of loading and boundary conditions in time and space. However, the steady-state-flow-based-transport model reflects the impact of different agricultural management practices on the resulting breakthrough curve concentrations, which is a critical piece of information for decision-making in water quality management.
LaBolle and Fogg (2001) observed the effects of transient flow systems on contaminant transport in a detailed modeling of plume-scale phenomena in an alluvial aquifer where the effects of the storage of contaminant mass in the finegrained materials mainly became important when the boundary conditions were changed from steady and forward plume movement, to the transient pump-and-treat remediation conditions. The current study does not resolve local-scale finegrained, low-permeability inclusions and their impact on the tail of breakthrough curves at discharging points. However, unlike most point-source-groundwater-pollution sources, nonpoint-source pollution will not be removed, rather its leaching rate will be gradually lowered. Hence, tailing behavior of plumes, as in remediation investigations, is not critical to decision support.

Central Valley groundwater quality management implications
For over a century, the primary landuse in the Central Valley of California (52,000 km 2 ) has been agricultural production, with increasing intensification. Groundwater nitrate frequently exceeds the regulatory standard of 10 mg N/L. The aquifer system of the Central Valley is mainly composed of alluvial sediments, with most extraction in the upper 300 m ). The climate is semiarid with transient patterns of dry and wet water years. Water users in the region rely on a combination of surface-water delivery and groundwater extraction. The hydrologic system has been simulated with 3D monthly based transient regional groundwater/ surface-water flow models using the finite-difference platform of MODFLOW (Faunt 2009) and the finite-element model IWFM (Brush et al. 2013), with spatial resolution of over one to several kilometers. Run-times for these models are on the order of hours. Development of a transport model for the Central Valley, especially for nonpoint source contaminants such as nitrate, would require additional spatial resolution in the flow model plus a transport model, multiplying computational expenses. Evaluating a single nitrate management practice scenario at the scale of the Central Valley would be an enormous computational expense not suitable for a decision-support tool. Therefore, using a flow model with upscaled time-resolution that provides accurate largescale trends of nitrate at significant computational savings may be a more promising decision support tool for evaluating impacts of different land use and water recharge scenarios in the Central Valley. Kourakos and Harter (2014) developed a transport modeling approach that generates a detailed 3D steady-state flow field at high (<100 m) resolution with one-dimensional (1D) transport equation solved along relevant flow streamlines (to wells and streams, considering longitudinal dispersivity, but neglecting transverse dispersivity). Since the flow field is solved only once and the 3D transport equation is decomposed into different 1D transfer function equations, it is a time-efficient nonpoint-source assessment tool for regional scale analysis. This study indicates that an approach based on a steady-state groundwater flow model with transient nitrate load is a promising approach to simulate the general pattern of breakthrough curves for management purposes at basin scales of the size of California's Central Valley. However, future research is needed to assess the loss of accuracy due to neglecting transverse dispersivity.

Conclusions
A main challenge in developing detailed nonpoint-source transport models is reliance of these models on highfrequency time-series of flow boundary conditions. This paper investigates the plausibility of developing accurate nonpoint-source-contaminant-transport models, albeit at coarser temporal resolution (or lower frequency of timeseries data), for flow stresses, which are, hence, computationally more efficient. Upscaling flow boundary conditions in time allows for computational savings also in transport modeling, thus paving the way for tools capable of scenario analysis and decision support with regionalscale analysis. In this study, various upscaled time discretizations were considered and the degradation of the solution accuracy relative to the reference simulation observed experimentally. For the upscaling analysis, a set of highly relevant decision-support scenarios were employed which were designed to assess the effects of crop management and winter recharge on groundwater quality across a large region. This work therefore provides insights both, on nonpoint source attenuation through management practices, and on computational algorithms to capture those future outcomes accurately and yet efficiently.
Among the management practices, "winter recharge (WR)" or "agricultural managed aquifer recharge" is not effective if it is applied without nutrient management improvements, across the regional scale. Its water quality benefits strongly depend on the area suitable for winter recharge (i.e., soil characteristics and crop type) and on the amount of available recharge water. The application of recharge on the most favorable sites (e.g., incisedvalley-fill deposits) might prove to have greater and faster improvements in water quality (Maples et al. 2019). On the other hand, improving nutrient management (e.g., changing crop type, LIC) alone or together with winter recharge (LICWR) leads to significant attenuation of nitrate across the regional agricultural groundwater basin. Management practices related to winter recharge are best used for targeting high priority polluted lands, with suitable winter recharge lands, rather than distributing winter recharge across the region. Focusing winter recharge on key areas may provide significant benefits to specific public supply wells that would otherwise need a treatment plant (Bastani and Harter 2019). Fig. 15 Second-order polynomial trend line of RMSE between "Reference" and "Scheme 6" nitrate concentrations at discharging wells in the BAU scenario for different subregions and for the study area Among the various computational time-upscaling schemes to develop these insights, even the coarsest temporal upscaling (a steady-state flow simulation) provides sufficient accuracy if the load of the nonpoint source contaminant is known as a time-dependent spatial function. Many studies have used the assumption of steady-state flow to apply transfer-functionbased approaches to analyze trends of nitrate and other nonpoint source solutes in a system (Green et al. 2018;Kennedy et al. 2009;Massoudieh et al. 2014;McCallum et al. 2014;Osenbruck et al. 2006;Jeffrey et al. 2014). Although not identical to the temporally fully resolved results, the model can predict sufficiently accurate long-term nitrate concentration changes. Aggregated results, such as the regionally averaged breakthrough curve of nitrate at public supply wells, are particularly well predicted even when transport simulations were based on steady-state flow (least time resolution); however, there are some limitations. Sudden changes in flow-system boundary conditions may lead to some rapid nitrate increases or decreases that are not captured when employing a steadystate flow simulation. For more exact predictions of the dynamic changes in contaminant concentrations at a particular well, a fully transient groundwater flow model, given the appropriate knowledge of boundary conditions, is a superior assessment tool, if at significantly higher cost.
Temporal variation of groundwater level, flow directions, and velocity is mainly due to the combination of temporal changes of boundary conditions in the flow system and the distribution of the spatial heterogeneity of the aquifer parameters, which eventually affects the path and the time that the solute spends in the porous medium. LaBolle and Fogg (2001) observed the effects of transient flow systems on contaminant transport in a detailed modeling of plume-scale phenomena in an alluvial aquifer where the effects of the storage of contaminant mass in the fine-grained materials mainly became important when the boundary conditions were changed from steady and forward plume movement, to the transient pump-and-treat remediation conditions. In the current study, it was not addressed that whether change of flow boundary condition from transient to the steady state pushes the nitrate that got historically stored in the typically abundant aquitard sediments in the porous media to come back out of those sediments in sufficient concentrations. Instead, the focus of the study is on the concentrations of the discharging points at large scale.
Finally, these findings are encouraging for further work on highly efficient decision support systems for use by stakeholders and policy makers. There are opportunities to build more efficient, yet accurate modeling frameworks for nonpoint-source assessment tools of nitrate and salinity at regional scales. However, future research is needed to understand the effects of upscaling time-resolution of temporally fluctuating groundwater flow on reactive transport of nitrate (i.e., denitrification process) in heterogeneous porous media with transient sources of nitrogen leaching into a nitrate-reactive aquifer.
This would also apply to other reactive nonpoint source pollutants such as pesticides, pharmaceutical compounds, and other emerging contaminants. Moreover, future investigations are needed to address the effects of uncertainty of model input parameters on the nonpoint-source-contaminant-concentration levels in a model with transient-flow boundary conditions compared to steady-state conditions.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. Fig. 16 Major changes in nitrate load used in the model due to land-use variation in the study area