^{1}

^{*}

^{2}

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-situleaching 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 giv ing 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.

Leaching of NO_{3}-N from the rootzone to the subsoil, groundwater, or surface water is a problem in high-input agriculture [_{3}-N leaching. A major problem encountered is the large spatial variability of NO_{3}-N concentrations in soil, soil solution, or groundwater samples [_{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 [

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

Elverdinge | Assenede | |
---|---|---|

Drain depth (cm) | 70 | 70 |

Drain spacing (cm) | 700 | 1200 |

Effective diameter (cm) | 2.5 | 2.5 |

Depth of soil profile (cm) | 400 | 400 |

q_{r} (cm^{3}∙cm^{−3}) | q_{s } (cm^{3}∙cm^{−3}) | a (cm^{−1}) | n (−) | K_{s } (cm∙d^{−1}) | l (−) | |
---|---|---|---|---|---|---|

Elverdinge | ||||||

0 - 30 cm | 0.0001 | 0.3811 | 0.0067 | 1.2726 | 5.76 | 0.5 |

30 - 60 cm | 0.0001 | 0.3891 | 0.0045 | 1.2381 | 22.62 | 0.5 |

60 - 400 cm | 0.0001 | 0.3756 | 0.0106 | 1.2457 | 36.42 | 0.5 |

Assenede | ||||||

0 - 30 cm | 0.0001 | 0.3746 | 0.0043 | 1.2837 | 425 | 0.5 |

30 - 60 cm | 0.0524 | 0.3554 | 0.0058 | 2.2185 | 393 | 0.5 |

60 - 400 cm | 0.0001 | 0.3565 | 0.0070 | 1.5984 | 25 | 0.5 |

on both water retention and multi-step outflow data, using the multi-step outflow program [

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 [

van Genuchten parameters were fitted to the measurement using a least- squares approach (RETC) [_{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 [_{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.

・

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

・ NO 3 − and NH 4 + on the soil samples are determined by a KCl (potassium chloride) extraction followed by spectrophotometric analysis [

HYDRUS-2D [_{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:

∂ θ ∂ t = ∇ ⋅ ( K ∇ h + K ∇ z ) − 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:

∂ ( θ c ) ∂ t = − ∇ . ( q c − θ D ∇ c ) − 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 [

θ D = θ D i j = D T | q | δ i j + ( D L − D T ) q i q j | q | + θ D d τ δ i j (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}], | q | 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).

DRAINMOD [_{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 [

∂ ( θ C ) ∂ t = ∂ ∂ z ( θ D ∂ C ∂ z ) − ∂ ( q C ) ∂ z + Γ (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 = λ | q θ | + τ 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 [

Γ = Γ + dep Γ fer + Γ mnl − Γ rnf − Γ upt − Γ 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}].

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 (

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

Nitrogen movement parameters in DRAINMOD-N

Nitrogen-related parameters required in DRAINMOD-N include standard

Elverdinge | Assenede | |
---|---|---|

Groundwater level | −60 | −77 |

NO_{3} (kg N/ha) | ||

0 - 30 cm | 70.0 | 54.8 |

30 - 60 cm | 77.8 | 40.7 |

60 - 90 cm | 20.1 | 10.7 |

90 - 120 cm | 60.2 | 10.7 |

120 - 150 cm | 100.5 | 10.7 |

150 - 400 cm | 65.3 | 10.7 |

Elverdinge | Assenede | |
---|---|---|

longitudinal dispersivity D_{L} (cm) | 10 | 10 |

transversal dispersivity D_{T} (cm) | 10 | 10 |

diffusion coefficient D_{d} (cm^{2}∙d^{−1}) | 4.32 | 4.32 |

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.

The qualitative judgement of when the model performance is good is a subjective matter [

Mean Absolute Error (MAE)

M A E = ∑ i = 1 n | ( O i − P i ) | 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)

R R M S E = 1 n ∑ i = 1 n ( P i − O i ) 2 O _ (8)

where 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)

E F = ∑ i = 1 n ( O i − O ¯ ) 2 − ∑ i = 1 n ( P i − O i ) 2 ∑ i = 1 n ( O i − O ¯ ) 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)

C R M = ∑ i = 1 n O i − ∑ i = 1 n P i ∑ 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)

C D = ∑ i = 1 n ( O i − O ¯ ) 2 ∑ i = 1 n ( P i − O ¯ ) 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 = [ ∑ i = 1 n ( O i − O ¯ ) ( P i − P ¯ ) ∑ i = 1 n ( O i − O ¯ ) 2 ∑ i = 1 n ( P i − P ¯ ) 2 ] 2 (12)

where 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 [

To ensure a good modeling of the nitrate leaching a good water table prediction is a necessity [

needed to carry nitrate within the soil profile, the results of the water quantity modeling as presented in

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

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

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

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

NO_{3}-N_{0}_{ }_{-}_{ }_{30 cm} | MAE | RRMSE | EF | CRM | R^{2} |
---|---|---|---|---|---|

Measured-HYDRUS_2D | 11.84 | 0.566 | 0.379 | +0.287 | 0.654 |

Measured-DRAINMOD | 13.10 | 0.616 | 0.266 | +0.351 | 0.698 |

NO_{3}-N_{0}_{ }_{-}_{ }_{90 cm} | |||||

Measured-HYDRUS_2D | 24.17 | 0.336 | 0.678 | −0.067 | 0.693 |

Measured-DRAINMOD | 26.59 | 0.348 | 0.654 | +0.069 | 0.712 |

NO_{3}-N in drainage discharge | |||||

Measured-HYDRUS_2D | 5.09 | 0.209 | 0.479 | −0.180 | 0.99 |

Measured-DRAINMOD | 0.29 | 0.014 | 0.998 | −0.002 | 1.00 |

NO_{3}-N_{0}_{ }_{-}_{ }_{30 cm} | MAE | RRMSE | EF | CRM | R^{2} |
---|---|---|---|---|---|

Measured-HYDRUS_2D | 4.47 | 0.289 | 0.850 | +0.200 | 0.939 |

Measured-DRAINMOD | 5.21 | 0.315 | 0.823 | −0.008 | 0.867 |

NO_{3}-N_{0}_{ }_{-}_{ }_{90 cm} | |||||

Measured-HYDRUS_2D | 22.19 | 0.404 | −0.089 | +0.289 | 0.499 |

Measured-DRAINMOD | 15.80 | 0.272 | 0.505 | +0.123 | 0.639 |

NO_{3}-N in drainage discharge | |||||

Measured-HYDRUS_2D | 1.313 | 0.139 | 6.143 | +0.134 | 0.877 |

Measured-DRAINMOD | 0.293 | 0.038 | 0.471 | −0.009 | 0.996 |

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.

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