SoDA project: A simulation of soil surface degradation by rainfall

ARTICLE IN PRESS

Computers & Graphics 30 (2006) 494–506 www.elsevier.com/locate/cag

Natural Phenomena

SoDA project: A simulation of soil surface degradation by rainfall Gilles Valettea, Ste´phanie Pre´vosta, Laurent Lucasa,, Joe¨l Le´onardb a

CReSTIC/LERI/MADS, Universite´ de Reims Champagne Ardenne, Reims, France b Unite´ d’Agronomie INRA de Laon-Reims-Mons, France

Abstract The main objective of the SoDA (Soil Degradation Assessment) project is to realize a simulator of soil surface degradation by rainfall at the meter scale and including visualization. Soil surface structure and morphology deeply influence a lot of processes of high agronomic and environmental relevance, such as mass and heat transfer through the soil–atmosphere interface, runoff and erosion, seed germination and seedling emergence. The soil surface structure of agricultural field is in continuous evolution: it is strongly affected by tillage, and in between tillage operations, erosion by rainfall and runoff causes a progressive degradation of the structure whose intensity and speed partly depend on the initial state associated to tillage modalities. A soil surface degradation model could allow one to predict this evolution of the soil surface structure, and even to help choosing adequate tillage practices and sowing dates. Erosion modeling has been addressed by soil scientists but also by computer graphic scientists in order to add realism to virtual landscapes. Mixing both of these points of view would be interesting to simulate and visualize the evolution of the soil surface of a cultivated soil. Based on a 3D cellular automata approach using the knowledge accumulated by soil scientists about the physical processes involved in erosion, the principles of our simulator and its first implementation are presented in this paper. r 2006 Elsevier Ltd. All rights reserved. Keywords: Simulation; Modeling; Cellular automata; Erosion; Soil surface degradation; Soil crusting; Direct volume rendering

1. Introduction Soil surface structure and morphology deeply influence a lot of processes of high agronomic and environmental relevance, such as mass and heat transfer through the soil–atmosphere interface [1,2], runoff and erosion [3,4], seed germination and seedling emergence [5,6]. The soil surface of arable soils is in constant evolution both because of farming operations and climate. Major aspects of this evolution are the formation of soil crusts and the development of cracks [7–9]: presence of aggregates and crust due to rainfalls can mechanically inhibit seedling emergence whereas presence of cracks can allow seedlings to break through the soil surface. A model predicting soil structure under different initial soil conditions and climatic scenarios would be an efficient way to select adequate tillage and sowing practices. A portion of land of 1 m2 Corresponding author. Laurent Lucas, CReSTIC LERI, Rue des Craye`res, B.P. 1035, 51687 REIMS Cedex 2. E-mail address: [email protected] (L. Lucas).

0097-8493/$ - see front matter r 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.cag.2006.03.016

seems to be specially adequate for such a study. This scale allows clods and cracks observation and a precise study of the solid particles transfer associated with rainfall and runoff. There has been many studies of erosion by soil scientists [10–12], at scales ranging from small plots to entire catchments, under various climatic conditions. Often, the main objective of the developed models is to allow soil loss evaluation rather than to analyze the evolution of the soil surface and its relief. Besides, during the last two decades considerable progress has been made towards developing efficient models for generating approximations to natural terrain [13–18]. Although these kinds of models generally focused on large scale relief evolution, they are potentially useful and should not be neglected. A popular type of model simulates erosion of the entire landscape using transport equations for the removal of solid material [19]. We think that it can be profitable to bring together both of these fields of research by blending the objectives of a pedological simulation and those of a realistic visualization, i.e., physical accuracy and visually plausible results.

ARTICLE IN PRESS G. Valette et al. / Computers & Graphics 30 (2006) 494–506

495

The objective of our project, called SoDA (Soil Degradation Assessment), is therefore to develop and validate a dynamic simulation of soil erosion at the meter scale, keeping in mind a constant care for visualization.1 To ensure the correctness of our model, we have a preoccupation with a validation by the images and numerical data that simulation will provide. For that purpose collaborations with experts in this field of research are essential and we do cooperate with an INRA2 laboratory (Agronomy Unit Laon-Reims-Mons, France). 2. Description of the processes involved in soil erosion and surface crusting In order to simulate the soil structure evolution we need to have a description of the different processes involved in this phenomenon. In our literature review we have encountered three main kinds of soil erosion: wind erosion [20], thermal erosion [21] and hydraulic erosion [17]. Although thermal weathering, due to thermal shocks and gravity, influences the topological aspect of the soil and has been modeled for several visual simulations [15,19,21,22], we will concentrate on hydraulic erosion which is the main process responsible for the evolution of the soil surface of agricultural soil. Hydraulic erosion can be caused by rainfall or running water (Fig. 1). First, raindrops cause disaggregation: continuous matter or aggregates break down into smaller fragments [23]. Then, these fragments can be mobilized and transported by splash effect or by runoff [24,25]. Rainfall and runoff cause structural reorganization by the formation of crusts. Crusts are thin soil surface layers more compact and hard, when dry, than the material directly beneath. Generally, two main types of crust are distinguished by their mode of formation: structural crusts and sedimentary crusts [26]. A structural crust develops in situ and is the result of gradual coalescing of aggregates caused both by particle translocation and by raindrop compaction, whereas a sedimentary crust is formed by deposition of the particles suspended in overland flow [27]. Crusts hamper seedling emergence, reduce infiltration and favor runoff, puddling and thus erosion. The photographs of Fig. 2 show the visual perception of soil evolution and degradation with the generation of both types of crust from the initial soil shown in the first photograph: structural crust, in the second photograph, and sedimentary crust, in the last one. This evolution can be characterized by four main points: (1) aggregates outlines become less sharp and can even disappear, (2) filled orifices become more and more numerous, 1 Informations and multimedia material are available at www.sodaproject.com 2 INRA stands for ‘‘Institut National de la Recherche Agronomique’’ (i.e., National Institute of Agronomic Research).

Fig. 1. The processes of erosion.

(3) surface roughness is decreasing, (4) the color of the soil and its reflectance properties are changing: the particle segregation results in either chemical (minerals) or physical (particle size) differentiation, which are both correlated with reflectance properties of soil [28]. It is worth noticing that these criteria are purely visual. For an agronomist the first and very important tool is direct visual observation. That is why in our project we have a constant care for visualization and we want to get images as well as numerical results from our simulator, in order to allow visual comparisons between simulation and reality. In particular, we aim to recreate the same visual evolution. 3. Related works We present in this section a brief review of models found in the field of Computer Graphics and Soil Science followed by a few comments. 3.1. Computer graphics models Musgrave et al. [19] demonstrate a new method for creating mountain fractal terrains with the use of height fields. They suggest two erosion algorithms which simulate hydraulic erosion by flowing water, ignoring evaporation and infiltration, and thermal weathering due to the thermal shocks which chips away steep inclines and forms talus slopes. Kelley et al. [13] produce images of realistic-looking terrain with an algorithm consisting of two distinct steps.

ARTICLE IN PRESS 496

G. Valette et al. / Computers & Graphics 30 (2006) 494–506

Fig. 2. Formation of crusts in the field (photos C. Du¨rr, INRA). (a) Initial soil, 11/25/97; (b) structural crust, 12/15/97; (c) sedimentary crust, 01/15/98.

The first steps represent the creation of a drainage network of rivers according to geologically known data and behavior. In the second step the terrain is modeled as a set of surfaces under tension fitting previously generated data. Benesˇ and Forsbach [21] introduce new data structure for visual simulation of 3D terrains. This representation is based on horizontal stratified layers consisting of one material. They run the previous algorithm simulating thermal erosion on this representation. The authors [17] divide the hydraulic erosion process into four independent steps that can be applied independently to achieve high level of realism: new water appears, water erodes the underlying terrain and captures the material, water and the suspended material are transported and finally, because of the evaporation, water deposits the material at another location. Benesˇ and Arriaga [29] present an efficient algorithm for visually modeling table mountains using two different types of erosion, one for hard material (rock) and the other for soft material (sand), corresponding to two regular height fields. In order to create some realistic landforms, Roudier et al. [14] simulate geologically contrasted terrains and apply deterministic erosion processes to them. The erosion on any point of the landsurface is therefore related to local geological parameters. The available erosion laws simulate mechanical erosion, chemical dissolution and alluvial deposition. Nagashima [15] creates eroded valley and mountain terrains with layer traces by using geological parameters and simulating the physical erosion of river flow, precipitation, and thermal weathering. Marak et al. [30] try to formalize some of the algorithms used for simulation of erosion under one algorithm that is based on rewriting matrices. A height field representing the terrain is considered as the matrix, and context sensitive rewriting of this matrix represents the erosion process. The authors define classes of matrices that are used for different erosion algorithms.

Chiba et al. [16] present a simple ‘‘quasi physically based’’ method for simulating the topography of eroded mountains based on velocity fields of water flow. In this method, infiltration and evaporation are ignored, the velocity fields of water flowing down the face of a mountain are calculated by simulation of the motion of ‘‘water particles’’, and these velocity fields are used for simulating erosion. Neidhold et al. [18] present a method that combines a non-expensive fluid simulation based on Newtonian physics approach, with an erosion algorithm, permitting the user to control parameters of several simulation steps such as rain fall, water sources, fluid transportation, material dissolution, material deposition and evaporation.

3.2. Soil erosion models In the field of Soil Science, prediction of soil erosion generally means prediction of soil loss rather than analysis of the evolution of the soil surface and its relief. Current models for soil erosion by water such as WEPP (Water Erosion Prediction Project) [31] and EUROSEM (European Soil Erosion Model) [32] consider an eroding hillslope to be conceptually partitioned into rill and interrill areas.3 Soil loss in rill and interrill areas is calculated separately and summed to yield a prediction of total soil erosion. Favis-Mortlock et al. [33] have developed the ‘‘RillGrow’’ model which applies simple rules at a millimeter scale, to govern the iterative interaction between microtopography, runoff routing and soil loss. In this model emphasis is put on the evolution of the soil surface rather than on sediment exportation. No distinction is made between interrill and rill areas. Soil is decomposed in a regular grid of cells and runoff is considered to be 3 A rill is a small channel as one formed by soil erosion. The interrill area is the area comprised between rills.

ARTICLE IN PRESS G. Valette et al. / Computers & Graphics 30 (2006) 494–506

discretized into packets which move from cell to cell. Splash erosion, deposition and infiltration are not modeled. D’Ambrosio et al. [34,35] and Avolio et al. [36] suggest a cellular automata (CA) model for soil erosion by water. This model involves a larger number of states, including altitude, water depth, total head, vegetation density, infiltration, erosion, sediment transport and deposition. A similar CA approach, based on the model of ‘‘precipitons’’ introduced by Chase [37], is used by Luo et al. [38] for the WILSIM (Web-based Interactive Landform Simulation Model) project. In order to simulate landscape erosion and deposition Haff [39] uses a CA model called ‘‘waterbot’’. Waterbots are discrete units of runoff that are moving autonomously and asynchronously across the landscape and are able to pick up and deposit sediment. Servat [40] suggests a description of flows in terms of heterogeneous agents that interact in a continuous space and whose various laws of interaction enable one to account for the coupling of concurrent processes. This research has led to the development of the RIVAGE4 simulator of runoff and infiltration, with some additional trials to incorporate erosion processes. 3.3. Comments It is worth noting that evaporation, splash and crusting are quasi absent of the simulations encountered in this review. This can be easily explained by the fact that these simulations generally aim to create virtual landscapes or to study soil erosion on a large scale for which these processes can be neglected, which is not the case of our project. Nevertheless some other processes like runoff or deposition have been modeled with success in different ways and these approaches are potentially useful. The most complete model is Servat’s multi-agent model RIVAGE. However, this approach does not seem to be adapted to our objective because it is centered on the water particles and their mobility and not on the structure of the soil. Techniques based on an extended CA model seem to be of great interest in order to obtain correct results in simulation of erosion. 4. Our approach Wolfram [41] has shown that CA with simple underlying rules can produce highly complex behavior and can mimic many features of what we see in nature. The principle of CA is not to describe a complex system with complex equations, but to let the complexity emerge by interaction of simple individuals following simple rules. On the one hand, we need a volumetric model because we must take into account the structure of the subsurface and its evolution, and on the other hand, soil erosion can be considered as a system evolving exclusively by means of 4

RIVAGE stands for ‘‘Ruissellement et Infiltration Vus par des AGEnts’’ (i.e., runoff and infiltration seen by agents).

497

local interactions. For both these reasons CA could be applied successfully and that is our basic approach. Furthermore, CA models have been employed in the field of Computer Graphics in order to simulate stone weathering [42], corrosion [43] or crack patterns [44,45]. Nevertheless, the classical CA model presents certain limits in simulating complex macroscopic phenomena: they are based upon elementary automata, with few states and a simple transition function [46]. Therefore, in order to deal with specific aspects of the processes we want to simulate, and to introduce some complexity and flexibility, we use an extended CA model allowing a large number of different states and a more complicated transition. 4.1. The cells CA are based on a regular division of space in cells, each one characterized by a state with a few possible values. In our extended model we consider a cell as a cube (2 mm side) with a large number of substates described afterwards. For each substate, permitted values must form a finite set. As they represent mostly physical quantities, they are continuous variable and therefore they have to be discretized with a sufficient number of digits. We have chosen to use 64 bits for each cell and Table 1 shows the distribution of these bits. The soil system is composed by three different kinds of components: soil, liquid and gas. We use a system where these three phases are additive. Such a system is provided by the use of specific volumes. We can write the total volume of a soil sample as: V cell ¼ V s þ V l þ V g ,

(1)

where V s is the volume of soil, V l is the volume of liquid and V g is the volume of gas. As our 3D-grid is regular, the height of each component in a cell can be considered equivalent to its corresponding volume, so we can do calculations on heights: H cell ¼ H s þ H l þ H g ,

(2)

where H s is the height of soil, H l is the height of liquid and H g is the height of gas. As H cell is a fixed constant (2 mm), H g is easily computed with the knowledge of H s and H l which are therefore the first substates of the cell. Table 1 Distribution of the 64 bits of the cell substates Bits

Number

1–6 7–12 13–15 16–21 22–30 31–32

6 6 3 6 9 2

33–45 46–61 62–64

13 16 3

Substate Height of matter Height of moisture # Fragments of class 1 # Fragments of class 2 # Fragments of class 3 Not used # Fragments of class 4 # Fragments of class 5 X-link, Y-link and Z-link

ARTICLE IN PRESS G. Valette et al. / Computers & Graphics 30 (2006) 494–506

498

During rainfall, the size of the fragments produced by the disaggregation processes varies with the initial size distribution of particles, the moisture state of the soil, the intensity of the stresses exerted by raindrops. In addition, sedimentary crusts are defined by the fact that particles are sorted according to their size during sedimentation. The size of soil fragments determines if runoff will be able to transport them or not. It is thus very important to have and to keep some information about particle size. So, in addition to these different heights we keep in a cell the numbers of particles of certain given sizes, the total volume of which is now given by: V cell ¼ Acell :ðH s þ H l þ H g Þ þ

NP X

N i :V i ,

(3)

i¼1

where Acell is the area of a cell (i.e., 4 mm2 ), NP the number of different particle sizes, V i the volume of a particle of size S i , N i the number of particles of that size. We do not need to keep the specific value H g because it can be obtained by a simple operation, even if it is less immediate than previously. In the current implementation we use five classes of cells corresponding to sizes of 1000, 500, 250, 100 and 50 mm which are common sizes for sieve analysis. The presence of aggregates, their location, their size and their shape are important features of the soil structure because they can make a direct link between the state of the soil and soil tillage, seedbed preparation and sowing techniques. In order to be able to represent aggregates we add three substates which correspond to links with adjacent cells (actually there are six possible links for one cell but they are bidirectional; in order to avoid redundant information and to save bits, when one link exists between two cells, we store this information in only one cell). We use also another important substate which defines if a cell is on the surface or not. That means if it undergoes to the influences of the external world and receives raindrops or soil particles transported by splash. For this type of cell we need other substates such as the heights of rainwater, infiltrated or running water and their porosity. These substates are meaningless for the subsoil cells so a 2D-grid is sufficient to store them. 4.2. The transition function In a classical CA model, at the time t ¼ 0, cells are in arbitrary states and the CA evolves changing the state of all cells simultaneously at discrete times, according to a transition function, the input of which is given by the states of the neighboring cells. In our model we replace this unique transition function by several transition functions, one for each simulated process, with particular neighborhoods. Moreover, we want an open model which offers a choice between different functions for each process and permits the inhibition of some processes. Therefore, we obtain a global transition function F which can be written

as a matrix, one column for 0 F 1;1 F 1;2 ... B F 2;1 F 2;2 ... B B . .. .. B . F ¼B . . . B BF @ k
1;1 F k
1;2 . . . F k;1

F k;2

...

each process: 1 F 1;m F 2;m C C .. C C , . C C C F k
1;m A F k;m

(4)

where m is the number of processes, k is the maximal number of available simulating functions for the processes. In order to be able to inhibit some processes we have 8i 2 N with 1pipm F 1;i ¼ Id. The input of each function F j;i is defined by its own neighborhood Gj;i . For example, we have for the first line 8i G1;i ¼ ð0; 0; 0Þ (this neighborhood defines internal transformations). Following Avolio et al. [36], we allow input from the external world to certain cells of the CA, in order to describe an external influence which cannot be represented by local rules. For our simulation these external influences are raindrops and particles projected by splash. It is worth noticing that for each process we can use various kinds of modeling (stochastic, physically based, phenomenological, empirical,. . .) and the next sections show some examples of different approaches. We formalize our CA model by the set ðC; Q; G; h; F ; s; PÞ:  C  Z3 is the 3D matrix representing the locations of the cubic cells in a discrete space,  Q ¼ ðQ1 ; . . . ; Qp Þ is a p-tuple of sets of substates,  the application h : N  C ! Q gives the state of C at the instant t; hðt; cÞ is denoted ct ; c0 defines the initial state of c,  F is the k  m matrix of transition functions,  G is also a k  m matrix of neighborhoods,  s is the application so that sðiÞ ¼ j if and only if we choose to use the function F j;i for the simulation of the process i,  P is the finite set of the global parameters which have an effect upon the transition functions. Fundamental global parameters are the cell size and the time step. The evolution of a cell c between the instants t and t þ 1 can now be written by this composition of functions: 8c 2 C;

ctþ1 ¼ ðF sðmÞ;m      F sð1Þ;1 Þðct Þ.

(5)

4.3. Current implementation We have based our main loop on the succession of the processes involved in the natural phenomenon we attempt to reproduce. Here is the simple algorithm of our simulation: Cells initialization Repeat Time increment Rain drops generation

ARTICLE IN PRESS G. Valette et al. / Computers & Graphics 30 (2006) 494–506

iteration and thus we have a better accuracy. The size of the raindrops can vary between a minimal and a maximal value but is kept under the size of one cell, which is compatible with the intensity of the rainfalls we want to simulate. The number Nr of raindrops which fall on the surface during an iteration is given by:

Detachment by splash Transport by splash (Evaporation) Infiltration Deposition or mobilization Runoff Transport by runoff Until simulation is completed The movement of water is due to runoff and infiltration. We do not take into account the effect of evaporation because the infiltration rate is always much greater than the evaporation rate, especially during a rainfall. The fragments are created by the splash detachment and all fragments can be transported by raindrops. Calculation of the quantity of fragments mobilized (or deposited) by the surface water is made before the runoff in order to compute the transported quantity by runoff from a correct value. 4.3.1. Soil initialization One first important step of our work is to generate an initial volume of soil whose properties reflect the state of the real soil we want to study. In the current implementation this stage is very simplified and we have only an accurate topography given by a digital elevation map obtained by laser profile metering. We assume that the initial soil is dry. The porosity is chosen between 25% and 40%, which are values measured in the field, and slight random variations are applied in each cell. 4.3.2. Time increment and rain generation A rainfall event can be defined by a hyetogram which represents rain intensity variations within a given duration (Fig. 3). We compute the time step in order to respect a given percentage of the surface cells which have to receive a raindrop (the current value is currently between 5% and 25%) and when we do the uniform random distribution of the raindrops (based on L’Ecuyer’s random number generator with Bays–Durham shuffle [47]), we insure that a cell cannot receive more than one raindrop until another given percentage (e.g., 75%) of the surface cells have received one raindrop. With these limitations we have not to deal with important quantities of rainwater in one

Intensity (mm / h)

120

Simulated hyetogram Real hyetogram

100 80 60 40 20 0 0

10

20

499

30 40 Time (min)

Fig. 3. Example of hyetogram.

50

60

70

Nr ¼ ðNxy  Acell Þ 

I  Dt , ðHw  Acell Þ

(6)

where I is the intensity of the rain (mm/h), Nxy is the number of surface cells, Dt is the time step (h), Acell is the area of a cell ðmm2 Þ, Hw is the height of water in a cell corresponding to the minimal size of a raindrop (mm). After a trivial simplification we obtain: Nr ¼ Nxy 

I  Dt . Hw

(7)

In order to water the terrain, we randomly pick a surface cell on the surface of the terrain. If this cell is not marked, we mark it, we pick a size between the minimal and maximal values and we add this height to the rainwater of this cell. This rainwater will run or infiltrate in the next steps of the algorithm. We repeat this process until we have reached Nr . When the number of marked cells represent the desired percentage, we delete all the marks. 4.3.3. Detachment Two distinct mechanisms for soil detachment are to be considered [27]. The first one is disaggregation by entrapped air compression when fast moistening of a dry aggregate occurs. The detached mass is given at a first approximation by the formula: dm ¼ a dh, with dh the unit of rain height and a an index of soil sensibility, function of its moisture and composition. The second one is the disaggregation without bursting, that is the one due to the impact of raindrops which is able to detach particles from the mass. The detached mass is given again at a first approximation by the formula: dm ¼ b ke dh, with ke the kinetic energy of the raindrop and b a sensibility factor, function of its wetness and compactness. As rain kinetic energy is not a commonly measured meteorological parameter, empirical laws relating the rain kinetic energy to the more easily available rain intensity have been proposed. We use the empirical linear-log formulation of Wischmeier and Smith [48,49]: ( 11:9 þ 8:7  logðIÞ if Ip76 mm=h; ke ¼ (8) 28:3 if I476 mm=h: With the assumptions that raindrops size is never greater than the cell size and that the products of detachment are particles smaller than a cell, disaggregation is considered as a local transformation. In each cell which is touched by a raindrop, particles are created. To decide the size of these particles, we use a fuzzy logic approach [50,51] based on simple rules which take a coefficient of compactness of the

ARTICLE IN PRESS G. Valette et al. / Computers & Graphics 30 (2006) 494–506

500

soil (representing the surface particle size distribution) and the size of the raindrop as inputs, and give an index in an array of sets of percentages as output. The result is shown in Fig. 4. Each of these sets contains one percentage for each class of particles. The index returned by the rules is actually a real number r which is comprised between two integer indexes i and j ¼ i þ 1. To get the percentage P0p we use for a class p we make a linear interpolation between the percentages in the two sets Pi and Piþ1 : 8p;

P0p

¼ ðr
iÞPi;p þ ½1
ðr
iÞPiþ1;p .

(9)

The coefficient of compactness C is defined as: C¼1


Ns X Mi

Mt

p

,

(10)

where Ns is the number of different particle sizes, p is the index of the first class of particles that are considered as small (i.e., smaller than 250 mm), Mi the mass of the particles of size i, Mt the total mass of the soil in the cell. To decide how many fragments are detached, two solutions are implemented: to use the empirical calculation of the detached mass (with arbitrarily fixed values a ¼ 0:5 and b ¼ 0:05) or to use other fuzzy rules which take the rain kinetic energy and the size of the raindrop as inputs.

Disaggregation index

4.3.4. Transport by splash Soil transport by rainsplash has been measured using a variety of approaches, including splash cups, trays, and boards. The most frequent assumption in the literature is that the spatial distribution of particles splashed from a point source can be described by a negative exponential function [24,52]. The intensity of the raindrops’ impact, the size of particles and the local topography determine projection distance. As there is a lack of knowledge about the size selectivity of the phenomenon, we have chosen to use a fuzzy logic approach to estimate the sizes and the number of the particles ejected from the cell where the raindrop falls. The distance is then randomly chosen with respect to an exponential distribution of given mean for each class. This mean depends on the soil and on the size of

5 4.5 4 3.5 3 2.5 2 1.5

5 4.5 4 3.5 3 2.5 2 1.5 1 0

-1

20

40 Dro 60 p si 80 ze

100 1

0

-0.5 s tnes

0.5

ac omp

C

Fig. 4. Index of detached fragments obtained by fuzzy logic.

the fragments and is given by an expert. After that one of these two rules is observed: (1) the point of destination is accepted if it is lower than the source point, (2) the distance is increased if the point of destination is lower than the source point. Assuming that our terrain patch is surrounded by patches that have the same behavior, we treat the horizontal boundaries as if the patch was a torus (both side edges are connected for both horizontal dimensions). 4.3.5. Infiltration Infiltration is modeled with the physically based Green– Ampt method [53] which makes the following assumptions: the wetting front advances at a constant rate, volumetric water contents remain constant above and below the wetting front as it moves and soil water suction below the wetting front remains constant. It needs four parameters: the saturated hydraulic conductivity Cs , the saturated and initial volumetric water contents ys and yi , and the soil water suction at the wetting front C. Values for these parameters can be found in the literature [54]. The potential infiltration rate is given by:
 Ws þ C Wi f ¼ Cs 1 þ , (11) with L ¼ L ys
yi where Ws is the local height of surface water and Wi is the height of infiltrated water. An example of infiltration is shown in Fig. 9. 4.3.6. Mobilization and deposition Mobilization and deposition are based on the principle of capacity of transport [17]. A quantity of water W is able to transport a quantity of matter S, which is calculated with a coefficient of saturation Ks (kg/l) which depends on the soil characteristics. Between two iterations, if the quantity of water has decreased because of infiltration or evaporation, a quantity of matter deposits. If the quantity of water has increased because of the rainfall, a quantity of matter is mobilized and can be transported by runoff. At the instant t the quantity of matter mobilized in the water of a surface cell is given by the formula: St ¼ Wt  Ks .

(12)

As we want stþ1 ¼ St , we compare this value with the quantity st really present in the water in order to calculate the quantity of matter deposited or mobilized Dm : ( deposition if st 4St ; Dm ¼ jst
St j with (13) mobilization if st oSt : As for the splash we have to determine the size and the number of the fragments mobilized or deposited by the flow. For that we use a simple rule: we try to get the quantity of matter St by picking first the coarsest

ARTICLE IN PRESS G. Valette et al. / Computers & Graphics 30 (2006) 494–506

fragments and we repeat this operation with the other classes of fragments, until we get the needed quantity.

4.3.7. Runoff Proving the versatility of our simulator the implementations of infiltration and runoff are very different. We do not use any physically based equation for runoff but we simply calculate the report of water between the surface cells, using the fact that water flows in the direction of the highest gradient [14,21] (Fig. 5). Considering the 2D-space SðCÞ formed by the surface cells, the neighborhood of a single cell can be the von Neumann neighborhood or the Moore neighborhood [35] as shown in Fig. 6. For each surface cell c we determine the lowest cell mðcÞ in its Moore neighborhood gðcÞ. mðcÞ ¼ arg minðHi þ Wi ; i 2 gðcÞÞ,

(14)

where Hi is the height of the terrain and Wi is the height of running water in the cell i. The height of water to redistribute from the cell c is:
 ðHc þ Wc Þ
ðHmðcÞ þ WmðcÞ Þ DwðcÞ ¼ min Wc ; . (15) 2 We compute the heights of water to redistribute for all the surface cells at the instant t and then we update all these

Fig. 5. Computation of runoff.

(a)

(b)

Fig. 6. Neighborhoods of a single sell c in a two-dimensional regular grid. (a) The von Neumann or ‘‘primary’’ neighborhood. (b) The Moore or ‘‘complete’’ neighborhood.

501

cells to obtain the heights of water at the instant t þ 1: ( tþ1 Wc ¼ Wtc
DwðcÞ; (16) 8c 2 SðCÞ t Wtþ1 mðcÞ ¼ WmðcÞ þ DwðcÞ: When a certain quantity of water moves because of runoff, it transports a proportional quantity of mobilized matter to the cell of destination. We calculate this proportion for all sizes of fragments and move this number of fragments to the cell of destination, treating in the same way all the sizes of fragments. 5. Visualization Our simulator has to produce numerical results (e.g., mass of displaced soil) as well as images or animations. Therefore, we want to visualize three-dimensional data sets and for that result we use direct volume rendering (DVR) [55]. DVR is based on the premise that the data values in the volume are themselves a sufficient basis for creating an informative image. What makes this possible is a mapping from the numbers which comprise the dataset to the optical properties that compose a rendered image, such as opacity and color. The most attractive aspect of DVR is actually its defining characteristic: it maps directly from the dataset to a rendered image without any intermediate geometric calculations. This is in contrast to traditional isosurface rendering, where the rendering step is proceeded by the calculation of an isosurface for a particular data value [56]. DVR can also avoid the time-consuming processes of segmentation and classification, because visualization can be done without any high-level information about the material content at each voxel.5 The basic benefit of DVR over isosurface rendering is that it provides much greater flexibility in determining how the voxels contribute to the final image. Voxels over a range of values can all contribute to the image, with varying amounts of importance, depending on the transfer function. The MC6 Slicing Algorithm we use, described by Benassarou et al. [57], is an accelerated slicing algorithm for interactive volume rendering of structured grids. It requires a small amount of memory and provides adaptive rendering for improved image accuracy as well as progressive rendering for rapid feedback at interaction time. This method, as other slicing methods, use the 3D texture mapping hardware. The rectilinear volume data are first converted to a 3D texture. Then, a number of planes perpendicular to the viewer’s line of sight are clipped against the volume bounding box. The texture coordinates in parametric object space are assigned to each vertex of the clipped polygons. During rasterization, fragments in the slice are trilinearly interpolated from 3D texture and projected onto the image plane using 5 Voxel stands for ‘‘VOlume piXEL’’ and represents a ‘‘three-dimensional pixel’’, i.e., the smallest unit of a cubic 3D-grid just as a pixel represents the smallest unit of a 2D-grid. 6 MC stands for ‘‘Marching Cubes’’: in the method we use [57], the original marching cubes algorithm [56] is used to efficiently polygonize an approximation of the intersection between a surface and a cube.

ARTICLE IN PRESS 502

G. Valette et al. / Computers & Graphics 30 (2006) 494–506

Fig. 7. Direct volume rendering. (a) View-aligned slice stacks with 3D texture mapping. (b) A sufficient number of slices is needed to create the visual impression of a continuous terrain.

Fig. 9. Examples of hybrid visualization during infiltration. (a) Beginning of the infiltration. The surface is shown as a triangulated surface, subsoil and infiltrated water are shown with direct volume rendering. A clipping plane is used on the right side of the volume. (b) Continuation of the infiltration. Another clipping plane is used horizontally in order to show the moisture inside the clods. Fig. 8. Two aspects of the same data volume. (a) Porosity. (b) Moisture.

adequate blending operations (Fig. 7). In the current implementation the 3D texture is a 8-bits texture and mapping values to optical values is made with 1D texture lookup which gives the 256 possible RGBA values. Thus, we have to reduce the amount of data contained in each cell (64 bits) to an appropriate 8-bits value. The difficulty to find a good transfer function for this operation is the main drawback of this method but it gives the opportunity of showing different aspects of the same data volume by retaining only certain data for the cell to voxel transformation (see Fig. 8). Moreover, it is possible to use hybrid visualization, i.e., to mix classical rendering techniques and DVR. Other visual tools like clipping planes may be useful in order to see for example the infiltrated water (see Fig. 9). 6. Results We have tested our model by comparing its predictions to the results of a rain-simulator experiment on a patch of

reconstituted soil (see Fig. 10). We use a very small area because we need a well controlled experiment. This is only possible in the lab, where we have severe constraints on the size of the studied area. In particular, the area covered with the rainfall simulator impose a size of 0:5 m  0:5 m to ensure that the rain is evenly distributed. However, we take care of the representativity of the soil surface studied, by sampling an unremoulded seedbed surface directly from the field. The real experiment duration was 80 min, with a 30 mm/h rainfall and raindrops of about 2 mm size. At different moments of this experiment, soil surface roughness was measured using a laser profile meter7 in order to compare results of the virtual and real simulations. The resulting images are close enough in their visual appearance (see Fig. 11): they show that the simple rules simulating the surface degradation on our soil model give 7 The total measurement process takes 4 h (62 500 point measurement for an area of 0:5 m  0:5 m). However, each point measurement is quasi instantaneous and does not cause any significant increase in heat nor evaporation.

ARTICLE IN PRESS G. Valette et al. / Computers & Graphics 30 (2006) 494–506

Fig. 10. Initial patch of soil, real and virtual. (a) A photograph of the initial soil. (b) The initial soil in our simulator ð40  40  10 cm3 Þ.

plausible results, comparable with those of the rainsimulator experiment. We observe on both images some characteristics of soil surface structure degradation: the decreasing of the aggregates outlines sharpness, the disappearing of orifices and the decreasing of surface roughness. Of course the changes in the color of the soil cannot be observed on the image of the real experiment which is only the visualization of a simple digital elevation map. On the contrary, on the image produced by our simulator, the choice of the transfer function permits to see a sort of crust formed by particles on the surface. Both the static results (visible on images) and the dynamics of the simulation (visible on recorded animations) are interesting. For example, at the beginning of the virtual simulation, water excess occurs and after a while one can see puddles, which is very close to reality. At the end of the virtual simulation rain stops and water can be infiltrated again. That shows that we have recreated a Hortonian runoff [58]: there was too much water for the infiltrability rate but the soil was not saturated. The results of this first experiment are qualitatively and visually satisfactory even if virtual simulation acts too fast (this default should be corrected if we take into account

503

Fig. 11. Comparison between the results of the virtual simulation and those of a rain-simulator experiment. (a) Result of our simulation after 600 iterations (corresponding to 10 min), crusts formed by particles are visible because we have chosen a transfer function which shows where the particles are concentrated. (b) Real terrain after 80 min of a 30 mm/h rain in a rain-simulator. (c) Our virtual simulation after 40 min with a 30 mm/h rain.

protection against splash which is constituted by surface water). To get some quantitative results, we make a comparison between the initial height and the final height of the terrain for each surface cell ði; jÞ: ðH final
H initial Þ i;j i;j 2 ½
1; 1, maxðjDH i;j jÞ ( ðDH i;j ; 0; 0Þ  red if DH i;j 40; colori;j ¼ ð0;
DH i;j ; 0Þ  green if DH i;j p0:

DH i;j ¼

The last formula simply indicates that increase of height appears in red and decrease appears in green in the visualization of this comparison. The virtual simulation data gives the result shown in Fig. 12(a). It seems normal that hollows and dips should be filled with particles and therefore that their height should increase. But comparison on the real simulation data gives another result: as Fig. 12(b) shows, there is no increase of the height. We

ARTICLE IN PRESS G. Valette et al. / Computers & Graphics 30 (2006) 494–506

504

simulation slower. The memory requirements depend on the size of the terrain and Fig. 13(b) shows that this variation is linear with the surface area and the depth. The slopes of these lines are different because of the data corresponding to the surface cells whose number is constant when the depth is varying. 7. Conclusions and future work

Fig. 12. Difference between virtual and real simulation. (a) Virtual simulation. (b) Real simulation.

Time for one iteration (s)

2

Rain intensity 30 mm/h 20 mm/h 10 mm/h

1.5

1

0.5

0 10x10 20x20

30x30

2

Memory occupation (Mo)

60

50x50

Terrain depth (cm) 4 6 8

10

12

Depth=10 cm area varying Area=50x50 cm2 depth varying

50 40 30 20 10 0 10x10 20x20

(b)

40x40

Surface area (cm2)

(a)

30x30

40x40

Surface area

50x50

(cm2)

Fig. 13. Run-time and memory requirements. (a) Time requested for one iteration. (b) Memory requested.

think that its due to a process we have not yet taken into account: the soil compaction due to the raindrops impact. The simulation times for the terrain of 40 cm side for a depth of 10 cm on a 3 GHz Intel CPU represent about 6 s per iteration of our algorithm for the parameters chosen in the virtual erosion detailed in this section (thus the 40 min required about 4 h for a time step of 1 s). Fig. 13(a) shows how the time requested for one iteration (with different parameters) of 1 s is linearly varying with the surface area. We can notice that a stronger rain intensity makes the

We are developing a simulator of the evolution of soil surface structure of cultivated soils under rainfall, at the meter scale. This is an important issue as it has both theoretical and practical interest. Knowledge has been accumulated by soil scientists about the various physical processes controlling the degradation and crusting of soil surfaces submitted to rainfall. However, there exist few tools that combine these processes to allow the simulation of the evolution of the structure of the soil surface. Thus, our simulator will first have a role in the aggregation of the different pieces of existing—and future—knowledge. This integration of the different processes in a model will also offer the possibility of an alternative validation of the description of the processes, by comparing simulated and observed soil surface evolution. In addition, as soil surface degradation and soil crusting control key variables such as hydraulic conductivity and mechanical resistance of the soil surface, we expect our simulator to have an important applied interest, to establish links between the state of the soil surface associated to tillage, sowing operations, and important processes such as runoff and erosion or seedling emergence. Our approach is original in several ways: (1) We have a 3D time dependent approach to soil surface structure evolution, that is we are able to analyze in time the spatial variations of this structure, both horizontally and vertically. (2) We try to add and keep information about the granulometry of the soil and its evolution. This may be useful for example to differentiate different types of crusts on the basis of their porosity, or on the fact that the material is sorted or not. (3) We use a cellular automaton framework which is homogeneous and open to different kinds of modeling, and in particular we have recourse to a fuzzy logic approach to compensate a deficiency of knowledge in certain natural processes. (4) Finally, we put an emphasis on the visualization of soil surface structure evolution, which will allow one to visually compare real and simulated dynamics. This is especially interesting because most field characterizations of soil crusting are based on visual observations and classifications. Our project is still in an early stage in its development, and some key issues have to be addressed, such as how to define in the context of the model the different types of crusts

ARTICLE IN PRESS G. Valette et al. / Computers & Graphics 30 (2006) 494–506

observed in the field. We have to complete and validate the different rules used in the transition functions of our model. We have to find a method to generate an initial soil with realistic topography and respect to the differences in aggregates shape, size and spatial organization in order to make a direct link between the state of the soil and soil tillage, seedbed preparation and sowing techniques. Nevertheless the first results obtained on disaggregation and transport by splash, despite the crude assumptions made in the description of the processes involved, are encouraging.

Acknowledgements This work is supported by the regions of ChampagneArdenne and Picardie. References [1] Kaune A, Horn R. The changes in the apparent temperature and heat conductivity of a loess due to soil structure. Mitteilungen der Deutschen Bodenkundlichen Gesellschaft 1991;66(1): 157–60. [2] Gue´rif J. Conse´quences de l’e´tat structural sur les proprie´te´s et les comportements physiques et me´caniques. In: Marin-Lafleche, JB et A, editor. La structure du sol et son e´volution: conse´quences agronomiques, maıˆ trise par l’agriculteur, INRA, Paris; 1990. p. 71–89. [3] Le Bissonnais Y, Singer MJ. Crusting, runoff, and erosion response to soil water content and successive rainfalls. Soil Science Society of America Journal 1992;56(6):1898–903. [4] Le Bissonnais Y, Bruand A. Crust micromorphology and runoff generation on silty soil materials during different seasons. Catena Supplement 1993;24:1–16. [5] Baumhardt RL, Unger PW, Dao TH. Seedbed surface geometry effects on soil crusting and seedling emergence. Agronomy Journal 2004;96(4):1112–7. [6] Tamet V, Boiffin J, Durr C, Souty N. Emergence and early growth of an epigeal seedling (Daucus Carota L.): influence of soil temperature, sowing depth, soil crusting and seed weight. Soil and Tillage Research 1996;40(1/2):25–38. [7] Hamblin AP. The influence of soil structure on water movement, crop root growth, and water uptake. Advances in Agronomy 1985;38:95–158. [8] Baranowski R, Bakowski B. The influence of differentiated soil structure on temperature dynamics. Roczniki Gleboznawcze 1977;28(1):37–44. [9] Boiffin J, Gue´rif J, Stengel P. Les processus d’e´volution de l’e´tat structural du sol: quelques exemples d’e´tudes expe´rimentales re´centes. In: Marin-Lafleche JB et A, editor. La structure du sol et son e´volution: conse´quences agronomiques, maıˆ trise par l’agriculteur. Paris: INRA; 1990. p. 37–69. [10] Toy TJ, Foster GR, Renard KG. Soil erosion: processes, prediction, measurement, and control. New York, USA: Wiley; 2002. [11] Boardman J, Favis-Mortlock D. Modelling soil erosion by water. In: Proceedings of the NATO advanced research workshop global change: modelling soil erosion by water. NATO ASI series, vol. 55, University of Oxford. Springer, Berlin, Germany; 1998. [12] Morgan R. Soil erosion and conservation. Harlow, UK: Longman; 1995. [13] Kelley AD, Malin MC, Nielson GM. Terrain simulation using a model of stream erosion. In: SIGGRAPH 1988, computer graphics proceedings, vol. 22. New York, NY, USA: ACM Press; 1988. p. 263–8.

505

[14] Roudier P, Peroche B, Perrin M. Landscapes synthesis achieved through erosion and deposition process simulation. Computer Graphics Forum 1993;12(3):375–83. [15] Nagashima K. Computer generation of eroded valley and mountain terrains. The Visual Computer 1997;13:456–64. [16] Chiba N, Muraoka K, Fujita K. An erosion model based on velocity fields for the visual simulation of mountain scenery. The Journal of Visualization and Computer Animation 1998;9:185–94. [17] Benesˇ B, Forsbach R. Visual simulation of hydraulic erosion. Journal of Winter School of Computer Graphics 2002;10:79–94. [18] Neidhold B, Wacker M, Deussen O. Interactive physically based fluid and erosion simulation. In: Galin E, Poulin P, editors. Eurographics workshop on natural phenomena, Dublin; 2005. p. 25–32. [19] Musgrave FK, Kolb CE, Mace RS. The synthesis and rendering of eroded fractal terrains. In: Proceedings of the 16th annual conference on computer graphics and interactive techniques. ACM Press; 1989. p. 41–50. [20] Shao Y. Physics and modelling of wind erosion, atmospheric and oceanographic sciences library. Dordrecht: Kluwer Academic Publishers; 2000. [21] Benesˇ B, Forsbach R. Layered data representation for visual simulation of terrain erosion. In: EEE proceedings of SCCG, vol. 25(4); 2001. p. 80–6. [22] Benesˇ B, Forsbach R. Parallel implementation of terrain erosion applied to the surface of Mars. In: Proceedings of the 1st international conference on computer graphics, virtual reality and visualisation. ACM Press; 2001. p. 53–7. [23] Boiffin J. La de´gradation structurale des couches superficielles sous l’action des pluies. The`se docteur inge´nieur, INA-PG; 1984. [24] Legout C, Legue´dois S, Le Bissonnais Y. Aggregate breakdown dynamics under rainfall compared with aggregate stability measurements. European Journal of Soil Science 2004;56(2):225–38. [25] Legout C, Legue´dois S, Le Bissonnais Y, Malam Issa O. Splash distance and size distributions for various soils. Geoderma 2005;124(3–4):279–92. [26] Valentin C, Bresson L-M. Morphology, genesis and classification of surface crusts in loamy and sandy soils. Geoderma 1992;55:225–45. [27] Bresson LM, Boiffin J. Morphological characterization of soil crust development stages on an experimental field. Geoderma 1990;47: 301–25. [28] Ben-Dor E, Goldlshleger N, Benyamini Y, Agassi M, Blumberg DG. The spectral reflectance properties of soil structural crusts in the 1.2to 2:5
mm spectral region. Soil Science Society of America Journal 2003;67(1):289–99. [29] Benesˇ B, Arriaga X. Table mountains by virtual erosion. In: Galin E, Poulin P, editors. Eurographics workshop on natural phenomena, Dublin; 2005. p. 33–40. [30] Marak I, Benesˇ B, Slavik P. Terrain erosion based on rewriting of matrices. In: Proceedings of the fifth international conference in Central Europe on computer graphics and visualization; 1997. p. 341–51. [31] Lane L, Nearing M. USDA—Water Erosion Prediction Project: Hillslope Model. Technical Report 2, USDA-ARS National Soil Erosion Research Laboratory, NSERL, West Lafayette, Indiana, USA; 1989. [32] Morgan R, Quinton J, Smith R, Govers G, Poesen J, Auerswald K, et al. The european soil erosion model (EUROSEM): a dynamic approach for predicting sediment transport from fields and small catchments. Earth Surface Processes and Landforms 1998;23(6): 527–44. [33] Favis-Mortlock D, Boardman J, Parsons A, Lascelles B. Emergence and erosion: a model for rill initiation and development. Hydrological Processes 2000;14(11–12):2173–205. [34] D’Ambrosio D, Gregorio SD, Gabriele S, Gaudio R. A cellular automata model for soil erosion by water. Physics and Chemistry of the Earth, Part B 2001;26(1):33–40.

ARTICLE IN PRESS 506

G. Valette et al. / Computers & Graphics 30 (2006) 494–506

[35] D’Ambrosio D, Gregorio SD, Iovine G. Simulating debris flows through a hexagonal cellular automata model: Sciddica S3-hex. Natural Hazards and Earth System Sciences 2003; 3:545–59. [36] Avolio MV, Crisci GM, D’Ambrosio D, Di Gregorio S, Iovine G, Rongo R, et al. An extended notion of cellular automata for surface flows modelling. WSEAS Transactions on Computers 2003;2:1080–5. [37] Chase CG. Fluvial landsculpting and the fractal dimension of topography. Geomorphology 1992;5:39–57. [38] Luo W, Duffin KL, Peronja E, Stravers JA, Henry GM. A web-based interactive landform simulation model (WILSIM). Computers and Geosciences 2004;30:215–20. [39] Haff PK. Waterbots. In: Harmon RS, Doe III WW, editors. Landscape erosion and evolution modeling. New York: Kluwer Academic/Plenum Publishers; 2001. [40] Servat D, Le´onard J, Perrier E, Treuil J. The RIVAGE-Project: a new approach for simulating runoff dynamics. In: Feyen J, Wiyo K, editors. Modelling of transport processes in soils at various scales in time and space. Wageningen, The Netherlands: Wageningen Pers; 1999. p. 592–601. [41] Wolfram S. Cellular automata and complexity: collected papers. New York: Addison-Wesley; 1994. [42] Dorsey J, Edelman A, Legakis J, Jensen HW, Pedersen HK. Modeling and rendering of weathered stone. In: Rockwood A, editor. SIGGRAPH 1999, computer graphics proceedings. Los Angeles: Addison-Wesley, Longman; 1999. p. 225–34. [43] Merillou S, Dischler J-M, Ghazanfarpour D. Corrosion: simulating and rendering. In: Graphics interface, Canadian information processing society, Toronto, Ont., Canada; 2001. p. 167–74. [44] Gobron S, Chiba N. Crack pattern simulation based on 3D surface cellular automata. The Visual Computer 2001;17:287–309. [45] Federl P, Prusinkiewicz P. Finite element model of fracture formation on growing surfaces. Lecture Notes in Computer Science 2004; 3037:138–45.

[46] Di Gregorio S, Serra R. An empirical method for modelling and simulating some complex macroscopic phenomena by cellular automata. Future Generation Computer Systems 1999;16(2–3): 259–71. [47] L’Ecuyer P. Efficient and portable combined random number generators. Communications of the ACM 1988;31(6):742–9. [48] Wischmeier W, Smith D. Rainfall energy and relationship to soil loss. Transactions—American Geophysical Union 1958;30(2):285–91. [49] Salles C, Poesen J, Sempere-Torres D. Kinetic energy of rain and its functional relationship with intensity. Journal of Hydrology 2002; 257:256–71. [50] Zadeh LA. Fuzzy sets. Information Control 1965;8:338–53. [51] McBratney A, Odeh I. Application of fuzzy sets in soil science: fuzzy logic, fuzzy measurements and fuzzy decisions. Geoderma 1997;77:85–113. [52] van Dijk AIJM, Meesters AGCA, Bruijnzeel LA. Exponential distribution theory and the interpretation of splash detachment and transport experiments. Science Society of America Journal 2002; 66:1466–74. [53] Green WH, Ampt GA. Studies on soil physics part I: the flow of air and water soils. Journal of Agricultural Science 1911;4(1):1–24. [54] Rawls WJ, Brakensiek DL, Miller N. Green–Ampt infiltration parameters from soils data. Journal of Hydraulic Engineering 1983;109(1):62–70. [55] Levoy M. Display of surfaces from volume data. IEEE Computer Graphics and Applications 1988;8(3):29–37. [56] Lorensen WE, Cline HE. Marching cubes: a high resolution 3D surface construction algorithm. In: SIGGRAPH 1987, computer graphics proceedings. New York, NY, USA: ACM Press; 1987. p. 163–9. [57] Benassarou A, Bittar E, John N, Lucas L. MC slicing for volume rendering applications. In: Computer graphics and geometric modeling. Atlanta, GA, USA: Springer; 2005. p. 314–21. [58] Horton RE. The role of infiltration in the hydrological cycle. Transaction—American Geophysical Union 1933;14:446–60.