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

© 2015 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and eproduction in any medium, provided the original work is properly cited.

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 onedimensional Boussinesq equation which is a result of the continuity equation, ∂(*υH* ) / ∂*t* + ∂(*Hq*)/ ∂ *x* =*Rw*, and the Darcy's law, *q* = − *Ks*∂*H* / ∂ *x*, namely:

$$
\mu\left(H\right)\frac{\partial H}{\partial t} = \frac{\partial}{\partial \mathbf{x}} \left[T\left(H\right)\frac{\partial H}{\partial \mathbf{x}}\right] + R\_{\mathbf{w}}\tag{1}
$$

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*), *<sup>T</sup>* (*<sup>H</sup>* ) is the transmissivity given by *<sup>T</sup>* (*<sup>H</sup>* )= *KsH <sup>L</sup>* <sup>2</sup> *<sup>T</sup>* <sup>−</sup><sup>1</sup> , *Rw* is the volume of recharge in the unit of time per unit of the aquifer *L* <sup>3</sup> , *υ* =*υ*(*H* ) is the drainable porosity as a head function, and *Ks* is the saturated hydraulic conductivity *<sup>L</sup> <sup>T</sup>* <sup>−</sup><sup>1</sup> .

The storage capacity, see [24], is: *μ*(*H* )=*θ<sup>s</sup>* −*θ*(*H* −*Hs*), where *θs* is the saturated volumetric water content *L* <sup>3</sup> *<sup>L</sup>* <sup>−</sup><sup>3</sup> , and *θ*(*<sup>H</sup>* <sup>−</sup>*Hs*) 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: *θ*(*ψ*)=*θ<sup>r</sup>* + (*θ<sup>s</sup>* −*θr*) 1 + (*ψ* / *ψ<sup>d</sup>* )*<sup>n</sup>* <sup>−</sup>*<sup>m</sup>*, where *ψ* is the soil water potential defined by *ψ*=(*H* − *z*) *L* , *ψd* is the pressure scale parameter *L* , *θs* is the saturated volumetric water content *L* <sup>3</sup> *<sup>L</sup>* <sup>−</sup><sup>3</sup> , *θr* is the residual volumetric content *<sup>L</sup>* <sup>3</sup> *L* <sup>−</sup><sup>3</sup> , *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* )=(*θ<sup>s</sup>* −*θr*){1− 1 + ((*Hs* −*H* ) / |*ψ<sup>d</sup>* |)*<sup>n</sup>* <sup>−</sup>*<sup>m</sup>*}.

The saturated volumetric water content can be assimilated to the soil porosity (*ϕ*), dimension‐ less, this is calculated with the formula *<sup>ϕ</sup>* =1−*ρ<sup>t</sup>* / *<sup>ρ</sup>o*, where *ρ<sup>t</sup>* is the bulk density *<sup>M</sup> <sup>L</sup>* <sup>−</sup>3 and *<sup>ρ</sup><sup>o</sup>* is the particles density *<sup>M</sup> <sup>L</sup>* <sup>−</sup><sup>3</sup> ; the residual volumetric water content (*θr*) is considered to be zero.

#### **2.3. Initial and boundary conditions**

substance which transported by the water flow and the substance which adsorbed and

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 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 onedimensional Boussinesq equation which is a result of the continuity equation,

> ( ) ( ) ¶¶ ¶ é ù = + ê ú ¶¶ ¶ ë û *<sup>w</sup> H H H TH R*

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

in the unit of time per unit of the aquifer *L* <sup>3</sup> , *υ* =*υ*(*H* ) is the drainable porosity as a head

The storage capacity, see [24], is: *μ*(*H* )=*θ<sup>s</sup>* −*θ*(*H* −*Hs*), where *θs* is the saturated volumetric

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

*z* =*Hs*, while the free surface decreases, and *z* is the elevation of ground surface *L* .

*<sup>L</sup>* <sup>−</sup><sup>3</sup> , and *θ*(*<sup>H</sup>* <sup>−</sup>*Hs*) represents the water content evolution in the position

*tx x* (1)

*<sup>T</sup>* <sup>−</sup><sup>1</sup> , *Rw* is the volume of recharge

∂(*υH* ) / ∂*t* + ∂(*Hq*)/ ∂ *x* =*Rw*, and the Darcy's law, *q* = − *Ks*∂*H* / ∂ *x*, namely:

m

time (*t*), *<sup>T</sup>* (*<sup>H</sup>* ) is the transmissivity given by *<sup>T</sup>* (*<sup>H</sup>* )= *KsH <sup>L</sup>* <sup>2</sup>

function, and *Ks* is the saturated hydraulic conductivity *<sup>L</sup> <sup>T</sup>* <sup>−</sup><sup>1</sup> .

in the liquid phase.

water content *L* <sup>3</sup>

**2.2. The drainable porosity**

**2.1. The Boussinesq equation**

**2. Theory**

102 Agroecology

exchanges in the soil solid structure is known as the adsorption isotherm [11,12,13].

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 propor‐ tional 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]:

$$-K\_s \frac{\partial h}{\partial \mathbf{x}} \pm q\_s \left(\frac{h}{h\_s}\right)^{2s} = 0; \; \mathbf{x} = \mathbf{0}, L\mathbf{x} = \mathbf{0}, L \tag{2}$$

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 *L T* <sup>−</sup><sup>1</sup> . 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−*ϕ*)*<sup>s</sup>* + *ϕ*2*<sup>s</sup>* =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* <sup>2</sup>*<sup>s</sup>* , and the cumulative drained depth is calculated by <sup>ℓ</sup>(*t*)= <sup>1</sup> *<sup>L</sup> ∫* 0 *t Qd* (*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 onedimensional form, is a result of the continuity equation, ∂(*H CT* ) / ∂*t* + ∂*Qs* / ∂ *x* =*Rs*, and the dynamic law given by *Qs* =*HqC* −*υHDa*(∂*C* / ∂ *x*), namely:

$$\frac{\partial \left( HC\_{\tau} \right)}{\partial t} + \frac{\partial \left( HqC \right)}{\partial \mathbf{x}} = \frac{\partial}{\partial \mathbf{x}} \left[ \nu H Da \frac{\partial C}{\partial \mathbf{x}} \right] + R\_s \tag{3}$$

where *Da* is the diffusion coefficient in the water *L* <sup>2</sup> *<sup>T</sup>* <sup>−</sup><sup>1</sup> ; *CT* is the total solute concentration in soil *<sup>M</sup> <sup>L</sup>* <sup>−</sup><sup>3</sup> ; *C* is the solute concentration in water *<sup>M</sup> <sup>L</sup>* <sup>−</sup><sup>3</sup> ; 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* / *υ L T* <sup>−</sup><sup>1</sup> .

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 concentra‐ tion of the compound (*CT* ), expressed in units of mass per volume of soil can be written as: *CT* <sup>=</sup>*υ<sup>C</sup>* <sup>+</sup> *<sup>ρ</sup>tCa*, where *Ca* is the concentration of the adsorbed compound *<sup>M</sup> <sup>L</sup>* <sup>−</sup><sup>3</sup> and is a function of the concentration of the solute in the mobile phase (*Cd* ) *<sup>M</sup> <sup>L</sup>* <sup>−</sup><sup>3</sup> 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*+*<sup>γ</sup>* − *xi* / *xi*+1 − *xi* ) and *ω* =(*t <sup>j</sup>*+*<sup>ω</sup>* −*tj* / *t <sup>j</sup>*+1 −*tj* ), 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:

$$\mathbf{d}\Phi\_{i\leftrightarrow\gamma}^{j} = \left(1-\chi\right)\Phi\_{i}^{j} + \chi\Phi\_{i\leftrightarrow 1}^{j} \tag{4}$$

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

$$
\Delta \Phi\_i^{/\* \alpha} = \left(1 - \alpha \nu \right) \Phi\_i^{/i} + \alpha \nu \Phi\_i^{/\* 1} \tag{5}
$$

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

$$\left.\frac{\partial \left(\nu H C\right)}{\partial t}\right|^{i+\theta} = \frac{\left(\nu H\right)\_i^{i+1} C\_i^{j+1} - \left(\nu H\right)\_i^j C\_i^j + \left(\rho\_i H\right)\_i^{j+1} C\_{\text{id}}^{j+1} - \left(\rho\_i H\right)\_i^j C\_{\text{id}}^j}{\Delta t\_j} = b\_2 C\_i^{j+1} - b\_1 C\_i^j + b\_0; \text{ } \Delta t\_j = t\_{j+1} - t\_j \tag{6}$$

where:

from the equation (1−*ϕ*)*<sup>s</sup>* + *ϕ*2*<sup>s</sup>* =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:

The advection-dispersion equation used to study the solute transport [28,29,30], in a onedimensional form, is a result of the continuity equation, ∂(*H CT* ) / ∂*t* + ∂*Qs* / ∂ *x* =*Rs*, and the

u

in soil *<sup>M</sup> <sup>L</sup>* <sup>−</sup><sup>3</sup> ; *C* is the solute concentration in water *<sup>M</sup> <sup>L</sup>* <sup>−</sup><sup>3</sup> ; 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

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 concentra‐ tion of the compound (*CT* ), expressed in units of mass per volume of soil can be written as: *CT* <sup>=</sup>*υ<sup>C</sup>* <sup>+</sup> *<sup>ρ</sup>tCa*, where *Ca* is the concentration of the adsorbed compound *<sup>M</sup> <sup>L</sup>* <sup>−</sup><sup>3</sup> and is a function of the concentration of the solute in the mobile phase (*Cd* ) *<sup>M</sup> <sup>L</sup>* <sup>−</sup><sup>3</sup> 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.

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

¶ ¶ ¶ ¶ é ù += + ê ú ¶ ¶¶ ¶ ë û

*HC HqC C HDa R*

¯ is the integration variable.

, and the cumulative drained depth is calculated by

*s*

*<sup>T</sup>* <sup>−</sup><sup>1</sup> ; *CT* is the total solute concentration

*t xx x* (3)

*Qd* (*t*)=2 *Do* + *h* (0, *t*) *qs h* (0, *t*)/ *hs* <sup>2</sup>*<sup>s</sup>*

**2.4. Solute transport equation**

calculated by *v* =*q* / *υ L T* <sup>−</sup><sup>1</sup> .

**2.5. Numerical scheme**

¯, where *t*

dynamic law given by *Qs* =*HqC* −*υHDa*(∂*C* / ∂ *x*), namely:

*T*

where *Da* is the diffusion coefficient in the water *L* <sup>2</sup>

( ) ( )

<sup>ℓ</sup>(*t*)= <sup>1</sup> *<sup>L</sup> ∫* 0

104 Agroecology

*t Qd* (*t* ¯)*dt*

$$b\_0 = \frac{\left(\rho\_r H\right)\_i^{j+1} C\_{di}^{j+1} - \left(\rho\_r H\right)\_i^j C\_{di}^j}{\Delta t\_j}; b\_1 = \frac{\left(\nu H\right)\_i^j}{\Delta t\_j}; b\_2 = \frac{\left(\nu H\right)\_i^{j+1}}{\Delta t\_j} \tag{7}$$

The spatial derivative discretization in the continuity equation is:

$$\left. \frac{\partial \mathcal{Q} \mathbf{s}}{\partial \mathbf{x}} \right|\_{i}^{l \ast a} = \frac{\mathcal{Q} \mathbf{s}\_{i \ast \mathbf{y}}^{l \ast a} - \mathcal{Q} \mathbf{s}\_{i \ast (l \ast \mathbf{y})}^{l \ast a}}{\Delta \mathbf{x}\_{i}} ; \ \Delta \mathbf{x}\_{i} = (1 - \boldsymbol{\gamma}) (\mathbf{x}\_{i} - \mathbf{x}\_{i \ast 1}) + \boldsymbol{\gamma} \left( \mathbf{x}\_{i \ast 1} - \mathbf{x}\_{i} \right) \tag{8}$$

According with the dynamic law:

$$\left[\left.Q\mathbf{s}\right\|\_{i\leftrightarrow\gamma}^{j\leftrightarrow a} = \left(Hq\right)\_{i\leftrightarrow\gamma}^{j\leftrightarrow a}C\_{i\leftrightarrow\gamma}^{j\leftrightarrow a} - \left(\nu H\right)\_{i\leftrightarrow\gamma}^{j\leftrightarrow a}\left(Da\right)\_{i\leftrightarrow\gamma}^{j\leftrightarrow a}\frac{C\_{i+1}^{j\leftrightarrow a} - C\_i^{j\leftrightarrow a}}{x\_{i+1} - x\_i} \tag{9}$$

$$\left[\left.Q\mathbf{s}\right\|\_{l-\left(1-\gamma\right)}^{l^{+\ast a}} = \left(Hq\right)\_{i-\left(1-\gamma\right)}^{l^{+\ast a}}C\_{l-\left(1-\gamma\right)}^{l^{+\ast a}} - \left(\nu H\right)\_{i-\left(1-\gamma\right)}^{l^{+\ast a}}\left(Da\right)\_{i-\left(1-\gamma\right)}^{l^{+\ast a}}\frac{C\_{l}^{j^{+\ast a}} - C\_{l-1}^{j^{+\ast a}}}{x\_{l} - x\_{l-1}}\tag{10}$$

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

$$C\_{i\*\gamma}^{j} = \left(1 - \gamma\right)C\_i^j + \gamma C\_{i\*1}^j; \; C\_{i-\{1-\gamma\}}^j = \left(1 - \gamma\right)C\_{i-1}^j + \gamma C\_i^j \tag{11}$$

and according with the equation (5) the temporal interpolation is *Ci <sup>j</sup>*+*<sup>ω</sup>* =(1−*ω*)*Ci <sup>j</sup>* <sup>+</sup> *<sup>ω</sup>Ci j*+1 .

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

$$C\_{i \leftrightarrow \gamma}^{\prime \leftrightarrow a} = \left(1 - \alpha \right) C\_{i \leftrightarrow \gamma}^{\prime} + \alpha C\_{i \leftrightarrow \gamma}^{\prime \leftrightarrow 1} = \left(1 - \alpha \right) \left[ \left(1 - \gamma \right) C\_i^{\prime} + \gamma C\_{i \leftrightarrow 1}^{\prime} \right] + \alpha \left[ \left(1 - \gamma \right) C\_i^{\prime + 1} + \gamma C\_{i \leftrightarrow 1}^{\prime + 1} \right] \tag{12}$$

$$C\_{i - (1 - \gamma)}^{j + \alpha} = (1 - \alpha)C\_{i - (1 - \gamma)}^{j} + \alpha C\_{i - (1 - \gamma)}^{j + 1} = (1 - \alpha) \left[ (1 - \gamma)C\_{i - 1}^{j} + \gamma C\_{i}^{j} \right] + \alpha \left[ (1 - \gamma)C\_{i - 1}^{j + 1} + \gamma C\_{i}^{j + 1} \right] \tag{13}$$

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

$$\mathbf{C}\_{i+1}^{\prime \ast a} = \left(1 - \alpha \right) \mathbf{C}\_{i+1}^{\prime} + \alpha \mathbf{C}\_{i+1}^{\prime \ast 1};\\\mathbf{C}\_{i}^{\prime \ast a} = \left(1 - \alpha \right) \mathbf{C}\_{i}^{\prime} + \alpha \mathbf{C}\_{i}^{\prime \ast 1};\\\mathbf{C}\_{i-1}^{\prime \ast a} = \left(1 - \alpha \right) \mathbf{C}\_{i-1}^{\prime} + \alpha \mathbf{C}\_{i-1}^{\prime \ast 1} \tag{14}$$

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

$$\left. \frac{\partial \mathcal{Q} \mathbf{s}}{\partial \mathbf{x}} \right|\_{i}^{i+a} = a\_1 \mathbf{C}\_{i+\mathbf{y}}^{j+a} - a\_2 \left( \mathbf{C}\_{i+1}^{j+a} - \mathbf{C}\_i^{j+a} \right) - a\_3 \mathbf{C}\_{i-\{1-\boldsymbol{\gamma}\}}^{j+a} + a\_4 \left( \mathbf{C}\_i^{j+a} - \mathbf{C}\_{i-1}^{j+a} \right) \tag{15}$$

where:

$$a\_1 = \frac{\left(Hq\right)\_{i+\mathcal{I}}^{i+o}}{\Delta \mathbf{x}\_i}; a\_2 = \frac{\left(\nu H\right)\_{i+\mathcal{I}}^{i+o} \left(Da\right)\_{i+\mathcal{I}}^{i+o}}{\Delta \mathbf{x}\_i \left(\mathbf{x}\_{i+1} - \mathbf{x}\_i\right)}; a\_3 = \frac{\left(Hq\right)\_{i-\left(1-\gamma\right)}^{i+o}}{\Delta \mathbf{x}\_i}; a\_4 = \frac{\left(\nu H\right)\_{i-\left(1-\gamma\right)}^{i+o} \left(Da\right)\_{i-\left(1-\gamma\right)}^{i+o}}{\Delta \mathbf{x}\_i \left(\mathbf{x}\_i - \mathbf{x}\_{i-1}\right)}\tag{16}$$

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

$$\begin{split} \left. \frac{\partial \mathcal{Q}s}{\partial \mathbf{x}} \right|\_{i}^{i+\alpha} &= -\alpha \left[ a\_4 + (1-\gamma)a\_3 \right] C\_{i-1}^{j+1} + \alpha \left[ (1-\gamma)a\_1 + a\_2 - \gamma a\_3 + a\_4 \right] C\_i^{j+1} + \alpha \left[ \gamma a\_1 - a\_2 \right] C\_{i+1}^{j+1} \\ &- (1-\alpha) \left[ a\_4 + (1-\gamma)a\_3 \right] C\_{i-1}^{j} + (1-\alpha) \left[ a\_4 - \gamma a\_3 + a\_2 + (1-\gamma)a\_1 \right] C\_i^{j} \\ &+ (1-\alpha) \left[ \gamma a\_1 - a\_2 \right] C\_{i+1}^{j} \end{split} \tag{17}$$

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

$$\rm{As}\_{i}\rm{C}\_{i-1}^{j+1} + \rm{Bs}\_{i}\rm{C}\_{i}^{j+1} + \rm{Ds}\_{i}\rm{C}\_{i+1}^{j+1} = \rm{Es}\_{i}; \; i = \rm{2,3}, \ldots, n-1 \tag{18}$$

where

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

gg

and according with the equation (5) the temporal interpolation is *Ci*

 w

g

 gg

 w

gg

 w

 w  ww

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

( )

1 2 3 4

1 1

+

*j i*

 ww

 gg

w

g

( ) ( ) ( )

u

w

w

w

w

+

*j*

*i*

w

DD

( )[ ]

w

w g

+- -

1

w

w

+

*j*

*i*

12 1

*a aC*

 g

 g

g

¶

¶

defined by:

106 Agroecology

w

g

w

defined by:

where:

ing:

g

w

( ) ( ) 1 1 ( ) <sup>1</sup> 1 1

g

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

( ) ( ) ( ) ( ) <sup>1</sup> 1 1 1 1 1 1 1 1

> g

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

( ) ( ) ( ) 11 1 1 11 1 11 1 11

 ww

 <sup>+</sup> + + + + <sup>+</sup> + ++ =- + =- + =- + - --

= - -- + -

ww

( ) ( ) 1 21 34 1 <sup>1</sup> ( )

 D

*a ;a ; a ; a x xx x <sup>x</sup> xx x* (16)

( ) ( ) [ ]

4 31 12 34 12 1

 g

*j j i i*

+ ++ + ++ + ++ - - -- --

Substituting equations (12)-(14) in equation (15) and associating similar terms allows obtain‐

4 31 4 32 1

*a aC a a a aC*

 wg

*j jj j jj i ii i ii i ii i i ii i*

*Hq H Da Hq H Da*

== = = - -

( ) ( ) ( ) ( )


*Qs a aC a a a a C a aC <sup>x</sup>*

11 1 1

¶ =- é + - ù + é - + - + ù + - ë ûë <sup>û</sup> ¶

 wg


*j jj j jj i ii i i i*

+ ++ + ++ + + - - -

 w

*Qs aC a C C aC a C C <sup>x</sup>* (15)

1 1

1 1 1

+ + + - +

*j j j i i i*

+ -

g

w

g

*j jj j j <sup>j</sup> j jj C C C ;C C C ;C C C i ii i i <sup>i</sup> i ii* (14)

 g

*<sup>i</sup> i i i i i i C CC CC C C* (12)

 g

( ) ( ) ( ) ( ) ( ) ( ) ( ) <sup>1</sup> 1 1 1 11 1 1 1 1 1 1

 <sup>+</sup> <sup>+</sup> + + - - -- -- - - =- + =- - + + - + <sup>é</sup> ù é <sup>ù</sup> <sup>ë</sup> û ë <sup>û</sup> *<sup>j</sup> j j j j j j i ii i i i i C CC CC C C* (13)

 w

 <sup>+</sup> <sup>+</sup> + + + ++ <sup>+</sup> <sup>+</sup> =- + =- - + + - + <sup>é</sup> ù é <sup>ù</sup> <sup>ë</sup> û ë <sup>û</sup> *<sup>j</sup> j j j j j j*

 g

++ - - - =- + =- + *<sup>j</sup> j jj j j C C C ;C i ii <sup>i</sup> C C i i* (11)

 wg

 w

 ww

( ) ( ) ( ) ( ) ( ) ( )

u

1 11

 D w g

 g ( )

ww

gg

(17)

 wg

 g  g

*<sup>j</sup>*+*<sup>ω</sup>* =(1−*ω*)*Ci*

 g

> g

 w *<sup>j</sup>* <sup>+</sup> *<sup>ω</sup>Ci j*+1 .

$$A\mathbf{s}\_i = -\alpha \mathbb{E}\left[a\_4 + \left(1 - \gamma\right)a\_3\right] \tag{19}$$

$$B\mathbf{s}\_{i} = a\overline{o}\left[ (1-\gamma)a\_{1} + a\_{2} - \gamma a\_{3} + a\_{4} \right] + b\_{2} \tag{20}$$

$$D\mathbf{s}\_{i} = o\left[\gamma a\_{1} - a\_{2}\right] \tag{21}$$

$$\begin{split} \mathbf{E}\mathbf{s}\_{i} &= \mathbf{R}\mathbf{s}\_{i}^{\prime + \alpha} + (\mathbf{l} - \alpha) \left[ \mathbf{a}\_{4} + (\mathbf{l} - \gamma) \mathbf{a}\_{3} \right] \mathbf{C}\_{i-1}^{\prime} \\ &- \left\langle \left( \mathbf{l} - \alpha \right) \left[ \mathbf{a}\_{4} - \gamma \mathbf{a}\_{3} + \mathbf{a}\_{2} + \left( \mathbf{l} - \gamma \right) \mathbf{a}\_{1} \right] - b\_{1} \right\rangle \mathbf{C}\_{i}^{\prime} - \left( \mathbf{l} - \alpha \right) \left[ \gamma \mathbf{a}\_{1} - \mathbf{a}\_{2} \right] \mathbf{C}\_{i+1}^{\prime} - b\_{0} \end{split} \tag{22}$$

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*(*C* −*Cext* / *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.
