- Email: [email protected]
ECOUL: an interactive computer tool to study hydraulic behavior of swelling and rigid soils
Computers & Geosciences 28 (2002) 1107–1118
ECOUL: an interactive computer tool to study hydraulic behavior of swelling and rigid soils$ Edith Perrier*, Patricia Garnier1, Christian Leclerc2 Unite de Recherches Geodes, Centre IRD Ile de France, 32 avenue Henri Varagnat, 93143 Bondy Cedex, France Received 15 June 2001; received in revised form 19 March 2002; accepted 25 March 2002
Abstract ECOUL is an interactive, didactic software package which simulates vertical water flow in unsaturated soils. Endusers are given an easily-used tool to predict the evolution of the soil water profile, with a large range of possible boundary conditions, through a classical numerical solution scheme for the Richards equation. Soils must be characterized by water retention curves and hydraulic conductivity curves, the form of which can be chosen among different analytical expressions from the literature. When the parameters are unknown, an inverse method is provided to estimate them from available experimental flow data. A significant original feature of the software is to include recent algorithms extending the water flow model to deal with deforming porous media: widespread swelling soils, the volume of which varies as a function of water content, must be described by a third hydraulic characteristic property, the deformation curve. Again, estimation of the parameters by means of inverse procedures and visualization facilities enable exploration, understanding and then prediction of soil hydraulic behavior under various experimental conditions. r 2002 Elsevier Science Ltd. All rights reserved. Keywords: Water flow modeling; Richards equation; Deforming soils; Unsaturated hydraulic properties; Inverse method
1. Introduction Water flow in soils has been described for nearly a century by extensions of Darcy’s Law (ca. 1856) to unsaturated media, which can be summarized by the widely known Richards equation (Richards, 1931). The theory and basic assumptions are fully described in any soil physics reference book (e.g. Bear, 1978). Whereas in the saturated situation a simple scalar K—or a matrix in multiple dimensions—characterizing the soil saturated conductivity, may explain proportionality between $
Code on server at http://www.iamg.org/CGEditor/index.htm *Corresponding author. Tel.: +33-1-48035689; fax: +33-148473088. E-mail address: [email protected] (E. Perrier). 1 Now at INRA, Laon, France. 2 Now at Ilog, Paris, France.
water flux and the gradient of applied potentials (for example gravitational potentials), the unsaturated situation is much more complicated, since the conductivity becomes a function of the water content KðyÞ and the variation of water content induces new capillary potentials described by the function hðyÞ; which superimpose gradients that also govern fluxes. This leads to the non-linear Richards equation—recalled later in the theory section—which is based on the knowledge of the two soil hydraulic characteristic functions KðyÞ and hðyÞ and which usually cannot be solved analytically. So soil physicists have been developing, then improving numerous numerical schemes to solve the Richards equation, with various conditions adapted to their research subjects. The codes may be found in general books (e.g. Campbell, 1985) as well as in specialized papers (e.g. Haverkamp et al., 1977; Celia et al., 1990; Gottardi and Venutelli, 1993; Ju and Kung, 1997). Most of them are not fully utilized and used only by researchers
0098-3004/02/$ - see front matterr 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 9 8 - 3 0 0 4 ( 0 2 ) 0 0 0 2 7 - 4
1108
E. Perrier et al. / Computers & Geosciences 28 (2002) 1107–1118
already familiar both with the model and with numerical coding. All of them deal with the numerical solution of the partial differential equation that rules the spatiotemporal redistribution of water y or associated pressure head h in the soil profile, and a major difference between these algorithms lies in their accuracy with respect to the mass-balance error generated by numerical approximations. In fact, using these models often requires not only specific understanding, but also additional coding to cope with input and output data, and to adapt them to any particular example (different types of input soil properties and different flow boundary conditions), which reveals the need for a complete user modeling environment. That is why some software packages have recently been developed which enable easier use of the basic flow model (e.g. Verburg et al., 1997; van Dam et al., 1997; Simunek et al., 1999). Such high-level packages are sold as operational tools for scientists working in conditions similar to the researchers who developed the programs, and successive versions increase the number of possible conditions. For example, Simunek et al. (1999) handles in the HYDRUS package only particular analytical expressions (Van Genuchten, 1980), whereas Verburg et al. (1997) includes in the new version of the SWIM package various analytical expressions to model the characteristic hydraulic properties. Simunek et al. (1999) added advanced procedures to cope with two-dimensional flow and to couple with solute transport. Up to now, even if most soils—as soon as they contain a small amount of clay—are subject to natural deformation with varying water content, deformation is still neglected. The new package presented allows simulation of onedimensional vertical flow in cores or in situ, which is the main application of this type of modeling in soil science, agronomy and hydrology. A classical numerical scheme is implemented, which is expected to behave correctly in most situations and will be only briefly referenced (it could be modified by advanced users if necessary). To estimate the soil hydraulic properties, the program includes an inverse process based on the simplex method, already used in Garnier et al. (1997b), the principle of which is extended and described in the theory section. The first contribution of this work is to add procedures to take into account soil deformation when necessary according to Garnier’s modeling algorithms which were described theoretically (Garnier, 1996), successfully used (Garnier et al., 1997a), but never published. The principle is given in the theory section and code sources can be downloaded on the web.3 The second contribution lies in the emphasis put upon didactical and open access to researchers unfamiliar with 3 The Ecoul Software. http://www.bondy.ird.fr/Bperrier/ data/ecoul/logecoul.html; 5pp.
both the modeling environment and the model: a helpful feature consisting of visualizing dynamically the soil moisture redistribution and the soil deformation throughout the simulation. Several years after the conception of a first version in rigid soils (Perrier, 1989), the present up-to-date software architecture was programmed (Leclerc, 1996) for internal use. Only one public demonstration of its applications was done at a workshop (Perrier et al., 1997). It has been developed using the object-oriented language C++ using the graphical and multi-windowing C++ procedures from the library Ilog Views. (The executable program can be downloaded freely by university researchers by sending an e-mail to the corresponding author but, for private companies, this would need run-time charges or purchase of the Ilog Views library.)
2. Theory 2.1. The water flow model The program simulates water flow in swelling soils and rigid soils (which may be considered as a special case of swelling soils). Two algorithms are proposed, one for rigid soils using the classical Richards equation following the resolution scheme of Touma (1987) and one for swelling soils using the modified Richards equation of Garnier et al. (1997a). 2.1.1. Rigid soils The Richards equation is frequently considered to govern water flow in partially saturated soils. In a onedimensional vertical system, this equation is given by qy q qh ¼ KðhÞ 1 ; ð1Þ qt qz qz where y (cm3/cm3) is the volumetric water content (water volume/total volume of soil), t is the time (h), z (cm) is the depth (positive downward from the soil surface), K is the hydraulic conductivity (cm/h), and h is the pressure head (cm). The numerical solution of Eq. (1) is carried out by using a classical finite-difference discretization scheme (e.g. Touma, 1987), where the unknowns are the pressure heads hi at each node and Ci is the capacity equal to the derivative qy=qh of the function yðhÞ: " ! kþ1 kþ1 hiþ1 hkþ1 hki 1 i k hi k Ci ¼ 1 ; K Dt Dz Dz iþ1=2 kþ1 hi hkþ1 k i1 1 ð2Þ Ki1=2 Dz where i refers to the depth increment and k refers to the time increment.
E. Perrier et al. / Computers & Geosciences 28 (2002) 1107–1118
2.1.2. Swelling soils In swelling soils, the water flow equation is usually expressed in a moving particle coordinate system – the Lagrangian coordinates. Another difference with the description of water flow in rigid soils is the inclusion of the overburden pressure due to the soil swelling. Garnier et al. (1997a) proposed recently a coordinate transformation based on a three-dimensional and anisotropic deformation of the soil considering a one-dimensional water flow. The water movement equation uses the variables moisture ratio W (water volume/solid particle volume) and void ratio e (void volume/solid particle volume) because these variables are expressed relative to the solid particle volume, which remains constant during the deformation. The equation that describes the water movement in such system is qW q qW ¼ Ið1 þ eÞ T1 T2 ; ð3Þ qt qZ qZ where Z (cm) is the Lagrangian coordinate. The terms I; T1 and T2 of Eq. (3) are defined, respectively, by 1 þ e0 1=rs I¼ ; ð4aÞ 1þe T1 ¼ KI
qh þ qW
Z 0
Z
qV% ; gI 1 dZ qW
ð4bÞ
% T2 ¼ Kð1 gVÞ;
ð4cÞ
where rs is the geometric factor, g is the apparent wet specific density, and V% is the slope of the deformation curve. The state of reference, denoted by the subscript ‘‘0’’, may be chosen arbitrarily. In the program, e0 (cm3/cm3) is taken as the void ratio at the shrinkage limit. The numerical solution of Eq. (3) may be carried out by using the same discretization method as for rigid soils (finite-difference scheme). Because of the depth-dependent lateral deformation of the soil, however, the discretization has to account explicitly for the crosssectional area of each spatial element so as to preserve mass balance. The discretization form of the water flow equation (Eq. (3)) is ! kþ1 Wiþ1 Wikþ1 Wkþ1 Wki k k k 1 i ¼ Ii ðl þ ei Þ ðT1 Þiþ1=2 dt dZ dZ sk i1=2 ðT2 Þkiþ1=2 k siþ1=2
ðT1 Þki1=2
Wkþ1 Wkþ1 i i1 dZ
!
!! ðT2 Þki1=2 ð5Þ
where i refers to the depth increment and k refers to the time increment.
1109
The cross-sectional area of the spatial element i; Sik ; is calculated by ðVs Þki ; ð6Þ dZ k where ðVs Þi is the solid particle volume within element i:
Sik ¼ Iik ð1 þ eki Þ
2.1.3. Hydraulic characteristics To solve the water flow equations for rigid and swelling soils (Eqs. (1) and (3), respectively), we need two hydraulic properties for rigid soil, the retention curve and the hydraulic conductivity curve and a third hydraulic property for swelling soils, the deformation curve. The retention curve expresses the relation between the water content (volumetric water content y or moisture ratio W) and the pressure head h: The hydraulic conductivity curve expresses the relation between the water content (volumetric water content y or moisture ratio W) and the hydraulic conductivity K (cm/h). The deformation curve expresses the relation between water content (volumetric water content y or moisture ratio W) and void ratio e: These three relations can be modeled by different mathematical expressions. The program gives the user a large choice among different analytical expressions from literature (Gardner, 1958; Brooks and Corey, 1964; Van Genuchten, 1980; Braudeau, 1988; Rieu and Sposito, 1991; Perrier et al., 1996) which can be expressed either as functions of volumetric water content y or moisture ratio W: 2.1.4. Boundary conditions The program allows different types of boundary conditions, flux, hydraulic head and porous plate conditions. 2.2. The inverse method The water flow (Eq. (1) or (3)) can be associated with an inverse process in order to estimate the hydraulic properties from experimental data (Garnier et al., 1997b). Here it is extended to any type of available experimental data and any type of hydraulic properties according to the large set of possibilities given in the package. The principle of the inverse method is to estimate parameters of the analytical expressions chosen to model the hydraulic properties (retention curve and hydraulic conductivity curve) so that the numerical simulation of the experiments with the water flow equations (Eq. (1) or (3)) gives the best representation of available measurements. The unknown parameters are determined by minimizing the least-squares objective function OðbÞ: OðbÞ ¼
N X % i ðbÞg2 ; fwi ½Mi M
ð7Þ
i¼1
where b ¼ ðb1 ; y; bp Þ is a vector of p components that are % i is ith the p parameters, Mi is the ith calculated value, M
1110
E. Perrier et al. / Computers & Geosciences 28 (2002) 1107–1118
measured value, wi is the weighting factor of the ith measurement, and N is the number of measurements. The minimization of OðbÞ is carried out using the simplex method (Nedler and Mead, 1965). The simplex method consists of finding an optimal vector bopt by means of a geometrical interpretation of the minimization process. One works in a space of p dimensions, calculating the value of the objective function OðbÞ of p þ 1 vectors, and oriented the search towards a minimal value of OðbÞ associated with an optimal value of b by successive geometrical transformations in the space. Most optimization methods are analytical and not geometrical which means that they calculate the spatial derivative of the objective function ðqOðbÞ=qbj Þ; where bj is the jth component of the vector b: Because of the mathematical form of Eqs. (1) and (3), analytical optimization methods determine the spatial derivative by complex numerical approximations that are avoided by the simplex method.
3. Program description 3.1. General presentation The software ECOUL has been developed using the object-oriented programming language C++ and the graphical and multi-windowing libraries Ilog Views distributed by the computer engineering firm Ilog. It runs either on PCs requiring the operating system Windows 95/98 or on Sun Workstations requiring the operating system Solaris 2/8. The code has been separated into two distinct parts of equal importance: *
The interface part which creates the different windows (editing spreadsheets or drawing panels), a
*
variety of menus and sub-menus, control buttons and so on. It provides an easy environment to handle the main input/output data associated with the specific modeling field. To modify the sources of this part, one has to use the Ilog Views product and a C++ compiler (Borland C++ has been used so far on PCS and Sun CC on Sun workstations). The modeling part which deals with the specific C++ callbacks coding the numerical solution of the flow equations, the mathematical expressions of the hydraulic properties, the algorithm to implement the inverse method, etc. So it comprises all the numerical calculations linked to the modeling choices. To modify this part, one needs only the C++ compiler.
3.1.1. Program architecture As shown in Fig. 1, the full installation of the software comprises 6 directories: 1. The first directory ‘‘bin’’ includes the executable file ecoul.exe. 2. The ‘‘data’’ directory includes a text file ecoul.dbm, which contains a list of all the French words used in the user interface and their translation in English. Other languages may be added using any word processor. (To switch to another language, one has only to type ecoul-l language in any command term instead of clicking on the executable icon.) The first sub-directory ‘‘Icones’’ gathers all the images used in the software interface as.gif files (about 70 images as in the example CndSurf.gif in Fig. 1 used to select the panel to choose the boundary conditions at the top of the domain). The second sub-directory ‘‘Panneau’’ includes text files in a special Ilog format to define the numerous windows, panels, buttons, menus, used to select parameters or simulation scenarios (about 30.ilv files, like nbretent.ilv in Fig. 1 used for the
Fig. 1. Overview of software architecture.
E. Perrier et al. / Computers & Geosciences 28 (2002) 1107–1118
3.
4.
5.
6.
menu where different analytical expressions for the retention curve are presented to the user). The directory ‘‘include’’ contains all the include file necessary in C++ programming to define classes of objects used either in the interface part or the modeling part. The directory ‘‘src’’ contains all the C++ source codes dealing with the core of the software and using objects of the type defined in the associated include files. It is also divided into two sub-directories ‘‘interface’’ and ‘‘modelization’’ with about 20 files each. For example, the nbretention.cpp files in ‘‘src/ interface’’ (Fig. 1) deal with all the graphics, the type of analytical model, the parameter choices, the associated soil layers if several layers are used with different retention curves, while retention.cpp in ‘‘src/ modelization’’ (Fig. 1) contains procedures to explicitly calculate the water content from pressure head or vice versa, numerical or mathematical derivatives used for each particular analytical expression of retention function: because we use an object-oriented programming style, these particular functions are defined as object classes derived from the general class retention. The directory ‘‘object’’ contains all the classical.o files created during the compilation process, plus .ide files (ecoul.ide in Fig. 1) which contain all the compilation commands and paths to libraries for a straightforward compilation if any changes in the source are made, so-called C++ projects to compile the whole software or only the modeling part. Finally, in the ‘‘work’’ directory, there are text files where all the simulation data and parameters can be recorded, examples (as demo.sim in Fig. 1) and files created by users that can be reopened when needed to carry out new simulations.
3.1.2. An interactive simulation tool Interactivity is first characterized by now classical integrated editing and graphical facilities. All the input data needed to run a simulation are registered into a main data text file using appropriate coding (the extension.sim is associated to such a file, demo.sim is used by default, and any mysimulation.sim can be created). Several specific spreadsheets have been prepared to avoid direct but uncomfortable editing tasks. Qualitative choices are made using menus (type of hydraulic soil properties, type of flow model, type of boundary conditions, etc.) while quantitative values are chosen using either the spreadsheets cells, or mouse actions on intensity sliders. Data can be also imported from already registered files. The soil is defined by one or several layers, hydraulic properties associated with each layer, and initial moisture or pressure head profile. All the data can be checked by visualizing associated graphics, e.g the
1111
curves representing the hydraulic properties (Fig. 2 or 6a and b). These various simulation scenarios associated with different types of boundary conditions can be chosen either according to a given experimental protocol (e.g. Fig. 3) or by ‘‘manually’’ clicking on specific choice menus and intensity gliders. These simulation conditions are also visualized at the top and bottoms of the soil column by schematic images (e.g. Figs. 4a and 5). The simulated soil is visualized as a column the height of which is divided into as many sub-parts as defined by the spatial discretization step chosen for the flow model numerical solution. At the beginning of the simulation, the initial distribution of the moisture content is calculated in each sub-part by interpolating available data. Colors are associated with a water content graduating in scale: it goes from yellow for low water content to brown for high water content (Figs. 4a, 5 and 6c). At each time step, colors are refreshed according to the calculated evolution of the moisture profile. When soil deformation occurs (e.g. Fig. 5), isotropic or anisotropic variations of volume are also visible on the simulated soil column. The size of each layer changes with the deformation. Mean values in the profile may be simultaneously displayed in another window where dynamical gauges indicate the evolution of the different terms of the water mass balance (e.g. Fig. 4b). The end-user can also decide which type of output data he/she wants to visualize or to register by choosing in available list menus (e.g. evolution of the moisture content in time and space, or dynamic evolution of the mean water content, runoff series or cumulative infiltration at the top of the soil) or associated suggestive buttons (right hand side of Fig. 4a). Output data can be visualized not only at the end of the simulation run, but they can also be displayed simultaneously. When no experimental data are available, the simulated temporal series are given for a set of successive times calculated from the total duration and the simulated water profile at three top/middle/bottom depths (Fig. 4c and d). Where observed data are available, simulated data are registered and displayed at the same location and time (Fig. 6). Comparison of simulated and observed data can be made numerically and graphically within the same software. These data can also be saved and then exported if needed to other statistical or graphical software packages. As far as the inverse method is concerned, one selects which parameters of the hydraulic properties have to be retro-estimated by a simple ‘‘yes’’ plus a modification of the default range of investigation if needed (Figs. 6a and b). The successive runs corresponding to successive estimates of new parameters (Figs. 6a0 , b0 , g and h) are clearly shown and the increasing quality of the match between observed and simulated flow data can also be visualized (e.g. Figs. 6d and f for initial and final
1112
E. Perrier et al. / Computers & Geosciences 28 (2002) 1107–1118
Fig. 2. Panels to select three soil hydraulic characteristics: (a) retention curve; (b) hydraulic conductivity curve; and (c) deformation curve, for deforming soils only.
comparisons at the beginning and at the end of the optimization process); there is also a specific control window showing the initial parameters, the values corresponding to the on-going run and the evolution of the optimization criterion (Fig. 6e). Interactivity goes beyond these features by introducing a new idea: graphics are used not only to visualize the input and output data of the flow model, but also to give as much information as possible about the model itself. One may define a type of soil as a column with one
or multiple layers. One clicks on the computer mouse to pour rain on it, using the glider to modify rain intensity without stopping the simulation, then making sunshine and choosing an evaporation intensity. The user can observe the spatio-temporal evolution of the water profile or any other indicator, stop the simulation at any time and change any boundary condition to understand better their effect on the results. The addition of such visual components may first appear as a game. Without neglecting the importance of
E. Perrier et al. / Computers & Geosciences 28 (2002) 1107–1118
1113
Fig. 3. Panels to select boundary conditions. Here, according to given table, first two top conditions are specified for simulation shown in Figs. 4 and 5.
educational games as far as students are concerned, this way to experiment with the basic concepts of a model may be useful also to many users. For advanced users, it is of course possible to de-activate for a while or systematically dynamical displays to spare simulation time. Anyway, this concept of interactivity is more than a computer game. It relies on a state of mind concerning mathematical modeling and computer simulations: computer software packages give nowadays the opportunity to run complex mathematical models (it could be any elaborated numerical schemes or statistical packages) without being a specialist in the underlying theory. Adding interactivity of the first category (graphical or editing facilities, plus on-line help associated with each interface item) may spare reading the practical part of the user’s manual. Adding interactivity of the second category (graphical features to make the modeling process as clear as possible) may not only spare some a priori reading of the theoretical part, it may be a way to learn about it and later to avoid misuse. 3.2. Test cases 3.2.1. Infiltration–evaporation cycles in a swelling soil This example registered as fissure.sim deals with a swelling soil surface (3 cm depth) submitted to infiltration–evaporation cycles. The swelling soil is a mixture of
loam and bentonite. Its hydraulic properties are described in Garnier et al. (1997a). The algorithms of transfer in swelling soil and direct method were selected. In this example, the moisture variable is the moisture ratio and there is one layer, characterized by 3 curves: the Van Genuchten (1980), Brooks and Corey (1964), and Braudeau (1988) models were selected for the retention, hydraulic conductivity and deformation curves, respectively. Parameters of these three curves are introduced on the panels (Fig. 2). The geometry factor Rs for the deformation curve is 10, which corresponds to a much larger variation of volume in the horizontal direction than in the vertical one. The swelling soil is subjected to different cycles of evaporation–infiltration depicted in Fig. 3 (3 evaporations and 2 rain events). Zero flux was imposed at the bottom. The main screen view (Fig. 4a) shows the distribution of water at a given time whereas Fig. 4b shows the cumulative value of some mean, global indicators. They are the evaporated water depth (cm), the soil water depth and the mass balance calculated from imposed fluxes and calculated water content. The mass error of 5% is acceptable. Fig. 4c shows the temporal changes of water content at 1, 2 and 3 cm depths while Fig. 4d gives the water content profiles at some intermediate times. We can see
1114
E. Perrier et al. / Computers & Geosciences 28 (2002) 1107–1118
Fig. 4. Simulation of infiltration–evaporation cycles in swelling soil Simulation of example fissure.sim (see test case A in text). Hydraulic properties are described in Fig. 2. Top boundary condition is described in Fig. 3. Bottom boundary condition is zero flux. (a) Main panel of software. Icons in main panel can be used to open other panels. (b) Panel to exhibit evolution of components (mean evaporation, infiltration, drainage, etc.) of mass balance. (c) Panel exhibiting at any time evolution of soil water content at three different depths. (d) Panel exhibiting at any depth evolution of soil water content at ten different times.
E. Perrier et al. / Computers & Geosciences 28 (2002) 1107–1118
1115
Fig. 5. Successive screen copies of main panel throughout simulation shown in Fig. 4.
the rapid changes of water content at 1 cm compared to the other depths. Fig. 5 shows the state of the soil at successive times throughout a total duration of about 10 days. When the soil is submitted to evaporation, the soil surface loses water with associated shrinkage so its shape becomes pyramidal with less water at the top (yellow) than at the bottom. This corresponds to the formation of cracks. When water infiltrates under a rain event, the soil swelling leads to crack closure at the surface. These processes simulated on a synthetic soil depict and explain rather well observed phenomena in real soils.
3.2.2. Parameter estimation by inverse method in a rigid soil during an infiltration experiment In a first step, before using the inverse method, an infiltration experiment was simulated with ECOUL using the direct method. Hydraulic properties of the Ariana silty clay loam soil obtained by Rieu and Sposito (1991) were introduced in the model. We used Rieu and Sposito’s model (1991) for the retention curve and Brooks and Corey’s model (1964) for the hydraulic conductivity curve. The soil (1 m depth) is submitted to an infiltration at its surface with a hydraulic head of 10 cm for 1 h. Initial conditions are those of a soil profile
1116
E. Perrier et al. / Computers & Geosciences 28 (2002) 1107–1118
E. Perrier et al. / Computers & Geosciences 28 (2002) 1107–1118
in equilibrium with an aquifer at the bottom. Simulated temporal evolution of water content at two depths and cumulative infiltrated water were registered. These data are called reference data. In a second step, we consider the previous soil as a real soil, submitted to the same conditions as in the direct simulation (as shown in Fig. 6c), and the reference data as ‘‘experimental’’ data that could have been easily measured on this soil using tensiometers at the two given depths and monitoring entering water. Then we test the inverse method in order to re-calculate some hydraulic parameters used in the direct simulation but here ‘‘forgotten’’. In this example, the file is registered as ariana.sim. When introducing hydraulic parameters, we consider that some parameters of the retention curve like residual and saturated water content and the minimal potential (Fig. 6a) are easily known and we will estimate by inverse method only the fractal parameter D linked to the slope of the curve. The forgotten value of D is 2.90 and we will try to investigate all possible values between 2.70 and 2.99 (Fig. 6a). As regards the hydraulic conductivity curve, we consider that the residual and saturated water contents are known as well as the conductivity value at saturation and that the parameter B is more difficult to measure. So the optimization process will deal with a search for B between 1 and 100, where the forgotten value was equal to 29 (Fig. 6b). The initial values for the search for D and B are set, respectively, to 2.85 and 10 (Figs. 6a and b) and correspond to the hydraulic characteristic curves displayed in Figs. 6a0 and b0 . The user can visualize the evolution of the objective function (Fig. 6e) which corresponds to the distance between here experimental data (red points) and simulated data (yellow points), presented in Fig. 6d at the beginning of the simulation and in Fig. 6f at the end. The weights given to each point of the surface infiltration curve are 1 (the default value) and 100 to each point of the sub-surface water content curves, to account for the different type of units and associated numerical values. The objective function decreased as expected when the number of iterations increased. Initial, current and optimal parameters are also shown in Fig. 6e after 59 iterations. The optimization process reached here 208 iterations and the final estimated values for D and B were, respectively, 2.899876 and
1117
29.010525. The estimated hydraulic properties at the end of the inverse method are shown in Figs. 6g and h which are to be compared to Figs. 6 a0 and b0 . This case study is given here as a test of the optimization algorithm. The inverse method is here totally satisfactory, of course because the experimental data were actually associated with a ‘‘perfect modeled soil’’, but in real situations, the exact solution is unknown, the number of estimated parameters as well as the number of iterations may be higher, it is even possible that the optimization process does not converge towards a unique solution.
4. Conclusion This paper aims first at informing a range of potential users about the existence of the software package ECOUL. ECOUL gives easy access to the modeling of deformation processes in swelling soils or other deforming porous media, which are widespread, but little studied as far as the modeling of their hydraulic behavior is concerned. We give an example of the way numerical modeling in a specific field could become available to a larger audience, or used for educational purposes by providing the means to implement interactivity and visualization concepts. The ECOUL software package is a new type of simulation tool that could later be enhanced by additional functionalities.
Acknowledgements The authors are grateful to Mrs. Chantal Bernard who designed all the iconography of the package, and to Mr. Nicolas Perrier who recently corrected some computer memory bugs.
References Bear, J., 1978. Hydraulics of Groundwater. McGraw-Hill International Book Co., New York, London, 569pp. Braudeau, E., 1988. Equation g!en!eralis!ee des courbes de retrait d’!echantillons de sol non structur!es (in French). Comptes Rendus de l’Acad!emie des Sciences S!erie 2 307, 1731–1734.
Fig. 6. Parameter estimation by inverse method in rigid soil during infiltration experiment. Simulation of example ariana.sim (see test case B in text). (a) and (a0 ) Initial parameters for water retention curve. (b) and (b0 ) Initial parameters for hydraulic conductivity curve. (c) Initial soil water profile as shown in main simulation panel. (d) Comparison between observed and simulated data (water content evolution at two different depth and amount of infiltrated water at top of soil) using (a) and (b) initial parameters for hydraulic properties. (e) Optimization process and evolution of difference between observed and simulated data using different sets of parameters ðD; BÞ for hydraulic properties. (f) Comparison between observed and simulated data using final, optimized parameters for hydraulic properties. (g) Optimized retention curve. (h) Optimized hydraulic conductivity curve.
1118
E. Perrier et al. / Computers & Geosciences 28 (2002) 1107–1118
Brooks, R.H., Corey, A.T., 1964. Hydraulic properties of porous media in hydrology Paper 3, Colorado State Univeristy, Fort Collins, pp. I–V, 1–27. Campbell, G.S., 1985. Soil Physics with BASIC—Transport Models for Soil–Plant Systems. Developments in Soil Science, Vol. 14, Elsevier, New York, 150pp. Celia, M.A., Bouloutas, E.T., Zarba, R.L., 1990. A general mass-conservative numerical solution for the unsaturated flow equation. Water Resources Research 26 (7), 1483– 1496. Gardner, W.R., 1958. Some steady-state solutions of the unsaturated flow equation with application to evaporation from a water table. Soil Science 85, 228–232. Garnier, P., 1996. D!etermination des caract!eristiques hydrodynamiques de sols d!eformables par la m!ethode inverse. Ph.D. Dissertation, Nancy University, France, 144pp. (in French). Garnier, P., Perrier, E., Angulo-Jaramillo, R., Baveye, P., 1997a. Numerical model of 3-dimensional anisotropic deformation and 1-dimensional water flow in swelling soils. Soil Science 162 (6), 410–420. Garnier, P., Rieu, M., Boivin, P., Vauclin, M., Baveye, P., 1997b. Determining the hydraulic properties of a swelling soil from a transient evaporation experiment. Soil Science Society American Journal 61 (6), 1555–1563. Gottardi, G., Venutelli, M., 1993. Richards: computer program for the numerical simulation of one dimensional infiltration into unsaturated soil. Computers & Geosciences 19 (9), 1239–1266. Haverkamp, R., Vauclin, M., Touma, J., Wierenga, P.J., Vachaud, G., 1977. A comparison of numerical simulation models for one dimensional infiltration. Soil Science Society American Journal 41 (2), 285–294. Ju, S.H., Kung, K.J.S., 1997. Mass types, element orders and solution schemes for the Richards equation. Computers & Geosciences 23 (2), 175–187. Leclerc, C., 1996. Le logiciel ECOUL, M!emoire de DESS Informatique. Universit!e Henri Poincar!e, Nancy, France, 50pp. (in French). Nedler, J.A., Mead, R., 1965. A simplex procedure for function minimization. The Computer Journal 7, 308–313. Perrier, E., 1989. ECOULv1.0: simulation num!erique et graphique des transferts hydriques dans un sol non satur!e. ORSTOM report FA41831/2, Institut de Recherches pour le D!eveloppement Documentation Services, France, 54pp. (in French).
Perrier, E., Garnier, P., Leclerc, C., 1997. Profile-scale simulation of water flow: a software package to visualize and estimate soil hydraulic properties effects, poster and demonstration, In: van Genuchten, M.Th., Leij, F.J., Wu, L. (Eds.), Proceedings of the International Workshop on Characterization and Measurement of the Hydraulic Properties of Unsaturated Porous Media, Riverside, California, October 22–24, 1997, University of California, pp. 349–354. Perrier, E., Rieu, M., Sposito, G., de Marsily, G., 1996. Models of the Water Retention Curve for soils with a fractal poresize distribution. Water Resources Research 32 (10), 3025– 3031. Richards, L.A., 1931. Capillary conduction of liquids through porous medium. Physics 1, 318–333. Rieu, M., Sposito, G., 1991. Fractal fragmentation, soil porosity, and soil water properties: I theory, II applications. Soil Science Society American Journal 55, 1231–1244. Simunek, J., van Genuchten, M.Th., Sejna, M., 1999. Using the HYDRUS-1D and HYDRUS-2D codes for estimating unsaturated soil hydraulic and solute. In: van Genuchten, M.Th., Leij, F.J., Wu, L. (Eds.), Proceedings of the International Workshop on Characterization and Measurement of the Hydraulic Properties of Unsaturated Porous Media, Riverside, California, October 22–24, 1997, University of California, pp. 1523–1536. Touma, J. 1987. TEST: mod"ele pour tester la representativit!e des caract!eristiques hydrodynamiques d’un sol non satur!e d!etermin!es in situ. ORSTOM report FA24686/2, Institut de Recherches pour le D!eveloppement Documentation Services, France, 26pp. (in French). Van Dam, J.C., Hugen J., Wesseling J.G., Feddes R.A., Kabat P., van Walsum, P.E.V., Groenendijk, P., van Diepen, C.A., 1997. Theory of Swap version 2.0. Simulation of water flow, solute transport and plant growth in the soil–water–atmosphere–plant environment. Report 71, Water Resources, Wageningen University, Technical Document 45, Alterra Green World Research, Wageningen, 167pp. Van Genuchten, M.T., 1980. A closed form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society American Journal 44, 892–898. Verburg, K., Ross, P.J., Bristow, K.L., 1997. SWIMv2.1. User manual. Division Report, CSIRO Division of Soil, CSIRO Information Services, East Melbourne, Australia, No. 130, 107pp.
ECOUL: an interactive computer tool to study hydraulic behavior of swelling and rigid soils$ Edith Perrier*, Patricia Garnier1, Christian Leclerc2 Unite de Recherches Geodes, Centre IRD Ile de France, 32 avenue Henri Varagnat, 93143 Bondy Cedex, France Received 15 June 2001; received in revised form 19 March 2002; accepted 25 March 2002
Abstract ECOUL is an interactive, didactic software package which simulates vertical water flow in unsaturated soils. Endusers are given an easily-used tool to predict the evolution of the soil water profile, with a large range of possible boundary conditions, through a classical numerical solution scheme for the Richards equation. Soils must be characterized by water retention curves and hydraulic conductivity curves, the form of which can be chosen among different analytical expressions from the literature. When the parameters are unknown, an inverse method is provided to estimate them from available experimental flow data. A significant original feature of the software is to include recent algorithms extending the water flow model to deal with deforming porous media: widespread swelling soils, the volume of which varies as a function of water content, must be described by a third hydraulic characteristic property, the deformation curve. Again, estimation of the parameters by means of inverse procedures and visualization facilities enable exploration, understanding and then prediction of soil hydraulic behavior under various experimental conditions. r 2002 Elsevier Science Ltd. All rights reserved. Keywords: Water flow modeling; Richards equation; Deforming soils; Unsaturated hydraulic properties; Inverse method
1. Introduction Water flow in soils has been described for nearly a century by extensions of Darcy’s Law (ca. 1856) to unsaturated media, which can be summarized by the widely known Richards equation (Richards, 1931). The theory and basic assumptions are fully described in any soil physics reference book (e.g. Bear, 1978). Whereas in the saturated situation a simple scalar K—or a matrix in multiple dimensions—characterizing the soil saturated conductivity, may explain proportionality between $
Code on server at http://www.iamg.org/CGEditor/index.htm *Corresponding author. Tel.: +33-1-48035689; fax: +33-148473088. E-mail address: [email protected] (E. Perrier). 1 Now at INRA, Laon, France. 2 Now at Ilog, Paris, France.
water flux and the gradient of applied potentials (for example gravitational potentials), the unsaturated situation is much more complicated, since the conductivity becomes a function of the water content KðyÞ and the variation of water content induces new capillary potentials described by the function hðyÞ; which superimpose gradients that also govern fluxes. This leads to the non-linear Richards equation—recalled later in the theory section—which is based on the knowledge of the two soil hydraulic characteristic functions KðyÞ and hðyÞ and which usually cannot be solved analytically. So soil physicists have been developing, then improving numerous numerical schemes to solve the Richards equation, with various conditions adapted to their research subjects. The codes may be found in general books (e.g. Campbell, 1985) as well as in specialized papers (e.g. Haverkamp et al., 1977; Celia et al., 1990; Gottardi and Venutelli, 1993; Ju and Kung, 1997). Most of them are not fully utilized and used only by researchers
0098-3004/02/$ - see front matterr 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 9 8 - 3 0 0 4 ( 0 2 ) 0 0 0 2 7 - 4
1108
E. Perrier et al. / Computers & Geosciences 28 (2002) 1107–1118
already familiar both with the model and with numerical coding. All of them deal with the numerical solution of the partial differential equation that rules the spatiotemporal redistribution of water y or associated pressure head h in the soil profile, and a major difference between these algorithms lies in their accuracy with respect to the mass-balance error generated by numerical approximations. In fact, using these models often requires not only specific understanding, but also additional coding to cope with input and output data, and to adapt them to any particular example (different types of input soil properties and different flow boundary conditions), which reveals the need for a complete user modeling environment. That is why some software packages have recently been developed which enable easier use of the basic flow model (e.g. Verburg et al., 1997; van Dam et al., 1997; Simunek et al., 1999). Such high-level packages are sold as operational tools for scientists working in conditions similar to the researchers who developed the programs, and successive versions increase the number of possible conditions. For example, Simunek et al. (1999) handles in the HYDRUS package only particular analytical expressions (Van Genuchten, 1980), whereas Verburg et al. (1997) includes in the new version of the SWIM package various analytical expressions to model the characteristic hydraulic properties. Simunek et al. (1999) added advanced procedures to cope with two-dimensional flow and to couple with solute transport. Up to now, even if most soils—as soon as they contain a small amount of clay—are subject to natural deformation with varying water content, deformation is still neglected. The new package presented allows simulation of onedimensional vertical flow in cores or in situ, which is the main application of this type of modeling in soil science, agronomy and hydrology. A classical numerical scheme is implemented, which is expected to behave correctly in most situations and will be only briefly referenced (it could be modified by advanced users if necessary). To estimate the soil hydraulic properties, the program includes an inverse process based on the simplex method, already used in Garnier et al. (1997b), the principle of which is extended and described in the theory section. The first contribution of this work is to add procedures to take into account soil deformation when necessary according to Garnier’s modeling algorithms which were described theoretically (Garnier, 1996), successfully used (Garnier et al., 1997a), but never published. The principle is given in the theory section and code sources can be downloaded on the web.3 The second contribution lies in the emphasis put upon didactical and open access to researchers unfamiliar with 3 The Ecoul Software. http://www.bondy.ird.fr/Bperrier/ data/ecoul/logecoul.html; 5pp.
both the modeling environment and the model: a helpful feature consisting of visualizing dynamically the soil moisture redistribution and the soil deformation throughout the simulation. Several years after the conception of a first version in rigid soils (Perrier, 1989), the present up-to-date software architecture was programmed (Leclerc, 1996) for internal use. Only one public demonstration of its applications was done at a workshop (Perrier et al., 1997). It has been developed using the object-oriented language C++ using the graphical and multi-windowing C++ procedures from the library Ilog Views. (The executable program can be downloaded freely by university researchers by sending an e-mail to the corresponding author but, for private companies, this would need run-time charges or purchase of the Ilog Views library.)
2. Theory 2.1. The water flow model The program simulates water flow in swelling soils and rigid soils (which may be considered as a special case of swelling soils). Two algorithms are proposed, one for rigid soils using the classical Richards equation following the resolution scheme of Touma (1987) and one for swelling soils using the modified Richards equation of Garnier et al. (1997a). 2.1.1. Rigid soils The Richards equation is frequently considered to govern water flow in partially saturated soils. In a onedimensional vertical system, this equation is given by qy q qh ¼ KðhÞ 1 ; ð1Þ qt qz qz where y (cm3/cm3) is the volumetric water content (water volume/total volume of soil), t is the time (h), z (cm) is the depth (positive downward from the soil surface), K is the hydraulic conductivity (cm/h), and h is the pressure head (cm). The numerical solution of Eq. (1) is carried out by using a classical finite-difference discretization scheme (e.g. Touma, 1987), where the unknowns are the pressure heads hi at each node and Ci is the capacity equal to the derivative qy=qh of the function yðhÞ: " ! kþ1 kþ1 hiþ1 hkþ1 hki 1 i k hi k Ci ¼ 1 ; K Dt Dz Dz iþ1=2 kþ1 hi hkþ1 k i1 1 ð2Þ Ki1=2 Dz where i refers to the depth increment and k refers to the time increment.
E. Perrier et al. / Computers & Geosciences 28 (2002) 1107–1118
2.1.2. Swelling soils In swelling soils, the water flow equation is usually expressed in a moving particle coordinate system – the Lagrangian coordinates. Another difference with the description of water flow in rigid soils is the inclusion of the overburden pressure due to the soil swelling. Garnier et al. (1997a) proposed recently a coordinate transformation based on a three-dimensional and anisotropic deformation of the soil considering a one-dimensional water flow. The water movement equation uses the variables moisture ratio W (water volume/solid particle volume) and void ratio e (void volume/solid particle volume) because these variables are expressed relative to the solid particle volume, which remains constant during the deformation. The equation that describes the water movement in such system is qW q qW ¼ Ið1 þ eÞ T1 T2 ; ð3Þ qt qZ qZ where Z (cm) is the Lagrangian coordinate. The terms I; T1 and T2 of Eq. (3) are defined, respectively, by 1 þ e0 1=rs I¼ ; ð4aÞ 1þe T1 ¼ KI
qh þ qW
Z 0
Z
qV% ; gI 1 dZ qW
ð4bÞ
% T2 ¼ Kð1 gVÞ;
ð4cÞ
where rs is the geometric factor, g is the apparent wet specific density, and V% is the slope of the deformation curve. The state of reference, denoted by the subscript ‘‘0’’, may be chosen arbitrarily. In the program, e0 (cm3/cm3) is taken as the void ratio at the shrinkage limit. The numerical solution of Eq. (3) may be carried out by using the same discretization method as for rigid soils (finite-difference scheme). Because of the depth-dependent lateral deformation of the soil, however, the discretization has to account explicitly for the crosssectional area of each spatial element so as to preserve mass balance. The discretization form of the water flow equation (Eq. (3)) is ! kþ1 Wiþ1 Wikþ1 Wkþ1 Wki k k k 1 i ¼ Ii ðl þ ei Þ ðT1 Þiþ1=2 dt dZ dZ sk i1=2 ðT2 Þkiþ1=2 k siþ1=2
ðT1 Þki1=2
Wkþ1 Wkþ1 i i1 dZ
!
!! ðT2 Þki1=2 ð5Þ
where i refers to the depth increment and k refers to the time increment.
1109
The cross-sectional area of the spatial element i; Sik ; is calculated by ðVs Þki ; ð6Þ dZ k where ðVs Þi is the solid particle volume within element i:
Sik ¼ Iik ð1 þ eki Þ
2.1.3. Hydraulic characteristics To solve the water flow equations for rigid and swelling soils (Eqs. (1) and (3), respectively), we need two hydraulic properties for rigid soil, the retention curve and the hydraulic conductivity curve and a third hydraulic property for swelling soils, the deformation curve. The retention curve expresses the relation between the water content (volumetric water content y or moisture ratio W) and the pressure head h: The hydraulic conductivity curve expresses the relation between the water content (volumetric water content y or moisture ratio W) and the hydraulic conductivity K (cm/h). The deformation curve expresses the relation between water content (volumetric water content y or moisture ratio W) and void ratio e: These three relations can be modeled by different mathematical expressions. The program gives the user a large choice among different analytical expressions from literature (Gardner, 1958; Brooks and Corey, 1964; Van Genuchten, 1980; Braudeau, 1988; Rieu and Sposito, 1991; Perrier et al., 1996) which can be expressed either as functions of volumetric water content y or moisture ratio W: 2.1.4. Boundary conditions The program allows different types of boundary conditions, flux, hydraulic head and porous plate conditions. 2.2. The inverse method The water flow (Eq. (1) or (3)) can be associated with an inverse process in order to estimate the hydraulic properties from experimental data (Garnier et al., 1997b). Here it is extended to any type of available experimental data and any type of hydraulic properties according to the large set of possibilities given in the package. The principle of the inverse method is to estimate parameters of the analytical expressions chosen to model the hydraulic properties (retention curve and hydraulic conductivity curve) so that the numerical simulation of the experiments with the water flow equations (Eq. (1) or (3)) gives the best representation of available measurements. The unknown parameters are determined by minimizing the least-squares objective function OðbÞ: OðbÞ ¼
N X % i ðbÞg2 ; fwi ½Mi M
ð7Þ
i¼1
where b ¼ ðb1 ; y; bp Þ is a vector of p components that are % i is ith the p parameters, Mi is the ith calculated value, M
1110
E. Perrier et al. / Computers & Geosciences 28 (2002) 1107–1118
measured value, wi is the weighting factor of the ith measurement, and N is the number of measurements. The minimization of OðbÞ is carried out using the simplex method (Nedler and Mead, 1965). The simplex method consists of finding an optimal vector bopt by means of a geometrical interpretation of the minimization process. One works in a space of p dimensions, calculating the value of the objective function OðbÞ of p þ 1 vectors, and oriented the search towards a minimal value of OðbÞ associated with an optimal value of b by successive geometrical transformations in the space. Most optimization methods are analytical and not geometrical which means that they calculate the spatial derivative of the objective function ðqOðbÞ=qbj Þ; where bj is the jth component of the vector b: Because of the mathematical form of Eqs. (1) and (3), analytical optimization methods determine the spatial derivative by complex numerical approximations that are avoided by the simplex method.
3. Program description 3.1. General presentation The software ECOUL has been developed using the object-oriented programming language C++ and the graphical and multi-windowing libraries Ilog Views distributed by the computer engineering firm Ilog. It runs either on PCs requiring the operating system Windows 95/98 or on Sun Workstations requiring the operating system Solaris 2/8. The code has been separated into two distinct parts of equal importance: *
The interface part which creates the different windows (editing spreadsheets or drawing panels), a
*
variety of menus and sub-menus, control buttons and so on. It provides an easy environment to handle the main input/output data associated with the specific modeling field. To modify the sources of this part, one has to use the Ilog Views product and a C++ compiler (Borland C++ has been used so far on PCS and Sun CC on Sun workstations). The modeling part which deals with the specific C++ callbacks coding the numerical solution of the flow equations, the mathematical expressions of the hydraulic properties, the algorithm to implement the inverse method, etc. So it comprises all the numerical calculations linked to the modeling choices. To modify this part, one needs only the C++ compiler.
3.1.1. Program architecture As shown in Fig. 1, the full installation of the software comprises 6 directories: 1. The first directory ‘‘bin’’ includes the executable file ecoul.exe. 2. The ‘‘data’’ directory includes a text file ecoul.dbm, which contains a list of all the French words used in the user interface and their translation in English. Other languages may be added using any word processor. (To switch to another language, one has only to type ecoul-l language in any command term instead of clicking on the executable icon.) The first sub-directory ‘‘Icones’’ gathers all the images used in the software interface as.gif files (about 70 images as in the example CndSurf.gif in Fig. 1 used to select the panel to choose the boundary conditions at the top of the domain). The second sub-directory ‘‘Panneau’’ includes text files in a special Ilog format to define the numerous windows, panels, buttons, menus, used to select parameters or simulation scenarios (about 30.ilv files, like nbretent.ilv in Fig. 1 used for the
Fig. 1. Overview of software architecture.
E. Perrier et al. / Computers & Geosciences 28 (2002) 1107–1118
3.
4.
5.
6.
menu where different analytical expressions for the retention curve are presented to the user). The directory ‘‘include’’ contains all the include file necessary in C++ programming to define classes of objects used either in the interface part or the modeling part. The directory ‘‘src’’ contains all the C++ source codes dealing with the core of the software and using objects of the type defined in the associated include files. It is also divided into two sub-directories ‘‘interface’’ and ‘‘modelization’’ with about 20 files each. For example, the nbretention.cpp files in ‘‘src/ interface’’ (Fig. 1) deal with all the graphics, the type of analytical model, the parameter choices, the associated soil layers if several layers are used with different retention curves, while retention.cpp in ‘‘src/ modelization’’ (Fig. 1) contains procedures to explicitly calculate the water content from pressure head or vice versa, numerical or mathematical derivatives used for each particular analytical expression of retention function: because we use an object-oriented programming style, these particular functions are defined as object classes derived from the general class retention. The directory ‘‘object’’ contains all the classical.o files created during the compilation process, plus .ide files (ecoul.ide in Fig. 1) which contain all the compilation commands and paths to libraries for a straightforward compilation if any changes in the source are made, so-called C++ projects to compile the whole software or only the modeling part. Finally, in the ‘‘work’’ directory, there are text files where all the simulation data and parameters can be recorded, examples (as demo.sim in Fig. 1) and files created by users that can be reopened when needed to carry out new simulations.
3.1.2. An interactive simulation tool Interactivity is first characterized by now classical integrated editing and graphical facilities. All the input data needed to run a simulation are registered into a main data text file using appropriate coding (the extension.sim is associated to such a file, demo.sim is used by default, and any mysimulation.sim can be created). Several specific spreadsheets have been prepared to avoid direct but uncomfortable editing tasks. Qualitative choices are made using menus (type of hydraulic soil properties, type of flow model, type of boundary conditions, etc.) while quantitative values are chosen using either the spreadsheets cells, or mouse actions on intensity sliders. Data can be also imported from already registered files. The soil is defined by one or several layers, hydraulic properties associated with each layer, and initial moisture or pressure head profile. All the data can be checked by visualizing associated graphics, e.g the
1111
curves representing the hydraulic properties (Fig. 2 or 6a and b). These various simulation scenarios associated with different types of boundary conditions can be chosen either according to a given experimental protocol (e.g. Fig. 3) or by ‘‘manually’’ clicking on specific choice menus and intensity gliders. These simulation conditions are also visualized at the top and bottoms of the soil column by schematic images (e.g. Figs. 4a and 5). The simulated soil is visualized as a column the height of which is divided into as many sub-parts as defined by the spatial discretization step chosen for the flow model numerical solution. At the beginning of the simulation, the initial distribution of the moisture content is calculated in each sub-part by interpolating available data. Colors are associated with a water content graduating in scale: it goes from yellow for low water content to brown for high water content (Figs. 4a, 5 and 6c). At each time step, colors are refreshed according to the calculated evolution of the moisture profile. When soil deformation occurs (e.g. Fig. 5), isotropic or anisotropic variations of volume are also visible on the simulated soil column. The size of each layer changes with the deformation. Mean values in the profile may be simultaneously displayed in another window where dynamical gauges indicate the evolution of the different terms of the water mass balance (e.g. Fig. 4b). The end-user can also decide which type of output data he/she wants to visualize or to register by choosing in available list menus (e.g. evolution of the moisture content in time and space, or dynamic evolution of the mean water content, runoff series or cumulative infiltration at the top of the soil) or associated suggestive buttons (right hand side of Fig. 4a). Output data can be visualized not only at the end of the simulation run, but they can also be displayed simultaneously. When no experimental data are available, the simulated temporal series are given for a set of successive times calculated from the total duration and the simulated water profile at three top/middle/bottom depths (Fig. 4c and d). Where observed data are available, simulated data are registered and displayed at the same location and time (Fig. 6). Comparison of simulated and observed data can be made numerically and graphically within the same software. These data can also be saved and then exported if needed to other statistical or graphical software packages. As far as the inverse method is concerned, one selects which parameters of the hydraulic properties have to be retro-estimated by a simple ‘‘yes’’ plus a modification of the default range of investigation if needed (Figs. 6a and b). The successive runs corresponding to successive estimates of new parameters (Figs. 6a0 , b0 , g and h) are clearly shown and the increasing quality of the match between observed and simulated flow data can also be visualized (e.g. Figs. 6d and f for initial and final
1112
E. Perrier et al. / Computers & Geosciences 28 (2002) 1107–1118
Fig. 2. Panels to select three soil hydraulic characteristics: (a) retention curve; (b) hydraulic conductivity curve; and (c) deformation curve, for deforming soils only.
comparisons at the beginning and at the end of the optimization process); there is also a specific control window showing the initial parameters, the values corresponding to the on-going run and the evolution of the optimization criterion (Fig. 6e). Interactivity goes beyond these features by introducing a new idea: graphics are used not only to visualize the input and output data of the flow model, but also to give as much information as possible about the model itself. One may define a type of soil as a column with one
or multiple layers. One clicks on the computer mouse to pour rain on it, using the glider to modify rain intensity without stopping the simulation, then making sunshine and choosing an evaporation intensity. The user can observe the spatio-temporal evolution of the water profile or any other indicator, stop the simulation at any time and change any boundary condition to understand better their effect on the results. The addition of such visual components may first appear as a game. Without neglecting the importance of
E. Perrier et al. / Computers & Geosciences 28 (2002) 1107–1118
1113
Fig. 3. Panels to select boundary conditions. Here, according to given table, first two top conditions are specified for simulation shown in Figs. 4 and 5.
educational games as far as students are concerned, this way to experiment with the basic concepts of a model may be useful also to many users. For advanced users, it is of course possible to de-activate for a while or systematically dynamical displays to spare simulation time. Anyway, this concept of interactivity is more than a computer game. It relies on a state of mind concerning mathematical modeling and computer simulations: computer software packages give nowadays the opportunity to run complex mathematical models (it could be any elaborated numerical schemes or statistical packages) without being a specialist in the underlying theory. Adding interactivity of the first category (graphical or editing facilities, plus on-line help associated with each interface item) may spare reading the practical part of the user’s manual. Adding interactivity of the second category (graphical features to make the modeling process as clear as possible) may not only spare some a priori reading of the theoretical part, it may be a way to learn about it and later to avoid misuse. 3.2. Test cases 3.2.1. Infiltration–evaporation cycles in a swelling soil This example registered as fissure.sim deals with a swelling soil surface (3 cm depth) submitted to infiltration–evaporation cycles. The swelling soil is a mixture of
loam and bentonite. Its hydraulic properties are described in Garnier et al. (1997a). The algorithms of transfer in swelling soil and direct method were selected. In this example, the moisture variable is the moisture ratio and there is one layer, characterized by 3 curves: the Van Genuchten (1980), Brooks and Corey (1964), and Braudeau (1988) models were selected for the retention, hydraulic conductivity and deformation curves, respectively. Parameters of these three curves are introduced on the panels (Fig. 2). The geometry factor Rs for the deformation curve is 10, which corresponds to a much larger variation of volume in the horizontal direction than in the vertical one. The swelling soil is subjected to different cycles of evaporation–infiltration depicted in Fig. 3 (3 evaporations and 2 rain events). Zero flux was imposed at the bottom. The main screen view (Fig. 4a) shows the distribution of water at a given time whereas Fig. 4b shows the cumulative value of some mean, global indicators. They are the evaporated water depth (cm), the soil water depth and the mass balance calculated from imposed fluxes and calculated water content. The mass error of 5% is acceptable. Fig. 4c shows the temporal changes of water content at 1, 2 and 3 cm depths while Fig. 4d gives the water content profiles at some intermediate times. We can see
1114
E. Perrier et al. / Computers & Geosciences 28 (2002) 1107–1118
Fig. 4. Simulation of infiltration–evaporation cycles in swelling soil Simulation of example fissure.sim (see test case A in text). Hydraulic properties are described in Fig. 2. Top boundary condition is described in Fig. 3. Bottom boundary condition is zero flux. (a) Main panel of software. Icons in main panel can be used to open other panels. (b) Panel to exhibit evolution of components (mean evaporation, infiltration, drainage, etc.) of mass balance. (c) Panel exhibiting at any time evolution of soil water content at three different depths. (d) Panel exhibiting at any depth evolution of soil water content at ten different times.
E. Perrier et al. / Computers & Geosciences 28 (2002) 1107–1118
1115
Fig. 5. Successive screen copies of main panel throughout simulation shown in Fig. 4.
the rapid changes of water content at 1 cm compared to the other depths. Fig. 5 shows the state of the soil at successive times throughout a total duration of about 10 days. When the soil is submitted to evaporation, the soil surface loses water with associated shrinkage so its shape becomes pyramidal with less water at the top (yellow) than at the bottom. This corresponds to the formation of cracks. When water infiltrates under a rain event, the soil swelling leads to crack closure at the surface. These processes simulated on a synthetic soil depict and explain rather well observed phenomena in real soils.
3.2.2. Parameter estimation by inverse method in a rigid soil during an infiltration experiment In a first step, before using the inverse method, an infiltration experiment was simulated with ECOUL using the direct method. Hydraulic properties of the Ariana silty clay loam soil obtained by Rieu and Sposito (1991) were introduced in the model. We used Rieu and Sposito’s model (1991) for the retention curve and Brooks and Corey’s model (1964) for the hydraulic conductivity curve. The soil (1 m depth) is submitted to an infiltration at its surface with a hydraulic head of 10 cm for 1 h. Initial conditions are those of a soil profile
1116
E. Perrier et al. / Computers & Geosciences 28 (2002) 1107–1118
E. Perrier et al. / Computers & Geosciences 28 (2002) 1107–1118
in equilibrium with an aquifer at the bottom. Simulated temporal evolution of water content at two depths and cumulative infiltrated water were registered. These data are called reference data. In a second step, we consider the previous soil as a real soil, submitted to the same conditions as in the direct simulation (as shown in Fig. 6c), and the reference data as ‘‘experimental’’ data that could have been easily measured on this soil using tensiometers at the two given depths and monitoring entering water. Then we test the inverse method in order to re-calculate some hydraulic parameters used in the direct simulation but here ‘‘forgotten’’. In this example, the file is registered as ariana.sim. When introducing hydraulic parameters, we consider that some parameters of the retention curve like residual and saturated water content and the minimal potential (Fig. 6a) are easily known and we will estimate by inverse method only the fractal parameter D linked to the slope of the curve. The forgotten value of D is 2.90 and we will try to investigate all possible values between 2.70 and 2.99 (Fig. 6a). As regards the hydraulic conductivity curve, we consider that the residual and saturated water contents are known as well as the conductivity value at saturation and that the parameter B is more difficult to measure. So the optimization process will deal with a search for B between 1 and 100, where the forgotten value was equal to 29 (Fig. 6b). The initial values for the search for D and B are set, respectively, to 2.85 and 10 (Figs. 6a and b) and correspond to the hydraulic characteristic curves displayed in Figs. 6a0 and b0 . The user can visualize the evolution of the objective function (Fig. 6e) which corresponds to the distance between here experimental data (red points) and simulated data (yellow points), presented in Fig. 6d at the beginning of the simulation and in Fig. 6f at the end. The weights given to each point of the surface infiltration curve are 1 (the default value) and 100 to each point of the sub-surface water content curves, to account for the different type of units and associated numerical values. The objective function decreased as expected when the number of iterations increased. Initial, current and optimal parameters are also shown in Fig. 6e after 59 iterations. The optimization process reached here 208 iterations and the final estimated values for D and B were, respectively, 2.899876 and
1117
29.010525. The estimated hydraulic properties at the end of the inverse method are shown in Figs. 6g and h which are to be compared to Figs. 6 a0 and b0 . This case study is given here as a test of the optimization algorithm. The inverse method is here totally satisfactory, of course because the experimental data were actually associated with a ‘‘perfect modeled soil’’, but in real situations, the exact solution is unknown, the number of estimated parameters as well as the number of iterations may be higher, it is even possible that the optimization process does not converge towards a unique solution.
4. Conclusion This paper aims first at informing a range of potential users about the existence of the software package ECOUL. ECOUL gives easy access to the modeling of deformation processes in swelling soils or other deforming porous media, which are widespread, but little studied as far as the modeling of their hydraulic behavior is concerned. We give an example of the way numerical modeling in a specific field could become available to a larger audience, or used for educational purposes by providing the means to implement interactivity and visualization concepts. The ECOUL software package is a new type of simulation tool that could later be enhanced by additional functionalities.
Acknowledgements The authors are grateful to Mrs. Chantal Bernard who designed all the iconography of the package, and to Mr. Nicolas Perrier who recently corrected some computer memory bugs.
References Bear, J., 1978. Hydraulics of Groundwater. McGraw-Hill International Book Co., New York, London, 569pp. Braudeau, E., 1988. Equation g!en!eralis!ee des courbes de retrait d’!echantillons de sol non structur!es (in French). Comptes Rendus de l’Acad!emie des Sciences S!erie 2 307, 1731–1734.
Fig. 6. Parameter estimation by inverse method in rigid soil during infiltration experiment. Simulation of example ariana.sim (see test case B in text). (a) and (a0 ) Initial parameters for water retention curve. (b) and (b0 ) Initial parameters for hydraulic conductivity curve. (c) Initial soil water profile as shown in main simulation panel. (d) Comparison between observed and simulated data (water content evolution at two different depth and amount of infiltrated water at top of soil) using (a) and (b) initial parameters for hydraulic properties. (e) Optimization process and evolution of difference between observed and simulated data using different sets of parameters ðD; BÞ for hydraulic properties. (f) Comparison between observed and simulated data using final, optimized parameters for hydraulic properties. (g) Optimized retention curve. (h) Optimized hydraulic conductivity curve.
1118
E. Perrier et al. / Computers & Geosciences 28 (2002) 1107–1118
Brooks, R.H., Corey, A.T., 1964. Hydraulic properties of porous media in hydrology Paper 3, Colorado State Univeristy, Fort Collins, pp. I–V, 1–27. Campbell, G.S., 1985. Soil Physics with BASIC—Transport Models for Soil–Plant Systems. Developments in Soil Science, Vol. 14, Elsevier, New York, 150pp. Celia, M.A., Bouloutas, E.T., Zarba, R.L., 1990. A general mass-conservative numerical solution for the unsaturated flow equation. Water Resources Research 26 (7), 1483– 1496. Gardner, W.R., 1958. Some steady-state solutions of the unsaturated flow equation with application to evaporation from a water table. Soil Science 85, 228–232. Garnier, P., 1996. D!etermination des caract!eristiques hydrodynamiques de sols d!eformables par la m!ethode inverse. Ph.D. Dissertation, Nancy University, France, 144pp. (in French). Garnier, P., Perrier, E., Angulo-Jaramillo, R., Baveye, P., 1997a. Numerical model of 3-dimensional anisotropic deformation and 1-dimensional water flow in swelling soils. Soil Science 162 (6), 410–420. Garnier, P., Rieu, M., Boivin, P., Vauclin, M., Baveye, P., 1997b. Determining the hydraulic properties of a swelling soil from a transient evaporation experiment. Soil Science Society American Journal 61 (6), 1555–1563. Gottardi, G., Venutelli, M., 1993. Richards: computer program for the numerical simulation of one dimensional infiltration into unsaturated soil. Computers & Geosciences 19 (9), 1239–1266. Haverkamp, R., Vauclin, M., Touma, J., Wierenga, P.J., Vachaud, G., 1977. A comparison of numerical simulation models for one dimensional infiltration. Soil Science Society American Journal 41 (2), 285–294. Ju, S.H., Kung, K.J.S., 1997. Mass types, element orders and solution schemes for the Richards equation. Computers & Geosciences 23 (2), 175–187. Leclerc, C., 1996. Le logiciel ECOUL, M!emoire de DESS Informatique. Universit!e Henri Poincar!e, Nancy, France, 50pp. (in French). Nedler, J.A., Mead, R., 1965. A simplex procedure for function minimization. The Computer Journal 7, 308–313. Perrier, E., 1989. ECOULv1.0: simulation num!erique et graphique des transferts hydriques dans un sol non satur!e. ORSTOM report FA41831/2, Institut de Recherches pour le D!eveloppement Documentation Services, France, 54pp. (in French).
Perrier, E., Garnier, P., Leclerc, C., 1997. Profile-scale simulation of water flow: a software package to visualize and estimate soil hydraulic properties effects, poster and demonstration, In: van Genuchten, M.Th., Leij, F.J., Wu, L. (Eds.), Proceedings of the International Workshop on Characterization and Measurement of the Hydraulic Properties of Unsaturated Porous Media, Riverside, California, October 22–24, 1997, University of California, pp. 349–354. Perrier, E., Rieu, M., Sposito, G., de Marsily, G., 1996. Models of the Water Retention Curve for soils with a fractal poresize distribution. Water Resources Research 32 (10), 3025– 3031. Richards, L.A., 1931. Capillary conduction of liquids through porous medium. Physics 1, 318–333. Rieu, M., Sposito, G., 1991. Fractal fragmentation, soil porosity, and soil water properties: I theory, II applications. Soil Science Society American Journal 55, 1231–1244. Simunek, J., van Genuchten, M.Th., Sejna, M., 1999. Using the HYDRUS-1D and HYDRUS-2D codes for estimating unsaturated soil hydraulic and solute. In: van Genuchten, M.Th., Leij, F.J., Wu, L. (Eds.), Proceedings of the International Workshop on Characterization and Measurement of the Hydraulic Properties of Unsaturated Porous Media, Riverside, California, October 22–24, 1997, University of California, pp. 1523–1536. Touma, J. 1987. TEST: mod"ele pour tester la representativit!e des caract!eristiques hydrodynamiques d’un sol non satur!e d!etermin!es in situ. ORSTOM report FA24686/2, Institut de Recherches pour le D!eveloppement Documentation Services, France, 26pp. (in French). Van Dam, J.C., Hugen J., Wesseling J.G., Feddes R.A., Kabat P., van Walsum, P.E.V., Groenendijk, P., van Diepen, C.A., 1997. Theory of Swap version 2.0. Simulation of water flow, solute transport and plant growth in the soil–water–atmosphere–plant environment. Report 71, Water Resources, Wageningen University, Technical Document 45, Alterra Green World Research, Wageningen, 167pp. Van Genuchten, M.T., 1980. A closed form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society American Journal 44, 892–898. Verburg, K., Ross, P.J., Bristow, K.L., 1997. SWIMv2.1. User manual. Division Report, CSIRO Division of Soil, CSIRO Information Services, East Melbourne, Australia, No. 130, 107pp.