Modeling of Water Flow and Nitrate Transport to Subsurface Drains

Show more

1. Introduction

Leaching of NO_{3}-N from the rootzone to the subsoil, groundwater, or surface water is a problem in high-input agriculture [1] [2] . Nitrate concentrations in these environmental compartments often exceed criteria set by the European Union [3] . In many field experiments, it has been difficult to quantify NO_{3}-N leaching. A major problem encountered is the large spatial variability of NO_{3}-N concentrations in soil, soil solution, or groundwater samples [4] . This variability makes it difficult to draw conclusions for an entire field [5] . More integrated information on the nitrate-nitrogen leaching can be derived from measuring, for example, the concentration of the drain water. However, most water and nitrate- nitrogen discharge from a drain originates from the catchment area of the drain [6] , which as a function of the dryness of the soil profile is often limited to the 2 m thick soil layer around the drain, and therefore not always representative for the entire field. Computer simulation models are useful tools to evaluate the complex mechanisms of nitrogen transport and transformation in agricultural fields [7] . In the fall-winter-spring period in Belgium, significant NO_{3}-N losses can occur due to leaching of nitrate that remains in the soil after harvest. Mineralization of organic nitrogen in soil, organic material, plant residue or manure, in combination with the rainfall excess increases in this period the leaching of NO_{3}-N. In order to meet the EU-norm of 11.3 mg/l of NO_{3}-N in surface and groundwater, the Flemish Government in Belgium states that the residual mineral nitrogen, measured in the soil profile of cropland (0 - 90 cm), may not exceed 90 kg N ha^{−1} between the 1^{st} of October and the 15^{th} of November [8] [9] . The objective of this article is to evaluate the performance of a simplified and detailed model approach to predict the nitrate transport to a subsurface drain, with application to the specific hydro-geological conditions of Elverdinge and Assenede, Belgium in the fall-winter period 2000-2001.

2. Field data

Elverdinge and Assenede are experimental fields situated in Flanders, Belgium. Soil physical properties were determined for each soil horizon, using undisturbed core samples. The soil types in the fields are classified as a sandy loam and clay soils respectively. The groundwater level fluctuates between 40 and 80 cm below surface in the study period. The fields were left fallow during the fall- winter period. Soil samples were taken with three weeks intervals in layers of 30 cm, up to a depth of 90 cm. Nitrate concentration of the soil water was determined using suction cups at a depth of 90 cm. Samples of drainage water and groundwater were also analyzed for nitrate. Also, the mineralization rate and denitrification capacity of the soils were measured to get a better estimation of the nitrate leaching. The soil moisture content was measured at several depths in the soil profile, taking with an auger soil samples in the layers 0 - 30 cm, 30 - 60 cm and 60 - 90 cm. From the wet and dry weight of the soil samples the water content was derived. The fields are equipped with a subsurface drainage system consisting of parallel, 10 cm diameter, corrugated plastic drains, placed at a depth of 70 cm below surface, and spacings of 7 m in Elverdinge and 12 m in Assenede. The soil profile was assumed to have a depth of 4 m, on top of an impermeable layer. The properties of the two fields are shown in Table 1 and Table 2. Soil physical properties were determined for each distinguishable soil horizon, using undisturbed soil samples taken with Kopecky rings. van Genuchten-Mualem parameters for describing the hydraulic functions [10] were fitted

Table 1. Drain properties of the fields in Elverdinge and Assenede.

Table 2. Soil hydraulic properties for the studied horizons.

on both water retention and multi-step outflow data, using the multi-step outflow program [11] . Basic water retention and hydraulic conductivity curves were established by averaging individual curves for each soil layer.

2.1. Hydraulic Properties

Hydraulic properties were determined on small undisturbed soil cores (5 cm diameter, 5.1 cm height). For both fields, three horizons could be distinguished in which three soil samples were taken. Saturated hydraulic conductivity was measured on the undisturbed soil cores using the constant head method [12] and five points of the retention curve were determined with a combination of hanging water columns and pressure cells [13] . The limited amount of samples (nine for each field) does not allow statements about the spatial variability of soil hydraulic properties within the fields, but the data does show a larger variability between the different soil layers than between the different replicates for each layer.

van Genuchten parameters were fitted to the measurement using a least- squares approach (RETC) [14] . In this approach, θ_{r} must essentially be seen as an empirical constant in the van Genuchten water retention functions rather than a true physical parameter representing the maximum amount of water that will not contribute to liquid flow [10] [15] . Therefore, no true physical meaning should be attributed to the low values for the residual moisture content for the Elverdinge field (θ_{r} = 0.0001). The rather low value for the saturated hydraulic conductivity of the top layer of the Elverdinge field is a measured value (the average of three samples) and may be ascribed to some form of compaction.

2.2. Measurement of N-Values

・ Figure 4 does not represent nitrate concentrations in the shallow groundwater but in the drainage water.

・ N-mineralisation and denitrification rate is measured for both fields based on incubation experiments [16] .

・ ${\text{NO}}_{3}^{-}$ and ${\text{NH}}_{4}^{+}$ on the soil samples are determined by a KCl (potassium chloride) extraction followed by spectrophotometric analysis [16] .

3. Theory

3.1. Simplified Model

HYDRUS-2D [17] was used in this study as a simplified model to describe NO_{3}-N transport and leaching without any production, transformation or exchange taking place. The HYDRUS-2D model numerically solves water flow and solute transport equations in a two-dimensional flow domain in the unsaturated-saturated zone using the Galerkin finite element method. Water flow is calculated solving the Richards’ equation, which can be written as:

$\frac{\partial \theta}{\partial t}=\nabla \cdot \left(K\nabla h+K\nabla z\right)-S$ (1)

where q is the volumetric water content [L^{3}∙L^{−}^{3}], t is time [T], Ñ a vector differential operator [L^{−}^{1}], K the hydraulic conductivity [L∙T^{−}^{1}], h is pressure head [L], z is gravitational head [L], and S a sink term accounting for root water uptake [L^{3}∙L^{−}^{3}∙T^{−}^{1}]. K and S can be functions of position, θ or h, and time. In this work, K is assumed to be isotropic. For a conservative solute, transport is described by the advective-dispersive equation as:

$\frac{\partial \left(\theta c\right)}{\partial t}=-\nabla .\left(qc-\theta D\nabla c\right)-{S}_{c}$ (2)

where c is solute concentration [M∙L^{−}^{3}], q is volumetric water flux density vector [L^{3}∙L^{−}^{2}∙T^{−}^{1}], D is the dispersion tensor [L^{2}∙T^{−}^{1}] and S_{c} a sink term accounting for root solute uptake [M∙L^{−}^{3}∙T^{−}^{1}]. The components of D are given by [18] :

$\theta D=\theta {D}_{ij}={D}_{T}\left|q\right|{\delta}_{ij}+\left({D}_{L}-{D}_{T}\right)\frac{{q}_{i}{q}_{j}}{\left|q\right|}+\theta {D}_{d}\tau {\delta}_{ij}$ (3)

with D_{ij} components of the dispersion tensor [L^{2}∙T^{−}^{1}], D_{L} the longitudinal dispersivity [L], D_{T} the transversal dispersivity [L], D_{d} the molecular diffusion coefficient in free water [L^{2}∙T^{−}^{1}], q_{i} the i-component of q [L^{3}∙L^{−}^{2}∙T^{−}^{1}],
$\left|q\right|$ the absolute value of the water flux density [L^{3}∙L^{−}^{2}∙T^{−}^{1}], t the tortuosity factor [−], and d_{ij} the Kronecker delta function. Equations (1) and (2) can be solved if initial and boundary conditions for a given problem are known. In this work, a fall-winter leaching period was studied so no root water uptake or plant NO_{3}-uptake were considered (S = 0; S_{c} = 0).

3.2. Detailed model

DRAINMOD [19] was used as a detailed model to simulate nitrogen transformation within the soil profile. DRAINMOD is a computer model that was developed to simulate the performance of drainage and related water table management systems. DRAINMOD-N [7] is an add-on module to DRAINMOD for simulating the nitrogen dynamics in artificially drained soils. Nitrate-nitrogen (NO_{3}-N) is the main N pool considered. The ammonium-nitrogen pool is ignored because in most soils ammonium nitrifies quickly or stays fixed to the soil; thus ammonium losses in subsurface drainage can be neglected. The controlling processes considered by the model [20] are rainfall deposition, fertilizer dissolution, net mineralization of organic nitrogen, denitrifcation, plant uptake, and surface runoff and subsurface drainage losses. Assuming one-dimensional (vertical) flow processes in the unsaturated zone the nitrogen cycle can be represented by the advective-dispersive-reactive (ADR) equation:

$\frac{\partial \left(\theta C\right)}{\partial t}=\frac{\partial}{\partial z}\left(\theta D\frac{\partial C}{\partial z}\right)-\frac{\partial \left(qC\right)}{\partial z}+\Gamma $ (4)

where: C is the NO_{3}-N concentration [M∙L^{−}^{3}], θ is the volumetric water content [L^{3}∙L^{−}^{3}], q is the vertical water flux [L∙T^{−}^{1}], D is the coefficient of hydrodynamic dispersion [L^{2}∙T^{-1}], Γ is a source/sink term [M∙L^{−}^{3}∙T^{−}^{1}] used to represent additional processes (plant uptake, transformations, etc.), z is the coordinate direction along the flow path [L], and t is the time [T]. The coefficient of hydrodynamic dispersion is defined as follows:

$D=\lambda \left|\frac{q}{\theta}\right|+\tau D*$ (5)

where λ is dispersivity [L], τ is a dimensionless tortuosity factor, and D* is the molecular diffusion coefficient [L^{2}∙M^{−}^{1}]. Assuming z is positive in the downward direction and water flows downward in the soil profile, in DRAINMOD-N, functional relationships are used to quantify processes other than NO_{3}-N transport (the source/sink term in Equation (4)), as follows [21] :

$\Gamma =\Gamma {}_{\text{dep}}+{\Gamma}_{\text{fer}}+{\Gamma}_{\text{mnl}}-{\Gamma}_{\text{rnf}}-{\Gamma}_{\text{upt}}-{\Gamma}_{\text{den}}$ (6)

where: Γ_{dep} stands for rainfall deposition [M∙L^{−}^{3}∙T^{−}^{1}], Γ_{fer} for fertilizer dissolution [M∙L^{−}^{3}∙T^{−}^{1}], Γ_{mnl} for net mineralization [M∙L^{−}^{3}∙T^{−}^{1}], Γ_{rnf} for loss [M∙L^{−}^{3}∙T^{−}^{1}] in surface runoff, Γ_{upt} for plant uptake [M∙L^{−}^{3}∙T^{−}^{1}], and Γ_{den} for denitrification [M∙L^{−}^{3}∙T^{−}^{1}].

4. Materials and Methods

4.1. Model input

Water flow and nitrate leaching are modeled in the flow domain of the drain spacing, and a depth of the soil profile of 4 m. The bottom of the soil profile was assumed to be impermeable. The drain was located at 70 cm depth, and was described as a half circular hole with real physical dimensions. The inner wall of the drain was described as a seepage face, implying that the drain is always practically empty. The models were applied to simulate the lateral subsurface drainage, groundwater level and nitrate leaching for the fall-winter period 2000 - 2001 in Elverdinge and Assenede. In the simulations no ponding at the soil surface was allowed.

As initial condition the measured nitrate concentrations were used (Table 3). For both fields, the flow domain is defined as a 400 cm deep soil profile with a width of one half drain spacing. Since flow on both sides of the drain is symmetrical, describing flow from the middle of the drain to one half of the drain spacing is sufficient. The lower, left and right boundaries are defined as no flow-boundaries, except for the position of the half drain tube at the left boundary. The drain tube is defined as a seepage face, through which water and solutes can leave the flow domain when saturated conditions occur [6] . The simulation period ranged from October 1, 2000 to March 31, 2001. Daily rainfall and evapotranspiration were used for the upper boundary.

4.2. Model Parameters for HYDRUS-2D

In HYDRUS-2D, the flow domain is divided into a network of triangular elements that are smaller close to the surface. The denser grid close to the surface will allow more rapid changes in water and solute content, which can be expected under natural meteorological conditions. The flow domain is divided in subdomains so separate water and NO_{3}-N mass balances could be calculated for the layers 0 - 30 cm, 30 - 60 cm and 60 - 90 cm below surface to obtain water content and NO_{3}-N content values that could be compared to the measurements. Groundwater levels were calculated at different distances of the drain and then averaged over the width of the flow domain to compare to measured levels. The transport parameters used in HYDRUS-2D are shown in Table 4.

4.3. Model parameters for DRAINMOD-N

Nitrogen movement parameters in DRAINMOD-N

Nitrogen-related parameters required in DRAINMOD-N include standard

Table 3. Initial conditions for groundwater and nitrate concentrations.

Table 4. Transport parameters used in HYDRUS-2D calculations.

rate coefficients for denitrification (K_{den}) and net mineralization (K_{min}), soil dispersivity (λ), the nitrogen content in rain and other inputs:

Dispersion parameters

- Soil dispersivity: 10 cm;

- Tortuosity factor: 1;

- Diffusion coefficient: 0.01.

Process rate constants

- Rate coefficient of net mineralization:

Elverdinge: 0.000480 d^{−}^{1};

Assenede: 0.000434 d^{−}^{1}.

- Rate coefficient of denitrification:

Elverdinge: 0.01 d^{−}^{1};^{ }

Assenede: 0.03 d^{−}^{1}.

- Continuous days of saturation for denitrification: 3 days.

Model input

- Damping depth: 50 cm;

- Nitrate concentration of rain: 0.8 mg∙l^{−}^{1};

There is no yield reduction by soil water stress.

4.4. Statistical Analysis

The qualitative judgement of when the model performance is good is a subjective matter [22] . Therefore statistical criteria are used for the quantitative judgement [23] . Statistical based criteria provide a more objective method for evaluation of the performance of the models [24] [25] . In this study the following statistical criteria were used to evaluate the performance of the models:

Mean Absolute Error (MAE)

$MAE=\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum}}\left|\left({O}_{i}-{P}_{i}\right)\right|}}{n}$ (7)

where O_{i} is the observation at time i, P_{i} is the prediction at time i. The MAE has a minimum value of 0.0.

Relative Root Mean Square Error (RRMSE)

$RRMSE=\frac{\sqrt{\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum}}{\left({P}_{i}-{O}_{i}\right)}^{2}}}}{\stackrel{\_}{O}}$ (8)

where $\stackrel{\xaf}{O}$ is the mean of the observed values over the time period (1 to n). The RRMSE has a minimum value of 0.0, with a better agreement close to 0.0.

Model Efficiency (EF)

$EF=\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum}}{\left({O}_{i}-\stackrel{\xaf}{O}\right)}^{2}}-{\displaystyle \underset{i=1}{\overset{n}{\sum}}{\left({P}_{i}-{O}_{i}\right)}^{2}}}{{\displaystyle \underset{i=1}{\overset{n}{\sum}}{\left({O}_{i}-\stackrel{\xaf}{O}\right)}^{2}}}$ (9)

EF ranges from minus infinity to 1.0, with higher values indicating better agreement. If EF is negative, the model prediction is worse than the mean observation.

Coefficient of Residual Mass (CRM)

$CRM=\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum}}{O}_{i}-{\displaystyle \underset{i=1}{\overset{n}{\sum}}{P}_{i}}}}{{\displaystyle \underset{i=1}{\overset{n}{\sum}}{O}_{i}}}$ (10)

The CRM has a maximum value is 1.0. If CRM is negative the model overestimates and vice versa.

Coefficient of Determination (CD)

$CD=\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum}}{\left({O}_{i}-\stackrel{\xaf}{O}\right)}^{2}}}{{\displaystyle \underset{i=1}{\overset{n}{\sum}}{\left({P}_{i}-\stackrel{\xaf}{O}\right)}^{2}}}$ (11)

The CD describes the ratio of the scatter of the simulated values and the observed values around the average of the observations. A CD value of one indicates to what extent the simulated and observed values match perfectly. It is positive defined without upper limit and with zero as a minimum.

Goodness of Fit (R^{2})

${R}^{2}={\left[\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum}}\left({O}_{i}-\stackrel{\xaf}{O}\right)}\left({P}_{i}-\stackrel{\xaf}{P}\right)}{\sqrt{{\displaystyle \underset{i=1}{\overset{n}{\sum}}{\left({O}_{i}-\stackrel{\xaf}{O}\right)}^{2}}}\sqrt{{\displaystyle \underset{i=1}{\overset{n}{\sum}}{\left({P}_{i}-\stackrel{\xaf}{P}\right)}^{2}}}}\right]}^{2}$ (12)

where
$\stackrel{\xaf}{P}$ is the mean of the predicted values over the time period (1 to n). R^{2} is ranging from 0.0 to 1.0 indicating a better agreement for values close to 1.0 and it is known as the goodness of fit [26] [27] .

5. Results and Discussion

5.1. Water Flow

To ensure a good modeling of the nitrate leaching a good water table prediction is a necessity [3] . Therefore in the first step of the analysis the subsurface drainage discharge and the related groundwater level were modeled for the simulation period October 1, 2000 to March 31, 2001, being in Belgium the leaching period. Comparison of measured and simulated results for water table level in Elverdinge field is shown in Figure 1. Assuming that the water is the vehicle

Figure 1. Simulated and measured groundwater level in Elverdinge.

needed to carry nitrate within the soil profile, the results of the water quantity modeling as presented in Figure 1 ensure a good nitrate transport and leaching as a second step. During the simulation period in Elverdinge field, the total amount of rainfall, subsurface drainage and evapotranspiration were 61.3, 50.6 and 11.3 cm, respectively. As expected, there is a high subsurface drainage and a low ET because the simulation period was the fall-winter period. Under extremely wet conditions, the major part of the drainage water originated from the topsoil. Under dry conditions with a relatively deep phreatic surface drainage water mainly originated from the zone close to the drain depth. In general both models (simplified and detailed) simulated quite accurately the hydrological variables.

5.2. Nitrate transport and Leaching

Groundwater levels simulated with both approaches are well. Simulated and measured drain discharge rates in time correspond well for both model approaches. Nitrate concentrations in the soil profiles as shown in Figure 2 and Figure 3 for Elverdinge and Assenede respectively match well with the measurements. For the Assenede field, the DRAINMOD-N model predicts higher nitrate concentrations in the soil profile. This is probably due to N-mineralisation, which is not calculated by the HYDRUS-2D model.

Nitrate concentrations in the drainage water simulated with HYDRUS-2D (simplified model to calculate the nitrate leaching) are initially higher for both fields. This could be an overestimation because denitrification is not taken into account. The changes in NO_{3}-N concentrations in the drainage water in time are simulated well. Nitrate concentrations in the drainage water simulated with DRAINMOD-N (detailed model to calculate the nitrate leaching) and measured correspond very well. At high infiltration rates the NO_{3}-N just above the drain is transported downwards rapidly. The measured, simplified simulated (HYDRUS-2D) and detailed simulated (DRAINMOD-N) are show in Figure 4 and Figure 5 for Elverdinge and Assenede respectively. The results indicate that

Figure 2. NO_{3}-N content in the soil profile in different layers (0 - 30, 30 - 60 and 60 - 90 from top) in Elverdinge.

to improve the description of N leaching, N mineralization, N deposition and denitrification have to be accounted for to allow a better description of N leaching to surface waters and groundwater.

The mean absolute error (MAE), the relative root mean square error (RRMSE) and the coefficient of determination (CD) have a minimum value of 0.0, while the model efficiency (EF) and the coefficient of residual mass (CRM) have a maximum value of 1.0. The goodness of fit (R^{2}) is ranging from 0.0 to 1.0 indicating a better agreement for values close to 1.0. Those statistical indices were calculated for the simulation periods 1991-1992 and 1994-1995, as to assess the performance of the simplified (HYDRUS-2D) and the more detailed (DRAINMOD) model with respect to measurements. According to the CRM the simplified model overestimates the NO_{3}-N content in the period 1991-1992, and underestimates the content in the period 1994-1995, whereas the detailed model overestimates the content in both periods. The results of the statistical analysis are presented in Table 5 and Table 6. The different performance indices clearly illustrate that on average the detailed (DRAINMOD-N) model performs better than the simplified (HYDRUS-2D) model, and this in the calibration and validation phase, as well as when the model is used in a predictive mode.

Figure 3. NO_{3}-N content in the soil profile in different layers (0 - 30, 30 - 60 and 60 - 90 from) in Assenede.

Figure 4. Simulated and measured nitrate-nitrogen (mg∙l^{−1}) in shallow groundwater during the fall-winter season in Elverdinge.

6. Conclusion

The simulation results indicate that both simplified and detailed approaches described the dynamics of the water balance well. The layering of the soil profile had a pronounced effect on the flow paths to the drain under different atmospheric

Figure 5. Simulated and measured nitrate-nitrogen (mg∙l^{−1}) in shallow groundwater during the fall-winter season in Assenede.

Table 5. Statistical performance analysers calculated between observed and simulated nitrate-nitrogen (kg∙ha^{−1}) in soil layers 0 - 30 cm and 0 - 90 cm and nitrate-nitrogen concentration (mg∙l^{−1}) in drainage water in Elverdinge field.

Table 6. Statistical performance analysers calculated between observed and simulated nitrate-nitrogen (kg∙ha^{−1}) in soil layers 0 - 30 cm and 0 - 90 cm and nitrate-nitrogen concentration (mg∙l^{−1}) in drainage water in Assenede field.

conditions. The nitrogen dynamics play an important role in determining actual nitrate concentrations in the soil profile. Considering NO_{3}-N as a conservative solute (as in HYDRUS-2D) leads to a higher estimate of NO_{3} concentrations in the drainage water. Field data of the N balance show that also during winter periods considerable increases in NO_{3}-N discharge occurred due to net mineralization and deposition especially during relatively dry periods. Other processes like immobilization and denitrification make the description of the N dynamics even more complex. The simplified approach is acceptable to give a general overview to predict the nitrate leaching. The detailed approach is required to give the decision maker the opportunity to take the right decision where nitrogen mineralization, N deposition and denitrification have to be accounted for to allow a better description of nitrate leaching. The results of the statistical analysis clearly illustrate that on average the detailed (DRAINMOD-N) model performs better than the simplified (HYDRUS-2D) model.

References

[1] Strebel, O., Duynisveld, W.H.M. and Bottcher, J. (1989) Nitrate Pollution of Groundwater in Western Europe. Agriculture, Ecosystems & Environment, 26, 189- 214.

https://doi.org/10.1016/0167-8809(89)90013-3

[2] El-Sadek, A. (2009) Comparison between Numerical and Analytical Solution of Solute Transport Models. Journal of African Earth Sciences, 55, 63-68.

https://doi.org/10.1016/j.jafrearsci.2008.12.007

[3] El-Sadek, A. (2010) Monte Carlo Approach to Developing a Water Quality Process- Factor. International Journal of Water Resources and Environmental Management, 1, 97-104.

[4] Addiscott, T.M., Whitmore, A.P. and Powlson, D.S. (1991) Farming, Fertilisers and the Nitrate Problem. Rothamsted Experimental Station, Harpenden, United Kingdom.

[5] Roth, K., Fluhler, H., Jury, W.A. and Parker, J.C. (Eds.) (1990) Field-Scale Water and Solute Flux in Soils. Monte Verita. Proceedings of the Centre Stefano Franscini, Ascona, Birkhauser Verlag, Basel, Switzerland.

https://doi.org/10.1007/978-3-0348-9264-3

[6] de Vos, J.A. (1997) Water Flow and Nutrient Transport in a Layered Silt Loam Soil. Ph.D. Thesis, Wageningen Agricultural University, Wageningen, the Netherlands, 287 p.

[7] Brevé, M.A., Skaggs, R.W., Parsons, J.E. and Gilliam, J.W. (1997) DRAINMOD-N: A Nitrogen Model for Artificially Drained Soils. Trans. American Society of Agricultural and Biological Engineers, 40, 1067-1075.

https://doi.org/10.13031/2013.21359

[8] Geypens, M., Feyen, J., Hofman, G., Merckx, R., Van Cleemput, O. and Van Orshoven, J. (2001) Mineral Nitrogen in the Soil as a Policy Instrument to Reduce N- Leaching from Agricultural Soils. Proceeding of 11th Nitrogen Workshop, 9-12 September 2001, Reims, France, 451-452.

[9] El-Sadek, A., Feyen, J., Skaggs, W. and Berlamont, J. (2002) Economics of Nitrate Loss from Drained Agricultural Land. Environmental Engineering, 128, 376-383.

https://doi.org/10.1061/(ASCE)0733-9372(2002)128:4(376)

[10] van Genuchten, M.Th. and Nielsen, D.R. (1985) On Describing and Predicting the Hydraulic Properties of Unsaturated Soils. Annales Geophysicae, 3, 615-628.

[11] van Dam, J.C., Stricker, J.N.M. and Droogers, P. (1990) From One-Step to Multi- Step. Determination of Soil Hydraulic Functions by Outflow Experiments. Rep. 7, Department of Water Resources, Agricultural University, Wageningen.

[12] Klute, A. and Dirksen, C. (1986) Hydraulic Conductivity and Diffusivity: Labora- tory Methods. In: Klute, A., Ed., Methods of Soil Analysis, Part 1 Physical and Mi- nerological Methods, 2nd Edition, Agronomy Series No. 9, American Society of Agronomy and Soil Science Society of America, Madison, 687-734.

[13] Klute, A. (1986) Water Retention: Laboratory Methods. In: Klute, A., Ed., Methods of Soil Analysis, Part 1 Physical and Minerological Methods, 2nd Edition, Agrono- my Series No. 9, American Society of Agronomy and Soil Science Society of Ameri- ca, Madison, 635-662.

[14] van Genuchten, M.Th., Leij, F.J. and Yates, S.R. (1991) The RETC Code for Quan- tifying the Hydraulic Functions of Unsaturated Soils, Version 1.0. EPA Report 600/ 2-91/065, U.S. Salinity Laboratory, USDA, ARS, Riverside.

[15] Luckner, L., van Genuchten, M.Th. and Nielsen, D.R. (1989) A Consistent Set of Parametric Models for the Two-Phase Flow of Immiscible Fluids in the Subsurface. Water Resources Research, 25, 2187-2193.

https://doi.org/10.1029/WR025i010p02187

[16] Geypens, M., Herelixka, E., Vogels, N., Vanongeval, L., Feyen, J., Oorts, K., Rombauts, S., Sammels, L., Verstraeten, W.W., El-Sadek, A., Merckx, R., Coppens, F., Hofman, G., Van Cleemput, O., D’Haene, K., Moreels, E., De Neve, S., Salomez, J., Boeckx, P., Van Orshoven, J., Librecht, I. and Wellens, J. (2002) Bepaling van de hoeveelheid minerale stikstof in de bodem als beleidsinstrument. Vlaamse Landmaatschappij, Brussels. (In Dutch)

[17] Simúnek, J., Sejna, M. and van Genuchten, M.Th. (1999) The HYDRUS-2D Software Package for Simulating Water Flow and Solute Transport in Two-Dimen- sional Variably Saturated Media (Version 2.0). U.S. Salinity Laboratory, Riverside.

[18] Bear, J. (1972) Dynamics of Fluid in Porous Media. Elsevier, New York.

[19] Skaggs, R.W. (1981) Methods for Design and Evaluation of Drainage Water Ma- nagement Systems for Soils with High Water Tables, DRAINMOD. North Carolina State University, Raleigh.

[20] Brevé, M.A., Skaggs, R.W., Kandil, H., Parsons, J.E. and Gilliam, J.W. (1997) DRAINMOD-N, a Nitrogen Model for Artificially Drained Soils. Transactions of the ASAE, 40, 1067-1075.

[21] Brevé, M.A., Skaggs, R.W., Parsons, J.E. and Gilliam, J.W. (1998) Using the DRAIN MOD-N Model to Study Effects of Drainage System Design and Management on Crop Productivity, Profitability and NO3-N Losses in Drainage Water. Agricultural Water Management, 35, 227-243.

https://doi.org/10.1016/S0378-3774(97)00035-8

[22] Anderson, M.P. and Woessner, W.W. (1992) Applied Groundwater Modeling Si- mulation of Flow and Advective Transport. University Press, Cambridge, 296 p.

[23] Vázquez, R.F., Feyen, L., Feyen, J. and Refsgaard, J.C. (2002) Effect of Grid-Size on Effective Parameters and Model Performance of the MIKE SHE Code Applied to a Medium Sized Catchment. Hydrological Processes, 16, 355-372.

https://doi.org/10.1002/hyp.334

[24] Ducheyne, S. (2000) Derivation of the Parameters of the WAVE Model using a Deterministic and a Stochastic Approach. PhD Thesis No. 434, Faculty of Agriculture and Applied Biological Sciences, KU Leuven, Belgium, 123 p.

[25] El-Sadek, A. (2007) Upscaling Field Scale Hydrology and Water Quality Modeling to Catchment Scale. Water Resources Management, 21, 149-169.

https://doi.org/10.1007/s11269-006-9046-y

[26] Shahin, M., Van Oorschot, H.J.L. and De Lange, S.J. (1993) Statistical Analysis in Water Resources Engineering. A. A. Balkema Publishers, Lisse, 394 p.

[27] Legates, D.R. and McCabe, G.J. (1999) Evaluating the Use of “Goodness-of-Fit” Measures in Hydrological and Hydroclimatic Model Validation. Water Resources Research, 35, 233-241.

https://doi.org/10.1029/1998WR900018