Base equations

In the computer program, the complete or simplified form of the one-dimensional Boussinesq equation of agricultural drainage is solved. In general it can be written as:

where H = H(x, t)= hydraulic head measured from an impervious layer or a reference level (L), and is a function of the horizontal coordinate x and time t; T(H) = aquifer transmissivity (L ² T ^-1), in the case of unconfined aquifers it is T(H) = K _s H; K _s = saturated hydraulic conductivity (LT ^-1); μ(H) = storage coefficient that can be a function of the hydraulic head H; R _w (t) = is the volume of recharge or discharge in the unit of time per unit area of the aquifer (LT ^-1).

The recharge or discharge R _w (t) is a difficult variable to determine since it depends on the flow conditions existing in the soil surface (rain, evaporation, water over soil), the water flow in the soil and the humidity regime in the vadose zone, the stratigraphy of the unsaturated zone of the soil, macroporosity, cracks, etc. In this new version of the numerical simulation model the term R _w (t) was incorporated representing it by a time-dependent polynomial:

where a _r , b _r , c _r and d _r = coefficients to be calculated considering experimental measurements or estimates of the recharge time evolution for a specific event.

The relationship between the storage coefficient and the soil-water retention curve established by ^{Fuentes et al. (2009)} was considered:

where θ _s = saturated volumetric water content (L ³ L ^-3); and H _s = elevation of the soil surface (L). It is possible to obtain an explicit analytical representation of the storage coefficient from equation (3) if the soil-water retention curve is known. The classic van Genuchten relationship (van Genuchten, 1980) was selected, because it is widely accepted in field and laboratory studies because of its descriptive flexibility:

where θ _r = residual volumetric water content (L ³ L ^-3); ψ _d < 0 = scale parameter of the soil-water pressure potential ψ, (L); m and n = dimensionless shape parameters. This new version of the computer program included mechanistic relations between m and n that satisfy the infiltration theory, the classic Burdine restriction (^{Burdine, 1953}) " m = 1 - 2/n " as well as the fractal relationships of ^{Fuentes, Brambila, Vauclin, Parlange & Haverkamp (2001)} that derive in their study of the hydraulic conductivity of unsaturated soils. Fractal constraints are the so-called neutral pore " m = (1 - 4s/n)/s ", geometric pore " m = (1 - 2s/n)/s " and big pore " m = (1 - 4s/n)/2s "; where s = quotient dimension of the fractal object, as the soil is considered in that study, "s" being the ratio between the fractal dimension of the object " D _f " (^{Mandelbrot, 1982}; ^{Falconer, 2014}) and the Euclidean space "E", s = D _f /E.

By entering equation (4) into equation (3) the following relationship was obtained:

The fractal radiation condition of ^{Zavala et al. (2007)} was retained in this second version of the computer program:

where q _d = drain flow; γ = dimensionless conductance coefficient; K _in = √(K _s K _d ) = soil-drain interface conductivity, K _s = saturated hydraulic conductivity of the soil, K _d = drain wall conductivity; P = depth of drains; H _d = hydraulic head in the drain position; D _o = aquifer thickness; s- = (1/2)(s ₁+s _d ) = quotient dimension of the soil-drain interface; s ₁= quotient dimension of the soil; and s _d = quotient dimension of the drain wall. The relationship between areal porosity (μ _areal ) and volumetric porosity (φ) is obtained according to the probabilistic idea as µ _areal = φ ^s φ^s= φ ^2s , it is necessary that volumetric porosity is:

and areal porosity:

The conductivity of the drain wall is obtained with a formula based on Poiseuille's law: K _d = (1/2)(g/ⱱ)µ _{areal_d} (R _HD )², where g = gravity acceleration; ⱱ = water kinematic viscosity; R _HD = hydraulic radius of the drain. The option to directly manage the linear radiation condition (s- = 0.5 in equation 6) was introduced in version 2.0 of DRENAS without using the s ₁ and s _d routines.

Simplified numerical modeling of agricultural drainage can be done by imposing Dirichlet type boundary conditions (hydraulic head in the drain position). In this new version of the computer program, this boundary condition was incorporated as an alternative to the radiation condition (6); the next function that depends on the time was selected:

where x _dren = horizontal coordinate where the drain is located; y a _d , b _d , c _d and d _d = coefficients that must be calculated from measurements taken of the evolution of the hydraulic head on the drain. Equation (9) has the flexibility to describe extreme behaviors of the hydraulic head above the drain (rise and fall water level).

Finally, to perform the numerical simulation with equation (1), the initial state of the hydraulic head in the system must be defined, in this work the option of using up to a third degree polynomial that is a function of the horizontal coordinate x is included:

where a _p , b _p , c _p and d _p = coefficients to be determined from measurements of the hydraulic head in the system at the initial moment of calculation or simulation.

In the numerical solution of the system (1)-(10) the Galerkin type finite element method was used for spatial discretization, a finite difference scheme for temporal discretization, Picard's iterative method for linearization of the resulting system and a Preconditioned conjugate gradient method for the solution of the system of algebraic equations (^{Zavala et al., 2014}; ^{Zienkiewicz,‎ Taylor, & Zhu, 2013}; ^{Noor & Peters, 1987}). The system of equations resulting from discretization was programmed in Visual Basic 2017.

It is known that simplified calculations of separation between drains and groundwater abatement can be performed by applying analytical solutions obtained for reduced forms of the Boussinesq equation (equation 1). For example, for the steady flow regime there is the Hooghoudt formula (^{Hooghoudt, 1940}) and in unsteady flow the relations of Glover-Dumm (^{Dumm, 1954}) and ^{Fuentes et al. (1997)}, which were already included in the original version of the computer program.

Graphical interface development

In the second version of the DRENAS computer program, all the information capture modules of the original version were reprogrammed in Visual Basic 2017, independent internal databases of Microsoft Access were developed, the processing of the graphics was disconnected from Microsoft Excel, and the executable file was generated to install DRENAS 2.0 on any computer that has a 64-bit Windows Operating System and even the option for 32-bit Windows Systems was also generated.

The DRENAS 2.0 program has two calculation modules, one called analytical solutions and another numerical model (Figure 1). The analytical solutions module has four sections, two to perform simulations for steady water flow conditions using the Hooghoudt solution (^{Hooghoudt, 1940}) and two for unsteady flow conditions, one to perform calculations with the Glover-Dumm solution (^{Dumm, 1954}) and the other with the solution of ^{Fuentes et al. (1997)}. In these four sections it is possible to calculate the separation between drains or the drainage module. Regarding the original version of the DRENAS model, in version 2.0 these sections were improved by incorporating the following alternatives: a) Storage of simulations performed in an internal database; b) Loading of previous examples; c) Generation of a simulation report and alternative to export it to a pdf file; d) The results graph can also be exported as an image. An example of a calculation section of the analytical solutions module is presented in Figure 2.

Figure 1 General screen of the DRENAS 2.0. 

Figure 2 A calculation section of the analytical solutions module. 

The numerical solution module was programmed in a Tab Index form that contains six sections: four are developed for the capture of data related to the drainage system, drain characteristics, soil properties and data necessary for numerical simulation such as those related time control, spatial discretization, boundary conditions and graphical display of results; The fifth tab serves to store the sample data in an internal database and the sixth section contains the button to execute the numerical simulation.

The new capabilities of the numerical simulation module that were programmed were established as follows:

First, from the UNSODA 2.0 database (^{Nemes, Schaap, Leij, & Wösten, 2001}), soil-water retention curves were selected (208 soils), were adjusted with Equation (4) considering the mechanistic restrictions between m y n described above (^{Zavala, Saucedo, & Fuentes, 2018}) and the results obtained were included in an internal database of the computer program (Figure 3). If the option to use data from this program base is selected, when choosing the type of soil required, the sections of physical, hydraulic properties and storage coefficient are automatically filled with their corresponding parameters. The improvement to highlight in this new version is that the relationships between m and n of Burdine models, neutral pore, geometric pore and big pore can be managed (Figure 4).

Figure 3 Numerical module: expanded UNSODA base. 

Figure 4 Numerical module: soil-water properties. 

In the numerical solution of the Boussinesq equation, the recharge term (R _w (t)) was incorporated as a function of time (equation 2), the necessary programming was carried out to be able to use the Dirichlet type boundary condition (equation 9) and the initial condition described by equation (10) was also programmed, and in Figure 5 the new section developed to include these three conditions is presented.

Figure 5 Numerical module: Limit conditions and recharge. 

A new section was developed to record the simulation data, load data to compare simulations as well as to select the type of graph to be displayed in real time that does not depend on external programs (Figure 6).

Figure 6 Numerical module: simulation graphs. 

Validation 1

An analytical solution of the Hooghoudt type presented in ^{Fragoza et al. (2003)} that considers linear radiation (s- = 0.5 in equation 6) in the drains to describe the evolution of the free surface in an subsurface drainage system with constant recharge (R _o ) and steady-state water flow:

where H _c = hydraulic head to center between drains. If defined h(x,t) = H(x,t) - D _o and consider the linear radiation condition:

where h _d = hydraulic head above the drain (x=0 and x=L); D _o = aquifer thickness; and γ = conductance coefficient of the linear radiation condition used in the drains. Analyzing equation (2) it is deduced for this case that a _r =b _r =c _r =0 and d _r =R _o .

The data of the drainage system are those reported by ^{Fragoza et al. (2003)} and correspond to a field experiment conducted in the irrigation district 076 Valle del Carrizo, Sinaloa. The drainage system evaluated has a separation between drains L = 50 m, drain depth P = 1.5 m, depth of the impervious layer D _o = 3.5 m, elevations of the impervious layer and the soil surface H _i =0.0 m and H _s =5.0 m. The soil is a clay with a saturation hydraulic conductivity of K _s = 0.557 m/d and a storage coefficient µ = 0.1087. It is assumed γ= 1.5 and h _c = 0.5 m, and with equation (12) allows to obtain h _d = 0.365 m and the recharge value obtained with the condition for steady-state flow that the drain discharge is equal to the recharge is R _o = 0.000944 m/d. With these values, equation (11) was applied to obtain the spatial distribution of the hydraulic head along the separation between drains.

The described scenario was reproduced with the computer program, modeling the unsteady process of water flow in the system until reaching the steady state, which corresponds to that described with the analytical solution. Initial condition of constant hydraulic head was assumed H(x, t=0) = 4.50 m (a _p = b _p = c _p = 0 and d _p = 4.50 m in the equation 10) and to have the same conditions of the analytical solution, the linear radiation boundary condition was used in the program, s-= 0.5 in the equation (6). The hypothesis of equality between the conductivity of the soil and the conductivity of the drain wall was considered K _s = K _d , which implies K _in = K _s = K _d . The conductance coefficient value reported by ^{Fragoza et al. (2003)}, to adapt it to the form of radiation that is handled in DRENAS 2.0, getting γ= 0.045. Spatial discretization was performed with the number of nodes N = 1000 (Δx = 0.05 m) and in the discretization of time was used Δt _ini = 1.157E-05 d, Δt _max = 6.94E-04 d and Δt _min = 1.157E-06 d. The simulation time was 720 days.

This comparison is evidence that the information capture modules of the numerical solution of the program and in particular the calculation module that handles the term of the vertical recharge, initial condition and boundary condition, operate correctly and are free of errors of programming; it is also observed that the simulated series are stable and free of numerical oscillations (Figure 7).

Figure 7. Simulation of an agricultural drainage scenario (constant storage coefficient and constant vertical recharge).  

Validation 2

^{Zavala, Fuentes & Saucedo (2004)} conducted an experiment in a drainage module in which two drains of 30 cm in length (l _D ) and 5 cm in diameter (D _D ) were installed, with N _o = 233 circular perforations of diameter d _o = 0.158 cm evenly distributed on the surface of each tube. The dimensions of each drain: L = 100 cm, P = 120 cm and D _o = 25 cm. The module was filled with sandy textured soil that was saturated by applying a constant water head on the soil surface until the air was removed. With the drains covered, excess water was removed on the surface and covered to eliminate evaporation. Finally, the drains were uncovered and the volume of drained water measured for 10 days (240 hours); it should be noted that the initial condition corresponds to H(x,0) = P + D _o (d _p = 145 cm in the equation 10 and the rest of its coefficients are null) and the recharge is null (R _w = 0 in the equation 2).

Soil porosity was Φ₁ = 0.5396 cm³/cm³ and with this value the equation (7) was solved obtaining s ₁ = 0.7026. The saturation hydraulic conductivity value is K _s = 18.3 cm/h. The aquifer storage coefficient is described with equation (5) for the parameters θ _s = Φ₁, θ _r = 0, ψ _d = -41.8 cm and n = 2/(1-m) = 3.19 (^{Burdine, 1953}).

The porous areal of the drain is µ _{areal_d} = A _o /A _d = 0.0097 cm²/cm² and when solving the relationship (8) s _d = 0.5688. With the hydraulic radius of the circular holes in the drain wall R _HD = 0.0395 cm, water kinematic viscosity ⱱ = 36 cm²/h, Poiseuille's law is applied to calculate the conductivity of the drain wall K_d = 2,670.8 cm/h. The conductivity values in the soil-drain interface and the relative dimension are K_in = 221.08 cm/h and s-= 0.6357. The coefficient γ of fractal radiation reported by ^{Zavala et al. (2004)} was reinterpreted to adapt it to the form of radiation that is handled in DRENAS 2.0 (equation 6), obtaining γ= 0.0749.

The variable storage coefficient and the initial condition of total saturation of the porous medium cause the simulation to be executed by starting small time steps at the beginning of the simulation, which minimizes problems of convergence and numerical stability. The time steps used are Δt _ini = 2.77E-04 h, Δt _ini = 2.77E-05 h and Δt _max = 0.2 h. The number of nodes to discretize the space is 200.

In Figure 8 (a and b), the simulation data capture process and the evolution of the drained depth are presented. In Figure 8b it can be seen that the computer program successfully reproduces the experimental data of the drained depth, having a stable and oscillation-free numerical solution, which is an indicator that the variable storage coeficient module works correctly. The time and space steps used in the simulation are adequate since the calculated total drained depth (23.915 cm) differs from the measured depth (23.92 cm) in only 0.03%.

Figure 8 Simulation of an agricultural drainage scenario (variable storage coefficient, null recharge and fractal radiation condition in the drains).