A Spatially Distributed, Physically-Based Modeling Approach for Estimating Agricultural Nitrate Leaching to Groundwater

: Nitrogen-nitrate, while being fundamental for crop production, is of particular concern in the agricultural sector, as it can easily leach to the water table, worsening groundwater quality. Numerical models and Geographic Information System may support the estimation of nitrate leaching rates in space and time, to support sustainable agricultural management practices. In this paper, we present a module for the simulation of the processes involved in the nitrogen cycle in the unsaturated zone, including nitrate leaching. This module was developed taking steps from the ANIMO and EPIC model frameworks and coupled to the hydrological models integrated within the FREEWAT platform. As such, the nitrogen cycle module was then included in the FREEWAT platform. The developed module and the coupling approach were tested using a simple synthetic application, where we simulated nitrate leaching through the unsaturated zone for a sunﬂower crop irrigated district during a dry year. The results of the simulation allow the estimation of daily nitrate concentration values at the water table. These spatially distributed values may then be further used as input concentration in models for simulating solute transport in aquifers.


Introduction
Nitrogen-nitrate, while being fundamental for crop production, is of particular concern in the agricultural sector, as it can easily leach to the water table, worsening groundwater quality.Nitrate pollution of water is a major issue worldwide (Mateo-Sagasta et al., 2018; Wang and Li 2019; Zhang et al., 2019) [1][2][3].It originates mostly from the use of inorganic fertilizers and animal manure in agricultural areas (Liao et al., 2012;Nakagawa et al., 2017;Ward et al., 2018;Zhang et al., 2018a) [4][5][6][7], and inefficient management of treated/untreated wastewater (Minnig et al., 2018) [8] and solid waste disposal (Wakida and Lerner 2005) [9].Focusing on groundwater, nitrate contamination is a worldwide problem (Re et [10][11][12][13].In Europe, policies to reduce NH 4 + and NO 3 − losses have been in place since nearly 30 years following the issue of the Nitrates Directive (ND; EC 1991) [14].The ND objectives have been to reduce water pollution caused or induced by nitrates, and to prevent further pollution.According to the ND, Member States are asked to identify Nitrate Vulnerable Zones (NVZs), i.e., zones with eutrophic water bodies or at risk of eutrophication, with nitrate concentration in ground-or surface-water exceeding 50 mg/L (NO 3 − regulatory limit for drinking water; Wild et al., 2018 [15]).Anyway, Wild et al., (2018) [15] argue that it will take several decades to significantly reduce nitrate concentrations in porous aquifers via denitrification, even if future nitrate inputs were considerably reduced.
The ND was followed in 2000 by the Water Framework Directive (WFD; EC, 2000 [16]), whose goal is to reach a good chemical and ecological status of water bodies by 2027.The purposes, objectives, and strategies of the WFD overlap with those of the ND, especially after the introduction of the Groundwater Directive (EC, 2006 [17]) in 2006.In the attempt to target the objectives of the above-mentioned directives, River Basin Management Plans envisage approaches to assess nitrogen occurrence and to evaluate the effectiveness of measures to reduce nitrogen content in groundwater.Such approaches include monitoring of groundwater quality in NVZs (e.g., Jørgensen et al., 2009 [18]) and evaluation of nitrate leaching through the unsaturated zone and nitrate fate and transport in surface-and ground-water by means of modeling (e.g., Pistocchi et al., 2012;Klement et al., 2018;Husic et al., 2019 [19-21]).
The nitrogen cycle (Figure 1) involves physical, chemical, and biological processes which determine an increase or decrease of nitrogen content in the soil and the vadose zone, as well as nitrogen transformation among organic nitrogen (NH 2 − ), ammonium (NH 4 + ) and nitrate (NO 3 − ) pools.Three sub-cycles may be distinguished: • the sub-cycle occurring in the atmosphere; • the sub-cycle occurring in the unsaturated zone; • the sub-cycle occurring in the plants, which involves NO 3 − uptake by roots.Sources/sinks contributing to the three nitrogen pools also play relevant roles (e.g., deriving from surface runoff, fertilizers, etc.).A detailed review on the N cycle in the vadose zone-groundwater system may be found in Xin et al., (2019) [22].
Several models are available to simulate watershed/field scale hydrological processes.Some of these models also include methods to simulate the nitrogen cycle in soils and nitrogen leaching, the downward transport of N passing the root zone by percolating water through the unsaturated zone (Golmohammadi et al., 2014 [23]).Appendix A reports the main characteristics of some of these models, in terms of the modeling approach adopted (lumped or spatially distributed), the scale of application, and space and time discretization.Details about the mathematical approach are reported as well, focusing on the hydrological and the nitrogen cycling components only.
Most of the codes listed in Appendix A include a hydrological component, except for ANIMO (Groenendijk and Kroes 1999 [24]), which can be coupled to external models (e.g., SIMGRO, Querner and Van Bakel 1989 [25]; SWAP, Van Dam et al., 1997 [26]) for the simulation of water fluxes and moisture contents.In ANIMO, such components are then algebraically summed to get the change in water volume over time, to be used as input in the nitrogen cycling component.In EPIC (Sharpley 1990 [27]), APEX (Williams et al., 2015 [28]) and SWAT (Neitsch et al., 2005 [29]), belonging to the CREAMS (Knisel 1980 [30]) family of models, the hydrological component is based on algebraic mass balance equations involving rainfall, runoff, and evapotranspiration and other hydrological processes.On the other hand, other codes (i.e., DAISY, Hansen et al. 1990 [31]); PRZM-3, Suárez 2005 [32]; and NTT-Watershed, Heng and Nikolaidis 1998 [33]) solve the Richards' equation for soil water dynamics in the unsaturated zone.
Regarding the nitrogen cycling component, all the codes reported in Appendix A integrate equations for the simulation of the processes involved in the nitrogen cycle.DAISY and PRZM-3 simulate NH 2 − mineralization in a very detailed way, including different NH 2 − sub-pools and resulting in large-parameterized models.EPIC, APEX, and SWAT use soil-dependent decay coefficients for the simulation of nitrogen leaching through the unsaturated zone, while ANIMO, DAISY, PRZM-3, and NTT-Watershed solve 1D advection-dispersion solute transport equations in the unsaturated zone.
Nitrate leaching research grew in the scientific community in the last decades (Padilla et al., 2018 [34]), focusing also on modeling methodologies (Canter 2019 [35]).Indeed, assessing nitrogen transport and turnover in the unsaturated zone and in groundwater is an expensive and complex task to be run at field scale (Sutton et al., 2011 [36]).At the same time, advances in distributed modeling efforts that incorporate spatial distribution of nitrogen inputs, soil and aquifer characteristics, land use data, and farming practices are promising (Xu et al., 2013;Ransom et al., 2017 [37,38]).While several modeling studies focused on simulating soil nitrogen balance or groundwater nitrogen transport independently (Wriedt 2004 [39]; Epelde et al., 2016 [40]), there is a clear trend to combine unsaturated zone and groundwater flow models for the simulation of nitrogen dynamics [41][42][43].
As such, in the following sections, we present a new modeling approach in order to compute nitrate leaching to the water table based on (a) the development of a numerical module for the simulation of the nitrogen cycle in the unsaturated zone, (b) the integration of such module in the Geographic Information System (GIS)-integrated FREEWAT platform (Foglia et al. 2018 [44]).The module was developed taking steps from the ANIMO and EPIC lumped models.Once developed, the module was then coupled to physically-based, spatially distributed codes for the simulation of unsaturated zone flow and agricultural water management and, finally, included into the GIS-integrated FREEWAT platform.

Overview of the Nitrogen Cycle Module
The nitrogen cycle module here developed allows to simulate the processes determining an increase or decrease of nitrogen content in the unsaturated zone, as well as nitrogen transformation among the NH 2 − , NH 4 + , and NO 3 − pools (Figure 1).The module developed is based on the conceptualization of the ANIMO model [24].ANIMO simulates the nitrogen budget for the three pools, by solving a partial differential equation for each of them (Table 1).All terms and symbols listed in Table 1 are detailed in Appendices B and C. Figure 2 reports a sketch of all the processes related to nitrogen mass flux in the unsaturated zone, occurring within a grid cell.By solving the equations reported in Table 1, the variation of nitrogen concentration for each pool in time and space, in the unsaturated zone, is calculated.
Table 1.Partial differential equations solved by the ANIMO model for the three nitrogen pools (modified after Groenendijk and Kroes 1999 [24]).

Crop Uptake Term Decomposition Term
Production Term The following terms appear in the equations presented in Table 1: • a conservation term, which expresses the variation of nitrogen concentration in time; • the conservation term for NH4+ includes also sorption by soil (ρ d ∂X e,NH4 (t) ∂t ); • a vertical, unidimensional transport term, which is driven by unsaturated zone flow; • a term related to lateral movement of nitrogen in the unsaturated zone (lateral outflow term).However, in our approach, we assume that vertical transport of nitrogen is the predominant process.As such, lateral movement of nitrogen is not simulated; • a crop uptake term for NO 3 − and NH 4 + .Such term is estimated for NO 3 − only, adopting the EPIC model approach (Sharpley 1990) [27], taking into account the phenological phases of crops.We assume, indeed, that NO 3 − is the predominant form of available nitrogen uptaken by roots (Nadelhoffer et al. 1984 1) and summed as production term for NO 3 − (R_(p,NO 3 ) in Table 1).
The equations reported in Table 1 are solved one by one, in the order presented above.Each of them can be rewritten as (subscripts NH 2 , NH 4 and NO 3 will be omitted in the following): By solving Equation (1), concentration in time (c(t)) for each pool is calculated in the unsaturated zone at the time t.
We decided to take steps from the ANIMO, among the models listed in the Appendix A, because we considered it a good compromise with respect to the following aspects: the mathematical approach adopted, based on the mass conservation and the solution of transport Equation (1); -the mathematical approach adopted is relatively frugal, in the sense that it requires few input parameters; -the integration between Equation (1) and the unsaturated zone flow term calculated by means of distributed models integrated within the FREEWAT platform was pretty straightforward with respect to the space and time dimension of the involved processes.

Overview of the FREEWAT Platform
FREEWAT [47] is a free and open source platform integrated in a GIS framework for water resource management, with a specific focus on groundwater, initially developed within the HORIZON 2020 FREEWAT project (FREE and open source software tools for WATer resource management; www.freewat.eu;[48]), a follow up of the SID & GRID platform [49].
FREEWAT is conceived as a QGIS [50] plugin, which integrates free and open source software tools for pre-and post-processing of spatial data, and numerical codes for the simulation of the hydrological cycle.On one side, GIS and spatial database management applications [51] are exploited for spatial data storage, visualization and analysis.Physicallybased and spatially distributed numerical models (i.e., MODFLOW, Harbaugh 2005 [52], and MODFLOW-related models) integrated within the FREEWAT platform allow to simulate several processes involved in the hydrological cycle (e.g., groundwater flow, densitydependent flow, conjunctive use of ground-and surface-water, solute transport, sensitivity analysis and model calibration).Dedicated tools for pre-processing of field data [53,54] and for post-processing of models' results are integrated as well.FREEWAT was applied to date to a variety of case studies [55].The present work may be considered as a continuation of the research activities performed to better describe agricultural water use in hydrological modelling.In particular, the Farm Process (FMP; embedded in the USGS's MODFLOW One-Water Hydrologic Model [56]) coupled to a Crop Growth Model (CGM [57]) for simulating crop yield based on solar radiation and biomass accumulation was previously integrated within the FREEWAT platform.

Coupling the Nitrogen Cycle and the Hydrological One
The nitrogen cycle module was coupled to the unsaturated flow process, in order to simulate advective transport of nitrogen compounds in the downward direction, across the unsaturated zone.According to the ANIMO approach, the nitrogen concentration in the outflowing water, moving from one unsaturated zone layer to another, is equal to the calculated nitrogen concentration within the ith unsaturated zone layer [24].In our approach, by activating the Unsaturated Zone Flow (UZF) MODFLOW package [58] in a groundwater flow model based on MODFLOW-2005 [52], we simulate a fraction of the total nitrogen mass within the unsaturated zone that reaches the water table.The nitrogen mass flow is driven by vertical water flow through the unsaturated zone.The UZF solves a kinematic wave approximation of the Richards' equation to simulate vertical flow and storage through the unsaturated zone [58].Conversely from most of the lumped available codes simulating the nitrogen cycle (e.g., ANIMO and EPIC), the codes integrated in the FREEWAT platform for the simulation of the unsaturated zone are spatially distributed.Hence, the partial differential equations are solved to get the involved variables (e.g., hydraulic head and solute concentration) also in a discretized space.
The coupling approach includes agricultural water management processes for the representation of inorganic nitrogen uptake by plants' roots.To this scope, the coupled FMP-CGM approach developed in FREEWAT for agricultural water management and crop yield simulation is adopted: in the nitrogen cycle module, input and output of an FMP-CGM scenario are used to calculate NO 3 − crop demand, as conceptualized in the EPIC model [27].Specifically, NO 3 − crop demand over the different phases of the cropping season is calculated in EPIC as the difference between the NO 3 − crop content, which is the nitrogen amount in plant biomass, and the NO 3 − availability in the soil.Crop transpiration fluxes simulated by FMP and cumulated biomass calculated by CGM are involved in the computation of these two variables.
All of the processes involved in the nitrogen cycle were simulated in the unsaturated zone layer at each cell of a finite-difference grid.This allows to obtain the spatial variation of nitrogen concentration for each nitrogen pool in the unsaturated zone.
The coupling approach just presented is drawn in Figure 3.The left side of Figure 3 (blue boxes) reports the codes that must be run before running the nitrogen cycle module.Two options may occur: if crop NO 3 − demand is simulated over the cropping season (Option 1), the following is needed:  if Option 1 occurs, all of the processes involved in the nitrogen cycle are simulated taking steps from the ANIMO model approach, except for the NO 3 − crop uptake process, which is simulated through a sub-routine based on the EPIC model approach; -if Option 2 occurs, all the processes involved in the nitrogen cycle are simulated taking steps from the ANIMO model approach.In such case, the sub-routine based on the EPIC model approach is not executed, and the NO 3 − crop uptake process is simulated by comparing the crop NO 3 − requirement defined by the User with the NO 3 − availability in the unsaturated zone.
In both options, input/output parameters involved in the execution of the UZF, FMP, and CGM components are exploited to run the nitrogen cycle module (i.e., to solve Equation ( 1)).As a result, the NO 3 − concentration at the water table is determined at each cell and for each stress period.This concentration can then be treated as a potential contaminant source term to the saturated zone in solute transport modelling of aquifers (e.g., through the application of the MT3DMS code [59]), as indicated in the green box on the right side of Figure 3.

Model Testing Using a Synthetic Case Study
The coupling approach discussed in this paper was tested using a synthetic case study presented in Rossetto et al. ( 2019) [57].There, the FMP-CGM methodology was applied for simulating daily conjunctive use of ground-and surface-water for the irrigation of sunflower.Sunflower yield at harvest was estimated based on crop phenology, climate conditions, and water availability.Here, we coupled the nitrogen module to that case study in order to exploit settings and results of the above-mentioned model.We focused our attention on the sunflower irrigated Water Demand Unit.Model settings are described in the following and reported in Figure 4.The whole model extends over an area of 115 km 2 and it is discretized by means of square cells with 500 m of side length.For the calculation of sunflower water root uptake, root depth values and capillary fringe for the three soil types (silt, silty clay, and sandy loam) are assigned.Groundwater is pumped from a sandy aquifer represented through a homogeneous, MODFLOW convertible model layer with variable thickness.Groundwater is delivered to the irrigation district by six pumping wells.Groundwater flow occurs from west to east.The model was run in transient conditions between 1 January and 31 October 2017, with daily stress periods.Irrigation of sunflower was simulated following sunflower seeding on 1 April and up to 31 August 2017 (when the crop was harvested).During the crop growth the cumulated rainfall depth was 99.8 mm, mainly distributed over four rainfall events occurring between April and August.Model results showed that sunflower water demand is primarily met by precipitation and root uptake, and secondly by irrigation water.Further information on the case study may be found in Rossetto et al., (2019) [57].
The nitrogen cycle module was applied to the present case study to simulate daily nitrate leaching to the water table.
In order to solve Equation ( 1) for the NH 2 − , NH 4 + and NO 3 − pools, initial concentration for each nitrogen compound was defined according to the total content of nitrogen in the soil and to the soil type.The following equations were used for such estimation (Masoni 2010) [60]: For the application of Equations ( 2)-( 4) we assumed: N T = 10 −3 kg/kg (Masoni 2010) [60], ρ a = 1360 kg/m 3 for silt, ρ a = 1440 kg/m 3 for sandy loam, and ρ a = 1200 kg/m 3 for silty clay (Yu et al., 1993) [61].
This resulted in an initial concentration in the irrigation district ranging between 1140 mg/L and 1368 mg/L for NH 2 − , between 48 mg/L and 58 mg/L for NH 4 + , and between 12 mg/L and 14 mg/L for NO 3 − .Nitrogen concentration in rain water varies according to the intensity of precipitation. Figure 5 reports the content of NH 4 + and NO 3 − in rain water over the simulation (Eriksson 1952) [62], compared to the cumulated rainfall.Peak values of the nitrogen content in rain water were input when rainfall events occurred.Decomposition processes for each nitrogen pool were simulated assuming that the first order decay rate constant (λ 1 ) was equal to 0.01 day −1 (Groot et al., 2012) [63].As a consequence, here we assumed that mineralization, nitrification, and denitrification contribute to a 1% daily decrease of NH 2 − , NH 4 + and NO 3 − content, respectively.Moreover, NH 4 + volatilization and sorption were simulated.We assumed that the volatilization fraction (f v ) was equal to 0.001.For simulating sorption, we used the same val-ues as the apparent density for ρ d , and the following values for the Cation Exchange Capacity (CEC): 148 meq/kg for silt, 110 meq/kg for sandy loam, 146 meq/kg for silty clay [64].
For simulating nitrate uptake by sunflower roots over the cropping season, we exploited results from the FMP-CGM model presented in [57].As such, we adopted the approach based on the EPIC conceptualization, by setting the following for sunflower [27]: bn 1 = 0.05, bn 2 = 0.023, bn 3 = 0.0146, PHU = 1600.
In order to define the N fertilizer rate to be applied for achieving the simulated sunflower yield, we applied a conventional N balance approach, as where N f is N fertilizer rate, N d is crop N demand, N r is N supplied by rainfall, N i is N supplied by irrigation, N m is mineralized N from soil organic matter, N l is N loss from leaching, N v is N loss from volatilization, N de is N loss from denitrification, and N p is N released (positive) or taken from soil (negative) for residue decomposition of the preceding wheat crop.All terms in the equation are expressed as kg N ha −1 .Sunflower N demand (N d ) for 5 Mg ha −1 grain yield with 3.3% grain N concentration, harvest index 0.4 and 1.0% residue N concentration is 240 kg N ha −1 .Nitrogen supply from rainfall and irrigation was 2 and 23.5 kg ha −1 , respectively, and from soil organic matter mineralization was 36.8 kg ha −1 from the silt soil and 43.6 kg ha −1 from the sandy-loam soil.Nitrogen loss from leaching, denitrification, and volatilization were assumed to be negligible.Indeed, the amount of irrigation was calculated to avoid percolation and since soils are well drained, they do not offer favorable anaerobic conditions for denitrification.Finally, we assume negligible NH 3 emissions occur with soil incorporation of urea (Rochette et al., 2013) [65].
Consequently, N f was estimated 185.6 kg N ha −1 for the silt soil and 178.8 kg N ha −1 for the sandy-loam soil.We then applied the following values for the dose of fertilization: D = 1.856 × 10 −2 kg/m 2 for silt and D = 1.788 × 10 −2 kg/m 2 for sandy loam.Such values were assigned the day before the beginning of the sunflower cropping season (31 March), and the whole amounts were assigned to the NH 2 − pool only ( f rac NH2 = 1; f rac NH4 = 0; f rac NO3 = 0).Table 2 presents the dataset used.

Results and Discussion
The simulation performed allowed us to estimate daily concentration values of NO 3 − arriving at the water table.Figure 6 presents the daily values of NO 3 − concentration simulated at the water table over the simulation period, with respect to the simulated effective rainfall infiltration through the unsaturated zone.We can notice that leaching occurs before the crop seeding (between early February and early April 2017) and then around mid-September 2017.These leaching events are mostly related to rainfall events concentrated in those periods (Figures 5 and 6).On February-April 2017, daily NO 3 − concentration reaching the water table may reach values up to about 1 mg/L, while values up to about 1.3 mg/L were simulated on mid-September 2017.From the results of the simulation, we explain the relatively high concentration arriving at the water table following rainfall events before the crop season due to the absence of the crop uptake term-similarly, such concentration increases at the end of the crop season, after crop harvesting (that is, no more crop nitrogen uptake).On the other hand, during the sunflower cropping season (April to end of August), the daily NO 3 − concentration values at the water table have very low values (in the order of 10 −3 mg/L), as a result of the drought period simulated and of crop uptake (Figure 6).Our results show a potential risk of NO 3 − leaching during the winter wet season.Moreover, substantial NO 3 − leaching losses may be due to high-intensity rainfalls, such as the one occurred in September 2017.Unfortunately, data coming from experimental activities about daily concentration of NO 3 − arriving at the water table are not published yet-and because of this it was not possible to validate the simulated values.− concentration arriving at the water table and the effective rainfall infiltration over the simulation period.
The time-varying simulated data in Figure 6 refers to one single point in space.Figure 7 shows the spatial distribution of NO 3 − concentrations arriving at the water table.At the selected time nitrogen crop demand is high.The available nitrogen is uptaken through the roots, and, as result of the simulation, we have very low values of NO 3 − concentrations arriving at the water table.These concentration values may then be used to define the time and spatially variant source term for simulation of nitrate transport in aquifers (e.g., using the MT3DMS code [59]).The coupling approach presented holds the following limitations: the unsaturated zone is discretized in only one layer.It is then considered as a unique layer extending from the ground surface to the water table.This representation may not be suitable in case of non-homogeneous soil sediments, with different physical and chemical properties.Furthermore, NO 3 − uptake at different depths according to the root distribution in the vertical domain may not be adequately reproduced; -purely advective transport of nitrogen through the unsaturated zone is simulated (dispersion is neglected); -the biological fixation process is not simulated.This can lead to underestimation of mineral nitrogen available in the soil; -volatilization of NH 4 + is taken into account through a volatilization coefficient (f v ) only, which is difficult to quantify; -NH 4 + uptake by plants' roots is not simulated (we assume, indeed, that NO 3 − is the predominant form of available nitrogen uptaken by roots); -taking steps from the UZF MODFLOW package approach, we assume that vertical movement of nitrogen is the predominant process.As such, lateral movement of nitrogen across neighboring cells is not simulated; -when the thickness of the unsaturated zone varies, according to water level fluctuations, the mass of nitrogen stored in the soil is simply redistributed over the whole unsaturated thickness.
Due to the number of required parameters, reduced compared to other simulation methods, but still very large, we were only able to validate this approach by means of a synthetic case study.Validation by means of real data is however needed, i.e., to verify in a real environment the processes that lead to the low concentration values in Figure 7. Indeed, coupling the nitrogen to the hydrological cycle means, when coming to model validation, to get data on the hydrology, the crop growth and nitrate cycle at the same time, which is more than 50 parameter's values required (see Appendix B for details) with related uncertainties evaluation.In this sense, experimental trials ranging from lysimeters to field scale pilot (as in [66][67][68]), even with high-resolution lysimeters [69] may allow to get the required data.Dedicated field scale experiments, where longer flowpaths through the unsaturated zone may be taken into account, using also state-of-the art ground and remote sensors will allow to complete the validation of the presented approach.

Conclusions
Agriculture-related activities very often cause leaching of high levels of nutrients in the soil and the connected aquatic ecosystems, including groundwater.Among nutrients, nitrate is of particular concern, as it is conservative and it may rapidly leach to the water table, causing high, or even exceeding, according to ongoing regulations, concentrations in aquifers and, consequently, potential negative effects on the environment and even human health, in case groundwater is used for drinking purposes.
In this paper, we described the development of a nitrogen module for the simulation of the nitrogen cycle in the unsaturated zone and the nitrate leaching process.The developed module takes steps from the mathematical framework of the ANIMO model, which solves a mass conservation and transport equation for dissolved organic nitrogen, nitrate, and ammonium to quantify the relation between fertilization, soil management and nutrients leaching occurring over time in a vertical soil column.The simulation of nitrate crop demand and root uptake, instead, is based on the more robust conceptualization of the EPIC model, which accounts for the different phases of the cropping season to calculate crop nitrate content and need.
We integrated the developed nitrogen cycle module within the GIS-based FREEWAT platform, in order to couple it with the spatially distributed codes (e.g., MODFLOW) for the simulation of unsaturated zone and groundwater flow processes.Thanks to the integration, we exploited inputs and outputs from the hydrological models to test the developed modules for simulating the processes involved in the nitrogen cycle.To do so, the nitrogen cycle module was developed to run at each cell and for the same stress periods of a MODFLOW model previously set.More specifically, the developed module exploits results from the UZF MODFLOW package for getting the water content and thickness of the unsaturated zone, and infiltration and runoff rates generated during rainfall events.Moreover, results from an agricultural water management scenario based on the application of FMP-CGM are used to determine the nitrate demand by crops over the growing season.
The coupling approach just summarized guarantees the possibility to implement nitrogen cycle scenarios requiring relatively few inputs.The simulate daily nitrate concentration at the water table distributed in space is the main relevant output derived from this simulation approach.
The presented effort is an initial stage of research in producing usable numerical tools for improving nitrogen management in agricultural fields.Further development of the code may include, among the others, the integration of lateral movement of nitrogen between neighboring cells, and the simulation of root uptake with a finer vertical discretization based on root depth and distribution.
Leaching rate [L 3 /T] Estimated by UZF package

Appendix C Conceptualization of the Nitrogen Cycle Module
In this appendix, we detail the terms involved in Equation (A1).To this scope, it is worth introducing some basic concepts which will be useful in the following.According to the coupling approach described in the paper, Equation (A1) is solved cell by cell.As such, the spatial dependence of all the involved variables will be omitted in the following for the sake of simplicity, and we will indicate the x and y dimensions of a cell with ∆x and ∆y, respectively.
In the following, c(t 0 ) and c(t) indicate the nitrogen concentration at the beginning and at the end of a stress period, respectively.We will also indicate with ∆t the length of a stress period.
Equation (A1) is solved for each nitrogen pool.As such, c(t 0 ) and c(t) must be read as NH 2 − , NH 4 + , or NO 3 − concentrations, unless otherwise explicated.Please notice that running the nitrogen cycle module requires defining initial concentration c 0 at each cell as initial condition.
Before running the nitrogen cycle module, setting up and running a MODFLOW-2005 model with the UZF MODFLOW package activated is needed, as nitrogen movement through the unsaturated zone is driven by vertical unsaturated flow.The following is calculated by the UZF MODFLOW package: the thickness of the unsaturated zone (∆z(t) in the following); -the water content (θ(t) in the following); -the runoff rate (q runoff (t) in the following); -the flow rate towards the saturated zone (q perc (t) in the following).
When running the MODFLOW-2005 model, the User has the possibility to store the above variables only for selected cells; that is, where the nitrogen cycle simulation will be run, in order to avoid long computational times.
We distinguish between input parameters, whose values must User-defined, and input parameters whose values are automatically retrieved from the groundwater model setup or derived as simulation results of the UZF package, CGM or FMP.
For the meaning of all symbols, please refer to Appendix B. Concentration change in the aqueous phase.This is expressed as: where ϕ is the water content change with time.Here, we assume a linear variation of the water content with time: where θ(t 0 ) and θ(t) are calculated by the UZF MODFLOW package.Equilibrium sorption of NH 4 + by soil grains.This is expressed as: where ρ d (User-defined) is the dry bulk density of a specific soil type and K d is the sorption coefficient.The latter is defined as the sorbed ammonium content change depending on the nitrogen concentration in the soil.Assuming a linear sorption isotherm, in a first estimate, K d may be related to the Cation Exchange Capacity (CEC; Groenendijk and Kroes 1999): Provided that CEC values are assigned by the User for each soil type, this equation is solved internally in the nitrogen cycle module.Vertical nitrogen flow through the unsaturated zone.This is expressed as: where q in (t) is the rainfall flux, retrieved among the inputs to the UZF MODFLOW package, while c rain (t) is the nitrogen concentration in the rain water (User-defined).NO 3 − uptake by plants' roots.This is expressed as: Two alternative approaches may be adopted to represent this process.
In the first approach, the EPIC model conceptualization is adopted to calculate NO 3 − crop demand.Such approach requires prior setting and running of an FMP-CGM scenario with daily stress periods.The following parameters are required, among the others: the air average temperature, the crop base temperature, the cumulated biomass, and the crop transpiration flux.The air average temperature and the crop base temperature are input variables to the CGM, while the cumulated biomass and the crop transpiration flux are calculated by the CGM and the FMP, respectively.As such, in this first approach, a coupling effort between FMP-CGM and the EPIC sub-routine for NO 3 − crop demand was made.Since the following equations are those integrated in the EPIC model, the involved variables are expressed in kg, ha and mm, as mass, area and length units of measurements, respectively.Conversions to model units were performed internally in the source code presented in this paper.
According to EPIC, the NO 3 − crop demand can be calculated as , where b n1 , b n2 and b n3 are crop parameters expressing NO 3 − concentration in the crop at emergence, at 0.5 maturity, at maturity, respectively.Values for such parameters are crop-dependent and may be found in the literature (e.g., Sharpley 1990).HUI(t) (the Heat Unit Index) involved in the calculation of c NB (t), in turn, is calculated as where HU sp (the daily Heat Unit accumulation in the sp-th stress period) is expressed in • C.This is calculated from air average temperature and crop base temperature, which are retrieved among the CGM inputs.PHU (the Potential Heat Units required for crop maturity) is expressed in • C too.This parameter is crop-dependent and values may be found in the literature (e.g., Sharpley 1990).B(t) is the cumulated biomass, expressed in kg/ha, as simulated by the CGM.The UN sp term is calculated according to the following equation, where the subscript sp has been replaced by t: where trp(t) (the crop transpiration flux) is calculated by FMP and it is converted, internally in the source code presented in this paper, into mm/day to fit units of measurements, as required in the EPIC equations.mass NO3 (t) (the NO 3 − mass in the unsaturated portion of the cell) is expressed in kg/ha.θ(t) (the water content) and ∆z(t) (the thickness of the unsaturated zone), are calculated by the UZF MODFLOW package.∆z(t) values are converted, internally in the source code presented in this paper, into mm to fit units of measurements, as required in the EPIC equations.
Once UND(t) is calculated, the NO 3 − mass that could be potentially uptaken by the crop (POT(t)) is calculated as following: POT(t) is then compared to the residual NO 3 − mass (RES(t)), available in the unsaturated portion of the cell as a result of all the other processes involved in the nitrogen cycle.In this way, the actual NO 3 − uptake (R u,NO3 (t)) is inferred: R u,NO3 (t) = RES(t) i f POT(t) ≥ RES(t) POT(t) i f POT(t) < RES(t) (A2) In the first case, an NO 3 − deficit is calculated, as the difference POT(t) − R u,NO3 (t).For stress periods other than the first one, such deficit is summed to the NO 3 − crop demand (UND(t)) in the following stress period.In the second case, no NO 3 − deficit occurs.After solving all the above equations in the units of measurements required by EPIC (i.e., kg for mass, ha for area and mm for length), conversions are made internally in the nitrogen cycle module, to express R u,NO3 (t) in model units.
In the second approach, the User has to define the optimal NO 3 − concentration for each crop type and for each stress period, not necessarily on a daily basis (i.e., running an FMP-CGM scenario is not needed).On the other hand, the concentration of NO 3 − available in the soil for crop uptake is calculated as: AWc(t), where AW is the plant available water.It is calculated as: where FC is the field capacity and WP is the wilting point.These are soil-dependent parameters and they have to be defined by the User for each soil type.In this second approach, two alternative cases may occur: if the optimal, User-defined, NO 3 − concentration is higher or equal to the concentration of NO 3 − available in the soil for crop uptake, then POT(t) equals the concentration of NO 3 − available in the soil for crop uptake; -if the optimal, User-defined, NO 3 − concentration is lower than the concentration of NO 3 − available in the soil for crop uptake, then POT(t) equals the optimal NO 3 − concentration.
POT(t) is then compared to RES(t), and R u,NO3 (t) is inferred from Equation (A2).In the first case of Equation (A2), an NO 3 − deficit is calculated, as the difference between the optimal NO 3 − concentration and R u,NO3 (t).For stress periods other than the first one, such deficit is summed to the optimal NO 3 − concentration in the following stress period.In the second case of Equation (A2), no NO 3 − deficit occurs.
Decomposition term.Decomposition processes are related to the mineralization of NH 2 − , nitrification of NH 4 + , denitrification of NO 3 − .This term is expressed as where λ 1 is the first order decay rate constant (User-defined).As highlighted above, each of these terms is subtracted from the related pool and summed to the R p term of the pool below; i.e., mineralization represents a decomposition term (sink, to be subtracted) for the NH 2 − pool and a production term (source, to be summed) for the NH 4 + pool, while nitrification represents a decomposition term (sink, to be subtracted) for the NH 4 + pool and a production term (source, to be summed) for the NO 3 − pool.Production term.This term is expressed as R p .It includes, for each nitrogen pool: the decomposition term from the above pool (source term to be summed to the current pool), as detailed in the above lines; -surface runoff (sink term to be subtracted), calculated as Other nitrogen point sources/sinks (sources or sinks of nitrogen to be summed or subtracted, respectively) may be defined at each cell and for each stress period, by assigning a specific concentration of NH 2 − , NH 4 + , or NO 3 − (a positive concentration must be assigned to simulate a source of nitrogen mass, a negative concentration must be assigned to simulate a sink of nitrogen mass).
The solution of Equation (A1) is expressed as where

Figure 1 .
Figure 1.Sketch of the main processes involved in the nitrogen cycle.

Figure 2 .
Figure 2. Sketch of all the processes related to nitrogen mass flux in the unsaturated zone, occurring within a grid cell.

•
running a MODFLOW-2005 model, including the UZF Package, to simulate water flow through the unsaturated zone.This is needed for the calculation of the thickness of the unsaturated zone, the water content, the runoff rate and the leaching rate in space and time; • running an FMP-CGM scenario.This is needed for the calculation of crop transpiration flux and cumulated above-ground biomass; -if crop NO 3 − demand is not simulated, but the User inputs crop NO 3 − requirement parameters over the cropping season (Option 2).In such case, only running a MODFLOW-2005 model, including the UZF Package, is needed to simulate water flow through the unsaturated zone as above-mentioned.

Figure 3 .
Figure 3. Overview of the coupling approach designed to integrate the nitrogen cycle module within the FREEWAT platform.

Figure 4 .
Figure 4. Model settings for the synthetic case study application.The dotted lines refer to the simulated groundwater level on 1 April 2017 (modified after Rossetto et al., 2019) [57].

Figure 5 .
Figure 5. Variation of nitrogen content in rain water and cumulated rainfall over the simulation period.

Figure 6 .
Figure 6.Daily simulated NO 3− concentration arriving at the water table and the effective rainfall

Figure 7 .
Figure 7. Example of simulated spatially distributed values of nitrate concentrations arriving at the water table (focus over the irrigation district).

3 −
content in the crop Derived variable D Dose of fertilizer application of the applied material [M/L 2 ] User-defined

flow term for NH 2 − 3 −
mass available in the grid cell [M/(L 2 *T)] Derived variable and converted into kg/(ha*day) trp Crop transpiration flux [L/T] Estimated by FMP and converted into mm/day

Table 2 .
Values set for the input parameters in the case study.
c(t) ∆z(t) ∆x ∆y -ammonium volatilization (sink term to be subtracted from the NH 4 + pool), calculated asR v,NH4 = f v c NH4 (t) ∆t ,where f v is the volatilization fraction for NH 4 + (User-defined).Please, notice that for the NH 4 + pool, this source term is actually calculated as is the dose of fertilizer application and frac N (t) is the fraction of NH 2 − , NH 4 + , or NO 3 − of the applied material, both defined by the User.