Hongtao Hao

Christian Hartmann
email: christian.hartmann@ird.fr

Jongskul Apichart

Siltacho Siwaporn

Sukanya Promsakha

Guy Richard

Ary Bruand

Anthony R Dexter

Slumping dynamics in tilled sandy soils under natural rainfall and experimental flooding

Keywords: Recompaction, Bulk density, Matric potential, Northeast Thailand, Tillage, Overburden pressure

forces developed by water bridges between the elementary grains. Models of wet granular material indicate that a drastic loss of cohesion can occur when the liquid phase becomes continuous. Indeed slumping was observed during the downward movement of the water front inside the profile which is consistent with this hypothesis. As the existence of a continuous liquid phase along the profile is the most determinant factor of slumping, all rainfall characteristics determining the existence of such a continuous liquid phase will affect slumping, i.e. rainfall intensity, duration and frequency. For given rainfall characteristics, slumping dynamic and intensity depend on the bulk density at the onset of the rainfall or irrigation event. The high variability recorded in earlier field result from the many possible interactions between these determining factors.

Abstract

Compaction of tilled layers under the single effect of rainfall or irrigation was called slumping.

Slumping affects strongly root development and plant biomass production and was observed in different soil types but sandy soils appear particularly prone to this physical degradation. Our objectives in this study were (i) to monitor in the field the changes in soil structure and water status simultaneously, (ii) to study the effects of rainfall and management practices on slumping and (iii) to propose a conceptual model for sandy soil slumping.

An experimental site was selected in Northeast Thailand and we tested the effect of tillage depth and initial water content on slumping dynamic. Plots (9 m×15 m) were tilled at (i) two depths (20 and 40 cm, called S and D respectively) in dry conditions, (ii) at 20 cm depth in dry or wet conditions (called Y and W respectively). These plots were submitted to natural rainfall for 20 or 61 days to get different total rainfall amounts (114 and 212 mm respectively). In addition, smaller plots (0.24 m 2 ) were used for experimental flooding irrigation (equivalent to measured rainfalls, i.e. 100 and 200 mm). Soil bulk density, soil surface elevation, soil water content and matric potential were recorded.

A decrease in soil elevation was observed in all treatments. In the absence of erosion it was interpreted as a loss of porosity which resulted from slumping. Bulk density increased in all layers of the tilled profile (from 1.38 to 1.57 g cm -3 ). In the surface layer (0-5 cm) this increase was systematically higher compared to deeper layer. No significant difference in final bulk density was found between the S and W treatments, and between the Y and W treatments. Bulk density increased more rapidly in the Y and W treatments than in the S and D treatments, even though the cumulative rainfall was lower. After flooding experiments, bulk density was higher than after natural rainfall despite similar water amount brought to the soil.

The existence of a layer with 50 % porosity straight after tillage was explained by the capillary

Introduction

Soil tillage aims at creating favorable physical conditions for crop growth by modifying the soil structure and associated properties of the tilled layer [START_REF] Guérif | A review of tillage effects on crop residue management, seedbed conditions and seedling establishment[END_REF]. However, the desired loose structure tends to be structurally unstable and is susceptible to structural collapse thus leading to a loss of the benefits for root development provided by soil (Hartmann et al., 2008a;[START_REF] Hamza | Soil compaction in cropping systems: A review of the nature, causes and possible solutions[END_REF]. Soil recompaction (i.e. bulk density increase to initial level) is a general phenomenon which has several causes related to external and internal forces. External force is defined as loading from outside the soil profile, such as tractor loading, rainfall kinetic energy; internal force is defined as force inside the soil profile, such as overburden pressure, capillary force.

Soil recompaction is commonly attributed to mechanized agriculture alone because the weight of the agricultural and forestry machinery has been increased three-to fourfold as well as the frequency of wheeling during recent decades (e.g. Soane et al., 1981a, b;[START_REF] Soane | Compaction by agricultural vehicles: A review, 3: Incidence and control of compaction in crop production[END_REF][START_REF] Soane | The role of organic matter in soil compactibility: A review of some practical aspects[END_REF][START_REF] Horn | Soil compaction processes and their effects on the structure of arable soils and the environment[END_REF][START_REF] Soane | Implications of soil compaction in crop production for the quality of the environment[END_REF][START_REF] Soane | Review of management and conduct of long-term tillage studies with special reference to a 25-yr experiment on barley in Scotland[END_REF][START_REF] Hamza | Soil compaction in cropping systems: A review of the nature, causes and possible solutions[END_REF]. The phenomenon of bulk density increase recorded after wetting and without application of any external load was termed 'slumping' [START_REF] Kemper | Soil cohesion as affected by time and water content[END_REF]Mullins et al.1990). While compaction, i.e. the bulk density increase resulting from mechanical loading was widely studied, much less attention has been paid to slumping.

For coarse-textured soils, [START_REF] Kemper | Soil cohesion as affected by time and water content[END_REF] suggested that slumping is related to soils saturation or matric potential ψ close to 0 hPa. Slumping is a ubiquitous phenomenon and was recorded under a wide range of conditions. It was recorded after a single rainfall event (20 mm in one hour) [START_REF] Mead | Effect of deep tillage and seedbed preparation on the growth and yield of wheat on a hard-setting soil[END_REF], two natural rainfall events (80 mm in 160 min) [START_REF] Hartmann | Soil compaction and plant growth in Northeast Thailand: The case of the Nam Pong soil serie[END_REF], eight weeks of measuring [START_REF] Osunbitan | Tillage effects on bulk density, hydraulic conductivity and strength of a loamy sand soil in Southwestern Nigeria[END_REF], and a cropping season [START_REF] Hamblin | Interactions between soil type and tillage level in a dryland situation[END_REF]. Moreover, slumping can affect the whole tilled layer (such as down to 40 cm depth in [START_REF] Moffat | The effectiveness of cultivation using the winged time on restored sand and gravel workings[END_REF], and can have similar adverse effects on crop development as well as compaction resulting from mechanical loading [START_REF] Kozlowski | Soil compaction and growth of woody plants[END_REF]. Despite its ubiquity and importance, no study has provided yet a detailed description of the dynamics of slumping and its main field characteristics. In the context of the extension of cultivated lands to marginal areas, which often include a high proportion of sandy soils [START_REF] Eswaran | Sandy soils of Asia: A new frontier for agricultural development[END_REF], it is a challenge to look for soil management that could minimize the occurrence of slumping. As a first step in this direction, the objectives of our study were: (i) to monitor the changes in soil structure and water status simultaneously at the field scale; (ii) to study the effects of rainfall and management practices on slumping and iii) to propose a conceptual model of sandy soil slumping.

Material and Methods

Field description

Northeast Thailand is a sandy alluvial plateau covered by an aeolian deposit (approximately 1 m thick) mainly made of quartz grains. More than 80% of the grains belong to the sand fraction (50-2000 µm) and about 10% to the silt fraction (2-50 µm) [START_REF] Lesturgez | Roots of Stylosanthes hamata create macropores in the compact layer of a sandy soil[END_REF][START_REF] Bruand | Composition, fabric, and porosity of an arenic haplustalf of Northeast Thailand: Relation to penetration resistance[END_REF]. The clay fraction (<2µm) is less than 10% and contains phylosilicates (mainly kaolinite and small amount of smectite) but also quartz grains [START_REF] Bruand | Composition, fabric, and porosity of an arenic haplustalf of Northeast Thailand: Relation to penetration resistance[END_REF]). An experimental site was selected in a village named Baan Nong Sang (WGS84: 16°10' N, 102°48' E), 30 km south of the city of Khon Kaen. The soil is representative of the region with a sandy texture (< 4 % clay), low organic matter content (< 5 g kg -1 ) and high bulk density (ρ b ≥ 1.6 g cm -3 ) (Table 1 and2). This field was planted with cassava in 2006 and harvested in February 2007. To minimize the effect of coarse crop residues fragments (cassava branches and leaves, weeds, etc.) on the different soil physical characteristics, they were removed from the plot before tillage operations. To prevent weed growth, an herbicide was applied three times during the experiment.

Control of the rainfall pattern and amount

Three experiments were conducted at different periods to get different rainfall distribution and rainfall amount:

-Experiment 1 (Exp. 1): plots were tilled on 25 May, i.e. at the beginning of the rainy season that is generally characterised by storms separated by several dry days. The experiment was stopped on 1 August; -Experiment 2 (Exp. 2): plots were tilled on 6 July, i.e. in the middle of the rainy season to get a different rainfall pattern (more regular distribution of the rainfalls and thus less frequent dry spells). The experiment was also stopped on 1 August; -Experiment 3 (Exp. 3): mini-plots (0.24 m 2 ) were artificially flooded; the amount of water added was similar to the rainfall amount recorded in Exp. 1 and Exp. 2. The experiment was done between 18 to 26 July.

Rainfall was recorded daily using the average readings from two rain gauges (±1 mm precision) that installed inside the experimental field.

Experiment 1

The objective of Exp. 1 was to study the effect of tillage depth on slumping characteristics. Soil preparation treatments consisted of two tillage depths: 20 and 40 cm, called shallow (S) and deep (D) treatments respectively. Each treatment was applied to five elementary plots (9 m×15 m each). Tillage was performed using a 120-horse-power tractor equipped with a set of disk plough. In order to increase the homogeneity of the resulting soil structure, the large clods left by mechanical tillage were broken into smaller pieces by labourers using rakes. To mimic the practice of farmers who want to avoid flooding of the seed rows, ridges were built by hand two days after tillage. The interval between ridges was 40 cm and their height ranged from 12 to 17 cm.

MonitorMeasuring days were identified by the number of days after ploughing (DAP, thus DAP 0 meaning the 25 May) and by accumulated rainfall (AR in mm) since DAP 0. Soil water potential was measured using a set of ceramic tensiometers each equipped with an electronic transducer (model 2150 and SMS-2500S, SDEC Company, France). They were installed at 5 cm intervals from 10 to 20 or 40 cm depth according to tillage depth in two different locations: underneath either the ridge or the furrow. The zero level (soil surface) was counted from the surface of the ridge or the furrow according to the location of the tensiometer. Two sets of tensiometers were installed for each treatment (S and D) as replicates. Water matric potential (ψ) was measured after each major rainfall event or after 3 to 7 days if no rainfall event occurred.

Changes in soil level were measured using a horizontal frame as a stable benchmark. It consisted of a 1×1 m 2 frame horizontally put on four metallic rods inserted 70 cm below the soil surface and fixed by cement. The distance between the frame and the soil surface was measured using a laser beam (Lasermeter Leica Disto 6A, Leica Geosystem, Switzerland) according to a 5×5 cm 2 grid (total of 361 data points). The vertical precision was 1 mm. For both treatments (D and S), replicates were installed in three of the five elementary plots. To measure changes in ρ b and soil water content (Wc), after each major rainfall event, undisturbed cylinders of soil (5 cm in height and diameter) were collected every 5 cm in depth under furrows. They were weighted before and after putting in an oven for 24 h at 105°C

Experiment 2

The objective of Experiment 2 was to study the effect on slumping of the water content at tillage. Five plots (9×15 m 2 ) were left under rainfall, while five other plots were protected from rainfall using greenhouses. These greenhouses were built with bamboo sticks and covered by plastic sheets (0.8 m high); to avoid overheating and water condensation, they were opened at both ends so that air could circulate. To prevent preferential water infiltration along the plots, rainwater falling on the surface of greenhouses was drained out of the field. The greenhouses were removed on the 6 July and water content was measured after putting the soil sample in an oven for 24 h at 105°C.

The soil was tilled on the 6 July (DAP 0), at one depth only (20 cm) using the same tractor as in Exp. 1. Therefore the treatments consisted in two different Wc during tillage; the treatments were named 'dry' (Y) (plots under greenhouse before tillage) and 'wet' (W) (plots that were submitted to rainfall before tillage). Five replicates were used for each treatment. Ridges and furrows were prepared using the same procedure as in Exp. 1.

In both treatments (Y and W), a set of four tensiometers was installed below the furrow at 5 cm intervals from 10 to 25 cm depth. Changes in soil level were measured using a 1 m long horizontal board as a stable benchmark, which was supported by metallic rods as "legs" inserted into soil at 70 cm depth as in Exp. 1. One replicate of this device was installed in each of the ten subplots. The distance from the benchmark to the soil surface was measured using the laser beam every 2.5 cm (total of 38 data points). The ρ b and Wc were measured using the same procedure as in Exp. 1.

Experiment 3

The objective was to mimic a continuous rainfall and to measure the consequences on slumping characteristics. Three supplementary plots were covered by greenhouses to keep the soil dry.

Metallic rings (55 cm in diameter) were inserted deeply in the soil (60 cm) to be used as stable benchmark for soil level recording. The soil inside the ring was hand tilled using a small paddle and rake to mimic the soil tractor tillage and hand raking as in Exp. 1 and 2. Soil preparation treatments consisted of two tillage depths: 20 and 40 cm, called "shallow-tillage flooding" (Sf) and "deep-tillage flooding" (Df) respectively. Three rings were installed for each treatment. Unlike Exp.

1 and 2, the soil surface was flat (no ridges and furrows were built). Tensiometers were not used because the time was too short to obtain a good equilibrium between tensiometer and soil.

The amount of added water was 100 mm and 200 mm for Sf and Df treatments respectively.

The soil surface was covered by a straw textile and water was added by increments of 25 mm water.

To avoid kinetic energy, the water was poured on a hand as a buffer so that it gently felt on the textile over the whole surface. Infiltration time was recorded and after each addition of 25 mm water, the soil surface level was measured every 5 cm along a horizontal "board" that was put on the top of the metallic ring (10 measurement points each time) using the lasermeter. After the last addition of water and measurement of soil level, a small pit was opened inside the ring. The bulk density and water content were measured with 5 cm interval from surface to 20 and 40 cm depth in Sf and Df respectively. Three replicates samples were taken inside each ring.

Statistical analysis

The data were analyzed statistically using ANOVA by SPSS 13.0 (SPSS Inc., 2004). The least significant difference (LSD) at P = 0.05 was used to establish the significance of differences between treatment means. When applicable, paired-samples T tests were done.

Results

Rainfall characteristics

Fig. 1 presents the daily rainfall events and the accumulated rainfall (AR). During 61 days for Exp. 1, AR was 212 mm, approximately twice as much as in Exp. 2 (114 mm) which lasted only 20 days. Exp. 2 was done during the middle of the rainy season and had higher average daily rainfall than Exp. 1 which started earlier: 5.7 and 3.4 mm d -1 , respectively. More than half of AR in Exp. 2 resulted from two successive big rainfall events (22 and 36 mm) which occurred at the end of the experiment. 1 and2 3.2.1. Soil matric potential (ψ) and water content (Wc) For Exp. 1 and Y treatment (dry protected plots) in Exp. 2, Wc was ranging 0.04 to 0.06 g g -1 ; or W (unprotected wet plots) in Exp. 2, Wc was significantly higher (p<0.05), ranging 0.08 to 0.10 g g -1 (data not shown).

Experiments

Fig. 2 presents the matric potential (ψ) measured under furrows and ridges in Exp. 1 (D treatment), under furrows in Exp. 2 (W treatment). During Exp. 1, minimum and maximum matric potential (i.e. driest and wettest soil respectively) were -120 and -40 hPa respectively. No significant difference was observed between tensiometers installed under a furrow or under a ridge, and no significant difference was observed between the D and S treatments (data not presented). At a given time, we recorded similar minimum and maximum ψ in Exp. 1 and Exp. 2. Moreover, no differences were observed with the W and Y treatments compared with the D treatment in the 0-20 cm layer (data not presented).

In the top layers (0-20 cm), a rapid increase in ψ was observed after rainfall events (e.g. DAP 41 and DAP 58), but a slow and regular ψ decrease was observed during the different dry spells (e.g. from DAP 47 to 51). In the 35 cm layer, changes were buffered: the ψ was significantly affected only by the two successive big rainfall events recorded at the end of the experiment (DAP 59 and 60, 22 and 36 mm respectively).

For Exp. 1 and 2, it is noteworthy that at the end of the experiments (DAP 59 and 60), the highest ψ at 10 cm depth was recorded after the first rainfall event (-40 and -50 hPa respectively in the top layer) while higher ψ was expected after the second rainfall event (approximately twice as much water input as the first event). This discrepancy is a consequence of manual monitoring: at DAP 59, ψ was recorded shortly (< 1 h) after the rainfall event; while at DAP 60, it was recorded more than 10 h after the rainfall event. For the latter, the water front had already drained deeper in the profile as indicated by the drastic increase in ψ in the 35 cm deep layer, an increase which was not observed after the first rainfall.

At the end of the experiment, in all plots, whatever the initial Wc at the tillage depth, Wc ranged from 0.08 to 0.13 g g -1 (data not shown). There was no significant difference between D and S, and between Y and W.

Soil surface elevation and bulk density

After each major rainfall event, a decrease in soil roughness was observed: the height of ridges decreased while furrows were filled with sandy material, and consequently the field became more flat (Fig. 3). This change in the topography of the soil surface indicated movements of solid particles from the top of ridges to the bottom of furrows, but we didn't observe any erosion (particles moving outside the field). The average soil level was calculated for the different treatments (Fig. 4). A level decrease was observed in all plots (ranging from 1.2 to 1.7 cm), but no significantly difference was found between the different treatments.

Average ρ b profiles under furrows are shown in Fig. 5. Immediately after tillage, all profiles were characterised by a regular increase with depth. For shallow tillage (S, Y, W treatments, 20 cm depth), the differences between ρ b at top and bottom of tilled layer ranged from 0.05 to 0.10 g cm -3 .

For deep tillage (D treatment, 40 cm depth), this increase was higher (0.15 g cm -3 ) indicating an effect of tillage depth. When the soil was tilled in dry conditions (S, D and Y treatments), a shift towards higher ρ b values was observed at all depths with increasing AR. Initial bulk density in the profile was already higher when ploughing was performed at wet conditions. The shift was higher in the top layer then in deeper layers. When the soil was tilled in wet conditions (W treatment), no change in ρ b was observed in relation with increasing AR, except in the top layer (0-5 cm) at the end of the experiment.

Experiment 3

During the first addition, the infiltration rate was 450 and 600 mm h -1 for Sf and Df respectively (Fig. 6a). During the following additions, IR decreased until it was stabilised around 100 mm h -1 .

The same value was observed in both treatments, suggesting that tillage depth did not influence the infiltration rate. For the first three water inputs, the sinkage rate increased at each input (-1.0, -1.2, -1.4 cm in Df treatment) (Fig. 6b). For all the following water inputs it decreased continuously (from 0.6 cm with 100 mm input to 0.1 cm with 200 mm input in Df treatment). Tillage depth affected the soil level decrease: at 100 mm water input, it was only 2.7 cm for Sf compared to 4.2 cm for Df treatment. Just after tillage, ρ b profile increased with depth: at the bottom of the tilled layer ρ b was about 0.3 g cm -3 higher than at the surface (Fig. 7). At the end of the experiment, for both Df and Sf, ρ b shifted to higher values on all the profile, ρ b being higher close to the surface than at the bottom (about 0.5 and 0.2 g cm -3 respectively).

Discussion and conclusion

Evidence of slumping

The soil surface reorganisation which was recorded during Exp. 1 and 2 and which consisted in a decrease in the soil roughness (Fig. 3) was commonly observed in coarse-textured soils and discussed as closely related to their low structural stability when submitted to rainfall [START_REF] Mwendera | Effects of tillage and rainfall on soil surface roughness and properties[END_REF][START_REF] Rudolph | Effect of antecedent soil water content and rainfall regime on microrelief changes[END_REF][START_REF] Twomlow | The influence of tillage on semi-arid soil-water regimes in zimbabwe[END_REF]. Our results show that soil reorganisation affected not only surface but also the whole tilled horizon. In the three experiments and whatever the treatment, a soil sinkage was recorded, i.e. a decrease in soil surface elevation (Fig. 4 and Fig. 6b). As no erosion was observed during the experiment (i.e. no mass movement out of the field), such an elevation decrease in soil surface necessarily resulted from a loss of pore volume, i.e. from an increase in soil compactness. Soil sinkage did not happen regularly over time but was related to the occurrence of rainfall events, and consequently can be considered as an indicator of soil slumping as defined by [START_REF] Mullins | Hardsetting soils: Behavior, occurrence, and management[END_REF]. Soil sinkage cannot be observed by the naked eye even if it is probably a common phenomenon in coarse textured soils, thus explaining the rarity of its description [START_REF] Wilton | Effect of cultivations on the level of the surface of a soil[END_REF][START_REF] Young | Hardsetting and structural regeneration in two unstable British sandy loams and their influence on crop growth[END_REF][START_REF] Moffat | The effectiveness of cultivation using the winged time on restored sand and gravel workings[END_REF][START_REF] Or | Modeling post-tillage soil structural dynamics: A review[END_REF]. Even if soil sinkage is a relevant indicator of slumping occurrence, it does not provide any information about how compactness evolves with depth in the soil.

Soil structure after tillage

After tillage, the ρ b profile presents two striking aspects. First, the ρ b increased regularly with depth (from 1.26 g cm -3 at surface to 1.44 g cm -3 at 40 cm depth), thus indicating a closer packing of the soil particles when going deeper in the profile (Fig. 5 and7). Second, the ρ b values (correspond to porosity ranging from 53 to 45%) were very low in sandy soils since it results mainly from the packing of coarse grains. This loose packing was created by tillage: the energy developed by the disks was enough to separate most elementary particles which were lifted up, increasing their potential energy. This energy was dissipated when the particles moved downward under the effect of gravity. During the piling process, each time a contact occurred between neighbouring particles, the downward movement was hampered due to (i) mechanical friction, and (ii) capillary bridges developing between the grains. At Wc <0.12 g cm -3 , as measured before tillage (data not presented), water can be located only at the contact points between two particles, forming bridges (Willet et al., 2000;[START_REF] Adams | Mapping the influence of gravity on pendular liquid bridges between rigid spheres[END_REF][START_REF] Herminghaus | Dynamics of wet granular matter[END_REF]. These bridges act as glue bonding the particles together in relation with the water interfacial tension and negative Laplace pressure [START_REF] Kemper | Soil cohesion as affected by time and water content[END_REF]. The development of water bridges hampered the flow of solid particles and blocked the sand grains into a loose packing. But during the piling process, the overburden pressure was gradually increasing with depth, forcing the grains to rearrange in a continuously closer packing, resulting in a regular ρ b increase with depth. This could be a general behaviour in tilled sandy soils [START_REF] Ampoorter | Impact of mechanized logging on compaction status of sandy forest soils[END_REF]. The stability of the tilled layer results from capillary bridges and is consequently low, but in the absence of mechanical loading, it is stable until occurrence of rain and subsequent slumping.

Water characteristics inducing slumping

No relation was found between rainfall amount and slumping characteristics (occurrence, intensity, dynamic): (i) slumping was triggered by a small but variable amount of water (34 mm in Exp. 2 but only 25 mm in Exp. 3), (ii) despite similar large amounts of water brought to the soil (200 mm), the intensity of slumping was much higher in Exp. 3 than in Exp.1 (soil sinkage of 4.6 and 1.2 cm respectively), (iii) when at the end of Exp.1 and 2, AR were very different (200 mm and 100 mm respectively), intensity of slumping was similar, i.e. same final ρ b (>1.4 g cm -3 below 10 cm), and (iv) slumping was faster in Exp. 3 compared to Exp 1 (nearly finished after 75 mm water input in Exp. 3 when still observed after 132 mm water input in Exp.1) (Fig. 5).

The lack of a global relation between slumping and rainfall characteristics, suggests that not all rainfall events have the same potential to induce slumping. A soil material is stable as long as it meets the Mohr-Coulomb criterions:

τ = μσ+c (1)
where τ is the shear stress, σ is the normal force, μ is the internal friction coefficient and c is cohesion. The cohesion is equivalent to the shear stress at zero normal force and it depends on (i) both the amount and type of clay and organic matter (i.e. permanent bonds), and (ii) the water content (Le [START_REF] Bissonnais | Interactions between soil properties and moisture content in crust formation, runoff and interrill erosion from tilled loess soils[END_REF]. Sandy soils have low clay and organic carbon content, consequently water becomes a major factor of cohesion between the elementary particles and globally of structural stability [START_REF] Panayiotopoulos | Packing of sands-a review[END_REF]Mullins, 1985, Panayiotopoulos, 1989). [START_REF] Kemper | Soil cohesion as affected by time and water content[END_REF] suggested that slumping was related to soil saturation and resulting lack of cohesion between the elementary particles. Recent development on the physics of wet granular material, like sandy soils, suggest that slumping would occur before reaching water saturation.

Indeed, in a sandy material, the relation between soil cohesion and water content is not linear, two successive stages need to be considered [START_REF] Pierrat | Tensile strength of wet granular materials[END_REF][START_REF] Iveson | The dynamic strength of partially saturated powder compacts: the effect of liquid properties[END_REF][START_REF] Mitarai | Wet granular materials[END_REF]: at the first stage, water is located only at the contact between grains and is filling the packing pores between neighbouring grains creating separated clusters of grains [START_REF] Kohonen | On capillary bridges in wet granular materials[END_REF][START_REF] Scheel | Morphological clues to wet granular pile stability[END_REF]. The second stage starts after all the packing pores were filled with water when the liquid phase becomes continuous [START_REF] Herminghaus | Dynamics of wet granular matter[END_REF]. The second stage is characterised by a drastic decrease in development of air/water interfaces and thus a drastic drecrease in macroscopic cohesion (Soulié et al., 2006;[START_REF] Grof | The strength of liquid bridges in random granular materials[END_REF]. Such a model suggests that slumping in sandy soils would be triggered before reaching saturation, when rainfall amount is high enough to obtain a continuous liquid phase draining down in the profile. Compared to natural tilled sandy soil, the experimental and mathematical models studied by physicists are simplified (the particles size distribution is limited and the grains are in close packing). Even so, it is still a challenge to determine Wc or ψ that corresponds to a continuous liquid phase in these models and it

is not yet possible for real soils. Anyway, such liquid phase continuity in sandy material can be observed during fast downward movement of the wetting front that correspond to a locally high Wc and ψ [START_REF] Kawamoto | Moisture structures of laterally expanding fingering flows in sandy soils[END_REF]

Slumping dynamics and water infiltration

In Exp.3 water supply was continuous, creating an infiltration and thus a downward movement of the wetting front. Soil and water characteristics were recorded during and immediately after (<1h) water supply. During the first water input, soil surface sinkage (indicating a rearrangement of elementary grains) was observed; sinkage intensity increased until the third addition of water, then it decreased again and even if it was small, sinkage was observed until the last water input (Fig. 6b).

This result confirms that slumping occurred during the downward movement of the wetting front, i.e. before saturation, contrary to the suggestion of [START_REF] Kemper | Soil cohesion as affected by time and water content[END_REF]. Similar fast rearrangement was already recorded in an independent field experiment during which the largest ρ b increase occurred within the first 10 min of rainfall [START_REF] Bedaiwy | Soil surface densification under simulated high intensity rainfall[END_REF].

Exp. 1 and 2 were conducted under natural rainfall events and were characterised by slower slumping dynamic. Soil and water characteristics (ρ b , ψ, Wc) were monitored manually at the interval of one to seven days; consequently our device could not record any fast water front movement. Anyway, structural changes (soil surface sinkage and bulk density increased) were observed only after the limited number of major rainfall events, events that occurred similarly during Exp. 1 and 2 (Fig. 1). This observation is consistent with the hypothesis that all rainfall events do not have the same capability to induce slumping; in the context of our experiment, only events > 15 mm induced fast wetting front infiltration and significant structural changes (surface level decrease, bulk density increase).

In Exp. 3, a steady state was reached earlier for infiltration rate than for soil elevation decrease (Fig. 6a and Fig. 6b). This suggests that infiltrability is not a good indicator of slumping (direct observations of structural changes have to be preferred). On the other hand, these data can also provide some indications on the characteristics of sand grains rearrangement at the onset of the first major rainfall. The steady state observed for infiltration when slumping was still on process indicates that the continuity of the porosity (a major factor of infiltration rate) decreased faster than total pore volume [START_REF] Nimmo | Hydraulic conductivity of a sandy soil at low water-content after compaction by various methods[END_REF][START_REF] Meek | Bulk density of a sandy loam: Traffic, tillage, and irrigation-method effects[END_REF]. In a pile of grains that results from a loose packing of grains linked together only by capillary water bridges, the first grains to be unbalanced are probably the biggest ones. They are indeed the most strongly pulled by gravity and less strongly hold by capillary bridges. Once unbalanced, those grains can 'fall' in the large pores (with 50% porosity, the volume occupied by porosity is the same as occupied by particles). This reorganisation would explain the fast infiltration decrease related to (i) a decrease in the volume of large pores (i.e. that allows fast water movement), and (ii) a decrease in the continuity of the remaining smaller packing pores.

Soil factors affecting slumping

Previous experiment demonstrated that under similar rainfall events, slumping could differ in relation with soil management [START_REF] Meek | Soil compaction and its effect on alfalfa in zone production systems[END_REF][START_REF] Bedaiwy | Soil surface densification under simulated high intensity rainfall[END_REF]Rolston 1993, Hartmann et al, 2008a,b). Our experiment confirmed that the relative change of ρ b was inversely correlated with the initial ρ b as already suggested by [START_REF] Bedaiwy | Soil surface densification under simulated high intensity rainfall[END_REF]. For exemple, compared to the shallow layers, the deep layers had the highest ρ b after tillage, and they had also the smallest ρ b increase during the experiment. Less slumping was also observed when the soil was initially wet compared to when it was initially dry (Fig. 6). As the wet sandy soils are more sensitive to mechanical stresses than dry soils (Panayiatopoulos et al., 1989;[START_REF] Meek | Bulk density of a sandy loam: Traffic, tillage, and irrigation-method effects[END_REF], during the tillage operations the W plots were probably more affected by mechanical constraints induced by the disk than the Y plots. This suggests that the factor hampering the slumping process in W compared to Y treatment was its higher initial ρ b and not its higher Wc.

While inside the profile, the relative ρ b change was inversely correlated to the initial ρ b , the top layer behaved differently: ρ b increase was constant and was finally huge compared to the underlying layers. In first estimate, a different behaviour of the top layer seemed to be related to a lateral sand movement from the ridge and accumulation on the surface of the furrow (Exp.1 and 2, Fig. 3). But similar increase was recorded in Exp. 3 which had a flat surface and where no lateral movement occurred. On a flat surface, the soil structure can also be degraded by vertical segregation of clay and sand fractions due to the kinetic energy of the raindrops and result in different surface crusts [START_REF] Casenave | A runoff capability classification system based on surface features criteria in semi-arid areas of West Africa[END_REF]. Such crusts are only some millimeters thick (Bielders and Bavey, 1995;[START_REF] Roth | Bulk density of surface crusts: depths functions and relationships to texture[END_REF][START_REF] Fohrer | Changing soil and surface conditions during rainfall. Single rainstorm/subsequent rainstorms[END_REF][START_REF] Ndiaye | Effect of rainfall and tillage direction on the evolution of surface crusts, soil hydraulic properties and runoff generation for a sandy loam soil[END_REF] and cannot explain changes at several centimeters depth and deeper in the tilled layer. Indeed, the surface layer is characterised by (i) a very small normal force σ (Eq (1)) due to quasi absence of overburden pressure, and (ii) a frequent occurrence of high Wc, even during minor rainfall events. These two specific characteristics can explain that the particle reorganisation at the surface is different from that observed below the surface layer.

Since surface heterogeneity can induce preferential water infiltration [START_REF] Twomlow | The influence of tillage on semi-arid soil-water regimes in zimbabwe[END_REF], the presence of ridges and furrows has to be mentionned: slumping was perhaps also influenced by that specific surface topography, but we have not collected any data that allows a estimation on the effect of this factor .

A conceptual model of soil slumping

When [START_REF] Kemper | Soil cohesion as affected by time and water content[END_REF] and [START_REF] Mullins | Hardsetting soils: Behavior, occurrence, and management[END_REF] suggested the existence of a specific degradation process named slumping, they did not provide any information about its kinetic, and the underlying processes and determinant factors. Here we suggest a scenario of slumping in sandy soils and a list of determinant factors based on our experiment. Slumping kinetic can be separated in three main steps (Fig. 8): (i) creation of an initial loose profile with a regular ρ b increase with depth after tillage, the cohesion and stability of these layers resulting from the strong capillary forces between elementary grains (Fig. 8-a), (ii) during major rainfall events or under irrigation, a fast increase in ρ b in all tilled layers, the loss of cohesion resulting from the development of a continuous (but transient) liquid phase along the profile (Fig. 8-b), and (iii) after each major rainfall, when the soil is gradually getting denser, the slumping intensity is decreasing because of steric constraints to reorganisation, except for the surface layer where bulk density seems to be continuously increasing (Fig. 8-c). As the existence of a continuous liquid phase along the profile is the most determinant factor of slumping, all rainfall characteristic determining the existence of such continuous liquid phase will affect slumping, i.e. rainfall intensity, duration and frequency. For given rainfall characteristics, slumping dynamic and intensity depends on the bulk density at the onset of the rainfall or irrigation event. The high variability recorded in earlier field works (for example: [START_REF] Mead | Effect of deep tillage and seedbed preparation on the growth and yield of wheat on a hard-setting soil[END_REF][START_REF] Meek | Soil compaction and its effect on alfalfa in zone production systems[END_REF][START_REF] Hartmann | Soil compaction and plant growth in Northeast Thailand: The case of the Nam Pong soil serie[END_REF][START_REF] Osunbitan | Tillage effects on bulk density, hydraulic conductivity and strength of a loamy sand soil in Southwestern Nigeria[END_REF] result from the many possible interactions between these determining factors.

Further experiments should include a better control of experimental conditions associated to more accurate recording at the time scale of the rainy events to obtain new data enabling to establish a quantitative model and consequently improve our knowledge and prevision of the elementary mechanisms of slumping.

and David Lees for helpful comments and improving the English. We thank the anonymous reviewers and the editors for their very useful suggestions to improve the manuscript. Hongtao Hao thanks IRD for the financial support that made his study of PhD in France possible. He also thanks the constant support from Prof. SUN Bo and Prof. ZHAO Qiguo (Institute of Soil Science, Chinese Academy of Sciences, China). Measurements were made along a regular grid 5x5 cm (i.e. 361 measurement points). Table 1.

Legend of the figures

Fig. 1 .

 1 Fig. 1. Rainfall distribution and accumulated rainfall (AR) in Exp.1 and Exp.2. DAP is day after tillage for Exp.1, which started from 25 May 2007. DAP' is day after tillage for Exp.2, which started from 6 July 2007.

Fig. 2 .

 2 Fig. 2. Matric potential of a) Top: D (40 cm tillage) under furrow ( ) and under ridge ( ); b) Bottom: S (under furrow). Average values (n=3) are shown.

Fig. 3 .

 3 Fig.3. Soil surface level recorded in Exp.1, example of one subplot; left side is DAP4 and right side is DAP60.

Fig. 4 .

 4 Fig.4. Mean soil height with cumulated rainfall. D and S are Deep (40 cm) and S shallow (20 cm) tillage respectively; W and Y indicate tillage made in wet and dry soil respectively. Each line is one replicate from the subplot. The numbers indicate mean values of level decrease with standard errors in the brackets. Note: The reference level was the board arm that support lasermeter. The arm was fixed with iron sticks fixed deep in the soil. The absolute value of soil height was not important, but the relative soil height that measured at different time was our aim.

Fig. 6

 6 Fig.6. a) Soil infiltration rate decreased with added water amount, b) Soil surface level change with added water amount. The soil was protected from rainfall under plastic before tillage. Error bars indicate standard error of means (n=9). D-f is deep flooding; S-f is shallow flooding.

Fig. 7 .

 7 Fig. 7. Bulk density after flooding. 200 mm water was added into "D-f" and 100 mm water into "S-f". Error bar indicates standard error of mean (SEM) (n=9). D-f is deep flooding; S-f is shallow flooding.

  Fig.3

Table 2 .

 2 Selected physical and chemical properties of soil measured one day before tillage. Resistance to penetration, measured by datalogger penetrometer. SS: Shear strength, measured by Hand Held Field Vane Shear Test(Helwany, 2007). Particle size distribution of soil (g kg -1 ) at 0-20 cm layer. It was measured by laser diffraction granulometer (Malvern Instruments).

	Depth	Sand (%) Silt (%)	Clay	Organic	Bulk	RP	SS
	(cm)	50-2000	2-50	(%) <2 µm	mater (g kg -1 )	density (g cm -3 )	(kPa)	(kPa)
		µm	µm
	0-20	90	9	1	4.5		1.60	13	4.6
	20-40	86	11	3	4.6		1.65	33	5.6
	40-60	86	11	3	1.9		1.62	32	5.1
	RP: Clay	Silt	Sand
	Size (μm)	0-2	2-50	50-100 100-150 150-200 200-250 250-500 500-2000
	Content	26	101	205	176	209	171	112	3

Acknowledgements

The authors thank the support from Ms Wanpen Wiriyakitnateekul and Ms Nopmanee Suvannang (LDD, Land Development Department, Thailand); and the help from Mr. Jean-Luc Maeght, Dr Claude Hammecker, Dr Christian Valentin (IRD, Institut de Recherche pour le Développement, France); Ms. Tum Nitjaporn Kunklang, Mr. Worraphan Chintachao. We thank Frédéric Darboux


Yohan D Fougerolle
email: yohan.fougerolle@u-bourgogne.fr

Frédéric Truchetet
email: frederic.truchetet@u-bourgogne.fr

Johan Gielis

A NEW POTENTIAL FUNCTION FOR SELF INTERSECTING GIELIS CURVES WITH RATIONAL SYMMETRIES

Keywords: Gielis curves and surfaces, implicit functions, parametric functions, R-functions, superquadrics, symmetry

We present a new potential field equation for self-intersecting Gielis curves with rational rotational symmetries. In the literature, potential field equations for these curves, and their extensions to surfaces, impose the rotational symmetries to be integers in order to guarantee the unicity of the intersection between the curve/surface and any ray starting from its center. Although the representation with natural symmetries has been applied to mechanical parts modeling and reconstruction, the lack of a potential function for Rational symmetry Gielis Curves (RGC) remains a major problem for natural object representation, such as flowers and phyllotaxis. We overcome this problem by combining the potential values associated with the multiple intersections using R-functions. With this technique, several differentiable potential fields can be defined for RGCs. Especially, by performing N-ary R-conjunction or R-disjunction, two specific potential fields can be generated: one corresponding to the inner curve, that is the curve inscribed within the whole curve, and the outer -or envelope-that is the curve from which self intersections have been removed.

INTRODUCTION

Describing and modeling nature is fascinating and, generally speaking, one of the most fundamental research activities: whether to model physical behaviors or geometric structures, to describe or to recognize natural shapes, every research community aims at representing nature as accurately as possible. Classical models are largely based on isotropic spaces with the Euclidean circle as (isotropic) unit circle. In nature however, anisotropy is the rule and different ways of measuring or geometrizing exist. In 2003, halfway between the fields of botany and computer graphics, Gielis et al. introduced the superformula [START_REF] Gielis | A generic geometric transformation that unifies a wide range of natural and abstract shapes[END_REF][START_REF] Gielis | Superquadrics with rational and irrational symmetry[END_REF], which can be seen as a parametric formulation for generalized circles or ellipses. Superellipses defining anisotropic unit circles led to notion of Minkowski distances and Minkowski geometry [START_REF] Thompson | Minkowski geometry, volume 63 of Encyclopedia of Mathematics and its applications[END_REF]. Superellipses have been extended to superquadrics in (Barr, 1981), which have found numerous applications due to the limited number of shape parameters and their ability to represent objects ranging from diamonds, cubes, spheres, and any intermediate shape. More interestingly, the superformula is now spreading to other fields. For instance, it has been used in recent work to study constant mean curvature surfaces for anisotropic energies [START_REF] Koiso | Equilibria for anisotropic surface energies and the Gielis formula[END_REF], in clustering and data mining in [START_REF] Morales | Clustering with an n-dimensional extension of Gielis superformula[END_REF], and in fluid dynamics in [START_REF] Wang | Investigation of trajectories of inviscid fluid particles in two-dimensional rotating boxes[END_REF]. Recently, Natalini et al. have presented a numerical algorithm to write down the explicit solution to the Dirichlet problem for the Laplace equation in a starlike domain [START_REF] Natalini | The Dirichlet problem for the Laplace equation in a starlike domain of a Riemann surface[END_REF], and presented closed form equations for various Gielis curves.

In 2005, based on the parametric formula proposed in [START_REF] Gielis | A generic geometric transformation that unifies a wide range of natural and abstract shapes[END_REF], potential fields for 3D Gielis surfaces have been proposed in [START_REF] Fougerolle | Boolean operations with implicit and parametric representation of primitives using R-functions[END_REF]. This representation has found applications for Gielis surface recovery for mechanical parts (Bokhabrine et al., 2007). Unfortunately, one of its major weakness is that the implicit field equations require the rotational symmetries to be integers. Such a restriction can be tolerated when manufactured objects are represented. Unfortunately, as initially remarked in [START_REF] Gielis | A generic geometric transformation that unifies a wide range of natural and abstract shapes[END_REF], some natural objects require the rotational symmetries to be rational numbers. Such symmetries are found in the phyllotaxy of plants. Leaves are arranged in a helical or spiral way around the stem. In rose, leaves are arranged around the stem in a 5/2 arrangement, meaning that the sixth leave will be precisely above the first one and the spacing between leaves is 144 • . In wild roses the five petals are arranged in a plane, but the sepals in the preceding (almost planar) whorl still show the helical arrangement as they are still 144 • apart. This arrangement can be seen using m = 5/2 in the superformula. It involves fusion of certain parts, while in the center an open structure is created, giving rise to the rose hip, as illustrated in Figure 1. Non-integer symmetries can be observed in biomolecules as well, e.g. DNA and proteins in which non-integer symmetries are observed frequently [START_REF] Janner | DNA enclosing forms from scaled growth forms of snow crystals[END_REF][START_REF] Janner | Strongly correlated structure of axialsymmetric proteins. ii. pentagonal, heptagonal, octagonal, nonagonal and ondecagonal symmetries[END_REF]. Complex objects that are defined as Boolean operations between multiple globally deformed Gielis surfaces can be modeled and reconstructed [START_REF] Fougerolle | Boolean operations with implicit and parametric representation of primitives using R-functions[END_REF]Bokhabrine et al., 2007). To transcribe the Boolean predicates between 3D Gielis surfaces into analytical equations, R-functions have been employed. The strategy adopted in this paper uses the same tools and concepts, i.e. we use radial distance functions to build a 2D potential field, and Rfunctions for the transcription of Boolean combinations into analytical equations. The difference is that now we perform R-functions not to combine implicit fields of several Gielis curves or surfaces, but the multiple implicit values of the same rational Gielis curve. Thus, we overcome the self-intersection issue through an auto-R-function operation and build 2D potential fields equation for RGCs, that can represent for instance the envelope of the curve or its "core". We present several advantages of this representation, ranging from flower modeling, from global shape to petals, and its potential for further research directions, such as parameter recovery and/or optimization.

The structure of the rest of paper is as follows: in section 2 we recall the initial parametric definition of Gielis curves, surfaces, and their associated potential fields. In section 3 we briefly present R-functions. Using R-function and the initial Gielis formula, two potential fields are presented for 2D RGCs in section 4. Several strategies about extension to Rational Gielis Surfaces are presented and discussed in section 5. We then present our future work and conclusions in section 6.

GIELIS CURVES

In polar coordinates, the radius r(θ) of a Gielis curve is defined by:

r (φ) = 1 n 1 1 a cos mφ 4 n 2 + 1 b sin mφ 4 n 3 , (1)
with n i ∈ R + , and a, b, and m ∈ R + * . Parameters a and b control the scale, m represents the number of rotational symmetries, n 1 , n 2 , and n 3 are the shape coefficients. Regular polygons and superellipses can be generated by setting the shape coefficients to specific values as shown in [START_REF] Gielis | A generic geometric transformation that unifies a wide range of natural and abstract shapes[END_REF].

Gielis only proposed the parametric formulation for 2D curves. In the case of closed non selfintersecting curve (m is positive integer), for a 2D point P(x, y) one can define the following potential field:

F 1 (x, y) = 1 - -→ OP -→ OI = 1 - x 2 + y 2 r 2 (θ) . ( 2
)
O is the center of the curve, and the point I = r (θ(x, y)) corresponds to the intersection between the curve and the half line [OP). If the symmetry parameter m is an integer, the intersection I is unique. If the curve is closed, the sign of the potential field F(x, y) generated by equation 2 can be used to define a partition of the 2D space. In this case, the set of points where F(x, y) is positive corresponds to the inside of the Gielis curve, the set of points where F(x, y) is negative corresponds to its outside, and the curve corresponds to the zero-set of the potential field.

In [START_REF] Fougerolle | Boolean operations with implicit and parametric representation of primitives using R-functions[END_REF], by setting a = b = 1 in equation 1 and by considering a 3D Gielis surface as the spherical product of two 2D Gielis curves, potential fields for 3D non self-intersecting Gielis surfaces have been introduced as:

F 2 (x, y, z) = 1 - 1 r 2 (φ) x 2 + y 2 + z 2 cos 2 φ r 2 1 (θ) -1 + 1 . (3)
Such representation of Gielis surfaces as iso-values of a potential field is crucial for Gielis surface reconstruction from 3D data, because it is used to build various cost functions to be optimized.

R-FUNCTIONS

R-functions find their origin in geometric algebra based on logic and have been introduced by Vladimir Logvinovich Rvachev in [START_REF] Rvachev | Geometric Applications of Logic Algebra[END_REF]. Since their introduction, R-functions have found direct applications in several fields, such as geometric modeling and boundary value problems. For conciseness purpose, we briefly present the most common Rfunctions in this section. The reader is invited to refer to the recent survey by Shapiro [START_REF] Shapiro | Semi analytic geometry with Rfunctions[END_REF] for more in depth presentation of R-functions and their applications.

The simplest R-function is R α (x 1 , x 2 ) and is defined by:

R α (x 1 , x 2 ) = 1 1 + α x 1 + x 2 ± x 2 1 + x 2 2 -2αx 1 x 2 , ( 4
) where α(x 1 , x 2 ) is an arbitrary symmetric function such that -1 < α(x 1 , x 2 ) ≤ 1. Setting α to 1 leads to the simplest and most popular R-functions: R-conjunction min(x 1 , x 2 ) and R-disjunction max(x 1 , x 2 ). Other useful R-functions with differential and normalization properties, namely R m 0 and R p , are studied in detail in [START_REF] Shapiro | Implicit functions with guaranteed differential properties[END_REF], and are respectively defined by:

R m 0 (x 1 , x 2 ) = x 1 + x 2 ± x 2 1 + x 2 2 x 2 1 + x 2 2 m 2 ,
(5) where m is any even positive integer and

R p (x 1 , x 2 ) = x 1 + x 2 ± x p 1 + x p 2 1 p , ( 6
)
for any even positive integer p. R α , R m 0 , and R p functions only handle two arguments. Rvachev introduced the N-ary R-conjunction and R-disjunction to handle more than two arguments, which are less restrictive than R α , R m 0 , and R p and more appropriate to RGCs. The parameter m is an integer and corresponds to the parameter used for R m functions.

i=n i=1 (m) x i ≡ n ∑ i=1 (-1) m x m i (x i -|x i |) + n ∏ i=1 x m i (x i + |x i |) . ( 7a
) i=n i=1 (m) x i ≡ n ∑ i=1 x m i (x i + |x i |) - n ∏ i=1 x m i (-1) m (|x i | -x i ) .
(7b) If no specific constraint about differentiability is required, m can be set to zero, which leads to the simplified version of N-ary R-conjunction and Rdisjunction defined in equations 8a and 8b.

i=n i=1 x i ≡ n ∑ i=1 (x i -|x i |) + n ∏ i=1 (x i + |x i |) . ( 8a
) (a) (b) (c) (d)
i=n i=1 x i ≡ n ∑ i=1 (x i + |x i |) - n ∏ i=1 (|x i | -x i ) . ( 8b
)
An R-function is a real-valued function characterized by some property that is completely determined by the corresponding property of its arguments. More specifically, the R-functions presented in this paper have the property that their sign is completely determined by the signs of their arguments. In the following section, we present how to use this property to build a signed potential field for RGCs.

POTENTIAL FIELD EQUATION FOR RATIONAL GIELIS CURVES

While RGC curves with multiple crossings and period k (with k even) can be studied in a XY graph without intersections, it is worthwhile to develop a potential field function which does take the intersections into account and in which multiple function values do occur. As observed in rose, such intersections could indeed give rise to a center which can develop into the rose hip. Alternatively, the multiple intersections can lead to the separation of specific, isolated sectors in which separate developmental processes may occur. Such sectors are indeed observed using R-conjunction of outer envelope and the complementary of the core (Figure 5). At the bud stage, flower development involves the sequential initiation of various whorls in well-defined yet separated sectors. In Greyia sutherlandii (Figures 2(b), 2(c), and 2(d)) for example, after the formation of sepals in the 'corners' petals are formed and inbetween the formation of stamens occurs.

We come back to the superformula initial formulation and introduce rational symmetry: the symmetry parameter m is no longer an integer and can be rational, and a = b = 1. By definition, RGCs are self inter- secting curves, i.e. the symmetry parameter m can be written as the ratio of two integers as m = p/q. The parametric formulation of a RGC is then written as:

r (φ) = 1 n 1 cos pφ 4q n 2 + sin pφ 4q n 3 , (9)
with p, q ∈ N + * . The parameter p is similar to m, i.e. still represents the rotational symmetry number. The parameter q corresponds to the maximum number of self intersections, and the angle θ now belongs to [0, 2qπ]. The condition for a RGC to be closed is to verify that the radius r(θ) is identical for angles θ = 0 and θ = 2qπ, which can be written as r(0) = r(2qπ).

Figure 3 shows an example of a RGC with p = 8 and q = 3. For a given RGC, there exist multiple intersections between a ray originated from its center and the curve, i.e. there exist multiple polar coordinates in [0, 2qπ] that generate points lying both on the curve and the ray, as illustrated in Figure 3, where three intersections, noted I 1 , I 2 and I 3 are detected. The number of intersections depends on q and double intersection points appear for angles θ = kp/q. For a given point P(x, y), a first angle θ 0 within [0, 2π] can be easily determined. The corresponding intersection point I 0 can be computed as r(θ 0 ). Other angles θ k and corresponding intersection points I k can easily be determined using θ k = θ 0 + 2kπ and I k = r(θ k ). We see that the several intersections can be simply computed and correspond to the multiple values of the radius for θ 0 + 2kπ values. For each intersection point I k , equation 2 can be applied to associate a potential. Now, the last problem to overcome is to build a continuous potential field from these k individual potentials. One of the simplest idea is to consider the maximum or minimum potential determined for intersection points I k using equation 4 with α = 1. This approach generates a potential field that is not smooth everywhere, as illustrated in Figures 4(b),4(d), 4(f), 4(h), ??, and ??, that have been obtained using min and max R-functions. Therefore, such a technique suffers a severe drawback, especially for reconstruction purposes, where differentiable functions are often preferred. To obtain a differentiable potential field, we can combine each potential using R-functions presented in equations 5 or 6. Using R p may be an advantage if normalization property is desired, but it also has a major drawback. With R p -functions, the generated potential field depends on the order the Rfunctions. Indeed, for instance, it is easy to verify that (x 1 ∧ x 2 ) ∧ x 3 = x 1 ∧ (x 2 ∧ x 3 ), except over its zero-set. Unicity of the potential field can be obtained using two techniques: sort individual values or restrict the q parameter to be set to 2 to keep binary R-functions, which is not satisfactory in both cases. Fortunately, for multiple self intersections, n potential values can be combined using n-ary R-functions as defined in equation 7. The justification for the introduction of equations 8a and 8b now clearly appears: using binary R-functions, such as R α , R p or R m , restricts the parameter q in equation 9 to be equal to 2, whereas with N-ary R-functions q ∈ N + * . The commutativity of N-ary R-functions is obvious since these functions are a sum of sums and products that are commutative.

Eventually, by replacing the arguments of equations 8a and 8b by the potential field defined in equation 2, we obtain the inner and outer potential fields for 2D-RGCs, respectively defined by:

i=q i=1 F 1 (x i , y i ) and i=q i=1 F 1 (x i , y i ), (10)
where F 1 (x i , y i ) corresponds to the potentials evaluated for the multiple intersection points (x i , y i ) using 2 combined through N-ary R-conjunction and R-disjunction, respectively. Figure 4 illustrates the relative intensity of the potential field for various RGCs and several R-functions.

Additionally, the inner and outer potential fields can be combined together, for example to represent developmental processes in flowers. For instance, one can want to perform the equivalent of the Boolean difference between the inner and outer curve, simply by using binary R-functions, as illustrated in Figure 5. Such a representation finds direct biological meaning: when a plant goes into flowering, the flower approximately becomes a planar structure, going from a spiral phyllotaxy (given by Fibonacci numbers) into a planar phyllotaxy or still with helical or spiral tendencies. But in a meristem of the flower, new entities have to develop sequentially (first sepals, followed by petals, stamen and finally pistil). This means that in a planar arrangement a clear separation of areas/sectors is necessary, as delineation of potential fields. And that is precisely what distinguishes Figures 4 and5. This illustrates the potential of our approach to build from bottom up (pure R-functions/logic, with RGC) a way of modeling flowers. This example also illustrates the strength of our approach, thanks to R-functions, because such sectors may not be represented using Natalini's approach.

(a) p = 5, q = 2 (b) p = 5, q = 3 (c) p = 5, q = 4 (d) p = 5, q = 2 (e) p = 5, q = 3 (f) p = 5, q = 4 Figure 5: R-conjunction between a RGC outer envelope and the complementary of its inner envelope to represent various petal shapes. First row: using R p -function. Second row: using R α .

EXTENSION TO RATIONAL GIELIS SURFACES

We present two approaches to extend RGCs to Rational Gielis Surfaces (RGSs). The first one is using spherical product, and the second is building revolution surfaces using a RGC and a profile. The ob- jective of 3D extension through spherical product is to generate a closed surface built upon two Gielis curves, which can be used in solid modeling for example. This operation requires one important restriction: the second generating curve cannot be rational. The main problems arising when considering a RGC as second generating curve are surface generation/continuity and potential field evaluation. In spherical product, half of the second curve, that corresponds to angles φ ∈ [-π 2 , π 2 ], is rotated around the rotation axis of the first generating curve. In this case, the 2qπ modulus to evaluate every intersection leads to curve discontinuities, which is a major problem for efficient and simple tesselation algorithms and potential field evaluation. Thus, one solution to build closed surfaces while preserving surface continuity and signed implicit field definition, is to perform the spherical product between one RGC with one GC. The implicit field generated can then be evaluated by:

i=q i=1 F 2 (x i , y i ) and i=q i=1 F 2 (x i , y i ), (11)
where F 2 (x, y) is the potential equation for 3D Gielis Surfaces presented in equation 3. An example of a mesh corresponding to such surface is illustrated in Figure 6(a).

The second approach consists in defining a profile, by any known techniques such as NURBS for instance, and to use this profile as an elevation profile. In this case, the surface created does not define a closed object, which makes impossible the expression of an implicit field. Nevertheless, such representation may be very useful for compact and efficient representation of elementary 3D flower patterns, especially in entertainment or video-game industry. An example of a revolution surface using RGC, using a cubic polynomial profile, is presented in 6(b).

We have presented new potential functions for closed Rational Gielis Curves and possible strategies for their extension to 3D surfaces. Our approach makes use of R-functions to overcome the self-intersections issue introduced by the rational rotational symmetry. With this technique, several differentiable potential fields can be defined for RGCs and, more specifically, one corresponding to the inner curve, that is the curve inscribed within the whole curve, and one corresponding to the envelope. Such representation offers promising perspectives, especially in botany with classification and morphology metrics: as illustrated in Figure 5, combining inner and outer potential fields through R-functions leads to the definition of sectors that are directly related to the flower developmental process.

Among the numerous other research directions, we consider the study of other possible 3D extensions and their applications to solid modeling, boundary value problems, and/or entertainment (fast 3D flower modeling and rendering, procedural flower field texture generation, etc). Another still highly challenging research concerns the shape and symmetry parameters recovery. To our knowledge, in the literature, there still exist very few papers dedicated to Gielis curves parameters recovery using integer symmetries, and none considering rational symmetries. The introduction of potential equations for such closed curves therefore opens new research perspectives in this field. Moreover, due to the complexity of the space parameters, deterministic methods, such as Levenberg-Marquardt method, can only be applied in restricted cases, with prior symmetry detections and strong assumptions. Our current and future works include the development of more suitable RGC potential functions for optimization processes combined with the study of appropriate stochastic algorithms for efficient RGCs parameters recovery and their application to classification, pattern recognition, and object segmentation both in 2D and 3D.

  Figure 1: a) Rose sepals. b) RGC with m = 5/2 and n 1 = n 2 = n 3 = 0.43.

Figure 2 :

 2 Figure 2: Examples of developing flower buds. a) Ochna atropurpurea: Development of ovary. b) Greyia sutherlandii 1: sepals protect the newly formed petals. c) Greyia sutherlandii stage 2: stamens develop with petals in the five corners of the pentagon (sepals removed). d) Greyia sutherlandii stage 3: stamens and pistils develop.

Figure 3 :

 3 Figure 3: Intersections between a RGC with m = 8/3 and a half ray.

Figure 4 :

 4 Figure 4: Color coding of relative intensity of the potential field generated by auto R-function of a RGC with a = b = 1, n 1 = 0.5, n 2 = n 3 = 3.5, p = 5. The RGC is in dark red. First row: m = 5/2. Second row m = 5/3. From left to right: relative potential field intensity using the N-ary R-disjunction, Maximum, N-ary R-conjunction, and Minimum.

Figure 6 :

 6 Figure 6: Rational Gielis Surfaces. a) Spherical product. b) Revolution surface.