- Email: [email protected]
Numerical routine for soil water dynamics from trickle irrigation
Journal Pre-proof
Numerical routine for soil water dynamics from trickle irrigation ´ Angel del Vigo , Sergio Zubelzu , Luis Juana PII: DOI: Reference:
S0307-904X(20)30066-4 https://doi.org/10.1016/j.apm.2020.01.058 APM 13271
To appear in:
Applied Mathematical Modelling
Received date: Revised date: Accepted date:
15 April 2019 13 December 2019 27 January 2020
´ Please cite this article as: Angel del Vigo , Sergio Zubelzu , Luis Juana , Numerical routine for soil water dynamics from trickle irrigation, Applied Mathematical Modelling (2020), doi: https://doi.org/10.1016/j.apm.2020.01.058
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier Inc.
Highlights
An explicit numerical routine for integrating Richards´ equation is proposed. The proposed numerical routine is validated via two different ways. Applications of proposed routine are shown. A discussion of advantage of this tool over other existing routines is given. A reference version of proposed routine is given at the end of the paper.
Numerical routine for soil water dynamics from trickle irrigation
Ángel del Vigo (a),(b)*, Sergio Zubelzu (c), Luis Juana (d).
*(a) Universidad Antonio de Nebrija. Dpto. Ingeniería industrial. C/ Pirineos 55. 28040. Madrid. Spain. Email: [email protected]. (b) Escuela Técnica Superior de Ingeniería Agronómica, Alimentaria y de Biosistemas, Universidad Politécnica de Madrid. Av. Puerta de hierro 2. 28040. Madrid. Spain. (c) Escuela Técnica Superior de Ingeniería Agronómica, Alimentaria y de Biosistemas, Universidad Politécnica de Madrid. Av. Puerta de hierro 2. 28040. Madrid. Spain. Email: [email protected] (d) Escuela Técnica Superior de Ingeniería Agronómica, Alimentaria y de Biosistemas, Universidad Politécnica de Madrid. Av. Puerta de hierro 2. 28040. Madrid. Spain. Email: [email protected]
Abstract An explicit finite differences routine has been developed to simulate the threedimensional water movement from a trickle irrigation source by direct numerical integration of Richards equation. Validation of the code was performed, on one hand, by comparison with existing analytical solutions, and on the other hand, with the Hydrus 2D/3D Software Package in order to check different conditions further than assumed analytical simplifications.
In the simulations carried out, the routine responded with accuracy and stability. Additionally, it is flexible to define different boundary conditions. Furthermore, some specific studies developed with this routine are presented in order to illustrate the applications of this tool. It is shown that the code is able to reproduce water infiltration features such as redistribution, unsteady surface puddle area, and time-dependent overpressure around a buried emitter, which are phenomena of interest difficult to simulate by the existing previous software.
Keywords 3D infiltration; Trickle irrigation; Soil water flux; Richards´ equation; Explicit finite differences method.
1. Introduction
Darcy [1] and Buckingham [2] established the foundations of water filtration for porous media, relating the velocity water field inside a soil with a potential gradient, where potential function includes the pressure (suction) and the gravitational effects. Furthermore, they introduced the concept of hydraulic conductivity, what is a function of moisture content at each point. These concepts, together with the continuity mass law, give rise to the filtration equation, also known as the Richards equation [3], what allows modelling the soil water flux.
Analytical or semi-analytical solutions to Richards equation may be obtained using a method of linearization, first proposed by Gardner [4]. It consists on defining an exponential relationship between hydraulic conductivity and matric pressure, what yields a linear form of Richards differential equation [4,5]. Despite the mathematical simplicity of this method, it presents the disadvantage of representing the soils with the aforementioned relationships, what not always fit with reality.
Otherwise, maybe the most frequently used expressions for the hydraulic conductivity and retention curves is the van Genuchten-Mualem model [6,7]. As linearization method is not applicable to these expressions, the solution of Richards equation in this case must be numerical [8], that is for example, the case of commercial package Hydrus 2D/3D.
The accurate modelling of the soil water dynamics becomes relevant for trickle irrigation systems. Trickle irrigation systems deliver water at discrete points either on the surface or, if buried, within the soil. In the first case, an area of ponded water may grow around the emitter. The shape of this area is influenced by the joint effect of the emitter discharge and the soil hydraulic characteristics. In the case of subsurface drip irrigation, soil saturation around the emitter occurs what causes a positive net pressure around the emitter that must be studied [9]. Studies which imply changing boundary conditions over time, are not easy to implement in closed programs, such as the aforementioned Hydrus 2D/3D. However, its implementation can be simpler in explicit open schemes, as the proposed in this work.
Studies based on analytical [10-12], semi-analytical [13], and numerical [14,15] solutions of Richards equation for trickle irrigation sources have been developed in the past fifty years. A detailed description of the most relevant models for
predicting soil water dynamics during trickle irrigation can be found in bibliography [16]. The applicability of analytical methods to real conditions are limited due to restrictive assumptions for source configuration, the error involved in the linearization procedure, and finally, because many of them lead to difficult formulations that include infinite series or integrals that must be numerically solve. Furthermore, analytical models do not allow the study of interesting phenomena in trickle irrigation such overpressure around the emitter, or surface puddle evolution. The most representative and established analytical formulations are the Green-Ampt [17], Warrick [10-11] and Philip [18-20] models. On the other hand, other problems are reported with existing numerical methods such as intensive computational effort, unrealistic boundary conditions or difficulties in operation for clay soils due to numerical instabilities.
Consequently, advances in this area are necessary to improve the predictability of water irrigation systems and water soil flow. In this way, the convergence for numerical approximations, soil characterization, and the evolution of puddle size in surface are examples of aspects that needs a better understanding.
In view of the above, the objectives of the work presented on this article are: (1) present an own-developed numerical routine able to simulate the soil water flux adapted to the particularities of the trickle irrigation systems, with independence of the soil characteristics, (2) validation of the proposed routine with existing analytical solutions and other numerical existing software (Hydrus 2D/3D) and (3) present some results obtained with the mentioned routine in order to show its capacity.
2. Theory
The Darcy-Buckingham equation is the more fundamental basis of water flow dynamics in a porous media [1,2]. This experimental equation express velocity field ⃗ as:
⃗
( )⃗⃗
( )⃗⃗( ( )
) ( )
) is a potential per unit of weight [m], and ( ) is a Where the function ( depending factor on medium known as hydraulic conductivity [m/s], what is a function of the water content at each point: ( ) [adim]. The potential function is usually expressed as the sum of gravity and suction potential (also ) ( ) known as matric potential): ( , either with length units, i.e., potential per unit of weight.
Using the continuity equation:
⃗⃗ ⃗ ( )
Over Darcy-Buckingham equation yields the known as Richards equation:
( ( )
)
( ( )
)
( ( )
)
( )
( )
This scalar differential equation describes the infiltration process in a medium with hydraulic conductivity ( ) and matric potential ( ). The solution of this equation provides the evolution of water content at each point of the soil: ( ).
In this work, a finite differences routine that solve this equation under a explicit scheme with adaptive boundary conditions is presented. Depending on the expression used for functions ( ) and ( ), equation (3) may be analytically solved. In this work, two different forms of these functions are used in order to check the validity of the code, first via analytical solutions, and secondly by other numerical approaches already tested.
2.1. Analytical solution Gardner supposed a hydraulic conductivity as exponential function of the matric potential in the next way:
( )
( )) ( )
(
Where [m/s] is the hydraulic conductivity under saturation conditions, and 1] is a depending of medium constant.
[m-
Hydraulic sub-saturation conductivity of equation (4), allows linearization of Richards equation, using Kirchhoff’s transformation (5). Thus, a matric flow potential, ϕ [m2/s] results [4]:
( )
( )
∫
( )
Assuming that: and , where is any constant with velocity units, the filtration equation simplifies in such a way:
( )
Equation (6) is a linear partial differential equation. Analytical solutions of this equation were provided by Carslaw & Jaeger [21] in the context of heat diffusion. Warrick [10], relying on the studies of Philip [22], derived a solution for an infiltration process since a superficial point source with a given constant flow ( ):
(
(
)) ∫
(
) [
]
( )
Equation (7) represents the dimensionless matric flow potential , where to the dimensionless buried matric flow potential calculated by Philip [22]:
( )
(
( )
(
√
√ )
(
)
(
√
√ ))
refers
( )
Equation (8) is given in terms of the Gauss complementary function error ‘erfc( x)’. Equations (7) and (8) are expressed in terms of dimensionless radius, depth and time ( ). Dimensionless variables are:
( )
In the case of a superficial saturated emitter disk [11], a matric flow potential is given, also with a constant flow (q = cte), by the next expression:
( ( ∫
(
)
( ( (
√
(
))
)[
√(
√
) )
√
√
(
))
√
) )
√(
√ (
)
(
√
√ ) ]
(
)
Where J0 and J1 are the Bessel functions and the dimensionless radius of the disk, . Again, [ ] are the dimensionless spatial-time coordinates expressed in (9). The analytical expressions given from (7) to (10) consider a z positive downwards.
Calculation of water content evolution is straightforward from matric flow potential since the aforementioned assumption: .
(
(
)) (
)
Equations (7)-(8)-(10) have been integrated with MATLAB under a specific condition in order to compare with results given for the numerical proposed routine.
2.2. Numerical solutions Numerical approximation of equation (3) mostly use the van Genuchten-Mualem relations for hydraulic conductivity and matric potential, since it is assumed that they describe the soil water movement phenomena better than the exponential functions of Gardner (4)-(5). These expressions are [6,7]:
(
(
Parameters
)
(
(
[
(
)) )
) ]
(
(
)
)
are constants of each medium. Normalized water content
is a function of residual variable (1 true, 0 false).
and saturation
water contents. (
,
) is a logical
The code was validated using commercial package Hydrus [8], that use spectral methods based on Galerkin-type linear finite element schemes, to solve equation (3) in a medium characterized by the van Genuchten-Mualem relations (12). Hydrus 2D/3D was selected because it is one of the most used and contrasted existing programs to simulate movement of water and solute in porous media [23], due to its intuitive graphics-based interface, efficiency, and success under comparison versus experimental data and/or other models [24-25].
3. Material and Methods
In this chapter, some features of proposed routine are presented. A version of the code is given in the Appendix.
3.1. Numeric scheme A finite differences explicit scheme is used to approximate Richards equation. Due to symmetry of trickle irrigation process a cylindrical system of coordinates ( ) was chosen to reduce the number of spatial variables and computational efforts. Thus equation (3) can be expressed as:
[
]
[
(
]
)
What have been approximated in terms of potential matric (see figure 1) as:
[
( (
)
((
) )(
)(
(( )
(
) )(
(
)( ))
)) ] (
)
All the terms at the right side of equation (13) are known for time . The unknown , is solved for each point of the grid with spatial resolution ( ), and time step . It may be notice that, once the matric potential is known, study of water content evolution is straightforward at each case from equations either (11) or (12.a).
Some tests were performed using, second order backward differences formula over time discretization, what allowed to use information of two previous time steps in order to calculate the unknown , obtaining a satisfactory approach.
Arithmetical media has been used on equation (14) to calculate hydraulic conductivity at each node.
1
4
2
z 3 r
Figure 1. Scheme of discretization
According to [26], an explicit integration of a parabolic differential equation is numerically stable if: ( ) . Following this idea, and the scheme shown in figure 1, the next restriction has been assumed that ensure numerical stability:
(
∑
)
(
(
)
)
Where is the spatial resolution (it has been considered for any case). Considering that a maximum value is given for saturation conditions, last equation transforms into:
(
)
(
)
Under van Genuchten-Mualem expressions ( ) tends to zero, what is not a realistic condition due to water compressibility and porosity of the media. To estimate a limit value for ( ) wells saturation flow specific capacity might be used. However, in this work it was assumed ( ) = -1 s/2000 cm , what is equivalent to approximate: near to the saturation conditions. Then, in the simulations developed, the spatial-time resolution is controlled by the next relation:
(
)
Nevertheless, in most cases was proved that numerical stability is guaranteed for a more relaxed .
It was observed that a high temporal resolution is convenient at the beginning of the simulations. For the sake of reduce computation effort, an adaptive spatial-time resolution scheme was developed, in a way that spatial-time resolution increase as simulation time progress.
3.2. Boundary conditions The emitter is placed at the origin of the r-z reference system. A summary of the boundary conditions is:
) saturation conditions are imposed: 1.- Drip region ( . Some studies were developed considering a variable source size under several sorts of geometries such as hemi-sphere or disk. In all these cases, the same boundary condition is taken for all the points involved. 2.- Rest of surface ( ; Thus, at this line.
3.- Symmetry axes (
);
); no vertical velocity upwards is considered.
is imposed (symmetry condition).
4.- Right vertical axes ( ); when a dividing line between two drippers is considered, such a symmetry axis. If a unique dripper is studied initial conditions are assumed.
5.- Bottom of the grid ( ); has been chosen far enough in order to water do not reach it. Thus initial conditions are also assumed for this line.
The size of saturation zone, , may be considered constant or variable. In the second case, when emitter discharge on a saturated area surface puddle raises with time if flow discharge is higher than infiltrated flow. The volume of ponded water is calculated by the difference between the discharged and infiltrated volumes. The stored volume on surface can be considered as a cylinder of radius rs and height hs. The roughness of the ground defines the maximum height (hr) of this cylinder, what is a necessary defined parameter in the code to get a time-dependent boundary condition for the water disc.
In the case of validation with analytical solutions, as flow rate is assumed constant (q = cte), matric potential is an unknown h around the emitter. In this case, h could be positive as a pressure potential, what affects to boundary condition flow.
3.3. Code structure
Figure 2 represents a flow chart for the developed routine. A reference version of the code is also presented in Appendix. Code is executed in MATLAB environment:
Figure 2. Flow chart. 4. Results
In this chapter validation results for the proposed routine are presented in two different forms: analytical and numerical solutions of Richards equation. Furthermore, some applications of interest are obtained to show the tool capacity.
4.1. Validation with analytical model The Sable DEK soil parameters, see table 1 [27], was used to compare analytical equations with numerical routine. Figure 3 presents the comparison between the water content profiles obtained with analytical expression (7) and numerical routine, from a superficial point source. Emitter flow used is q = 50 cm3/min (3 L/h) and time of irrigation: t = 60 min.
Integrals (7)-(8) were solved analytically for these soil parameters with MATLAB, what is the same environment where numerical routine is executed.
Table 1. Parameters used for code validation with analytical solutions. Initial water content (0) was deduced from equation (11) assuming at large distances of the emitter.
s
0
min-1)
k (cm·min-1)
ks (cm·min-1)
q (L/h)
0.3
0.125
0.073
0.557
0.0975
3
Numerical (Inf-3d.)
Analytical 0
0.14
-15 -20
-25
-25
z (cm)
-20
-30
0.14 -35
0.15
-30
0.14
0.1 3
-35 -40
-45
0. 13
-40
-45
-50
-50 0.13
-55 -60
0.13
z (cm)
-10 0.13
-15
265 000....222432 00.221 00..2 9 0.1 8 0.1 7 0.1 6 . 01
-5
0.15
-10
0.15
3 0. 254 00...22232 00. 21 00..2 9 0.1 8 0.1 7 0.1 6 0.1
-5
0.14
0
0.13
-55
0
5
10
15
20
r (cm)
(a)
25
30
35
-60
0
5
10
15
20
r (cm)
(b)
25
30
35
40
Figure 3. Irrigation simulation for a point source. Water content distribution in a Sable DEK soil irrigated from a point water source (analytical 3a, numerical 3b). q = 50 cm3/min (3 L/h) at t = 60 min.
Following the same idea, analytical and numerical results were also compared for a saturated disk source, equation (10). Figure 4 shows the comparison in this case for a disk radius of . The same parameters presented in table 1 were used:
Analytical.
Numerical (Inf-3d).
0
0
0.15
-10
7 0.1
-20
0.14
-15
-20
0.15
-25
z (cm)
0.13 0.14
-30
0.14
-35
-35
0. 13
z (cm)
-25 -30
6 0.1
-15
0.13
6 0.1
0.14
-10
2 0.2 0.21 0.2 9 0.1 8 0.1 7 0.1
-5
0.15
-5
3 0.1
3 0.2 0.221 0.2 .2 0 0.19 0.18
-40
-40
13 0.
-45
-45
-50
-50
-55
-55
0.13
-60
-60
0
5
10
15
20
25
30
35
40
0
5
10
15
20
r (cm)
r (cm)
(a)
(b)
25
30
35
40
Figure 4. Irrigation simulation for a disc source. Water content distribution in a Sable DEK soil irrigated from a disc source of r = 10 cm. Analytical (a), and numerical (b) are obtained with an introduced flow rate of: qi = 50 cm3/min (3 L/h). Time of simulation: t = 60 min.
Once the emitter stops and infiltration phase ends, the water flow keep filtering in the soil in a process that is usually known as redistribution. Since the partial differential equation (6) is linear, principles of proportionality and superposition, can be applied to study the redistribution process via analytical equations. Figure 5 shows the water content at time t = 90 min, i.e, 30 minutes after the situation of figure 3 once irrigation process had stopped (q = 3 L/h for 60 minutes). Numerical routine is able to control time of irrigation with the goal of study redistribution process.
On the other hand, to get analytical results was necessary to subtract the matric flow potential (ϕ) of 30 minutes of irrigation, from ϕ after 90 minutes of irrigation with the same flow rate, according to the principles mentioned before. A representation of the water content obtained is shown in next figure:
Analytical.
Numerical (Inf-3d).
0
0 -10
0. 12 6
0.128 -80
-80
0.127
6 12 0.
0.127 -90
-90
0.126
-100
0. 12 5
-100
-110
0.129
0.129
-70
0.128
5 0.12
0.125
-60
0.1 28
z (cm)
2 0.13 31 0.1 3 0.1
0.1 27
0.129 0.1 28 0.1 27
0. 13 3
z (cm)
-50
25 0.1
0.129
-70
-40
27 0.1
-60
0.126
2 0.13 31 0.1 3 0.1
-30
0. 12 8
0.13
5 0.12
0.13
-50
1 13 0.
-20
-30 -40
0.1 32
0.126
0. 12 8
27 0.1
-20
1 13 0.
0.1 32
0. 13 3
-10
0.12 9
6 12 0.
0.12 9
0
10
20
30
40
50
(a)
60
70
-110
0
10
20
30
40
50
60
70
r (cm)
r (cm) 25 0.1
(b)
Figure 5. Redistribution at time 90 min. Picture represents water content 30 minutes after irrigation stopped. Rate flow q = 50 cm3/min (3 L/h), during a time of 60 min. Comparison between analytical (a), and numerical (b) solutions. Sable DEK soil. 0.125
It should be noticed how in redistribution process the wetted region sinks under the surface due to the gravitational effect.
Finally, routine was tested for sub-superficial irrigation, equation (8). Figure 6 represents the water content profiles for subsurface point sources with the same rate flow of cm3/min (3 L/h) after one hour of irrigation. Again, Sable DEK soil was used in this example.
Analitical.
Numerical (Inf-3d)
0
0
0.13
-10
0.13
-10
-20
-20
z (cm)
-50
4 0.1
-50
4 0.1
-60
-70
0.1 3
-70
-80
-80
-90
-90
-100
-100
0.1 3
-60
0.23 212 00.00.02..2 0. 0.1.71819 16 0. 15
.23 0000..2 0 00.1.2212 0. .1.718 9 1 0. 15 6
z (cm)
-40
3 0.1
0.13
-40
-30
4 0.1
4 0.1
-30
0.13
0
10
20
30
r (cm)
(a)
40
50
0
10
20
30
40
50
r (cm)
(b)
Figure 6. Subsurface irrigation in Sable DEK soil. Comparison between analytical solution (a), and numerical approximation (b) after 60 minutes of irrigation. Flux rate: cm3/min (3 L/h).
After all these tests, it can be concluded that the proposed numerical routine reproduces properly the soil water dynamics according to analytical equations of Warrick and Lomen, in the case of infiltration, but also in redistribution phases. It is also proved this concordance for superficial and subsurface irrigation cases.
Analytical solutions consider an infiltrated constant flow. This hypothesis may not be true in reality. Actually, for superficial irrigation, when the emitter flow is excessive, part of the water is not infiltrated, and it keep stored over the surface. Thus, a water puddle appears and grows until the infiltrated rate flow became equal to the emitter rate. This obviously affects to the water distribution. Routine is able to describe this phenomenon what is a better approach to reality, taking in account the roughness of the floor, what also, affect to the manner in which puddle arises.
However, the greatest limitation of the analytical solutions is its definition itself, since an exponential function for the hydraulic conductivity with the matric potential, and a linear variation of the hydraulic conductivity with the water content, are two aspects difficult to assume in many realistic situations.
4.2. Numerical validation with Hydrus 2D/3D
Hydrus 2D/3D [8] was used to solve the general 3D infiltration equation over conditions presented in table 2. The soil characterization is defined by the van Genutchen-Mualem model, equations (12). Parameters used in these two equations are retrieved from Carsel and Parrish soil catalogue [28]. A summary of them are shown in the next table:
Table 2. Parameters used for validation with Hydrus Package 2D/3D. Catalogue of Carsel and Parrish [28]. Initial water content (assumed for the simulations) is .
Type soil Loam
( 0.43
0.078
0.2
0.036
)
( 0.0173
) 0.359
1.56
0.5
Figure 7 compares results provided by Hydrus-2D/3D (lines) and Inf3d (points) using a saturated constant radius disk as a water source. Several disk constants radii were studied for the soil parameters presented in table 2.
Figure 7. Infiltrated volume with time from a saturated disk. Saturated areas increase from 6 to 15 cm of radius, at the surface. Hydrus 2D/3D (lines) and the proposed numerical routine (points).
It may be note, how infiltration increases with the size of the water puddle.
Figure 7 reports a good agreement between the results obtained by Hydrus 2D/3D and the proposed numerical routine. The small observed deviations in those cases could be reduced with a minor temporal step in Hydrus 2D/3D.
4.3. Applications
The aim of this paragraph is to introduce some of the features of soil flow dynamics that can be studied with the proposed routine once its validity has been proved, such as the evolution of the puddle area, emitter overpressure or redistribution research. Some of these phenomena are difficult to analyze with other existing software. All the results presented in this section were obtained for the soil described in table 2, i.e., Loam Soil from the Carsel and Parrish Catalogue [28].
Firstly, the evolution of the source area was estimated with the proposed numerical routine. The results for different discharged flows are shown in figure 8. Parameter that defines the average roughness of the floor surface are introduced as an input in the code ( ). When emitter rate flow is higher than infiltrated flow, the water is stored over a surface cylinder of radius and height . When height of the stored water overcomes , the radius of the cylinder increases as , where is the spatial resolution over radius, assuming that water has reached the roughness limit. Then, the stored volume of the saturated cylinder spread over a bigger base cylinder that reduce in a same proportion the height of stored water . The process recurrences when stored water reach again the limit .
Figure 8. Evolution of the puddle radius as a function of infiltrated volume. Different simulations are developed under the same conditions with proposed routine, increasing the rate flow from 0.75L/h to 6L/h. Loam soil.
As expected, puddle radius arises more quickly in the case of a higher rate flow. It can be noticed, from Figure 8, how the stationary flow is reached before for the simulations where lower rate is applied.
Figure 9 presents results for redistribution, i.e., water content profile once irrigation has stopped. On this case, the same volume (6L) has been applied with two different flow rates: ( ), and ( ). It is interesting, how the simulation developed with a higher flow rate (blue) generated a bigger puddle on surface what produce a more spread redistribution of water content.
REDISTRIBUCION.
0
6 L/h 1 day after 0.75 L/h 1 day after
z (cm)
-10
-20
-30
-40
-50 0
10
20 30 r (cm)
40
50
Figure 9. Redistribution study. Wetting front is shown for . Irrigation applied volume was for either cases. Two emitters with different rate flow are presented: ( ) -red-, and ( ) -blue-. Dashed lines represent wetting front after redistribution 1 day after irrigation began.
Moreover, the proposed routine is able to calculate overpressure around the emitter, i.e., the increment of pressure around the drip due to flow rate when water infiltrates inside the ground. Overpressure increments with flow rate. Type of soil also may affect to this parameter. Figure 10 represents the overpressure, at 2 cm from a point source, as a function of time for different flow rates.
Figure 10. Overpressure (ha), at 2 cm from the point source (emitter), as a function of time. Several emitter discharges from 0.75L/h to 6L/h are analyzed. Loam soil with θ0 = 0.20. Soils of catalogue of Carsel and Parrish [28].
5. Conclusions
A novel numerical routine for simulating soil water dynamics by direct numerical integration of Richards equation is presented. It is specially designed to trickle irrigation. Validation showed a satisfactory similitude for both tests, i.e., comparison with analytical model and other numerical schemes (Hydrus-3D) for the same soil and irrigation event. In the studied examples, the proposed routine showed a good stability and convergence with a reasonable computational effort and calculation times. We have not observed relevant divergence problems within the spatial-time limit of resolution imposed.
The analytical solutions, based on the linearized form, require conductivity and retention curves that cannot represent most soils. Further, analytical solutions assume boundary conditions, as constant infiltrate flow, which do not exactly represent real situations. Actually, for superficial irrigation, when the emitter flow is excessive, part of the water is not infiltrated and stored over the surface. Thus, a water puddle appears and grows until the infiltrated rate flow equals the emitter rate. This obviously affects to the water distribution in the soil. The developed routine is able to describe this phenomenon what is a better approach to reality,
taking in account the roughness of the floor what also affect to the manner in which puddle arises.
Specific studies such as redistribution, unsteady surface puddle area, buried emitters with variable discharged flow and time-dependent overpressure around the emitter, may be performed by this routine as have been shown. On the contrary, difficulties arises for study some of these features with existing commercial programs as Hydrus 2D/3D.
The routine presented in this article is potentially useful to study puddle evolution with time from surface drip irrigation. In this way, the routine is currently used to infer the evolution of superficial water ponded area and its maximum radius (in the cases that exists), for the soils included in the three main catalogues, i.e. van Genuchten-Mualem [28], Clapp & Hornberger [29] and Gardner [30] type of soils, as a function of applied water rate. Comparison with Wooding model [13] is also aimed. Moreover, the program is being used to validate a set of analytical expressions that have been directly deduced from Darcy’s law under particular conditions and simplifications [31], similar to the Green-Ampt model [17].
As other existing software it might be used to predict water and solute front advance in a known soil, improving efficiency of trickle irrigation design.
Acknowledgment.
Author’s desire to express their gratitude to Universidad Politécnica de Madrid and Universidad Antonio de Nebrija, for economical and institutional support. Moreover, the authors wish to thank all the revision suggested, and the appreciable work of editors and reviewers in order to improve this article.
References
[1]
Darcy, H. Les fontaines publiques de la ville de Dijon. Dalmont, Paris, 1856.
[2] Buckingham, E. Studies on the movement of soil moisture. Bull. 38. U. S. Dept. of Agr. Bureau of Soils, Washington D. C., 1907. [3] Richards, L.A. Capillary conduction of liquids in porous mediums. Physics. 1(1931), 318–333. [4] Gardner, W.R. Some steady-state solutions for the unsaturated flow equation with application to evaporation from a water table. Soil Science. 85(1958), 228-232. [5] Philip, J.R. Theory of Infiltration. Advances in hydro-science. Academic Press, New York, NY. 9(1969). Pag.215. [6] van Genuchten M.T. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal. 44(5) (1980), 892-898. [7] Mualem Y. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resource journal. 12(1976), 513-522. [8] Šimůnek, J. van Genuchten M. and Šejna M. The HYDRUS Software Package for Simulating the Two- and Three-Dimensional Movement of Water, Heat, and Multiple Solutes in Variably-Saturated Media, Technical Manual Version 1.0, University of California Riverside, Riverside, CA, 3PC Progress, Prague, Czech Republic, 2006. [9] Ben-Gal, A. Lazarovitch, N. Shani, U. Subsurface drip irrigation in gravel filled cavities. Vadose Zone J. 3(2004.), 1407-1413. [10] Warrick, A.W. Time-dependent linearized infiltration I. Point sources. Soil Sci Soc Am Pro. 38(1974), 382-386. [11] Warrick, A. W. Lomen, D.O. Time-dependent Linearized Infiltration III: Strip and Disc Sources. Soil Sci Soc Am J. 40(1976), 639-643. [12] Pullan, A.J. Boundary element solutions of quasilinearized time dependent infiltration. App Math model. 12(2) (1988), 9-17. [13] Wooding, R.A. Steady infiltration from a shallow circular pond. Water resources research. 4(6) (1968), 1259-1273. [14] Bhatnagar, P.R. Chauhan, H.S. Srivastava, V.K. Unsteady unsaturated flow from a surface disc source. Journal of Hydrology. 203(1997), 154-161. [15] Revol, P. Clthier, B.E. Mailhol, J-C. Vachaud, G. Vauclin, M. Infiltration from a Surface Point Source and Drip Irrigation 2. An approximate time-dependet solution. Water Resources Research. 33(1997), 1869-1874
[16] Subaiah, R. A review of models for predicting soil water dynamics during trickle irrigation. Irrigation Science. 31(2013), 225-258. [17] Green, W. H., Ampt, G. A. Studies in soil physic, I: the flow of air and water through soils. J. Agr. Sci. 4 (1911),10 1-24. [18] Philip, J.R. Travel times from buried and surface infiltration point sources. Water Resources Research. 20 (1984a), 990-994. [19] Philip, J.R. Steady infiltration from spherical cavities. Soil Science Society American Journal. 48 (1984b), 724-729. [20] Philip, J.R. Steady infiltration from circular and cylindrical cavities. Soil Science Society American Journal. 48 (1984c), 270-278. [21] Carslaw, H.S. Jaeger, J.C. 1959. Conduction of Heat in Solids. 2nd Ed. Clarendon Press, Oxford, 1959. [22] Philip J.R. General theorem on steady infiltration from surfaces sources, with application to point and line sources. Soil Sc Soc Am Proc. 35(1971), 867-871. [23] Kandelous, M. & Simunek, J. Comparison of numerical analytical and empirical models to estimate wetting patterns to surface and subsurface drip irrigation. Irrigation Science 28 (2010a) 435-444. [24] El-Nesr, M., Alazba, A., Simunek, J. HYDRUS simulations of the effects of dual-drip subsurface irrigation and a physical barrier on the water movement and solute transport in soils. Irrigation Science. 32 (2014) 111-125. [25] Provenzano, G. Using HYDRUS-2D simulation model to evaluate wetted soil volume in subsurface drip irrigation systems. J. Irrig. Sci. Eng. 133 (2007) 342-349. [26] Carnahan, B. Luther, H.A. Wilkes, J.O. 1979. Cálculo numérico: métodos y aplicaciones. Copyright John Wliley & Sons. Editorial Rueda, Madrid, 1979. [27] Elmaloglou S, Diamantopoulos E. Simulation of soil water dynamics under subsurface drip irrigation from line sources. Agricultural Water Management 96 (2009) 1587–1595. [28] Carsel, R. Parrish, R. Developing joint probability of soil water retention characteristics. Water Resources Research. 24(5) (1988), 755-769. [29] Clapp, R.B. & Hornberger, G.M. Empirical Equations for Some Soil Hydraulic Properties. Water Resources Research. 14:4 (1978), 601-604. [30] Mualem, Y. A catalogue of Hydraulic properties of unsaturated soils. Hydrodynamics and Hydraulic Laboratory, Technion. Israel Institute of Technology, (1976). [31] del Vigo, A. et al. Soluciones analíticas aproximadas bajo hipótesis de GreenAmpt desde fuentes semiesférica y circular en superficie. JIA (2019). Toledo.
Appendix clc clear tic %Soils from Carsel & Parrish (1988). %qr,qs,a,n,ks y qo, q0= Cc -2/3(Cc-Pm),Cc(h=-1/3 atm); Pm(h=-15 atm) %Estimated with: h = 1/3 y 15 atm %Cc := soil field capacity; Pm := ground wilting point. DSH=[0.045 0.43
0.057 0.065 0.1
0.078 0.067 0.095 0.034 0.068 0.1
0.089
0.07
0.41
0.43
0.43
0.36
0.41
0.39
0.45
0.41
0.46
0.38
0.38
0.145 0.124 0.075 0.059 0.036 0.02
0.019 0.016 0.008 0.027 0.01
0.005
2.68
1.31
1.09
2.28
1.89
1.48
1.56
1.41
1.37
1.09
1.23
1.23
0.4950 0.24319444444 0.07368055556 0.02183333333 0.01733333333 0.00750 0.00433333333 0.00416666667 0.00333333333 0.0020 0.00116666667 0.00033333333 0.0452 0.0580 0.0718 0.1300 0.1133 0.1481 0.1886 0.1446 0.2955 0.2018 0.2427 0.2894];
%Soil selection j = 5; qr=DSH(1,j); qs=DSH(2,j); a= DSH(3,j); n=DSH(4,j);
m=1-1/n; ks=DSH(5,j); la=0.5;
%3D flux, axis symmetrical %Spatial resolution. hsupmax := assumed ground roughness dr = .25; dz = dr; Nr=64; Nz = 72; r = dr/2:dr:Nr*dr;z =-dz/2:-dz:-Nz*dz; mr=ones(Nz,1)*r;mz=z'*ones(1,Nr); hsupmax = 0.5; hsup = 0;
%Flow rate: qgot = 100;
%Initial and boundary condition. rmoj := initial puddle radius. qo =0.20; seo = (qo - qr) / (qs - qr); ho = -((seo.^ (-1/m) - 1) .^(1 / n) / a);
ci = 0; mh = ho-ci*mz; rmoj = dr; mh(1,:)=mh(1,:).*(r>rmoj); mha = mh; mqa = (1+(-a*mh.*(mh<0)).^n).^-m; mq = mqa*(qs-qr)+qr;mq0 = mq; v0b=sum(sum(2*pi*dz*dr*mq0.*mr));
%Integration control under stability criteria. dtmax = dz^2*qs/ks/8000; dtmin = 1e-7; rt = 1.0002; dhmax = 5;hs=-1; cmax = m*n*a*(qs-qr)*(1+(-a*hs).^n).^(-m-1).*(-a*hs)^(n-1);
%Irrigation and studying times tar = 10; tf = 10000;
vt=10^-3*[1,1.2,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,6.5,7,7.5,8,9]; vto = vt; for i =1:6; vt=[vt, vto*10^i]; end vt=[vt, 10000]; vti = zeros(length(vt),6); qsol=qo*ones(Nz, Nr, length(vt));
t = 0; ta = 0; jt = 1; qi = 0; while t=hs);
mh4(:,2:Nr) = mh(:,1:Nr-1);mh4(:,1)= mh(:,1); mqa = (1+(abs(a*mh4)).^n).^-m.*(mh4<0)+(mh4>=0); mk4 = ks*mqa.^la.*(1-(1-mqa.^(1/m)).^ m).^ 2;
mh2(:,1:Nr-1) = mh(:,2:Nr);mh2(:,Nr)= mh(:,Nr); mqa = (1+(abs(a*mh2)).^n).^-m.*(mh2<0)+(mh2>=0); mk2 = ks*mqa.^la.*(1-(1-mqa.^(1/m)).^ m).^ 2;
mh1(2:Nz,:) = mh(1:Nz-1,:);mh1(1,:)= mh(1,:); mqa = (1+(abs(a*mh1)).^n).^-m.*(mh1<0)+(mh1>=0); mk1 = ks*mqa.^la.*(1-(1-mqa.^(1/m)).^ m).^ 2;
mh3(1:Nz-1,:) = mh(2:Nz,:);mh3(Nz,:)= mh(Nz,:); mqa = (1+(abs(a*mh3)).^n).^-m.*(mh3<0)+(mh3>=0); mk3 = ks*mqa.^la.*(1-(1-mqa.^(1/m)).^ m).^ 2;
su1=(mk+mk2)/2.*(mh2-mh)/dr^2+(mk+mk4)/2.*(mh4-mh)/dr^2+mk./mr.*(mh2-mh4)/2/dr; su2=(mk+mk1)/2.*(mh1-mh)/dz^2+(mk+mk3)/2.*(mh3-mh)/dz^2+(mk1-mk3)/2/dz;
su2(1,:) = -(mk3(1,:)+mk(1,:))/dz.*((mh(1,:)-mh3(1,:))/dz+1); su = su1+su2; dht = su./mc; dhta = dht.*(dht>0)+0*(dht<0); dhtmax = max(dhta(:)); dt= max(min(dhmax/dhtmax, dtmax), dtmin); dtmin = min(dtmax*(1+10*(t>1.1*tar)),dtmin*rt);
mh = mh + dht*dt;
if trmoj); end mh(mh
vb=sum(sum(2*pi*dz*dr*mq.*mr)); vi = vb-v0b;
vsup = qgot*min(t,tar)-vi; if rmoj>0 hsup = vsup/(pi*rmoj^2); end if hsup>hsupmax rmoj = rmoj+dr;
elseif hsup<0 & t>tar rmoj = 0; end
mhaa=mha; mha = mh;
if t>vt(jt) vti(jt,1)=t; vti(jt,2)=vi; vti(jt,3)=dt; vti(jt,4)=dr; vti(jt,6)=rmoj; if jt>1 qi = (vi-vti(jt-1,2))/(t-vti(jt-1,1)); vti(jt,5)= qi; [t vi/1000 dt*1000 dr qi rmoj] end qsol(:,:,jt) = mq; jt = jt+1; end
if sum(mq(Nz/2+1+Nz/4:Nz,1))>1.1*sum(mq0(Nz/2+1+Nz/4:Nz,1)) | … sum(mq(1,Nr/2+1+Nr/4:Nr))>1.1*sum(mq0(1, Nr/2+1+Nr/4:Nr))
dr = 2*dr; dz = 2*dz; r = 2*r; z = 2*z; mr = mr*2;mz = mz*2; mh0 = ho-ci*mz; mh0(1,:)=mh0(1,:).*(r>rmoj); mqa0 = (1+(-a*mh0).^n).^-m; mq0 = mqa0*(qs-qr)+qr;
v0b=sum(sum(2*pi*dz*dr*mq0.*mr)); mh(1:Nz/2,1:Nr/2) = mh(1:2:Nz-1, 1:2:Nr-1); mh(Nz/2+1:Nz,:) = mh0(Nz/2+1:Nz,:); mh(:,Nr/2+1:Nr) = mh0(:,Nr/2+1:Nr); mqa = (1+(-a*mh.*(mh<0)).^n).^-m;mq = mqa*(qs-qr)+qr; dtmax = dtmax*3;dtmin=dtmin*3; end
end clear mk mk1 mk2 mk3 mk4 mq0 mq1 mq2 mq3 mq4 mh1 mh2 mh3 mh4 su su1 su2 su3 mc mhaa mha save inf3dr1RefRmojq100
toc
Numerical routine for soil water dynamics from trickle irrigation ´ Angel del Vigo , Sergio Zubelzu , Luis Juana PII: DOI: Reference:
S0307-904X(20)30066-4 https://doi.org/10.1016/j.apm.2020.01.058 APM 13271
To appear in:
Applied Mathematical Modelling
Received date: Revised date: Accepted date:
15 April 2019 13 December 2019 27 January 2020
´ Please cite this article as: Angel del Vigo , Sergio Zubelzu , Luis Juana , Numerical routine for soil water dynamics from trickle irrigation, Applied Mathematical Modelling (2020), doi: https://doi.org/10.1016/j.apm.2020.01.058
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier Inc.
Highlights
An explicit numerical routine for integrating Richards´ equation is proposed. The proposed numerical routine is validated via two different ways. Applications of proposed routine are shown. A discussion of advantage of this tool over other existing routines is given. A reference version of proposed routine is given at the end of the paper.
Numerical routine for soil water dynamics from trickle irrigation
Ángel del Vigo (a),(b)*, Sergio Zubelzu (c), Luis Juana (d).
*(a) Universidad Antonio de Nebrija. Dpto. Ingeniería industrial. C/ Pirineos 55. 28040. Madrid. Spain. Email: [email protected]. (b) Escuela Técnica Superior de Ingeniería Agronómica, Alimentaria y de Biosistemas, Universidad Politécnica de Madrid. Av. Puerta de hierro 2. 28040. Madrid. Spain. (c) Escuela Técnica Superior de Ingeniería Agronómica, Alimentaria y de Biosistemas, Universidad Politécnica de Madrid. Av. Puerta de hierro 2. 28040. Madrid. Spain. Email: [email protected] (d) Escuela Técnica Superior de Ingeniería Agronómica, Alimentaria y de Biosistemas, Universidad Politécnica de Madrid. Av. Puerta de hierro 2. 28040. Madrid. Spain. Email: [email protected]
Abstract An explicit finite differences routine has been developed to simulate the threedimensional water movement from a trickle irrigation source by direct numerical integration of Richards equation. Validation of the code was performed, on one hand, by comparison with existing analytical solutions, and on the other hand, with the Hydrus 2D/3D Software Package in order to check different conditions further than assumed analytical simplifications.
In the simulations carried out, the routine responded with accuracy and stability. Additionally, it is flexible to define different boundary conditions. Furthermore, some specific studies developed with this routine are presented in order to illustrate the applications of this tool. It is shown that the code is able to reproduce water infiltration features such as redistribution, unsteady surface puddle area, and time-dependent overpressure around a buried emitter, which are phenomena of interest difficult to simulate by the existing previous software.
Keywords 3D infiltration; Trickle irrigation; Soil water flux; Richards´ equation; Explicit finite differences method.
1. Introduction
Darcy [1] and Buckingham [2] established the foundations of water filtration for porous media, relating the velocity water field inside a soil with a potential gradient, where potential function includes the pressure (suction) and the gravitational effects. Furthermore, they introduced the concept of hydraulic conductivity, what is a function of moisture content at each point. These concepts, together with the continuity mass law, give rise to the filtration equation, also known as the Richards equation [3], what allows modelling the soil water flux.
Analytical or semi-analytical solutions to Richards equation may be obtained using a method of linearization, first proposed by Gardner [4]. It consists on defining an exponential relationship between hydraulic conductivity and matric pressure, what yields a linear form of Richards differential equation [4,5]. Despite the mathematical simplicity of this method, it presents the disadvantage of representing the soils with the aforementioned relationships, what not always fit with reality.
Otherwise, maybe the most frequently used expressions for the hydraulic conductivity and retention curves is the van Genuchten-Mualem model [6,7]. As linearization method is not applicable to these expressions, the solution of Richards equation in this case must be numerical [8], that is for example, the case of commercial package Hydrus 2D/3D.
The accurate modelling of the soil water dynamics becomes relevant for trickle irrigation systems. Trickle irrigation systems deliver water at discrete points either on the surface or, if buried, within the soil. In the first case, an area of ponded water may grow around the emitter. The shape of this area is influenced by the joint effect of the emitter discharge and the soil hydraulic characteristics. In the case of subsurface drip irrigation, soil saturation around the emitter occurs what causes a positive net pressure around the emitter that must be studied [9]. Studies which imply changing boundary conditions over time, are not easy to implement in closed programs, such as the aforementioned Hydrus 2D/3D. However, its implementation can be simpler in explicit open schemes, as the proposed in this work.
Studies based on analytical [10-12], semi-analytical [13], and numerical [14,15] solutions of Richards equation for trickle irrigation sources have been developed in the past fifty years. A detailed description of the most relevant models for
predicting soil water dynamics during trickle irrigation can be found in bibliography [16]. The applicability of analytical methods to real conditions are limited due to restrictive assumptions for source configuration, the error involved in the linearization procedure, and finally, because many of them lead to difficult formulations that include infinite series or integrals that must be numerically solve. Furthermore, analytical models do not allow the study of interesting phenomena in trickle irrigation such overpressure around the emitter, or surface puddle evolution. The most representative and established analytical formulations are the Green-Ampt [17], Warrick [10-11] and Philip [18-20] models. On the other hand, other problems are reported with existing numerical methods such as intensive computational effort, unrealistic boundary conditions or difficulties in operation for clay soils due to numerical instabilities.
Consequently, advances in this area are necessary to improve the predictability of water irrigation systems and water soil flow. In this way, the convergence for numerical approximations, soil characterization, and the evolution of puddle size in surface are examples of aspects that needs a better understanding.
In view of the above, the objectives of the work presented on this article are: (1) present an own-developed numerical routine able to simulate the soil water flux adapted to the particularities of the trickle irrigation systems, with independence of the soil characteristics, (2) validation of the proposed routine with existing analytical solutions and other numerical existing software (Hydrus 2D/3D) and (3) present some results obtained with the mentioned routine in order to show its capacity.
2. Theory
The Darcy-Buckingham equation is the more fundamental basis of water flow dynamics in a porous media [1,2]. This experimental equation express velocity field ⃗ as:
⃗
( )⃗⃗
( )⃗⃗( ( )
) ( )
) is a potential per unit of weight [m], and ( ) is a Where the function ( depending factor on medium known as hydraulic conductivity [m/s], what is a function of the water content at each point: ( ) [adim]. The potential function is usually expressed as the sum of gravity and suction potential (also ) ( ) known as matric potential): ( , either with length units, i.e., potential per unit of weight.
Using the continuity equation:
⃗⃗ ⃗ ( )
Over Darcy-Buckingham equation yields the known as Richards equation:
( ( )
)
( ( )
)
( ( )
)
( )
( )
This scalar differential equation describes the infiltration process in a medium with hydraulic conductivity ( ) and matric potential ( ). The solution of this equation provides the evolution of water content at each point of the soil: ( ).
In this work, a finite differences routine that solve this equation under a explicit scheme with adaptive boundary conditions is presented. Depending on the expression used for functions ( ) and ( ), equation (3) may be analytically solved. In this work, two different forms of these functions are used in order to check the validity of the code, first via analytical solutions, and secondly by other numerical approaches already tested.
2.1. Analytical solution Gardner supposed a hydraulic conductivity as exponential function of the matric potential in the next way:
( )
( )) ( )
(
Where [m/s] is the hydraulic conductivity under saturation conditions, and 1] is a depending of medium constant.
[m-
Hydraulic sub-saturation conductivity of equation (4), allows linearization of Richards equation, using Kirchhoff’s transformation (5). Thus, a matric flow potential, ϕ [m2/s] results [4]:
( )
( )
∫
( )
Assuming that: and , where is any constant with velocity units, the filtration equation simplifies in such a way:
( )
Equation (6) is a linear partial differential equation. Analytical solutions of this equation were provided by Carslaw & Jaeger [21] in the context of heat diffusion. Warrick [10], relying on the studies of Philip [22], derived a solution for an infiltration process since a superficial point source with a given constant flow ( ):
(
(
)) ∫
(
) [
]
( )
Equation (7) represents the dimensionless matric flow potential , where to the dimensionless buried matric flow potential calculated by Philip [22]:
( )
(
( )
(
√
√ )
(
)
(
√
√ ))
refers
( )
Equation (8) is given in terms of the Gauss complementary function error ‘erfc( x)’. Equations (7) and (8) are expressed in terms of dimensionless radius, depth and time ( ). Dimensionless variables are:
( )
In the case of a superficial saturated emitter disk [11], a matric flow potential is given, also with a constant flow (q = cte), by the next expression:
( ( ∫
(
)
( ( (
√
(
))
)[
√(
√
) )
√
√
(
))
√
) )
√(
√ (
)
(
√
√ ) ]
(
)
Where J0 and J1 are the Bessel functions and the dimensionless radius of the disk, . Again, [ ] are the dimensionless spatial-time coordinates expressed in (9). The analytical expressions given from (7) to (10) consider a z positive downwards.
Calculation of water content evolution is straightforward from matric flow potential since the aforementioned assumption: .
(
(
)) (
)
Equations (7)-(8)-(10) have been integrated with MATLAB under a specific condition in order to compare with results given for the numerical proposed routine.
2.2. Numerical solutions Numerical approximation of equation (3) mostly use the van Genuchten-Mualem relations for hydraulic conductivity and matric potential, since it is assumed that they describe the soil water movement phenomena better than the exponential functions of Gardner (4)-(5). These expressions are [6,7]:
(
(
Parameters
)
(
(
[
(
)) )
) ]
(
(
)
)
are constants of each medium. Normalized water content
is a function of residual variable (1 true, 0 false).
and saturation
water contents. (
,
) is a logical
The code was validated using commercial package Hydrus [8], that use spectral methods based on Galerkin-type linear finite element schemes, to solve equation (3) in a medium characterized by the van Genuchten-Mualem relations (12). Hydrus 2D/3D was selected because it is one of the most used and contrasted existing programs to simulate movement of water and solute in porous media [23], due to its intuitive graphics-based interface, efficiency, and success under comparison versus experimental data and/or other models [24-25].
3. Material and Methods
In this chapter, some features of proposed routine are presented. A version of the code is given in the Appendix.
3.1. Numeric scheme A finite differences explicit scheme is used to approximate Richards equation. Due to symmetry of trickle irrigation process a cylindrical system of coordinates ( ) was chosen to reduce the number of spatial variables and computational efforts. Thus equation (3) can be expressed as:
[
]
[
(
]
)
What have been approximated in terms of potential matric (see figure 1) as:
[
( (
)
((
) )(
)(
(( )
(
) )(
(
)( ))
)) ] (
)
All the terms at the right side of equation (13) are known for time . The unknown , is solved for each point of the grid with spatial resolution ( ), and time step . It may be notice that, once the matric potential is known, study of water content evolution is straightforward at each case from equations either (11) or (12.a).
Some tests were performed using, second order backward differences formula over time discretization, what allowed to use information of two previous time steps in order to calculate the unknown , obtaining a satisfactory approach.
Arithmetical media has been used on equation (14) to calculate hydraulic conductivity at each node.
1
4
2
z 3 r
Figure 1. Scheme of discretization
According to [26], an explicit integration of a parabolic differential equation is numerically stable if: ( ) . Following this idea, and the scheme shown in figure 1, the next restriction has been assumed that ensure numerical stability:
(
∑
)
(
(
)
)
Where is the spatial resolution (it has been considered for any case). Considering that a maximum value is given for saturation conditions, last equation transforms into:
(
)
(
)
Under van Genuchten-Mualem expressions ( ) tends to zero, what is not a realistic condition due to water compressibility and porosity of the media. To estimate a limit value for ( ) wells saturation flow specific capacity might be used. However, in this work it was assumed ( ) = -1 s/2000 cm , what is equivalent to approximate: near to the saturation conditions. Then, in the simulations developed, the spatial-time resolution is controlled by the next relation:
(
)
Nevertheless, in most cases was proved that numerical stability is guaranteed for a more relaxed .
It was observed that a high temporal resolution is convenient at the beginning of the simulations. For the sake of reduce computation effort, an adaptive spatial-time resolution scheme was developed, in a way that spatial-time resolution increase as simulation time progress.
3.2. Boundary conditions The emitter is placed at the origin of the r-z reference system. A summary of the boundary conditions is:
) saturation conditions are imposed: 1.- Drip region ( . Some studies were developed considering a variable source size under several sorts of geometries such as hemi-sphere or disk. In all these cases, the same boundary condition is taken for all the points involved. 2.- Rest of surface ( ; Thus, at this line.
3.- Symmetry axes (
);
); no vertical velocity upwards is considered.
is imposed (symmetry condition).
4.- Right vertical axes ( ); when a dividing line between two drippers is considered, such a symmetry axis. If a unique dripper is studied initial conditions are assumed.
5.- Bottom of the grid ( ); has been chosen far enough in order to water do not reach it. Thus initial conditions are also assumed for this line.
The size of saturation zone, , may be considered constant or variable. In the second case, when emitter discharge on a saturated area surface puddle raises with time if flow discharge is higher than infiltrated flow. The volume of ponded water is calculated by the difference between the discharged and infiltrated volumes. The stored volume on surface can be considered as a cylinder of radius rs and height hs. The roughness of the ground defines the maximum height (hr) of this cylinder, what is a necessary defined parameter in the code to get a time-dependent boundary condition for the water disc.
In the case of validation with analytical solutions, as flow rate is assumed constant (q = cte), matric potential is an unknown h around the emitter. In this case, h could be positive as a pressure potential, what affects to boundary condition flow.
3.3. Code structure
Figure 2 represents a flow chart for the developed routine. A reference version of the code is also presented in Appendix. Code is executed in MATLAB environment:
Figure 2. Flow chart. 4. Results
In this chapter validation results for the proposed routine are presented in two different forms: analytical and numerical solutions of Richards equation. Furthermore, some applications of interest are obtained to show the tool capacity.
4.1. Validation with analytical model The Sable DEK soil parameters, see table 1 [27], was used to compare analytical equations with numerical routine. Figure 3 presents the comparison between the water content profiles obtained with analytical expression (7) and numerical routine, from a superficial point source. Emitter flow used is q = 50 cm3/min (3 L/h) and time of irrigation: t = 60 min.
Integrals (7)-(8) were solved analytically for these soil parameters with MATLAB, what is the same environment where numerical routine is executed.
Table 1. Parameters used for code validation with analytical solutions. Initial water content (0) was deduced from equation (11) assuming at large distances of the emitter.
s
0
min-1)
k (cm·min-1)
ks (cm·min-1)
q (L/h)
0.3
0.125
0.073
0.557
0.0975
3
Numerical (Inf-3d.)
Analytical 0
0.14
-15 -20
-25
-25
z (cm)
-20
-30
0.14 -35
0.15
-30
0.14
0.1 3
-35 -40
-45
0. 13
-40
-45
-50
-50 0.13
-55 -60
0.13
z (cm)
-10 0.13
-15
265 000....222432 00.221 00..2 9 0.1 8 0.1 7 0.1 6 . 01
-5
0.15
-10
0.15
3 0. 254 00...22232 00. 21 00..2 9 0.1 8 0.1 7 0.1 6 0.1
-5
0.14
0
0.13
-55
0
5
10
15
20
r (cm)
(a)
25
30
35
-60
0
5
10
15
20
r (cm)
(b)
25
30
35
40
Figure 3. Irrigation simulation for a point source. Water content distribution in a Sable DEK soil irrigated from a point water source (analytical 3a, numerical 3b). q = 50 cm3/min (3 L/h) at t = 60 min.
Following the same idea, analytical and numerical results were also compared for a saturated disk source, equation (10). Figure 4 shows the comparison in this case for a disk radius of . The same parameters presented in table 1 were used:
Analytical.
Numerical (Inf-3d).
0
0
0.15
-10
7 0.1
-20
0.14
-15
-20
0.15
-25
z (cm)
0.13 0.14
-30
0.14
-35
-35
0. 13
z (cm)
-25 -30
6 0.1
-15
0.13
6 0.1
0.14
-10
2 0.2 0.21 0.2 9 0.1 8 0.1 7 0.1
-5
0.15
-5
3 0.1
3 0.2 0.221 0.2 .2 0 0.19 0.18
-40
-40
13 0.
-45
-45
-50
-50
-55
-55
0.13
-60
-60
0
5
10
15
20
25
30
35
40
0
5
10
15
20
r (cm)
r (cm)
(a)
(b)
25
30
35
40
Figure 4. Irrigation simulation for a disc source. Water content distribution in a Sable DEK soil irrigated from a disc source of r = 10 cm. Analytical (a), and numerical (b) are obtained with an introduced flow rate of: qi = 50 cm3/min (3 L/h). Time of simulation: t = 60 min.
Once the emitter stops and infiltration phase ends, the water flow keep filtering in the soil in a process that is usually known as redistribution. Since the partial differential equation (6) is linear, principles of proportionality and superposition, can be applied to study the redistribution process via analytical equations. Figure 5 shows the water content at time t = 90 min, i.e, 30 minutes after the situation of figure 3 once irrigation process had stopped (q = 3 L/h for 60 minutes). Numerical routine is able to control time of irrigation with the goal of study redistribution process.
On the other hand, to get analytical results was necessary to subtract the matric flow potential (ϕ) of 30 minutes of irrigation, from ϕ after 90 minutes of irrigation with the same flow rate, according to the principles mentioned before. A representation of the water content obtained is shown in next figure:
Analytical.
Numerical (Inf-3d).
0
0 -10
0. 12 6
0.128 -80
-80
0.127
6 12 0.
0.127 -90
-90
0.126
-100
0. 12 5
-100
-110
0.129
0.129
-70
0.128
5 0.12
0.125
-60
0.1 28
z (cm)
2 0.13 31 0.1 3 0.1
0.1 27
0.129 0.1 28 0.1 27
0. 13 3
z (cm)
-50
25 0.1
0.129
-70
-40
27 0.1
-60
0.126
2 0.13 31 0.1 3 0.1
-30
0. 12 8
0.13
5 0.12
0.13
-50
1 13 0.
-20
-30 -40
0.1 32
0.126
0. 12 8
27 0.1
-20
1 13 0.
0.1 32
0. 13 3
-10
0.12 9
6 12 0.
0.12 9
0
10
20
30
40
50
(a)
60
70
-110
0
10
20
30
40
50
60
70
r (cm)
r (cm) 25 0.1
(b)
Figure 5. Redistribution at time 90 min. Picture represents water content 30 minutes after irrigation stopped. Rate flow q = 50 cm3/min (3 L/h), during a time of 60 min. Comparison between analytical (a), and numerical (b) solutions. Sable DEK soil. 0.125
It should be noticed how in redistribution process the wetted region sinks under the surface due to the gravitational effect.
Finally, routine was tested for sub-superficial irrigation, equation (8). Figure 6 represents the water content profiles for subsurface point sources with the same rate flow of cm3/min (3 L/h) after one hour of irrigation. Again, Sable DEK soil was used in this example.
Analitical.
Numerical (Inf-3d)
0
0
0.13
-10
0.13
-10
-20
-20
z (cm)
-50
4 0.1
-50
4 0.1
-60
-70
0.1 3
-70
-80
-80
-90
-90
-100
-100
0.1 3
-60
0.23 212 00.00.02..2 0. 0.1.71819 16 0. 15
.23 0000..2 0 00.1.2212 0. .1.718 9 1 0. 15 6
z (cm)
-40
3 0.1
0.13
-40
-30
4 0.1
4 0.1
-30
0.13
0
10
20
30
r (cm)
(a)
40
50
0
10
20
30
40
50
r (cm)
(b)
Figure 6. Subsurface irrigation in Sable DEK soil. Comparison between analytical solution (a), and numerical approximation (b) after 60 minutes of irrigation. Flux rate: cm3/min (3 L/h).
After all these tests, it can be concluded that the proposed numerical routine reproduces properly the soil water dynamics according to analytical equations of Warrick and Lomen, in the case of infiltration, but also in redistribution phases. It is also proved this concordance for superficial and subsurface irrigation cases.
Analytical solutions consider an infiltrated constant flow. This hypothesis may not be true in reality. Actually, for superficial irrigation, when the emitter flow is excessive, part of the water is not infiltrated, and it keep stored over the surface. Thus, a water puddle appears and grows until the infiltrated rate flow became equal to the emitter rate. This obviously affects to the water distribution. Routine is able to describe this phenomenon what is a better approach to reality, taking in account the roughness of the floor, what also, affect to the manner in which puddle arises.
However, the greatest limitation of the analytical solutions is its definition itself, since an exponential function for the hydraulic conductivity with the matric potential, and a linear variation of the hydraulic conductivity with the water content, are two aspects difficult to assume in many realistic situations.
4.2. Numerical validation with Hydrus 2D/3D
Hydrus 2D/3D [8] was used to solve the general 3D infiltration equation over conditions presented in table 2. The soil characterization is defined by the van Genutchen-Mualem model, equations (12). Parameters used in these two equations are retrieved from Carsel and Parrish soil catalogue [28]. A summary of them are shown in the next table:
Table 2. Parameters used for validation with Hydrus Package 2D/3D. Catalogue of Carsel and Parrish [28]. Initial water content (assumed for the simulations) is .
Type soil Loam
( 0.43
0.078
0.2
0.036
)
( 0.0173
) 0.359
1.56
0.5
Figure 7 compares results provided by Hydrus-2D/3D (lines) and Inf3d (points) using a saturated constant radius disk as a water source. Several disk constants radii were studied for the soil parameters presented in table 2.
Figure 7. Infiltrated volume with time from a saturated disk. Saturated areas increase from 6 to 15 cm of radius, at the surface. Hydrus 2D/3D (lines) and the proposed numerical routine (points).
It may be note, how infiltration increases with the size of the water puddle.
Figure 7 reports a good agreement between the results obtained by Hydrus 2D/3D and the proposed numerical routine. The small observed deviations in those cases could be reduced with a minor temporal step in Hydrus 2D/3D.
4.3. Applications
The aim of this paragraph is to introduce some of the features of soil flow dynamics that can be studied with the proposed routine once its validity has been proved, such as the evolution of the puddle area, emitter overpressure or redistribution research. Some of these phenomena are difficult to analyze with other existing software. All the results presented in this section were obtained for the soil described in table 2, i.e., Loam Soil from the Carsel and Parrish Catalogue [28].
Firstly, the evolution of the source area was estimated with the proposed numerical routine. The results for different discharged flows are shown in figure 8. Parameter that defines the average roughness of the floor surface are introduced as an input in the code ( ). When emitter rate flow is higher than infiltrated flow, the water is stored over a surface cylinder of radius and height . When height of the stored water overcomes , the radius of the cylinder increases as , where is the spatial resolution over radius, assuming that water has reached the roughness limit. Then, the stored volume of the saturated cylinder spread over a bigger base cylinder that reduce in a same proportion the height of stored water . The process recurrences when stored water reach again the limit .
Figure 8. Evolution of the puddle radius as a function of infiltrated volume. Different simulations are developed under the same conditions with proposed routine, increasing the rate flow from 0.75L/h to 6L/h. Loam soil.
As expected, puddle radius arises more quickly in the case of a higher rate flow. It can be noticed, from Figure 8, how the stationary flow is reached before for the simulations where lower rate is applied.
Figure 9 presents results for redistribution, i.e., water content profile once irrigation has stopped. On this case, the same volume (6L) has been applied with two different flow rates: ( ), and ( ). It is interesting, how the simulation developed with a higher flow rate (blue) generated a bigger puddle on surface what produce a more spread redistribution of water content.
REDISTRIBUCION.
0
6 L/h 1 day after 0.75 L/h 1 day after
z (cm)
-10
-20
-30
-40
-50 0
10
20 30 r (cm)
40
50
Figure 9. Redistribution study. Wetting front is shown for . Irrigation applied volume was for either cases. Two emitters with different rate flow are presented: ( ) -red-, and ( ) -blue-. Dashed lines represent wetting front after redistribution 1 day after irrigation began.
Moreover, the proposed routine is able to calculate overpressure around the emitter, i.e., the increment of pressure around the drip due to flow rate when water infiltrates inside the ground. Overpressure increments with flow rate. Type of soil also may affect to this parameter. Figure 10 represents the overpressure, at 2 cm from a point source, as a function of time for different flow rates.
Figure 10. Overpressure (ha), at 2 cm from the point source (emitter), as a function of time. Several emitter discharges from 0.75L/h to 6L/h are analyzed. Loam soil with θ0 = 0.20. Soils of catalogue of Carsel and Parrish [28].
5. Conclusions
A novel numerical routine for simulating soil water dynamics by direct numerical integration of Richards equation is presented. It is specially designed to trickle irrigation. Validation showed a satisfactory similitude for both tests, i.e., comparison with analytical model and other numerical schemes (Hydrus-3D) for the same soil and irrigation event. In the studied examples, the proposed routine showed a good stability and convergence with a reasonable computational effort and calculation times. We have not observed relevant divergence problems within the spatial-time limit of resolution imposed.
The analytical solutions, based on the linearized form, require conductivity and retention curves that cannot represent most soils. Further, analytical solutions assume boundary conditions, as constant infiltrate flow, which do not exactly represent real situations. Actually, for superficial irrigation, when the emitter flow is excessive, part of the water is not infiltrated and stored over the surface. Thus, a water puddle appears and grows until the infiltrated rate flow equals the emitter rate. This obviously affects to the water distribution in the soil. The developed routine is able to describe this phenomenon what is a better approach to reality,
taking in account the roughness of the floor what also affect to the manner in which puddle arises.
Specific studies such as redistribution, unsteady surface puddle area, buried emitters with variable discharged flow and time-dependent overpressure around the emitter, may be performed by this routine as have been shown. On the contrary, difficulties arises for study some of these features with existing commercial programs as Hydrus 2D/3D.
The routine presented in this article is potentially useful to study puddle evolution with time from surface drip irrigation. In this way, the routine is currently used to infer the evolution of superficial water ponded area and its maximum radius (in the cases that exists), for the soils included in the three main catalogues, i.e. van Genuchten-Mualem [28], Clapp & Hornberger [29] and Gardner [30] type of soils, as a function of applied water rate. Comparison with Wooding model [13] is also aimed. Moreover, the program is being used to validate a set of analytical expressions that have been directly deduced from Darcy’s law under particular conditions and simplifications [31], similar to the Green-Ampt model [17].
As other existing software it might be used to predict water and solute front advance in a known soil, improving efficiency of trickle irrigation design.
Acknowledgment.
Author’s desire to express their gratitude to Universidad Politécnica de Madrid and Universidad Antonio de Nebrija, for economical and institutional support. Moreover, the authors wish to thank all the revision suggested, and the appreciable work of editors and reviewers in order to improve this article.
References
[1]
Darcy, H. Les fontaines publiques de la ville de Dijon. Dalmont, Paris, 1856.
[2] Buckingham, E. Studies on the movement of soil moisture. Bull. 38. U. S. Dept. of Agr. Bureau of Soils, Washington D. C., 1907. [3] Richards, L.A. Capillary conduction of liquids in porous mediums. Physics. 1(1931), 318–333. [4] Gardner, W.R. Some steady-state solutions for the unsaturated flow equation with application to evaporation from a water table. Soil Science. 85(1958), 228-232. [5] Philip, J.R. Theory of Infiltration. Advances in hydro-science. Academic Press, New York, NY. 9(1969). Pag.215. [6] van Genuchten M.T. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal. 44(5) (1980), 892-898. [7] Mualem Y. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resource journal. 12(1976), 513-522. [8] Šimůnek, J. van Genuchten M. and Šejna M. The HYDRUS Software Package for Simulating the Two- and Three-Dimensional Movement of Water, Heat, and Multiple Solutes in Variably-Saturated Media, Technical Manual Version 1.0, University of California Riverside, Riverside, CA, 3PC Progress, Prague, Czech Republic, 2006. [9] Ben-Gal, A. Lazarovitch, N. Shani, U. Subsurface drip irrigation in gravel filled cavities. Vadose Zone J. 3(2004.), 1407-1413. [10] Warrick, A.W. Time-dependent linearized infiltration I. Point sources. Soil Sci Soc Am Pro. 38(1974), 382-386. [11] Warrick, A. W. Lomen, D.O. Time-dependent Linearized Infiltration III: Strip and Disc Sources. Soil Sci Soc Am J. 40(1976), 639-643. [12] Pullan, A.J. Boundary element solutions of quasilinearized time dependent infiltration. App Math model. 12(2) (1988), 9-17. [13] Wooding, R.A. Steady infiltration from a shallow circular pond. Water resources research. 4(6) (1968), 1259-1273. [14] Bhatnagar, P.R. Chauhan, H.S. Srivastava, V.K. Unsteady unsaturated flow from a surface disc source. Journal of Hydrology. 203(1997), 154-161. [15] Revol, P. Clthier, B.E. Mailhol, J-C. Vachaud, G. Vauclin, M. Infiltration from a Surface Point Source and Drip Irrigation 2. An approximate time-dependet solution. Water Resources Research. 33(1997), 1869-1874
[16] Subaiah, R. A review of models for predicting soil water dynamics during trickle irrigation. Irrigation Science. 31(2013), 225-258. [17] Green, W. H., Ampt, G. A. Studies in soil physic, I: the flow of air and water through soils. J. Agr. Sci. 4 (1911),10 1-24. [18] Philip, J.R. Travel times from buried and surface infiltration point sources. Water Resources Research. 20 (1984a), 990-994. [19] Philip, J.R. Steady infiltration from spherical cavities. Soil Science Society American Journal. 48 (1984b), 724-729. [20] Philip, J.R. Steady infiltration from circular and cylindrical cavities. Soil Science Society American Journal. 48 (1984c), 270-278. [21] Carslaw, H.S. Jaeger, J.C. 1959. Conduction of Heat in Solids. 2nd Ed. Clarendon Press, Oxford, 1959. [22] Philip J.R. General theorem on steady infiltration from surfaces sources, with application to point and line sources. Soil Sc Soc Am Proc. 35(1971), 867-871. [23] Kandelous, M. & Simunek, J. Comparison of numerical analytical and empirical models to estimate wetting patterns to surface and subsurface drip irrigation. Irrigation Science 28 (2010a) 435-444. [24] El-Nesr, M., Alazba, A., Simunek, J. HYDRUS simulations of the effects of dual-drip subsurface irrigation and a physical barrier on the water movement and solute transport in soils. Irrigation Science. 32 (2014) 111-125. [25] Provenzano, G. Using HYDRUS-2D simulation model to evaluate wetted soil volume in subsurface drip irrigation systems. J. Irrig. Sci. Eng. 133 (2007) 342-349. [26] Carnahan, B. Luther, H.A. Wilkes, J.O. 1979. Cálculo numérico: métodos y aplicaciones. Copyright John Wliley & Sons. Editorial Rueda, Madrid, 1979. [27] Elmaloglou S, Diamantopoulos E. Simulation of soil water dynamics under subsurface drip irrigation from line sources. Agricultural Water Management 96 (2009) 1587–1595. [28] Carsel, R. Parrish, R. Developing joint probability of soil water retention characteristics. Water Resources Research. 24(5) (1988), 755-769. [29] Clapp, R.B. & Hornberger, G.M. Empirical Equations for Some Soil Hydraulic Properties. Water Resources Research. 14:4 (1978), 601-604. [30] Mualem, Y. A catalogue of Hydraulic properties of unsaturated soils. Hydrodynamics and Hydraulic Laboratory, Technion. Israel Institute of Technology, (1976). [31] del Vigo, A. et al. Soluciones analíticas aproximadas bajo hipótesis de GreenAmpt desde fuentes semiesférica y circular en superficie. JIA (2019). Toledo.
Appendix clc clear tic %Soils from Carsel & Parrish (1988). %qr,qs,a,n,ks y qo, q0= Cc -2/3(Cc-Pm),Cc(h=-1/3 atm); Pm(h=-15 atm) %Estimated with: h = 1/3 y 15 atm %Cc := soil field capacity; Pm := ground wilting point. DSH=[0.045 0.43
0.057 0.065 0.1
0.078 0.067 0.095 0.034 0.068 0.1
0.089
0.07
0.41
0.43
0.43
0.36
0.41
0.39
0.45
0.41
0.46
0.38
0.38
0.145 0.124 0.075 0.059 0.036 0.02
0.019 0.016 0.008 0.027 0.01
0.005
2.68
1.31
1.09
2.28
1.89
1.48
1.56
1.41
1.37
1.09
1.23
1.23
0.4950 0.24319444444 0.07368055556 0.02183333333 0.01733333333 0.00750 0.00433333333 0.00416666667 0.00333333333 0.0020 0.00116666667 0.00033333333 0.0452 0.0580 0.0718 0.1300 0.1133 0.1481 0.1886 0.1446 0.2955 0.2018 0.2427 0.2894];
%Soil selection j = 5; qr=DSH(1,j); qs=DSH(2,j); a= DSH(3,j); n=DSH(4,j);
m=1-1/n; ks=DSH(5,j); la=0.5;
%3D flux, axis symmetrical %Spatial resolution. hsupmax := assumed ground roughness dr = .25; dz = dr; Nr=64; Nz = 72; r = dr/2:dr:Nr*dr;z =-dz/2:-dz:-Nz*dz; mr=ones(Nz,1)*r;mz=z'*ones(1,Nr); hsupmax = 0.5; hsup = 0;
%Flow rate: qgot = 100;
%Initial and boundary condition. rmoj := initial puddle radius. qo =0.20; seo = (qo - qr) / (qs - qr); ho = -((seo.^ (-1/m) - 1) .^(1 / n) / a);
ci = 0; mh = ho-ci*mz; rmoj = dr; mh(1,:)=mh(1,:).*(r>rmoj); mha = mh; mqa = (1+(-a*mh.*(mh<0)).^n).^-m; mq = mqa*(qs-qr)+qr;mq0 = mq; v0b=sum(sum(2*pi*dz*dr*mq0.*mr));
%Integration control under stability criteria. dtmax = dz^2*qs/ks/8000; dtmin = 1e-7; rt = 1.0002; dhmax = 5;hs=-1; cmax = m*n*a*(qs-qr)*(1+(-a*hs).^n).^(-m-1).*(-a*hs)^(n-1);
%Irrigation and studying times tar = 10; tf = 10000;
vt=10^-3*[1,1.2,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,6.5,7,7.5,8,9]; vto = vt; for i =1:6; vt=[vt, vto*10^i]; end vt=[vt, 10000]; vti = zeros(length(vt),6); qsol=qo*ones(Nz, Nr, length(vt));
t = 0; ta = 0; jt = 1; qi = 0; while t
mh4(:,2:Nr) = mh(:,1:Nr-1);mh4(:,1)= mh(:,1); mqa = (1+(abs(a*mh4)).^n).^-m.*(mh4<0)+(mh4>=0); mk4 = ks*mqa.^la.*(1-(1-mqa.^(1/m)).^ m).^ 2;
mh2(:,1:Nr-1) = mh(:,2:Nr);mh2(:,Nr)= mh(:,Nr); mqa = (1+(abs(a*mh2)).^n).^-m.*(mh2<0)+(mh2>=0); mk2 = ks*mqa.^la.*(1-(1-mqa.^(1/m)).^ m).^ 2;
mh1(2:Nz,:) = mh(1:Nz-1,:);mh1(1,:)= mh(1,:); mqa = (1+(abs(a*mh1)).^n).^-m.*(mh1<0)+(mh1>=0); mk1 = ks*mqa.^la.*(1-(1-mqa.^(1/m)).^ m).^ 2;
mh3(1:Nz-1,:) = mh(2:Nz,:);mh3(Nz,:)= mh(Nz,:); mqa = (1+(abs(a*mh3)).^n).^-m.*(mh3<0)+(mh3>=0); mk3 = ks*mqa.^la.*(1-(1-mqa.^(1/m)).^ m).^ 2;
su1=(mk+mk2)/2.*(mh2-mh)/dr^2+(mk+mk4)/2.*(mh4-mh)/dr^2+mk./mr.*(mh2-mh4)/2/dr; su2=(mk+mk1)/2.*(mh1-mh)/dz^2+(mk+mk3)/2.*(mh3-mh)/dz^2+(mk1-mk3)/2/dz;
su2(1,:) = -(mk3(1,:)+mk(1,:))/dz.*((mh(1,:)-mh3(1,:))/dz+1); su = su1+su2; dht = su./mc; dhta = dht.*(dht>0)+0*(dht<0); dhtmax = max(dhta(:)); dt= max(min(dhmax/dhtmax, dtmax), dtmin); dtmin = min(dtmax*(1+10*(t>1.1*tar)),dtmin*rt);
mh = mh + dht*dt;
if t
vb=sum(sum(2*pi*dz*dr*mq.*mr)); vi = vb-v0b;
vsup = qgot*min(t,tar)-vi; if rmoj>0 hsup = vsup/(pi*rmoj^2); end if hsup>hsupmax rmoj = rmoj+dr;
elseif hsup<0 & t>tar rmoj = 0; end
mhaa=mha; mha = mh;
if t>vt(jt) vti(jt,1)=t; vti(jt,2)=vi; vti(jt,3)=dt; vti(jt,4)=dr; vti(jt,6)=rmoj; if jt>1 qi = (vi-vti(jt-1,2))/(t-vti(jt-1,1)); vti(jt,5)= qi; [t vi/1000 dt*1000 dr qi rmoj] end qsol(:,:,jt) = mq; jt = jt+1; end
if sum(mq(Nz/2+1+Nz/4:Nz,1))>1.1*sum(mq0(Nz/2+1+Nz/4:Nz,1)) | … sum(mq(1,Nr/2+1+Nr/4:Nr))>1.1*sum(mq0(1, Nr/2+1+Nr/4:Nr))
dr = 2*dr; dz = 2*dz; r = 2*r; z = 2*z; mr = mr*2;mz = mz*2; mh0 = ho-ci*mz; mh0(1,:)=mh0(1,:).*(r>rmoj); mqa0 = (1+(-a*mh0).^n).^-m; mq0 = mqa0*(qs-qr)+qr;
v0b=sum(sum(2*pi*dz*dr*mq0.*mr)); mh(1:Nz/2,1:Nr/2) = mh(1:2:Nz-1, 1:2:Nr-1); mh(Nz/2+1:Nz,:) = mh0(Nz/2+1:Nz,:); mh(:,Nr/2+1:Nr) = mh0(:,Nr/2+1:Nr); mqa = (1+(-a*mh.*(mh<0)).^n).^-m;mq = mqa*(qs-qr)+qr; dtmax = dtmax*3;dtmin=dtmin*3; end
end clear mk mk1 mk2 mk3 mk4 mq0 mq1 mq2 mq3 mq4 mh1 mh2 mh3 mh4 su su1 su2 su3 mc mhaa mha save inf3dr1RefRmojq100
toc