Open access

Soil Salinity Control in Irrigated Land with Agricultural Drainage Systems

Written By

Carlos Chávez, Carlos Fuentes, Fernando Brambila and Nami Morales-Durán

Submitted: June 19th, 2014 Published: June 10th, 2015

DOI: 10.5772/59370

Chapter metrics overview

1,635 Chapter Downloads

View Full Metrics

1. Introduction

Soil salinity is a worldwide problem as well as in Central and Northern Mexico. Nearly 8.4 million ha worldwide are affected by soil salinity and alkalinity, of which about 5.5 million ha are waterlogged [1]. The problem worsens in arid and semiarid areas, in soils with insufficient drainage [2] and high evaporation [3]. In Mexico, there are 6.46 million ha irrigated mainly in the center and north territory areas [4]; 10-30 % of irrigated land is affected by salinity and nearly two thirds of this area is located in the North [5].

The salinization of these irrigated areas is an increasing problem and the lands are abandoned; therefore, a technical and economic alternative to recover this land is needed. Agricultural subsurface drainage is a solution which takes into account the technology by environment interaction, as well as lowering the water table levels along with the salt concentration in the soil profile [1].

The dynamics of water drainage systems has been studied by applying the Boussinesq equation [6] for unconfined aquifers using the finite element technique [7,8] and the finite difference [9,10], and the solutes dynamics has been studied by applying the Fick’s law [11,12,13]. These results in the advection-dispersion equation, namely by gravity and Fick’s law.

The solutes are also found in the gas phase and adsorbed by soil in the solid phase, the first phase is disregarded for purposes of transport modeling in water, but it is really important in terms of the amount of fertilizer transferred into the atmosphere at a given time [14,15,16], and incorporating the adsorbed substance in the solid phase. The relationship between the substance which transported by the water flow and the substance which adsorbed and exchanges in the soil solid structure is known as the adsorption isotherm [11,12,13].

A large number of models for simulating solute transport in the unsaturated zone are now increasingly being used for a wide range of applications in both research and management [17], some of the more popular models include SWAP [18], HYDRUS-1D [19], STANMOD (STudio of ANalytical MODels) [20], UNSATH [21] and COUP [22] but the majority of applications for water flow in the vadose zone requires a numerical solution of the Richards equation [23], also requires more calculation time in order to find the equation solution.

This study aims to solve the one-dimensional advection-dispersion equation using the technique of finite differences, coupled with the Boussinesq equation in order to model the transport of solutes in subsurface drainage systems, assuming that the solute is concentrated in the liquid phase.


2. Theory

2.1. The Boussinesq equation

In the study of the water dynamics in agricultural subsurface drainage systems using the Boussinesq equation, the variations in hydraulic head along the drain pipes (direction y) are negligible with respect to head variations in the cross section (direction x). It is the one-dimensional Boussinesq equation which is a result of the continuity equation, (υH)/t +(Hq)/x = Rw, and the Darcy’s law, q=KsH/x, namely:


where μ(H) is the storage capacity, H=H(x,t) is the elevations of the free surface or hydraulic head above the impervious layer [L], and is a function the horizontal coordinate (x) and the time (t), T(H) is the transmissivity given by T(H)=KsH [L2T1], Rw is the volume of recharge in the unit of time per unit of the aquifer [L3], υ=υ(H) is the drainable porosity as a head function, and Ks is the saturated hydraulic conductivity [LT1].

The storage capacity, see [24], is: μ(H)=θsθ(HHs), where θs is the saturated volumetric water content [L3L3], and θ(HHs) represents the water content evolution in the position z=Hs, while the free surface decreases, and z is the elevation of ground surface [L].

2.2. The drainable porosity

To calculate the storage capacity and the drainable porosity it is necessary to provide the soil water retention curve. The model of van Genuchten [25] was accepted in field and laboratory studies: θ(ψ)=θr+(θsθr)[1+(ψ/ψd)n]m, where ψ is the soil water potential defined by ψ=(Hz) [L], ψd is the pressure scale parameter [L], θs is the saturated volumetric water content [L3L3], θr is the residual volumetric content [L3L3], m and n are parameters (dimensionless) that determine the shape of the soil water retention curve. The introduction of this equation in the storage capacity results in the following expression for storage capacity: μ(H)=(θsθr){1[1+((HsH)/|ψd|)n]m}.

The saturated volumetric water content can be assimilated to the soil porosity (ϕ), dimensionless, this is calculated with the formula ϕ=1ρt/ρo, where ρt is the bulk density [ML3] and ρo is the particles density [ML3]; the residual volumetric water content (θr) is considered to be zero.

2.3. Initial and boundary conditions

To study the agricultural drainage with equation (1), the initial and boundary conditions should be defined at the domain. The initial condition is established from the water table position at the initial time. Dirichlet and Neumann boundary type conditions can be used on drains to solve equation (1), the pressure head on the drains is required in the first condition whereas the drainage flux is required in the second one [8]. A third type of boundary condition is a linear combination of the precedent conditions; this condition includes a resistance parameter to the flow at the soil–drain interface. Null resistance corresponds to the Dirichlet condition and infinite resistance corresponds to Neumann condition. The third condition is a radiation type condition [26]. In the case of drainage, the radiation condition establishes that drainage flux is directly proportional to the pressure head on the drain and inversely proportional to the resistance in the interface between soil and the drainpipe wall in concordance to the Ohm law.

The hydraulic head measured above the impermeable barrier H(x,t) is associated with the head h(x,t) measured from above the drains using: H(x,t)=Do+h(x,t), where Do is the distance from the impermeable barrier to the drains [L]. Transversal variation of h at the beginning is considered as the initial condition h(x,0)=hs(x), where hs is the head on the drain in the initial time [L]. The fractal radiation condition for the Boussinesq equations is given by [8]:


where the positive sign corresponds to x=0 and the negative sign to x=L. L is the distance between drains; qs is the corresponding flux to hs and it is function of the soil-drain interface characteristic [LT1]. For the s parameter, the authors argued that it is defined by s=D/E, where D is the effective fractal dimension to the soil-drain interface, and E=3 is the Euclidean dimension of physical space. The relation of the s parameter and effective porosity is obtained from the equation (1ϕ)s+ϕ2s=1 given by [27]. Equation (2) contains as particular cases the lineal radiation condition when s=1/2 and the quadratic radiation condition when s=1. In a system of parallel drains, the drained water flows by length unit at each drain is: Qd(t)=2[Do+h(0,t)]qs[h(0,t)/hs]2s, and the cumulative drained depth is calculated by (t)=1L0tQd(t¯)dt¯, where t¯ is the integration variable.

2.4. Solute transport equation

The advection-dispersion equation used to study the solute transport [28,29,30], in a one-dimensional form, is a result of the continuity equation, (HCT)/t +Qs/x =Rs, and the dynamic law given by Qs=HqC υHDa(C/x), namely:

(HCT)t +(HqC)x =x[υHDaCx] +RsE3

where Da is the diffusion coefficient in the water [L2T1]; CT is the total solute concentration in soil [ML3]; C is the solute concentration in water [ML3]; and Rs is the term which includes gains or losses of the solute due to chemical reactions and the extraction plant [M]. Note that q and υ are obtained from the water flow model. The diffusion coefficient in the water is calculated by Da=λv, where λ is the dispersivity [L] and v the interstitial velocity of water calculated by v=q/υ [LT1].

The water soluble compounds which have a negligible vapor pressure can exist in three phases of soil: 1) dissolved in water, 2) as vapor in the soil atmosphere and 3) as stationary phase adsorbed to soil organic matter or in the clay mineral surfaces [11,12,13]. The total concentration of the compound (CT), expressed in units of mass per volume of soil can be written as: CT=υC+ρtCa, where Ca is the concentration of the adsorbed compound [ML3] and is a function of the concentration of the solute in the mobile phase (Cd) [ML3] and the adsorption constant of the solute to the stationary phase surface (κ), Ca=κCd, namely linear isotherm. Thus, the concentration of the substance compared to the volume of the porous medium (CT) will be the result of a part that is in the water, air and the dynamic equilibrium with the phase that generates it. Generally, in studies in small time scales, such as irrigation and drainage in a porous medium, the gas phase is not considered [29]. Thus, in this work, the concentration in the adsorbed phase, the concentration in the gas phase and the term Rs are ignored.

2.5. Numerical scheme

The numerical scheme used is based on the assumption that the solute is concentrated mainly in the liquid phase. Thus, the advection-dispersion equation in one-dimensional is given by equation (3). To solve this equation, we use the same discretization scheme to transfer water in the Boussinesq equation [10], for which two interpolation parameters are introduced: γ=(xi+γxi/xi+1xi) and ω=(tj+ωtj/tj+1tj), where 0γ1 and 0ω1; i=1,2,... and j=1,2,... are the space and time indices, respectively.

The dependent variable (Φ) in an intermediate node i+γ for all j is estimated as:


while the intermediate time j+ω for all i is estimated as:


The discretization of the temporal derivative in the equation (3) is:




The spatial derivative discretization in the continuity equation is:


According with the dynamic law:


According with the equation (4), the spatial interpolation is:


and according with the equation (5) the temporal interpolation is Cij+ω=(1ω)Cij+ωCij+1.

The dependent variables involved in the advective term of the equations (9) and (10) are defined by:


while the dependent variables involved in the dispersive term of the same equations are defined by:


Equation (8) considering equations (9) and (10) can be written as:




Substituting equations (12)-(14) in equation (15) and associating similar terms allows obtaining:


Substituting equations (6) and (17) in the continuity equation, the following algebraic equations system is obtained:

AsiCi1j+1+BsiCij+1+DsiCi+1j+1 = Esi;i=2,3,...,n1E18


Asi= ω[a4+(1γ)a3]E19
Esi = Rsij+ω+(1ω)[a4+(1γ)a3]Ci1j{(1ω)[a4γa3+a2+(1γ)a1]b1}Cij(1ω)[γa1a2]Ci+1j b0E22

The water flow and the head are obtained from the Boussinesq equation solution, so that they should be included in the system (18). To find the solution of the water transfer equation, it is necessary to specify the initial and boundary conditions, equation (18) can be solved with the Thomas Algorithm, see [31,10].

The Thomas algorithm, also known as the tridiagonal matrix algorithm (TDMA), is a simplified form of Gaussian elimination that can be used to solve tridiagonal matrix systems (equation 18) [32]. It is based on LU decomposition in wich the matrix system Mx=r where L is a lower triangular matrix and U is an upper triangular matrix. The system can be efficiently solved by setting Ux=p and then solving first Lp=r for p and then Ux=p for x. The Thomas algorithm consist of two steps. In the first step decomposing the matrix into M=LU and solving Lp=r are accomplished in a single downwards sweep, taking us straight from Mx=r to Ux=p. In the second step the equation Ux=p is solved for x in an upwards sweep [33].

2.6. Linear radiation condition

The radiation boundary condition, or mixed condition, is used to accept a linear variation between the dispersive flux and concentration difference with the external medium (Cext) and the border, for all time. The linear radiation condition is due originally to Newton, who postulated that the heat flow at the border of a body is proportional to the temperature difference between the body and the medium that surrounds it; the result is equivalent to a Ohm law into electricity. To linearize these conditions, we introduce a generalization of the dimensionless conductance coefficient (κs), as follows: (C/x)+κs(CCext/L)=0. If we observe the one-dimensional equation of solute transport, the dimensionless conductance coefficient (κs) must be zero by the advective component, however, the solution is allowed only for purposes of illustration to derive the boundary conditions.

2.7. Selection of the space (Δx) and time (Δt) increments

In reference [10] the authors discuss the selection of spatial and temporal increments pointing out a comparison of the depletion of the free surface for all time between the results obtained with the finite difference solution of the Boussinesq equation and the results obtained with an analytical solution reported in the literature. The same authors [10] concluded that the optimal interpolation that minimizes the sum of the squares errors are γ=0.5Δx (cm) and ω=0.98Δt (h), for space and time respectively.


3. Application

3.1. Laboratory experiment

To evaluate the descriptive capacity of the numerical solution, a drainage experiment was conducted in a laboratory. The drainage module (see Figure 1) is the one used by [8] and [10]. The module dimensions are: L=100 cm, Hs=120 cm and Do=25 cm. The drain diameter is d=5 cm and the drain length is =30 cm.

Figure 1.

Drainage module

The module was filled with altered sample of salty soil of Celaya, Guanajuato, México (see Figure 2). Soil passed through a 2 mm sieve and was disposed on 5 cm thick layers, in order to maintain the bulk density at a constant value. The soil was saturated by applying a constant water head (no salt) 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 volume; the initial condition was equivalent to h(x,0)=hs and the recharge was null Rw=0 during the drainage phase. Soil porosity (ϕ) was calculated with the formula ϕ=1ρt/ρo (the bulk density was determined by the weight and volume of the soil of drainage module ρt=1.14 g/cm3 and the particles density ρs=2.65 g/cm3, ϕ=0.5695 cm3/cm3 was obtained). The soil fractal dimension obtained was equal to 0.7026.

Figure 2.

Study site

3.2. Analysis of the salt content

During the module drainage process (154 h), measurements of pH, temperature and electrical conductivity of water samples were made at defined time intervals (each hour the first 20 hours and subsequently increased to the range 2, 4, 6 and 8 h). The sensor used for measurement is a CONDUCTRONIC PC 18 sensor. The electrical conductivity at room temperature was recorded with it. However, in order to accurately quantify conductivity, it is important to consider a standard value of 25° C, which can be used to correct the values obtained. The correction factor used in accordance with [34] is 2-3 % for every Celsius degree that is measured under standard temperature. According with [34] (1964), the relationship between electrical conductivity and concentration is:


where C is the concentration given in mg/l and EC the electrical conductivity given in dS/m or mmhos/cm.

3.3. The hydrodynamic characteristic

To solve the Boussinesq equation, the van Genuchten model [25] for the water retention curve was used, along with a model of hydraulic conductivity of Fuentes [27] namely geometric mean model {K(Θ)=Ks[1(1Θ1/m)sm]2} with the restriction 0<sm=12s/n<1; where Θ is the effective saturation defined by Θ=(θθr)/(θsθr).

3.4. The granulometric curve

The m and n form parameters from the water retention curve are obtained from the granulometric curve [35] adjusted with the equation F(D)=[1+(Dg/D)N]M, where F(D) is the cumulative frequency, based on the weight of the particles whose diameters are less than or equal to D; Dg is a characteristic parameter of particle size, M and N are two form empirical parameters. These parameters are rewritten as follows: M=m and N=[1/2(1s)]n.

3.5. Inverse problem

To evaluate the capacity of the numerical solution of the Advection-Dispersion Equation, the experimental information presented by [36] is used. The characteristics of the drainage module and the soil parameters used in the simulation are: hs=120 cm, D0=25 cm, L=100 cm, ϕ=0.5695 cm3/cm3, and s=0.7026. The hydrodynamic characteristics used are those of van Genuchten and Fuentes [25,27]. The scale parameters (ψd,Ks) are obtained from the inverse problem, using the experimental drained depth and the drained depth calculated with the numerical solution of the Boussinesq equation [10], given an error criterion between the previous and the new estimator (1x1012cm), using a constant head test and fractal radiation condition with variable storage capacity and a nonlinear optimization algorithm [37]. The calculations were performed on a dual-core AMD Opteron machine with 2.6 GHz CPU and 8 GB RAM. The computational time required to solve the inverse problem was 5 h.

In order to model the salt concentration in the soil profile, with the numerical solution of the solute transport, the hydraulic parameters obtained from the previous analysis were used. In the numerical solution, the unknown parameter is the dispersivity coefficient (λ), which is estimated by minimizing the sum of squares errors between the salt concentration measured and the salt concentration calculated with the numerical solution over time, using a Levenberg-Marquardt [37], given an error criterion between the previous and the new estimator (1x109g/l). The initial condition is the sample initial, taken as a constant in all the system and radiation as the boundary condition applied in the drains.


4. Results and discussion

4.1. The granulometric curve

The adjusted parameters are shown in Table 1. Figure 3 shows the experimental granulometric curve and best fit is obtained with Dg=36.2993 μm and m=0.3410 with a root mean square error RMSE=0.1477.

Model Ajusted parameters
Ks ψd κ RMSE
(cm/h) (cm) (non-dimensional) (cm)
Geometric mean model 1.5458 143.87 0.0616 0.2195

Table 1.

Values of the adjusted parameters from the granulometric curve and the drained depth.

Figure 3.

The experimental granulometric curve and adjusted with the model

4.2. The hydrodynamic characteristic

In order to obtain the values of ψd and Ks, the spatial and temporal increments used in all the simulation are Δz=0.0010 cm and Δt=5x105 h. Figure 4 shows the experimental drained depth and the drained depth calculated with the finite difference solution [10], using a storage capacity variable, fractal radiation condition in the drains and the geometric mean model. To linearize the boundary condition, one generalization of the conductance coefficient is optimized (κ) [8,10]. The residual volumetric water content is considered to be zero (θr=0.0 cm3/cm3) [38].

Figure 4.

Comparison between the experimental drained depth and the calculated drained depth.

4.3. Analysis of the salt content

The EC data are shown in Figure 5 using a 2.5% like correction factor. Applying equation (23) to the data shown in Figure 5, we obtain the concentration in grams per liter (see Figure 6). The initial condition using in the numerical solution is the sample initial (Cini=2.4g/l), taken as a constant in all the system and radiation as the boundary condition applied in the drains. The dispersivity value obtained is λ=91.80 cm, with RMSE = 0.1063g/l between the experimental values and the values obtained from the numerical solution. The computational time required to solve the advection-dispersion model was 2.7 h. The dispersivity value found is only for this soil, because this value change with depth [39], increase with the flow rate and is a soil type function. This increase was explained by the activation of large pores at higher flow rates [40]. Figure 6 shows the experimental salt concentration evolution and the concentration obtained with the numerical solution.

Comparison shows that the salt concentration obtained with the numerical solution, according to RMSE, reproduce the experimental salt concentration. Figure 7 shows that in the short time, when the water flow increased, the salt concentration increases sharply, and in the long time tends to an asymptote, indicating that the system could not continue removing salts from the system. However, the value of the dispersivity obtained (λ=91.80 cm) overestimates the measured data in the long time. Second simulation was performed with the accumulated mass. To obtain the accumulated mass, it is necessary to obtain the solute flow, which is estimated by multiplying the water flow by the measured salt concentration in the time interval (see Figure 7). The cumulative solute mass is obtained by multiplying the solute flow by the time interval (see Figure 8).

Figure 5.

Evolution of the electrical conductivity of drainage water

Figure 6.

Comparison between the experimental and the calculated drainage water salt concentration with numerical solution

The results obtained with the numerical solution, the solute flow and cumulative mass evolution are shown in Figures 7 and 8, respectively, which demonstrate that the reproductions of the data were acceptable. The solute flow decreases rapidly, as seen in Figure 7 the concentration decreased 3.5 g/l after 20 hours. In the long time, the theoretical water flow and experimental water flow tends to be constant. Comparison showed that the solute flow and the cumulative mass evolution obtained with the numerical solution, according to RMSE, reproduce the experimental salt concentration. The RMSE values for estimating the solute flow and cumulative solute mass were 0.1842g/l and 0.1104g, respectively. The dispersivity value obtained is λ=98.03 cm, with RMSE=0.1010g between the experimental values and the values obtained from the numerical solution. The dispersivity value for this new optimization (cumulative mass evolution) compared to the previous (salt concentration evolution) increases 6.2 cm.

Figure 7.

Solute flow (g/h) in the drainage system.

Figure 8.

Cumulative mass evolution: experimental and obtained with the numerical solution.

4.4. Using the solution to simulate the leaching of saline soils

To recover saline soils it is necessary to apply irrigation so that the salts are transported to deeper horizons without harming the roots and are evacuated to other areas through the drainage channel. For purposes of illustrating the leaching of salts in the soil by applying the finite difference solution, we assumed a soil with hydraulic and hydrodynamic characteristics previously found. The initial soil concentration was 10 dS/m and the problem was reduced to finding the number of irrigation that must be applied to carry a given concentration.

The final average concentration obtained in the profile at the end of the first simulation was the initial concentration in the system for the next simulation, and so on. Figure 9 shows the reduced concentration of salts in the soil profile based on an initial concentration. The values shown are an average concentration in the soil profile at 1-m depth. Depth of drains was assumed to be 2.0 m.

Figure 9.

Evolution of the salt concentration in the soil by applying the leaching.

The simulations were performed with 5, 10, 15, 20 and 25 m of drains distances. It can be seen that the decrease in the concentration of salts in the soil profile is similar in all the separations between drains after applying 6 leachings. However, the time of drainage in each system was different. For example, with 5 days and 5 m of separation a decrease of the water table profile was more than one meter, while in the system with separation of 25 m decreased gives only a few centimeters (see Figure 10), other simulations was realized with 5, 10, 15, 20 and 25 m of drains distances, but the depth of drains was assumed to be 1.5 m (see Figure 11), therefore the time of drainage of the soil was a function of the distance between drains

Figure 10.

Decrease of the midpoint water table at different separations between drains under a drain depth of 2.00 m.

Figure 11.

Decrease of the midpoint water table at different separations between drains under a drain depth of 1.50 m.


5. Conclusions

The irrigation in the arid and semi-arid regions to sustain agricultural production against the unpredictable of the rainfall have resulted in the double problem of salinity in many hectares of good agricultural land. Subsurface drainage systems are used to control the depth of the water table and to reduce or prevent soil salinity.

The advection-dispersion equation was solved in order to model the temporal evolution of the concentration of salts removed through an agricultural drainage system with the method of finite differences. The solution requires the values of the flow of water previously obtained from the solution of the Boussinesq equation. The hydrodynamic characteristic were obtained by the inverse problem from the depth drained.

The optimization of the accumulated mass which gave better results in terms of mean square error criterion between the theoretical and experimental values, since it is a property integrated in the time and concentration observed at specific levels. The solution presented coupled to the Boussinesq equation, satisfactorily reproduced the measured data, both in the short time where the change in concentration was high, as in long times where the concentration values tended to an asymptote. This asymptotic value of the concentration depended on the distance between drains of the drainage system.

Finally, the solution of differential equations of transfer processes of water and solute transport, and hydrodynamic characterization of the soil in an agricultural drainage system, will be a useful tool for designing new systems for the optimal development of crops according to their water needs and the degree of tolerance to salinity. In addition, this study can be help us for quantify crop yield reductions due to salinity on irrigation areas, in order to prevent future problems such as food shortages.


  1. 1. Ritzema HP., Satyanarayama TV., Raman S., Boonstra J. Subsurface Drainage to Combat Waterlogging and Salinity in Irrigated Lands in India: Lessons Learned in Farmers’ Fields. Agricultural Water Management 2008; 95 (3) 179-189.
  2. 2. Mousavi1 SF., Mostafazadeh-Fard B., Farkhondeh A., Feizi M. Effects of Deficit Irrigation with Saline Water on Yield, Fruit Quality and Water Use Efficiency of Cantaloupe in an Arid Region. Journal of Agricultural Science Technology 2009; 11 469-479.
  3. 3. Ruiz-Cerda E., Aldaco NR., Montemayor TJ., Fortis HM., Olague RJ., Villagómez GJ. 2007. Aprovechamiento y mejoramiento de un suelo salino mediante vermicomposta. Tecnologías Pecuarias en Mexico 2007; 45 (1) 19-24.
  4. 4. CONAGUA. Statistics on Water in Mexico (in Spanish). México. D.F; 2010.
  5. 5. IMTA. Manual de Diseño e Instalación de Drenaje Parcelario en Zonas áridas y Semiáridas bajo Riego. México; 1998.
  6. 6. Boussinesq J. 1904. Recherches The´oriques sur L’e´coulement des Nappes d’eau Infiltre´es Dans le Sol et Sur le De´bit des Sources. J. Math. Pure. Appl., 1904; 10 5-78.
  7. 7. Verhoest N., Pauwels V., Troch P., de Troch F. Analytical Solution for Transient Water Table Heights and Outflows from Inclined Ditch-Drained Terrains. J. Irrig. Drain Eng., 2002; 128(6) 358–364.
  8. 8. Zavala M., Fuentes C., Saucedo H. Nonlinear Radiation in the Boussinesq Equation of Agricultural Drainage. J. Hydrol. 2007; 332(3) 374-380.
  9. 9. Singh S., Ghosh NC., Pandey RP., GalkateRV., Thomas T., Jaiswal RK. Numerical Solution of 1D Boussinesq Equation for Water Table Fluctuation between Drains in Response to Recharge and ET in A Sloping Aquifer. Int. J. Eco. Econ. Stat. 2009; 14(9) 45-54.
  10. 10. Chávez C., Fuentes C., Zataráin F., Zavala M. Finite Difference Solution of the Boussinesq Equation with Variable Drainable Porosity and Fractal Radiation Boundary Condition. Agrociencia 2011; 45(8) 911-927.
  11. 11. Taylor GI. The Disperson of Matter in Turbulent Flow Through a Pipe. Proc. R. Soc. London, Ser. A. 1954; 223 446–48.
  12. 12. Elder JW. The Dispersion of Marked Fluid in Turbulent Shear Flow. J. Fluid Mech., 1959; 5 544-560.
  13. 13. Fischer HB. The Mechanics of Dispersion in Natural Streams. J. Hydraul. Div., Am. Soc. Civ. Eng., 1967; 93(6) 187-216.
  14. 14. Holly FM. Two Dimensional Mass Dispersion in Rivers. Hydrologic papers, Colorado State University Press, Fort Collins, Colorado, 1975; pp. 78.
  15. 15. Holly FM. Dispersion in Rivers and Coastal Waters, 1, Physical Principles and Dispersion Equations. In: Novak, P. (ed). Developments in Hydraulic Engineering, 3. Elsevier, New York. 1985; Chap. 1: 1-38.
  16. 16. Rutherford JC. River mixing, Wiley, New York; 1994.
  17. 17. Mirabzadeh M., Mohammadi K. A Dynamic Programming Solution to Solute Transport and Dispersion Equations in Groundwater. J. Agric. Sci. Technol. 2006; 8 233-241.
  18. 18. van Dam JC., Huygen J., Wesseling JG., Feddes RA., Kabat P., van Valsum PEV., Groenendijk P., van Diepen CA. Theory of SWAP, Version 2.0. Simulation of Water Flow, Solute Transport and Plant Growth in the Soil- Water-Atmosphere- Plant Environment, Department Water Resources, WAU, Report 71, DLO Winand Staring Centre, Technical Document 45, Wageningen; 1997.
  19. 19. Simunek J., Sejna M., van Genuchten MTh. The HYDRUS-1D software package for simulating the movement of water, heat, and multiple solutes in variably saturated media, version 2.0, United States Salinity Laboratory, USDA-ARS, Riverside, Calif; 1998.
  20. 20. Simunek J., van Genuchten MTh., Sejna M., Toride N., Leij FJ. The STANMOD Computer Software for Evaluating Solute Transport in Porous Media Using Analytical Solutions of Convection-Dispersion Equation, Versions 1.0 and 2.0, IGWMC – TPS-71, International Ground Water Modeling Center, Colorado School of Mines: Golden; 1999.
  21. 21. Fayer MJ. UNSAT-H Version 3.0: Unsaturated Soil Water and Heat Flow Model. Theory, User Manual, and Examples. Pacific Northwest National Laboratory 13249; 2000, USA, 184 pp.
  22. 22. Jansson PE., and Karlberg L. Coupled Heat and Mass Transfer Model for Soil-Plant-Atmosphere Systems, Royal Institute of Technology, Department of Civil and Environmental Engineering: Stockholm; 2001.
  23. 23. Richards LA. Capillary conduction of liquids through porous mediums. Physics 1931; 1 318-333.
  24. 24. Fuentes C., Zavala M., Saucedo H. Relationship between the Storage Coefficient and the Soil-Water Retention Curve in Subsurface Agricultural Drainage Systems: Water Table Drawdown. J. Irrig. Drain. Eng., 2009; 135(3) 279-285.
  25. 25. van Genuchten MTh. A Closed-Form Equation for Predicting the Hydraulic Conductivity of the Unsaturated Soils. Soil Sci. Soc. Amer. J. 1980; 44 892-898.
  26. 26. Carslaw HS., Jaeger JC. Conduction of Heat in Solids. Oxford University Press, Oxford; 1959.
  27. 27. Fuentes C., Brambila F., Vauclin M., Parlange J-Y., Haverkamp R. Fractal modeling of Hydraulic Conductivity in Non-Saturated soils. Hydraul. Eng. México 2001; 16(2) 119-137.
  28. 28. Abassi F., Simunek J., van Genuchten MTh., Feyen J., Adamsen FJ., Hunsaker DJ., Strelkoff TS., Shouse P. Overland Flow and Solute Transport: Model Development and Field-Data Analysis. J. Irrig. Drain. Eng. 2003; 129(2) 71–81.
  29. 29. Zerihun D., Furman A., Warrick AW,. Sánchez CA. Coupled Surface-Subsurface Solute Transport Model for Irrigation Borders and Basin. I. Model Development. J. Irri. Drain. Eng. 2005; 131(3) 396-406.
  30. 30. Simunek J. Models of Water Flow and Solute Transport in the Unsaturated Zone. Encyclopedia of Hydrological Sciences. Edited by M G Anderson; John Wiley & Sons, Ltd., Chichester, England, 2005; 1171-1180
  31. 31. Zataráin F., Fuentes C., Palacios VOL., Mercado E., Brambila F., Villanueva N. Modelación del Transporte de Agua y Solutos en el Suelo (In Spanish). Agrociencia 1998; 32(4) 373-383.
  32. 32. Freund RW., Hoppe RW. Stoer/Bulirsch: Numerische Mathematik 1. Springer-Lehrbuch, Germany; 2007.
  33. 33. Conte SD., De Boor C. Elementary Numerical Analysis: An Algorithmic Approach. McGraw-Hill, New York; 1980.
  34. 34. Villareal E., Bello S. The concentration and electrical conductivity in aqueous solutions of electrolytes. Rev. Mex. Fis. 1964; 13(2) 55-74.
  35. 35. Fuentes C. Approche Fractale des Transferts Hydriques Dans les Sols No-saturés. Tesis de Doctorado, Universidad Joseph Fourier de Grenoble, Francia; 1992.
  36. 36. Chávez C. Solución numérica de las ecuaciones de transferencia de agua y solutos en riego y drenaje. Dr. in Eng. Thesis, Universidad Autónoma de Querétaro (In spanish). México; 2010.
  37. 37. Marquardt DW. An Algorithm for Least-Squares Estimation of Nonlinear Parameters. SIAM J. Appl. Math. 1963; 11 431-441.
  38. 38. Haverkamp R., Leij FJ., Fuentes C., Sciortino A., Ross PJ. Soil Water Retention: I. Introduction of a Shape Index. Soil Sci. Soc. Am. J. 2005; 69 1881-1890.
  39. 39. Simunek J., van Genuchten MTh. Using the HYDRUS-1D and HYDRUS-2D Codes for Estimating Unsaturated Soil Hydraulic and Solute Transport Parameters. pp. 1523–1536. In M.Th. van Genuchten. Simunek, J. and "Sejna, M. (ed.) Characterization and measurement of the hydraulic properties of unsaturated porous media. University of California, Riverside, CA; 1999.
  40. 40. Feyen J., Jacques D., Timmerman A., Vanderborght J. Modelling Water Flow and Solute Transport in Heterogeneous Soils: A Review of Recent Approaches. J. Agric. Eng. Res. 1998; 70 231-256.

Written By

Carlos Chávez, Carlos Fuentes, Fernando Brambila and Nami Morales-Durán

Submitted: June 19th, 2014 Published: June 10th, 2015