ABSTRACT

In this study, the water flow and nitrate transport to a subsurface drain, using a simplified and detailed model, are simulated for the specific hydro-geologi- cal conditions of Elverdinge and Assenede, Belgium. Previously, the DRAIN MOD-N model proved to be able to simulate nitrate concentrations and drainage well for an in-situ leaching experiment, the “Hooibeekhoeve” in the community of Geel (north-eastern part of Belgium), conducted in 1992-1995. In this study, the calibrated model is used to simulate the nitrate leaching for the winter period 2000-2001 in Elverdinge and Assenede and is compared to a model with a simplified nitrate transport description. The comparative analysis between both model approaches reveals that the simplified model is able to predict sufficiently accurate the observed nitrate leaching. The detailed approach however has the advantage of giving a more accurate estimate of the nitrogen mineralization, N deposition and denitrification, resulting in a more precise modeling of the nitrate leaching to surface waters and groundwater.

Keywords:

non-Point Source, subsurface Drainage, Nitrogen cycle, nitrateleaching, water Quality

1. Introduction

Leaching of NO₃-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₃-N leaching. A major problem encountered is the large spatial variability of NO₃-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₃-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₃-N. In order to meet the EU-norm of 11.3 mg/l of NO₃-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⁻¹ 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

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₃-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³∙L⁻³], t is time [T], Ñ a vector differential operator [L⁻¹], K the hydraulic conductivity [L∙T⁻¹], h is pressure head [L], z is gravitational head [L], and S a sink term accounting for root water uptake [L³∙L⁻³∙T⁻¹]. 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⁻³], q is volumetric water flux density vector [L³∙L⁻²∙T⁻¹], D is the dispersion tensor [L²∙T⁻¹] and S_c a sink term accounting for root solute uptake [M∙L⁻³∙T⁻¹]. 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²∙T⁻¹], D_L the longitudinal dispersivity [L], D_T the transversal dispersivity [L], D_d the molecular diffusion coefficient in free water [L²∙T⁻¹], q_i the i-component of q [L³∙L⁻²∙T⁻¹], $\left|q\right|$ the absolute value of the water flux density [L³∙L⁻²∙T⁻¹], 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₃-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₃-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₃-N concentration [M∙L⁻³], θ is the volumetric water content [L³∙L⁻³], q is the vertical water flux [L∙T⁻¹], D is the coefficient of hydrodynamic dispersion [L²∙T^-1], Γ is a source/sink term [M∙L⁻³∙T⁻¹] 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²∙M⁻¹]. 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₃-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⁻³∙T⁻¹], Γ_fer for fertilizer dissolution [M∙L⁻³∙T⁻¹], Γ_mnl for net mineralization [M∙L⁻³∙T⁻¹], Γ_rnf for loss [M∙L⁻³∙T⁻¹] in surface runoff, Γ_upt for plant uptake [M∙L⁻³∙T⁻¹], and Γ_den for denitrification [M∙L⁻³∙T⁻¹].

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₃-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₃-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

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⁻¹;

Assenede: 0.000434 d⁻¹.

- Rate coefficient of denitrification:

Elverdinge: 0.01 d⁻¹;

Assenede: 0.03 d⁻¹.

- Continuous days of saturation for denitrification: 3 days.

Model input

- Damping depth: 50 cm;

- Nitrate concentration of rain: 0.8 mg∙l⁻¹;

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 \sum _{i=1}^n\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 \sum _{i=1}^n{\left(P_i-O_i\right)}^2}}}{\stackrel{\_}{O}}$ (8)

where $\overline{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 \sum _{i=1}^n{\left(O_i-\overline{O}\right)}^2}-{\displaystyle \sum _{i=1}^n{\left(P_i-O_i\right)}^2}}{{\displaystyle \sum _{i=1}^n{\left(O_i-\overline{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 \sum _{i=1}^n{O}_i-{\displaystyle \sum _{i=1}^n{P}_i}}}{{\displaystyle \sum _{i=1}^n{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 \sum _{i=1}^n{\left(O_i-\overline{O}\right)}^2}}{{\displaystyle \sum _{i=1}^n{\left(P_i-\overline{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²)

$R^2={\left[\frac{{\displaystyle \sum _{i=1}^n\left(O_i-\overline{O}\right)}\left(P_i-\overline{P}\right)}{\sqrt{{\displaystyle \sum _{i=1}^n{\left(O_i-\overline{O}\right)}^2}}\sqrt{{\displaystyle \sum _{i=1}^n{\left(P_i-\overline{P}\right)}^2}}}\right]}^2$ (12)

where $\overline{P}$ is the mean of the predicted values over the time period (1 to n). R² 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

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₃-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₃-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

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²) 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₃-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.

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

conditions. The nitrogen dynamics play an important role in determining actual nitrate concentrations in the soil profile. Considering NO₃-N as a conservative solute (as in HYDRUS-2D) leads to a higher estimate of NO₃ concentrations in the drainage water. Field data of the N balance show that also during winter periods considerable increases in NO₃-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.

Cite this paper

El-Sadek, A. and Radwan, M. (2018) Modeling of Water Flow and Nitrate Transport to Subsurface Drains. Journal of Agricultural Chemistry and Environment, 7, 45-59. https://doi.org/10.4236/jacen.2018.71005

References