Comparison of Two Nonlinear Radiation Models for Agricultural Subsurface Drainage

The subject of 'drainage: draining the water off' is as important as 'irrigation: application of water', if not more. 'Drainage' has a deep impact on food security, agricultural activity, hygiene and sanitation, municipal usage, land reclamation and usage, flood and debris flow control, hydrological disaster management, ecological and environmental balance, and water resource management. 'Drainage Systems' provides the reader with a tri-dimensional expose of drainage in terms of sustainable systems, surface drainage and subsurface drainage. varied technical backgrounds and experiences from around the world extensive range of issues concerning the drainage phenomenon. Field engineers, hydrologists,


Introduction
One of the basic functions of subsurface agricultural drainage systems is to avoid the establishment of a moisture regime adverse to crop development by means of the depletion of shallow groundwater and the timely evacuation of water excess stemming from overirrigation, rainfall, losses due to infiltration in channels and contributions from underground streams. In order to make an efficient water evacuation, it is crucial to know the way in which it moves in the soil so as to determine the water flow that can be removed from the porous medium by means of drainage with the corresponding variation of the aquifer level.
The determination of the variables of an agricultural drainage system requires the analysis of the mass and energy transfers occurring in the soil. The study of these transfers, despite it being highly complex -since it is about analysing processes that are basically nonlinear and occurring in a medium in which properties vary over time and space -may be executed by considering some of the following analysis scales: Microscopic scale: In this analysis scale, and corresponding to each soil pore, the mean velocity or microscopic water flow is estimated with the Poiseuille law, which comes from the Navier-Stokes equations, and the pressure in each pore is estimated with a Laplace equation. This analysis scale is recommended for a fine understanding of the fundamental mechanisms of water transfer processes in the soil.
Macroscopic scale: The complexity outlined by the specific definition of the geometric shape of the pore space means that the microscopic description cannot be implemented without a change of the scale, whose essential stage consists in the introduction of the concept of the representative volume element (RVE) which allows for the establishment of an equivalence between the real porous medium (dispersed) and a fictitious porous medium (continuous). In this analysis scale, which concurs with a set of pores of a size's range, the mean velocity in the pores filled with water or macroscopic flow is estimated with the Darcy -Buckingham law (1907), and the water pressure associated to the set is estimated with their own Laplace law applied to a larger pore size. The corresponding transfer equation is known as Richards equation (1931).
Megascopic Scale: In this analysis scale, which corresponds to a set of soils, the mean velocity -or megascopic flow -is estimated with the Darcy law and averaged through consideration of the Dupuit-Forcheimer hypothesis, relevant to a hydrostatic pressure distribution; moreover, water pressure is provided by the piezometers. The relevant transfer equation is known as the Boussinesq equation of agriculture drainage.
As a practical matter, the analysis of agricultural drainage is made either with the Richards equation (Zaradny & Feddes, 1979;Fipps & Skaags, 1986;Saucedo et al. 2002) or else with the Boussinesq equation of agriculture drainage (Dumm, 1954;Pandey, et al., 1992;Gupta et al., 1994;Samani et al., 2007). The Richards equation allows the executing of descriptions of the transfer processes occurring in the saturated and unsaturated zones of the soil; however, its application to the scale of an irrigation district and even that of the a farm field, is limited due to the difficulty and cost of the experim e n t a l w o r k r e q u i r e d t o d e p i c t t h e s o i l ' s hydrodynamic characteristics (the moisture retention curve and the hydraulic conductivity curve) as well as the necessary effort of calculating a three-dimensional water movement in the soil. These limitations have led to the analysis of agricultural drainage being mainly performed with the Boussinesq equation, an approach that considers in a simplified manner the transfers occurring in the soil's unsaturated zone, but with a smaller amount of data requirements than the Richards equation, adequately describes water dynamics in the saturated layer of the soil.
Recently, two mechanistic models have been developed for agricultural drainage that improve the traditional hypotheses of the models reported in the literature. On the one hand, Zavala et al. (2005) have developed a model for agricultural subsurface drainage based on the two-dimensional Richards equation. This differential equation is subjected in the drains boundary to a nonlinear radiation condition, and in this form the mass and energy transfers in a drainage system are better represented. On the other hand, Zavala et al. (2004) and Fuentes et al. (2009) have studied agricultural drainage with the Boussinesq equation and have deduced, respectively, the boundary condition to be used in agricultural drains by this equation and the relation between the moisture retention curve and the storage coefficient in shallow unconfined aquifers.
The aim of this chapter is to present the two models just described to develop their numerical solutions and compare the mass and energy transfers obtained with the Richards equation and the Boussinesq equation, both of which are subject to nonlinear radiation conditions in the drains.

Macroscopic scale
Soil water movement in a subsurface drainage system is a three-dimensional phenomenon, for which reason its description should be made using the Richards equation in three dimensions; nonetheless, due to the effort of calculation that the resolution of this form of Richards equation entails, it is convenient to accept the hypothesis that the phenomenon is basically two-dimensional (that is to say, it is made according to planes perpendicular to the direction of the drain). If, in addition, it is assumed that the water uptake by plant roots is negligible, the two-dimensional form of the Richards equation may be written as follows: www.intechopen.com where  is the pressure head   L ;     If a drainage system with equidistant parallel pipes installed at the same depth is considered, it is possible to define a domain for equation (1) as with the one shown in Fig. 1. The description of agricultural subsurface drainage by equation (1) requires the definition of the initial status of the pressure head in the porous media as well as its boundary conditions. The initial condition of the water pressure in the soil is specified as a known space function: By the flow symmetry, it is known that the Darcy flow in a perpendicular direction to segments AF, DE and BC is null (Neumann boundary condition), and a similar situation occurs in the boundary segment CD due to the impervious layer.
where E P is the depth of the impervious layer, measured as from the soil surface; and L is the space between drains.
In the soil surface (segment AB) when the rain or evaporation intensity (i) is known, a Neumann type boundary condition can be implemented: Consistent with Zavala et al. (2005), the soil water transfer to the drain shall be described with the following nonlinear radiation condition: where n  is the normal derivative; o q is a particular value of the water flow in the soil 1 LT    ;  and  are dimensionless shape parameters; P is the depth of the drain   L; a n d t h is the pressure inside the drain   L , equal to the atmospheric pressure   t h0  in the segment of the drain's internal perimeter in contact with the air and equal to the water depth had at every point of the drain's internal perimeter in contact with water ( Fig. 2).
The application of relation (7) requires knowledge of the evolution of mean depth of the water in the drain. If the variations in the longitudinal direction are negligible, the depth of the water in the drain may be supposed to be uniform in space but variable in time. In this situation, the evolution in time of the depth of water may be calculated as from an equation that relates the flow velocity with the energy loss in the movement direction. Because of its generality, the fractal resistance law proposed by Fuentes et al. (2004) is used in this work: where V is the water's mean velocity in the drain 1 LT      ; d is a dimensionless parameter that varies between 12 d 1   in terms of the type of flow (turbulent or laminar);  is a dimensionless coefficient; g is the gravitational acceleration R is the hydraulic radius  L ; and J is the friction slope The combination of relation (8) and the continuity equation for steady flow -which indicates that the flow Q is the product of the hydraulic area and the mean velocity ( QV A  ) -allows the obtaining of the relation between the mean depth of the water in the drain and the water flow that leads to However, the application of this relation displays a limitation: there are two unknown variables (flow and water depth) and only one equation. This problem is resolved by outlining a second equation that is obtained when the nonlinear radiation condition is integrated with (7): where  is the perimeter of the drain semi-circumference and  is its length.
In order to model agricultural drainage with the system of equations (1-9), it is crucial to have the analytical representations of the soil hydrodynamic characteristics   and  K  . In field and in laboratory applications, Fuentes et al. (1992) recommend using the van Genuchten model for the moisture retention curve (van Genuchten, 1980), subject to the Burdine restriction (Burdine, 1953): where s  is the saturated volumetric water content; r  is the residual volumetric water content; d 0  is a pressure scale parameter; m and n are the shape parameters.
As for the hydraulic conductivity curve, they suggest using the Brooks & Corey model (1964): where s K is the saturated hydraulic conductivity 1 LT      ; and  is a positive dimensionless shape parameter.

Megascopic scale
Rough descriptions of the mass and energy transfers of subsurface agricultural drainage systems can be obtained with the Boussinesq equation for unconfined aquifers. As per the hypothesis that variations in the direction of the drain are negligible and that the null www.intechopen.com Drainage Systems 170 recharge -the dynamic of the water in the saturated thickness of the soil -can be described with the Boussinesq equation of agricultural drainage: where H and i H are, respectively, elevations of the free surface or hydraulic head and of the impervious layer, measured from the same reference level   L , when the impervious layer is approximately horizontal it may be supposed as the marker level and take i H0  ;  , which, in a shallow unconfined aquifer is a function of the hydraulic head (Hilberts et al., 2005;Fuentes et al., 2009).
Taking into account the van Genuchten model for the moisture retention curve, subject to the Burdine restriction and the hydrostatic pressure distribution hypothesis, this allows the obtaining of the following analytical representation for the storage coefficient (Fuentes et al., 2009): where H s is the soil surface elevation.
To resolve equation (12) on the domain shown in Fig. 3, it is necessary to define the initial conditions and the boundary conditions. The specification of these limit conditions is more convenient if the free surface position is counted as from the impervious layer: where   h x, t is the hydraulic head counted as from the position of the drains; and o D is the depth of the impervious layer measured as from the drain's position (see Fig. 3).
In general terms, the pressure's initial condition shall be specified as the elevation of the free surface throughout the horizontal coordinate x: Zavala et al. (2004) have shown that the relation that will subject the Boussinesq equation in the boundary of the drains is the following fractal radiation condition: where  is a dimensionless conductance coefficient; s is the quotient of the soil fractal dimension f D and the dimensional Euclidean space ( f sD3  ); s K is the hydraulic conductivity of the soil-drain interface. The positive sign in equation (16) is taken for the drain located in coordinate x 0  and the negative sign for x L  . As per Zavala et al. (2007), it is convenient to express equation (16)

Application
The comparison of the mass and energy transfers provided by the systems of equations (1-11) and (12-18) is executed considering the drainage experimental information of Zavala et al. (2004). The experiment was carried out in a drainage module made with acrylic sheets in which two PVC drains were installed (Fig. 4) . The module was filled with an altered sample of sandy soil of the Mexican region of Tezoyuca, Morelos, passed through a 2 mm sieve; the soil was disposed at 0.20 m thick layers, seeking to maintain a constant bulk density. The soil was saturated by applying a constant water head on its surface until the entrapped air was virtually removed. Once the drains were closed, the water head was removed from the soil surface; the surface of the module was then covered with a plastic in order to avoid evaporation. Finally, the drains were opened to measure the drained water v o l u m e ( t e n d a y s ) ; i t i s w o r t h n o t i c i n g that the initial condition corresponded to a hydrostatic pressure distribution and the recharge was null (R = 0) during the drainage phase.
The hydrodynamic characterisation of the soil, executed independently of the transient drainage test, allowed the determination of the parameters of the van Genuchten model and the Brooks & Corey model (Zavala et al., 2004  To compare the transfers described, with two flow models, it is necessary to resolve the systems (1-11) and (12)(13)(14)(15)(16)(17)(18). Both systems of equations are numerically solved following the process employed by Zavala et al. (2004Zavala et al. ( ) & (2005. The spatial discretisation is carried out by using the Galerkin finite-element method; temporal discretisation is performed with a finitedifference implicit method. The resulting system becomes lineal using the Picard iterative method; the algebraic equation system is solved using a preconditioned conjugated gradient method. These methods are well documented, for example in Zienkiewicz et al. (2005).
The solution domain discretisation of the two-dimensional Richards equation is carried out by applying the Argus-One program, with which a finite element mesh of 10,795 nodes was generated and distributed in 21,082 elements (Fig. 5)   By applying the numerical solutions of the systems (1-11) and (12-18), the drainage experiment is simulated to determine and to compare the evolution in time of the water depth evacuated by the drain and the corresponding variation of the free surface to half of the space between the drains. The results obtained for the mass transfer are presented in Fig.  6a and 6b and, for the energy transfer, the results are presented in Fig. 7a and 7b; in each case, as the drainage time increases, the calculated evolutions trend to a limit value because the recharge is null.
A good agreement between the evolution of the drained depth described with the Richards equation and the evolution obtained with the Boussinesq equation ( 2 R 0.9847  ) can be appreciated in Fig. 6a and 6b;  product of the estimation of the radiation conditions' parameters (7) and (16.1), as from the experimental data of the drained depth. However, the evolution of the water table was not considered to optimise parameters, for which reason this variable is a good comparison element. When observing the results shown in Fig. 7a and 7b, an important discrepancy can be appreciated between the evolution described by the model based on the Richards equation as to the evolution obtained with the Boussinesq model; during all the drainage time, the water table drawdown calculated with the Boussinesq model is slower than the one obtained with the Richards model.
Taking into account that the Boussinesq equation may be obtained by integrating the Richards equation in the vertical with the hypothesis of a hydrostatic pressure distribution (Bear, 1972) -circumstances that are not met in the regions closer to the drain -we have the most accurate description of the hydraulic variables in a subsurface drainage system corresponding to the one provided by the Richards equation. Considering, in addition, the simulation results, it can be seen that the model based on the Boussinesq equation (equations 12-18) cannot simultaneously reproduce the evolutions of the mass and energy transfers that are described by the model based on the Richards equation; that is to say, if it reproduces the evolution of the mass, it is not feasible that it reproduces the energy evolution or, if it reproduces the energy evolution, it cannot reproduce the mass evolution.
The limitations of an accurate simultaneous description of mass and energy transfers with the Boussinesq equation should be had in mind when it is used to estimate soil's hydraulic properties or when it is applied to the design of drainage systems. On one hand, it is traditional to consider this equation in estimating the saturated hydraulic conductivity as with the lowering or recovery measurements of the groundwater; if this were the case, the determined value would be higher than the real value of the hydraulic conductivity of the porous medium, because -as per the results of this work -the Boussinesq equation describes a minor lowering of the free surface than the one occurring in the drainage system. On the other hand, if the saturated hydraulic conductivity has been estimated according to field or laboratory tests that consider relations that are more accurate than the Boussinesq equation, and this value is used together with the Boussinesq equation to calculate the space between drains, it is possible to obtain separations shorter than that which is really needed to satisfy the water-table drawdown.
To illustrate both situations, the drainage results obtained with the Richards model are regarded as benchmarks, and numerical simulations with the Boussinesq model are carried out. The first case involves determining the saturated hydraulic conductivity value which allows for approaching with the Boussinesq model, with the water-table drawdown to half of space between drains as described with the Richards equation. In the second case, the saturated hydraulic conductivity value determined in the laboratory by Zavala et al. (2004) is taken up again, and the space between drains required by the Boussinesq model is calculated in order to approximate the water-table drawdown obtained with the Richards model.
The results obtained for the first case are shown in Fig. 8a and 8b. The first Figure shows the better approach of the Boussinesq equation compared to the Richards model, as to the evolution of the free surface to half of the space between drains ( 2 R 0.9761  ), obtained with a saturated hydraulic conductivity value for the Boussinesq model of s K 0.500 m h  ; this value is 2.78 times higher than the value determined in the laboratory by Zavala et al. (2004). With this, it is shown that the hypotheses considered in the derivation of the Boussinesq equation noticeably affect their s K estimation capacity. The overestimation of the s K value necessarily results in an overestimation of the evolution of the drained depth, as shown in Fig. 8b ). The difference between the real separation of the experimental drainage system considered by the Richards model and the theoretical separation required by the Boussinesq model to reproduce the lowering settled is 143%.

Conclusions
Two mechanistic models to simulate the mass and energy transfers in agricultural subsurface drainage systems have been analysed: the first model resolves the Richards equation on a two-dimensional domain by using a nonlinear radiation boundary condition in the drain's perimeter; the second model considers the Boussinesq equation with a variable storage coefficient on a one-dimensional domain using a fractal radiation condition in the drain's.
Drawing upon experimental drainage information in an unconfined aquifer with null vertical recharge, the description capacity of both simulation models has been evaluated, obtaining that the Boussinesq equation cannot simultaneously reproduce the mass and energy transfers that the Richards equation provides. On one hand, if the Boussinesq equation is used to reproduce the mass evolution described with the Richards equation, the Boussinesq necessarily describes an energy evolution slower than the one provided by Richards. On the other hand, if the energy evolution that is described by the Richards equation is reproduced with the Boussinesq equation, it over-predicts the mass evolution associated with the Richards equation.
It has been shown that the description limits of the Boussinesq equation give rise to the overestimation of the saturated hydraulic conductivity when this equation is considered in the hydrodynamic characterisation of the soils, or else the overestimation of the space between drains if the saturated hydraulic conductivity is a value estimated as by the Richards equation.
The problem in the simultaneous description the mass and energy transfers with the Boussinesq equation is attributable to the hypothesis of a hydrostatic pressure distribution considered in its derivation; this hypothesis is not satisfied in the vicinity of the drains since, in this zone, the stream lines show an important curvature.
Once the usefulness and advantages the use of the Boussinesq equation in the study of the agricultural drainage has been informed and its description limits have become known, it is recommended that the determination of the parameters of a drainage system and the estimation of the hydraulic properties of the porous media with the Boussinesq equation be executed simultaneously (considering the optimisation procedure, the drained depth evolution as well as the water table variations) in order to proportionally distribute in the adjustment parameters the effect of the hypothesis considered in its derivation and so obtain a more appropriate description of the transferences of mass and energy in a subsurface drainage system.
The numerical models have been applied by considering one laboratory drainage test. However, both models can be applied equally to the description of drainage systems installed in the field (a real farm environment), the only requirement being to carry out an adequate hydrodynamic characterisation of the soil, using either direct methods or indirect methods (inverse problems). If the hydrodynamic characteristics of the soil (moisture retention curve and hydraulic conductivity curve) are well-identified in the field and the application of both models is performed as described in this study, the results will be similar to those presented in this paper.
If the soil hydrodynamic characterisation is carried out by an inverse method, it is recommended that the procedure developed by Fuentes is applied (see Saucedo et al., 2002;Zavala et al., 2003), based on the volumetric porosity of the soil, the granulometric curve and the drainage test (local or global). This methodology takes into account the Laplace law, Stoke's law and concepts of fractal geometry. The methodology is very precise for the representation of laboratory conditions and field conditions.
The two numerical models presented in this study considers the classical hypothesis of a deterministic model, accordingly a good description of real farm conditions can be to carry out with an adequately represented of the spatial variability of the hydrodynamic properties of the soil.