Jacques Bibaï

Pierre Savéant

Marc Schoenauer
email: marc.schoenauer@inria.fr

Vincent Vidal
email: vidal@cril.univ-artois.fr

Learning Divide-and-Evolve Parameter Configurations with Racing

. For a given domain, the learning phase uses a racing procedure to choose the rates of the different variation operators used in DAE, and processes the results obtained during this process to specify the predicates that will be later used to describe the intermediate goals.

Introduction

Divide-and-Evolve (DAE) is a generic hybrid approach to solve Temporal Planning Problems (TPPs), originally introduced in [START_REF] Kaufmann | Divide-and-Evolve: a New Memetic Scheme for Domain-Independent Temporal Planning[END_REF]2007). It uses an evolutionary algorithm to evolve an ordered sequence of subgoals; the resulting TPPs (going from one subgoal to the next) is then passed on to an embedded planner; if all those TPPs are solved, the concatenation of all corresponding subplans (after some compression step) is a solution of the initial TPP. The makespan (or number of actions) of this solution defines the fitness of the sequence of subgoals used by the evolutionary algorithm.

A general issue in Evolutionary Computation (EC), that somewhat hinders its wide use in spite of some highly successful applications, lies in the number of parameters the programmer has to tune (from population size to selection operators to rates of applications of variation operators), and the lack of theoretical guidance to help him. Experimental statistical procedures have been proposed, that build on standard Design of Experiments methods and use the specificities of the EC domain to reduce the amount of computation. Among those, the racing approach [START_REF] Yuan | Statistical Racing Techniques for Improved Empirical Evaluation of Evolutionary Algorithms[END_REF]) has been chosen here, and is used to learn, for a given domain, the best rates of application of variation operators.

However, another issue for DAE lies in the choice of the atoms that are used to describe each subgoal: the goal of a TPP can be (and usually is) described only partially, i.e. by requiring only the value of a small number of atoms (instanciated predicates from the domain). Furthermore, searching the space of complete states would result in a rapid explosion of the size of the search space. It thus seems more practical to search only sequences of partial states, and to limit the choice of possible atoms used within such partial states. However, previous experiments on different domains of TPP from the IPC benchmark series [START_REF] Bibai | Divide-And-Evolve Facing State-of-the-art Temporal Planners during the 6 th International Planning Competition[END_REF] has demonstrated the need for a very careful choice of such atoms. The approach proposed here builds on the results of the racing procedure to select a few atoms that will be later used to construct sequences of partial goals for their evolutionary optimization.

Next section briefly introduces the Divide-and-Evolve planner and details the different components of the evolutionary algorithm, as well as the way it interacts with the embedded planner to compute the fitness of potential solutions. After having described the learning procedure, the last section presents and discusses preliminary results of DAE using this learning procedure on the IPC-6 learning benchmarks.

The Divide-and-Evolve Planner

In order to solve a planning problem P D (I, G) (D for the domain), the basic idea of DAE is to find a sequence of states S 1 , . . . , S n , and to use some embedded planner to solve the series of planning problems P D (S k , S k+1 ), for k ∈ [0, n] (with the convention that S 0 = I and S n+1 = G). The generation and optimization of the sequence of states (S i ) is driven by an evolutionary algorithm, and we will now describe its main components: representation, variation operators, and fitness.

Representation

An individual, possible solution of the TPP at hand, is a (variable length) list of subgoals, or partial states of the given domain. In STRIPS representation model [START_REF] Fikes | STRIPS: A New Approach to the Application of Theorem Proving to Problem Solving[END_REF], a state is a list of boolean atoms. However, searching the space of complete states would result in a very fast combinatorial explosion of the search space size. Moreover, goals of TPP need only to be defined as partial states. It thus seemed practical to search only sequences of partial states. However, this raises the issue of the choice of the atoms to be used to represent individuals, among all possible atoms.

The result of the experiments in many kinds of TPP on the IPC benchmarks [START_REF] Bibai | Divide-And-Evolve Facing State-of-the-art Temporal Planners during the 6 th International Planning Competition[END_REF] demonstrates the need for a very careful choice of the atoms that are used to build the partial states. This lead to propose a new method to build the partial states, based on the earliest time from which an atom can become true. Such date can be estimated by a classical heuristic function (e.g h 1 , h 2 ... [START_REF] Haslum | Admissible Heuristics for Optimal Planning[END_REF]). The dates from which all atoms become true can then be discretized into the number of restrictions of the set of all possibly atoms, and a partial state is built at each date by choosing randomly among several atoms that are possible true at this date. The sequence of states is then built by preserving the estimated chronology between atoms.

Although these restrictions can contain a large number of atoms, they can be reduced by choosing to keep only atoms built with a set of allowed predicates. We expect that this set can be learned by analyzing several optimized sequences of states given by the algorithm on several problems of the same domain.

Nevertheless, even when restricted to specific choices of atoms, the choice of random atoms can lead to inconsistent partial states, because some sets of atoms can be mutually exclusive (mutex in short). Whereas it could be possible to allow mutex atoms in the partial states generated by DAE, and to let evolution discard them, it seems more efficient to a priori forbid them, as much as possible. In practice, it is difficult to decide if two atoms are mutex. Nevertheless, it can be estimated with h 2 heuristic function [START_REF] Haslum | Admissible Heuristics for Optimal Planning[END_REF] in order to build mutex-free states.

Initialization and Variation Operators

The initialization phase and the variation operators of the DAE algorithm respectively build the initial sequences of states and randomly modify some sequences during its evolutionary run.

The initialization of an individual is the following: first, the number of states is uniformly drawn between one and the number of estimated dates; For every chosen date, the number of atoms per state is chosen uniformly between 1 and the number of atoms of the corresponding restriction. Atoms are then chosen one by one, uniformly in the allowed set of atoms, and added to the individual if not mutex with any other atom already there.

A 1-point crossover is used, adapted to variable-length representation in that both crossover points are uniformly independently chosen in both parents.

Because an individual is a variable length list of states, and a state is a variable length list of atoms, a mutation operator can act here at two levels: at the individual level by adding (addStation) or removing (delStation) a state; or at the state level by changing (addAtom) or removing (delAtom) some atoms in the given state.

Note that the initialization process and these variation operators maintain the estimated chronology between atoms in a sequence of states and the local consistency of a state, i.e. estimated mutual exclusion relations between atoms.

Applying Variation Operators Several parameters control the application of the variation operators. During an evolutionary run, two parents are chosen according to the selection procedure. With probability p cross , they are recombined using the crossover operator. Each one then undergoes mutation with probability p mut . When an individual must undergo mutation, four additional user-defined relative weights (w addStation , w delStation , w addAtom , w delAtom ) are used to choose among the four mutation operators defined above: each operator has a probability proportional to its weight of being applied. At most one mutation operator is thus applied to each individual.

Fitness

The fitness of a list of partial states S 1 , . . . , S n is computed by repeatedly calling a embedded planner to solve the sequence of problems P D (S k , S k+1 ) (k = 0, . . . , n). Any existing planner could be used here, in this paper DAE uses YAHSP [START_REF] Vidal | A Lookahead Strategy for Heuristic Search Planning[END_REF]), a lookahead strategy planning system for non-optimal STRIPS planning which uses the actions in the relaxed plan to compute reachable states in order to speed up the search process.

For any given k, if the chosen embedded planner succeeds in solving P D (S k , S k+1 ), the final complete state is computed, and becomes the initial state of the next problem: initial states need to be completed (and denoting each problem as P D (S k , S k+1 ) is indeed an abusive notation). If all problems, P D (S k , S k+1 ) are solved by the chosen embedded planner, the individual is called feasible, and the concatenation of all solutions plans for all P D (S k , S k+1 ) is a global solution plan for P D (S 0 = I, S n+1 = G). However, this plan can in general be optimised by parallelising some of its actions, in a step call compression (see [START_REF] Schoenauer | Divideand-Evolve: a Sequential Hybridization Strategy using Evolutionary Algorithms[END_REF] for detailed discussion). The fitness of a feasible individual is the makespan (or the number of actions) of the compressed plan.

However, when the chosen embedded planner fails to solve one P D (S k , S k+1 ) problem, the following problem P D (S k+1 , S k+2 ) cannot be even tackled by the chosen embedded planner, as its complete initial state is in fact unknown, and no makespan can be given to that individual. All such plans receive a fixed penalty cost such that the fitness of any infeasible individual is higher than that of any feasible individual. In order to nevertheless give some selection pressure toward feasible individuals, the relative rank of the first problem that the chosen embedded planner fails to solve is added to the fixed penalty, so infeasible individuals which solve the more subproblems are favoured by selection.

Finally, because the initial population contains randomly generated individuals, some of them might contain some subproblems that are in fact more difficult than the original global problems. It was necessary to limit the chosen embedded planner by adding some constraints in order to discard those subproblems. And because, ultimately, it is hoped that all subproblems will be easy to solve, such limitation should not harm the search for solutions.

We have constrained YAHSP with a maximal number of nodes that it is allowed to use to solve any of the subproblems. We fixed the maximum number of nodes in two steps: first to solve the initial population (at the initialization) we allowed a large number of node (e.g. 100000); and in the second step, for the rest of the algorithm, we allowed the median of nodes used in the solutions found at the initialization.

Parameter Learning

Because there are too many parameters to tune, some of them were fixed after preliminary experiments reported in the previous work [START_REF] Kaufmann | Divide-and-Evolve: a New Memetic Scheme for Domain-Independent Temporal Planning[END_REF]2007), i.e. evolution strategy and stopping criteria. So the two steps learning process only involves choosing the probability and weights of each of the variation operators being used (the best domain-dependent search strategy) and choosing predicates (the representation domain-dependent knowledge) for the intermediate goals. It first finds the best domain-dependent search strategy and then for this best domain-dependent search strategy finds the best set of allowed predicates.

The Best Domain-Dependent Search Strategy

The naive way of tuning the parameters of evolutionary algorithms in order to solve a class of problems is to try exhaustively all the parameter configurations over problems and, by means of statistical analysis over the results, the best configuration is extracted. The problem of doing so is that a lot of computational processing time might be thrown away while extensively evaluating very bad candidates. Originally proposed for solving the model selection problem in Machine Learning [START_REF] Maron | Hoeffding Races: Accelerating Model Selection Search for Classification and Function Approximation[END_REF], racing technique was introduced [START_REF] Birattari | A racing algorithm for configuring metaheuristics[END_REF] in order to focus the search in the most performing parameter configurations.

The general idea is that, while performing all runs on each parameter configuration, as soon as there is enough statiticals evidence that the current parameter configuration is worst than the best parameter configuration found so far, there is no reason to keep performing experiments with such configuration and so it can be discarded. Such cycle execution-comparison-elimination is repeated until there is just one parameter configuration left or if the maximum number of allowed experiments has been reached.

However, the efficiency of such technique totally depends on the selection of the statistical test to be applied in the comparison. Because no assumption can be made about the distribution of the results (e.g. normality), we have chosen to use the nonparametric Friedman's two-way analysis of variances by ranks. The application of racing using the Friedman's test, the so-called F-RACE, was deeply examined in [START_REF] Birattari | A racing algorithm for configuring metaheuristics[END_REF][START_REF] Yuan | Statistical Racing Techniques for Improved Empirical Evaluation of Evolutionary Algorithms[END_REF]2007).

Moreover, two parameters are inherent to whatever statistical test chosen: the number of initial runs before starting the comparison; and the confidence level. In order to choose the best parameter configuration, we make 11 runs (the lowest significant number for the statitical test) before starting the comparison for each problem tested during the learning and parameter configuration. We used 0.025 confidence level (strong constraint for the acceptation of equality hypothesis between two parameter configurations) to select the best set of parameters in terms of the lowest number of actions and the lowest execution time.

Because the test of a parameter configuration can take too much time, we decided to choose only some problems for the racing. These problems were selected by using YAHSP: we tried to solve each bootstrap (example) problem with YAHSP with a limited time (e.g. 10 minutes). After this step we selected two problems for the learning procedure, one of which has been solved by YAHSP (unless all problems are unsolved). The solved problem chosen was that with the longest execution time; the unsolved problem was that with the lowest memory.

In order to reduce the space of the variation operator parameters we selected twenty sets of parameter configurations after preliminary experiments. The racing process was stopped after at most 50 runs.

The Set of Allowed Predicates

At the end of the racing process we analyze all the results obtained by the best parameter configuration to derive a set of useful predicates of the domain. We chose to keep the predicates that appear in at least 50% of all best solutions found by DAE on 11 first instances. However, we could also choose predicates according to the proportion (e.g. more than 20%) of their occurrence with regard to all the atoms contained in the solutions found by DAE.

Detailed DAE Results

Divide-and-Evolve has been implemented within the Evolving Objects framework1 , an open source, template-based, ANSI C++-compliant evolutionary computation library. In order to compare the solution quality of DAE with the optimal results found by CPT, we used the gold-miner domain of the sixth International Planning Competition (IPC6) learning track for the preliminary tests.

The chosen evolution engine is a (10+70)-ES: 10 parents generate 70 offspring using variation operators, and the best of those 80 individuals become the parents of the next generation. The same stopping criterion was used for all experiments: after a minimum number of 10 generations, evolution is stopped if no improvement of the best fitness in the population is made during 20 generations, with a maximum of 100 generations altogether.

After the racing, the DAE planner is run with the best parameter configuration found on each target problem of the gold-miner domain in at most 15min of CPU time.

Figures 1, 3 and6 show for all 3 algorithms, the makespan of all target instances of the gold-miner domain, each column corresponding to an instance (number on the X axis). For the deterministic YAHSP and CPT, symbols ('@' and '#' respectively) indicate the makespan found. For the stochastic DAE, standard boxplots sketch the distribution of the 11 makespans.

Figure 2 shows, for all results found by DAE, the distribution of predicate frequency (i.e. the number of solution plans containing at least one occurrence of a given predicate divided by the total number of solution plans).

Figure 5 shows the occurence proportion of a predicate (i.e. the number of occurrences of a given predicate divided by the total number of atoms).

Significantly, after the racing, the solutions found by DAE are very close to the optimal solution found by CPT (see Figure 1). Solution plans found by allowing the predicates for which the distribution frequency is greater than 50% are on average better than those using more than 20% of proportions of predicate occurrences (see Figure 3 and6). We note that the distribution frequency of a predicate is not correlated to its occurence proportion (see Figures 2 and5).

The other observation concerns the quality of the results (see Figure 4), for problem 27 for instance, the choice of predicates (with the distribution frequency) often reduces the variance of the solutions plans found by DAE. However, for the problem 11, these choices can also increase the variance of the solutions plans found by DAE. Nevertheless, for the problem 30, the racing and predicates selection (with the distribution frequency) improve more often the average of results found by DAE on each problem.

We can also notice that at least one run of DAE with racing and predicates selection found 13 times the optimal value when DAE after the racing has found just 7 (see Figure 4). However, DAE racing only found some optima solutions whereas DAE with racing and predicates selection cannot. The choice of the allowed set of predicates could explain this behavior. Indeed, a set of predicates can be relevant for one problem and not for another one. Gold-miner instances makespan @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ 1 2 3 4 5 6 7 8

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

Discussion and Conclusion

It is well known that parameter tuning is one of the weaknesses of evolutionary algorithms in general. Divide-and-Evolve is not an exception. In this paper we introduced a two step learning approach in order to enhance DAE performance on a specific domain. Gold-miner instances makespan @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ 1 2 3 4 5 6 7 8 # CPT-optimal @ YAHSP Figure 3: Optimal number of actions for CPT (#), YAHSP value (@) and DAE (standard boxplots sketch the distribution of the 11 makespans after the racing step and the predicate selection (more than 50% of frequency distribution of all solutions -see Figure 2)) on the gold-miner IPC-6 domain. Gold-miner instances makespan @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ 1 2 3 4 5 6 7 8 IPC-6 gold-miner domain showed that our approach can improve the average quality of solutions obtained by DAE after the racing step. But there is still room for improvement in tuning DAE's parameters. First, the choice of the set of parameter configurations for the racing step is still an open issue. Although we obtained good results with twenty parameter configurations, we might miss a parameter configuration that would improve these results.

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 50
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 10
According to the results, the selection of allowed predicates for the second step of our learning approach is still open. We plan to combine the frequency distribution with the occurrence proportions of predicates in order to choose more efficiently the most relevant set of allowed predicates of a specific domain. These directions will be pursued during further research.

Figure 1 :

 1 Figure1: Optimal number of actions for CPT (#), YAHSP value (@) and DAE (standard boxplots sketch the distribution of the 11 makespans after the racing step) on the gold-miner IPC-6 domain.

Figure 4 :Figure 5 :

 45 Figure4: Comparison between DAE after the racing step and the predicate selection (blue boxplots) (more than 50% of frequency distribution of all solutions -see Figure2)) and DAE after the racing step (black boxplots) on the gold-miner IPC-6 domain. Optimal number of actions for CPT (#).

Figure 6 :

 6 Figure 6: Optimal number of actions for CPT (#), YAHSP value (@) and DAE (standard boxplots sketch the distribution of the 11 makespans after the racing step and predicate selection (more than 20% of occurences proportion-see Figure 5)) on the gold-miner IPC-6 domain.

http://eodev.sourceforge.net/


Francesco Dolce
email: dolcefra@fit.cvut.cz

Pierre-Adrien Tahay
email: pierre.adrien.tahay@cvut.cz

Column representation of Sturmian words in cellular automata

Keywords: Sturmian words, cellular automata, quadratic numbers, continued fraction expansion

We prove that, given a Sturmian word w with quadratic slope, it is possible to construct a one-dimensional cellular automaton such that w is represented in a chosen column in its space-time diagram. Our proof is constructive and use the continued fraction expansion of the slope of the word.

Introduction

Sturmian words are infinite words over a binary alphabet that have exactly n + 1 factors of length n for each non-negative n. Their origin can be traced back to the astronomer J. Bernoulli III. Their first in-depth study is by Morse and Hedlund [START_REF] Morse | Symbolic dynamics II. Sturmian trajectories[END_REF]. Many combinatorial properties were described in the paper by Coven and Hedlund [START_REF] Coven | Sequences with minimal block growth[END_REF]. Sturmian words are one of the most studied topics in combinatorics on words. They can be defined in different ways and have various interpretations in several domains, including combinatorics and discrete geometry. A possible way to describe them is by using mechanical words.

In this paper we consider characteristic Sturmian words with quadratic slope (see Section 2 for precise definitions). In particular, given such a Sturmian word w with continued fraction of its slope α = [0, 1 + b 1 , b 2 , . . . , b m , a 1 , a 2 , . . . , a k ] and corresponding directive sequence

∆ = (b 1 , b 2 , • • • , b m , (a 1 , a 2 , • • • , a k ) ω ),
we have w = lim n→∞ w n , where w 0 = a, w 1 = b, w n = w bn n-1 w n-2 for 1 ≤ n ≤ m and w m+n = w a (n mod k) m+n-1 w m+n-2 for n > 0 (where a 0 := a k ). We use such a characterisation to construct a machine, called cellular automaton, that will "print" us the infinite word w. A cellular automaton is a dynamical system defined by an infinite string of symbols over an alphabet and a map, called local rule, that transforms every symbol of the string according to its neighbourhood (see Section 3 for a more precise definition). A classical example of a 2-dimensional

The research received funding from the Ministry of Education, Youth and Sports of the Czech Republic through the projects CZ.02.1.01/0.0/0.0/16_019/0000765 and CZ.02.1.01/0.0/0.0/16_019/0000778. cellular automaton is given by 1970 Conway's Game of Life. In spite of its very simple definition, the Game of Life has some quite remarkable properties. Indeed, Rendell proved that starting from it, it is possible to simulate any Turing machine [START_REF] Rendell | Turing machine universality of the game of life[END_REF]. In this paper we consider one-dimensional cellular automata. The initial configuration will be given by an infinite string over a (finite) alphabet.

Our main result is the following one.

Theorem 1. A Sturmian word with quadratic slope can be represented as a column in the space-time diagram of a one-dimensional cellular automaton.

The task of representing a sequence over a finite alphabet in the space-time diagram of a cellular automaton is a non-trivial and still not entirely explored topic. One of the first results on the subject is the construction of the characteristic sequence of primes numbers, done by Fischer in 1965, using a cellular automaton with more than 30,000 states [START_REF] Fischer | Generation of primes by a one-dimensional real-time iterative array[END_REF]. In 1997, Korec gives another construction with only 11 states [START_REF] Korec | Real-time generation of primes by a one-dimensional cellular automaton with 11 states[END_REF]. In 1999, Mazoyer and Terrier establish several geometric constructions of increasing functions, which they call Fischer constructible [START_REF] Mazoyer | Signals in one-dimensional cellular automata[END_REF]. A very interesting result is given by Rowland and Yassawi in 2015 [START_REF] Rowland | A characterization of p-automatic sequences as columns of linear cellular automata[END_REF]: they give a complete characterisation of the construction of q-automatic sequence, with q a power of a prime number, in the columns of the space-time diagram of linear cellular automata. Finally, in 2018, Marcovici, Stoll and Tahay construct different non-automatic sequences, such as the characteristic sequences of polynomials (squares, cubes, etc.) and the Fibonacci word [START_REF] Marcovici | Construction of some nonautomatic sequences by cellular automata[END_REF]. Our main result in this paper can be seen as an extension of the construction obtained by Marcovici, Stoll and Tahay for this last infinite word. While in [START_REF] Marcovici | Construction of some nonautomatic sequences by cellular automata[END_REF] the authors use ad hoc properties of Fibonacci numbers and Fibonacci finite words, in this article we consider the development of the continued fraction associated a Sturmian word with quadratic slope to define a new algorithm. Such an algorithm could even be generalised to larger families of infinite words (see Section 6).

Preliminaries

In this section we recall some basic definitions on finite and infinite words. For all undefined terms we refer to [START_REF] Lothaire | Algebraic combinatorics on words[END_REF]. We denote the set of integers, of positive integers and of non-negative integers respectively by Z, Z + and N, while Q and R denote the set of rational and of real numbers.

Words

An alphabet Σ is a (finite) set of symbols called letters. The set of finite words Σ * over Σ is the free monoid having neutral element the empty word ε. We also denote by Σ + the free semigroup over Σ, e.g., Σ + = Σ * \ {ε}. The length of a word w = a 0 a 1 • • • a n-1 , with a i ∈ Σ, is the non-negative integer |w| = n. The length of ε is considered to be 0. When it is possible to write w = pus, with p, u, s ∈ Σ * , we call p (resp., s, u) a prefix (resp., suffix, factor ) of w.

An infinite word over Σ is a sequence w = a 0 a 1 a 2 • • • , with a i ∈ Σ for all i. Similarly to finite words we can define the set of infinite words Σ N and extend in a natural way to Σ * ∪ Σ N the notions of prefix, suffix and factor. An infinite word w is eventually periodic if w = uv ω = uvvv • • • . When u = ε we say that w is purely periodic. An infinite word that is not eventually periodic is called aperiodic. The factor complexity of an infinite word w is the mapping C w : N → N defined by C w (n) = #{u | u is a factor of w and |u| = n}.

Aperiodic infinite words with the lowest possible factor complexity, i.e., such that C w (n) = n + 1 for all n ∈ N, are called Sturmian words (for other equivalent definitions see [START_REF] Balková | Sturmian jungle (or garden?) on multiliteral alphabets[END_REF]). It follows from the definition that all Sturmian words are defined over a binary alphabet, e.g., {a, b}. If both sequences aw and bw are Sturmian, we call w a characteristic Sturmian word. The family of Sturmian words coincide with the family of irrational mechanical words, as well as the family of binary balanced aperiodic words (for more on balanced words see, for instance, [START_REF] Balková | Sturmian jungle (or garden?) on multiliteral alphabets[END_REF][START_REF] Dolce | On balanced sequences and their asymptotic critical exponent[END_REF][START_REF] Hubert | Suites équilibrées[END_REF]). In particular, characteristic Sturmian words correspond to balanced irrational mechanical words with intercept equal to the slope [START_REF] Lothaire | Algebraic combinatorics on words[END_REF]. We denote by c α the unique characteristic Sturmian word with slope (and intercept) α. Since two Sturmian words with the same slope (but different intercept) share the same set of factors, we will focus only on the study of characteristic Sturmian words. Note that every mechanical word with rational slope is a purely periodic word.

Continued fraction expansion

Let θ be a real number. A continued fraction expansion of θ is defined as

[c 0 , c 1 , c 2 , . . .] whenever θ = c 0 + 1 c 1 + 1 c2+ . . .
with c 0 ∈ Z and c i ∈ Z + for every positive i. It is known that if θ ∈ Q then there exist exactly two continued fraction expansion of θ and they are both finite. On the other hand every positive irrational number corresponds to a unique infinite continued fraction expansion with c 0 ∈ N and

c i ∈ Z + for all i ≥ 1. Note that, if c 0 = 0 then 0 ≤ θ ≤ 1.
If θ is a quadratic irrational, then its continued fraction expansion is eventually periodic, that is it will be of the form

θ = [b 0 , . . . , b m , a 1 , . . . , a k ] = [b 0 , . . . , b m , a 1 , . . . , a k , a 1 , . . . , a k , . . .] . Example 2. The golden ratio ϕ = 1+ √ 5 2
is a quadratic irrational number. Its continued fraction expansion is ϕ = 1 = [1, 1, 1, . . .] The continued fraction expansions of e and π are [2, 1, 2, 1, 1, 4, 1, 1, 6, . . .] and [3, 7, 15, 1, 292, 1, 1, 1, 2, . . .] respectively (sequences A003417 and A001203 in the OEIS [START_REF] Sloane | On-line Encyclopedia of Integer Sequences[END_REF]).

Standard sequences and directed sequences

Let ∆ = (d n ) n≥1 be an integer sequence with d 1 ∈ N and d n ∈ Z + for every positive integer n. The standard sequence associated to ∆ is the sequence of finite words (w n ) n≥-1 defined by

w -1 = b, w 0 = a, w n = w dn n-1 w n-2 for every n ≥ 1.
The sequence ∆ is also called the directive sequence of (w n ) n≥-1 . Note that if d 1 > 0, every w n starts with a. Otherwise, w n starts with b for every n = 0. Let us consider the infinite word w = lim n→∞ w n . Such an infinite word is well defined and w n is a prefix of w for every positive n. We say that ∆ (resp. (w n ) n≥-1 ) is the directive sequence (resp. the standard sequence) of w.

Example 3. Let us consider the directive sequence ∆ = (1 ω ) = (1, 1, 1, . . .). The associated standard sequence is f -1 = b, f 0 = a, and f n = f n-1 f n-2 for every n ≥ 1.
It is known that the Fibonacci word f defined in Example 1 can be obtained as the limit f = lim n→∞ f n .

Proposition 1 ([14]

). Let α be an irrational number with 0 < α < 1 having continued fraction expansion α = [0,

d 1 + 1,
(v n ) n≥-1 is defined by v -1 = b, v 0 = a, v 2n+1 = v 2n v 2n-1 and v 2n = v 2 2n-1 v 2n-2 for every n ∈ N. Since 3- √ 3 3 = 0, 2, 2, 1 , we have, according to Proposition 1, v = c (3- √ 3)/3 .
In the next sections we will also need the lengths of the prefixes of a characteristic Sturmian sequence. Given a standard sequence (w n ) n≥-1 we define for every integer n ≥ -1 the number W n = |w n |.

F n = F n-1 + F n-2 , V 2n+1 = V 2n + V 2n-1 and V 2n = 2V 2n-1 + V 2n-2 .

Cellular automata

In the following we use the terminology developed by Mazoyer and Terrier in [START_REF] Mazoyer | Signals in one-dimensional cellular automata[END_REF] and Marcovici, Stoll and Tahay in [START_REF] Marcovici | Construction of some nonautomatic sequences by cellular automata[END_REF]. Definition 1. A one-dimensional cellular automaton (CA) is a dynamical system (A Z , T ), where A is a finite set, and where the map T : A Z → A Z is defined by a local rule which acts uniformly and synchronously on the configuration space. More precisely, there exists an integer r ∈ N called the radius of the CA, and a local rule τ : A 2r+1 → A such that for every x = (x k ) k∈Z and for every k ∈ Z, we have T (x) k = τ ((x k+i ) -r≤i≤r ).

When the set A is understood, we will call cellular automaton just the map T . The elements of A Z are called configurations. By the Curtis-Hedlund-Lyndon Theorem [START_REF] Hedlund | Endomorphisms and automorphisms of the shift dynamical system[END_REF], a map T : A Z → A Z is a CA if and only if it is continuous with respect to the product topology and it commutes with the shift map σ defined by σ(x) k = x k-1 , for every configuration x = (x k ) k∈Z and every k ∈ Z. Let 0 ∈ A and T : A Z → A Z a cellular automaton. We say that

T is 0-quiescent if T (0 Z ) = 0 Z = • • • 000 • • • . A configuration x = (x k ) k∈Z is called finite if the set {k ∈ Z : x k = 0} is finite.
A cellular automaton can be visualized by using a space-time diagram which is a 2-dimensional grid where each cell contains an element of the set A and is represented by a space coordinate and a time coordinate.

Let us consider the set S = (T n (x) 0 ) n≥0 ∈ A N : T is a 0-quiescent CA on A Z and x is finite .

In other words, S is the set of sequences of A N that can occur as the first column (and thus as any column) in the space-time diagram of some onedimensional 0-quiescent CA, starting from a finite initial configuration. This set corresponds to the set of Fischer's produced sequences in [START_REF] Mazoyer | Signals in one-dimensional cellular automata[END_REF].

In a space-time diagram it is also possible to "transmit information" through signals, that is to connect two cells (m, n) and (m , n + t) through a monotonous path; we call slope of the signal the number t m -m (see [16, Definitions 3 and 4] for a formal definition). When m = m , we call such a signal a vertical signal or a signal of infinite slope. For the sake of simplicity, we usually represent a signal as a straight line between the cells (m, n) and (m , n + t). Signals are usually "porous", i.e., they do not interact between each other. In some case, however, we also need to consider "concrete" signals. In particular, let us define two distinct kinds of walls. We say that a wall is of type (i) whenever a given signal hitting the wall bounces from the cell just above, i.e., when a given signal of slope d arrives in a cell ( , t), then such a signal dies and a new signal of slope -d starts from the cell ( , t + 1). We say that a wall is of type (ii) whenever a given signal hitting the wall bounces from the same cell, i.e., when a given signal of slope d arrives in a cell ( , t), then such a signal dies and a new signal of slope -d starts from the same cell ( , t). We usually represent a wall of type (i) as a line inside the column + 1 when the signals comes from the left (resp. -1 when the signal comes from the right), and a wall of type (ii) as a rectangle containing the cells in the column (see Figure 1). When two signals meet, we can mark the cell at the intersection, i.e., assign to it a value from the set A, and define new signals starting from it.

Construction of numbers

To prove our main result we proceed in two steps. First let us construct a CA recognising the lengths W n of the prefixes w n of our Sturmian word w.

Let X ⊂ Z. Let us denote by 1 X the characteristic function of X, that is the map 1 X : Z → {0, 1} defined by 1 X (x) = 1 iff x ∈ X.

Proposition 2 ([16]

). Let (S n ) n≥0 be an integer sequence defined by S n+p = p-1 i=0 a i S n+i , where p, a i ∈ N. Then 1 {Sn} n≥0 ∈ S.

Mazoyer and Terrier give an explicit method to build this sequences in a column of a CA. We propose here a different construction for a particular case that will be necessary for representing a Sturmian word of quadratic slope. Proposition 3. Let (d n ) n≥1 be an eventually periodic integer sequence with d 1 ∈ N and d i ∈ Z + for every i ≥ 2. Let (S n ) n≥0 be the integer sequence defined by S n = d n S n-1 + S n-2 for every n ≥ 0, with S -1 , S 0 ∈ Z + . Then

1 {Sn} n≥0 ∈ S. Proof (Sketch). Since the sequence (d n ) n≥1 is eventually periodic, then there exist m ∈ N, k ∈ Z + and b 1 , . . . , b m , a 1 , . . . a k ∈ N such that (d n ) n≥1 = (b 1 , . . . , b m , (a 1 , . . . , a k ) ω
) . Note that we can consider the first S m rows of the cellular automaton as initial conditions, i.e., we can start the construction from the row of rank S m . Let n ≥ 1 be an integer. We are going to consider two distinct cases according to the value of d n+1 .

Let us first suppose that d n+1 = 1. Assume that we have already marked the cells (0, S n ), (S n-2 , S n ), (S n-1 , S n ), (S n-1 + S n-2 , S n ) and (S n , S n ). We claim that we can mark the cells (0, S n+1 ), (S n-1 , S n+1 ), (S n , S n+1 ), (S n +S n-1 , S n+1 ) and (S n+1 , S n+1 ) In order to do that, we use the relation

S n+1 = S n + d n (d n+1 -1)S n-1 + (d n+1 -1)S n-2 + S n-1 .
The idea is to consider intermediate rows, such that their distance is given by the addends in the previous sum. When two signals meet, they died and we can use the cell on the intersection to define other signals. The slope of each signal is determined by the ratio between the difference between the time coordinates and the difference between the space coordinates. For example, the slope of the signal between the cell (0, S n ) and the cell

(S n-1 , S n + d n (d n+1 -1)S n-1 ) is (S n + d n (d n+1 -1)S n-1 ) -S n S n-1 -0 = d n (d n+1 -1)
When d n+1 = 1 the construction is different. In this case the three rows S n , S n +d n (d n+1 -1)S n-1 and d n+1 S n coincide, as well as the two rows S n +S n-1 and S n+1 (resp. the two columns S n + S n-1 and S n+1 ). We start with the cells (0, S n ), (S n-2 , S n ), (S n-1 , S n ) and (S n , S n ) and we claim that we can mark the cells (0, S n+1 ), (S n-1 , S n+1 ), (S n , S n+1 ) and (S n+1 , S n+1 ).

The construction of the sequence (S n ) n≥0 is illustrated in Figures 2 and3. Using these figures it is not hard to recover exact definitions of the signals. For sake of simplicity we use the same colour when two signals have the same slope and we represent only the cells on the intersections between two signals. In particular, in both cases we are able to mark the cell (0, S n ). Hence we can mark the sequence {S n } n≥0 on the column 0. To complete the proof it is enough to put the letter 1 in the cells (0, S n ), for every n ≥ 0 and the letter 0 in all other cells in the column 0.

Sn Sn + S n-1 Sn + dn(d n+1 -1)S n-1 d n+1 Sn S n+1 = d n+1 Sn + S n-1 0 S n-2 S n-1 S n-1 + S n-2 Sn Sn + S n-1 S n+1 slope 1 slope dn(d n+1 -1) slope d n+1 -1 slope dn (d n+1 -1) dn +1 slope 1 -d n+1 slope d n+1 -1 d n+1 slope -1 slope 1 dn vertical signal
Note that the hypothesis of eventual periodicity of the sequence (d n ) n≥1 in the previous proof is essential to guarantee that the cellular automaton is defined over a finite set A. Indeed, since the signals (and their slope) are periodically repeated, we have Card(A) = O(k), where k is the length of the maximum between the pre-period and the eventual period of (d n ) n≥1 .

Example 7. Let us consider the word f defined in Example 1 and the associated numerical sequence (F n ) n≥0 . According to Proposition 3 we have 1 {Fn} n≥0 ∈ S.

Construction of prefixes

In this section we prove our main theorem. In order to do that, we need some preliminary results. Proposition 4. Let ∆ = (d n ) n≥1 be an eventually periodic integer sequence with d 1 ≥ 0 and d n ≥ 1 for every n > 1; let (w n ) n≥-1 be the standard sequence associated to ∆ defined by w -1 = b, w 0 = a, and w n = w dn n-1 w n-2 for every n ≥ 1. Then w = lim n→∞ w n ∈ S In order to prove Proposition 4, we need the following result stating that a letter in a cell of a CA can be recopied in the same column and in any row above. Moreover, we can do it by using only walls and signals of slope 1 and -1.

Lemma 1. Let T : A Z → A Z be a CA, a ∈ A, and n, m, t, t ∈ Z with m > 0, and t > t ≥ 0. Suppose that the cells (n, t ), (n + m, t ) and (n, t) are marked, the last one with a. Then it is possible to mark the cell (n, t + m) with a.

Proof. Without loss of generality, suppose that n = 0 and t = 0. and so a is in the cell (0, t). Let us prove that we can recopy a into the cell (0, t + m). To mark the wall we are going to use the cell (m, 0). We consider two distinct cases, according to the parity of m.

Let m be odd. We consider a wall of type (i) to the right of the column m 2 . Such a wall can be defined, for instance, using two signals of slope 1 and -1 starting respectively at (0, 0) and at (m, 0) (see left of Figure 4).

We send a signal of slope 1 from the cell (0, t); when this signal touches the wall we define a new signal of slope -1 (starting from the cell above); when this new signal meets the column 0, we write the letter a. Since m 2 = m-1 2 , we have that the new a is exactly on the row t + m-1 2 + m-1 2 + 1 = t + m. as wanted. Suppose now that m is even. We consider this time a wall of type (ii) in the column m 2 . Such a wall can be defined, for instance, using two signals of slope 1 and -1 starting respectively at (0, 0) and at (m, 0) (see right of Figure 4). We send a signal of slope 1 from the cell (0, t); when this signal hit the wall, we define a new signal of slope -1 (starting from the same cell); when this new signal meets the column 0, we write the letter a. The new a is exactly on the row t + m 2 + m 2 -1 + 1 = t + m, as wanted. We can now prove Proposition 4.

Proof (of Proposition 4). Let us denote S n = |w n | for all n ≥ -1. For every n ≥ 1 we have S n = d n S n-1 + S n-2 . Since, for all n ≥ 2, w n-3 is a suffix of w n-1 , the word w n-3 w n-2 is a suffix of w n . Suppose the word w n constructed for a given n ≥ 2 and suppose that the last letter of w n is in the cell (0, S n ).

Let us first suppose that d n = 1. We will show that it is possible to construct the word w n+1 with its last letter in the cell (0, S n+1 ). In order to do that, we will use the relation

w n+1 = w n ((w dn-1 n-2 w n-3 ) dn w n-2 ) dn+1-1 w dn-1 n-2 w n-3 .
Let us take up the construction of Proposition 3 until the number S n . Moreover, in the column 0 we mark the cell (0, S n-1 + d n-1 (d n -1)S n-2 ). For this, we define a signal of slope -(d n-1 (d n -1)) from the cell (S n-2 , S n-1 ). This signal meets the vertical signal defined from (0, S n-1 ) in the required cell (0, S n-1 + d n-1 (d n -1)S n-2 ). From this cell we define a signal P (n-1) 1 of slope 1 and from the cells n-2 w n-3 = w n+1 in another way. We have to copy the letters of the word w n-3 w n-2 in the column 0 above at distance S n-3 + S n-2 + d n-1 S n-2 = S n-1 + S n-2 = S n . In order to do that, we define a wall M (n-1) 5 of type (i) or (ii) (according to the parity) by intersecting a signal P (n-1) 3 of slope 1 from (S n-2 , S n-1 ) and a signal N (n-1) 2 of slope -1 from (S n-1 , S n-1 ) defined as in the previous case. Such a wall is in the column

(S n-2 , S n-1 + d n-1 (d n -1)S n-2 ), (S n-1 , S n-1 + d n-1 (d n -1)S n-2 ) and (S n-1 + S n-2 , S n-1 + d n-1 (d n -1)S n-2 )
S n-2 + Sn-1-Sn-2 2 = S n -S n-1 + S n-1 -Sn 2 = Sn 2 .
We send a signal of slope 1 from each letter of w n-3 w n-2 to the wall M (n-1) 5

and from there we send a signal of slope -1 until the column 0. Therefore, we have constructed the word w n w dn-1 n-2 w n-3 w n-2 . The rest of the construction is the same as before, with M (n) 4

playing the role of the wall M (n-1) 3 for the rows above the word w n-3 w n-2 .

Note that the signals P , the last one is destroyed and the signal P

Conclusions

To prove our results we used the continued fraction expansion associated with w. However, to use Proposition 4 it is enough to know how to decompose each prefix w n in compositions of powers of smaller prefixes. A different approach could be to use the morphisms Indeed, it is known that for every standard Sturmian sequence w there exist a unique sequence of words w (i) i and an infinite sequence (ψ i ) i ∈ {G, D} N of morphisms such that w = lim n→∞ ψ 0 ψ 1 . . . ψ n (w (n) ) (see, for instance, [START_REF] Justin | Episturmian words and episturmian morphisms[END_REF]). Using these notions, and the strictly related notion of S-adicity, we think it is possible to generalise our results to larger classes of words and languages, namely Arnoux-Rauzy words and dendric words (see, e.g., [START_REF] Berthé | Acyclic, connected and tree sets[END_REF][START_REF] Berthé | Rigidity and substitutive tree words[END_REF][START_REF] Dolce | Eventually dendric shift spaces[END_REF]).

Blanchard and Kůrka [START_REF] Blanchard | Language complexity of rotations and Sturmian sequences[END_REF] consider a larger family of languages that can be recognised by a non-deterministic Turing machine. This family contains the languages corresponding to quadratic numbers but also the ones corresponding to Hurwitz numbers, such as e (a Hurwitz number is an irrational such that its continued fraction expansion is a polynomial mixture [START_REF] Hurwitz | Über die Bedingungen, unter welchen eine Gleichung nur Wurzeln mit negativen reellen Theilen besitzt[END_REF]). An interesting question is whether it is possible to generalise our construction to such a family as well.

Example 1 .

 1 Let us consider the well-known Fibonacci word f = abaababaab • • • . It is defined as the fixed point of the morphism sending a → ab and b → a. The word f is the characteristic Sturmian word c 1/ϕ 2 , where ϕ = 1+

  d 2 , d 3 , . . .], and let (w n ) n≥-1 be the standard sequence associated to (d 1 , d 2 , d 3 , . . .). Then c α = lim n→∞ w n . Note that mechanical words are defined for 0 ≤ α ≤ 1. It is possible to generalise such definition to every α ∈ R [14, Remark 2.1.12]: the fraction expansion in the statement of Proposition 1 would be α = [c 0 , d 1 + 1, d 2 , d 3 , . . .] with c 0 ∈ Z, but the standard associated sequence would not change. Example 4. As seen in Example 3, the directive sequence of the Fibonacci word f is (1 ω ). The continued fraction expansion of the corresponding irrational slope (see also Example 1) is 1 ϕ 2 = 0, 2, 1 . Example 5. Let us consider the characteristic Sturmian word v having associated directive sequence ∆ = ((1, 2) ω ). We have v = lim n→∞ v n , where the standard sequence

Example 6 .

 6 Let (f n ) n≥-1 , and (v n ) n≥-1 be the standard sequences defined in Examples 3, and 5. Let F n = |f n | and V n = |v n |. For every n ∈ N we have the relations

t t + 1 tFig. 1 .

 11 Fig. 1. Walls of type (i) (on the left) and of type (ii) (on the right).

Fig. 2 .Fig. 3 .

 23 Fig. 2. Construction of the number sequence (Sn) n≥0 when dn+1 = 1.

Fig. 4 .

 4 Fig. 4. Recopying of letters.

3 of slope - 1 . 3 in

 313 At the intersection of these with the signal P the columns S n-2 , S n-1 and to construct the word w n-3 w n-2 in w n (w dn-1 n-2 w n-3 w n-2 )

  can also be used to destroy the walls M (n-1) i previously constructed. Formally, when a signal P (n) i meets a wall M (n-1) i

Example 8 .

 8 The Fibonacci word f defined in Example 1 is in S.We can now prove our main result.Proof (of Theorem 1). Let w be a Sturmian word of quadratic slope, α = [0, 1 + b 1 , b 2 , . . . , b m , a 1 , a 2 , . . . , a k ] the continued fraction expansion of its slope, and (w n ) n≥-1 the standard sequence associated to the eventually periodic in-teger sequence ∆ = (b 1 , • • • , b m , (a 1 , • • • , a k ) ω) so that ∆ is an eventually periodic integer sequence. Following Proposition 1, we have w = lim n→∞ w n . Using Propositions 3 and 4, it is clear that w ∈ S.

S n-1 + S n-2 respectively. These walls will be of type (i) or of type (ii) according to the parity of the columns where N (n-1) i originate. Suppose we have already recopied the words w n-3 and w n-2 in the suffix of w n , i.e., we have found the two corresponding heights in the column 0. Now, from these two words, we are going to define signals of slope 1 that will stop against one of the three walls previously defined and send back signals of slope -1 that will recopy the same word in the column 0. Using Lemma 1 we can recopy one by one the letters of w n-2 and w n-3 . First, we send a signal of slope 1 from each letter of w n-2 . When this signal meets the wall M (n-1) 1

we send a signal of slope -1 until the column 0. Since N (n-1) 1 is generated from the column S n-2 the letters of w n-2 are recopied in the same column but above at distance S n-2 . This means that we have recopied a second word w n-2 above the first. We repeat this procedure d n-1 times and we get the word w dn-1 n-2 above the w n already constructed. Next, we send a signal of slope 1 from each letter of w n-3 to the wall M (n-1) 3 and, in a similar way, from there we send a signal of slope -1 to the column 0. The letters of w n-3 are recopied in the same column but above at distance

Thus, we have recopied a word w n-3 above the word w n-3 w n-2 w dn-1 n-2 , where the subword w n-3 w n-2 corresponds to a suffix of w n . So far we have constructed the word w n (w dn-1 n-2 w n-3 ) on the column 0. For the next step we need to use also the word w n-2 in the suffix of w n . This time, we send signals of slope 1 from cells in column 0 to the wall M (n-1) 3

and from there, signals of slope -1 to the column 0. Therefore, the letters of w n-2 are recopied in the same column but above at distance S n-1 + S n-2 = (d n-1 + 1)S n-2 + S n-3 . Hence we obtain the word w n (w dn-1 n-2 w n-3 )w n-2 in the column 0. By using the wall M and, from there, signals of slope -1 back to column 0. The letters of w n-3 are recopied in the same column but above at distance S n-1 = d n-1 S n-2 + S n-3 . Hence, we obtain the word w n (w dn-1 n-2 w n-3 )w dn-1 n-2 w n-3 . Similarly, it is easy to obtain the word w n (w

w n-3 ). Following the same idea, we use the wall M (n-1) 2 to recopy w n-2 in the column 0 but S n-1 = (d n-1 -1)S n-2 + S n-3 + S n-2 cells above and so, we obtain the word w n (w dn-1 n-2 w n-3 ) dn w n-2 . Since the suffix of this word is w n-3 w n-2 the previous steps can be applied again to obtain the word w n ((w dn-1 n-2 w n-3 ) dn w n-2 ) dn+1-1 which also has the suffix w n-3 w n-2 . Thus we easily obtain the word w n ((w

Let us now consider the case d n = 1. Here we have S n = S n-1 + S n-2 and the cell (S n-1 + S n-2 , S n-1 + d n-1 (d n -1)S n-2 ) = (S n , S n-1 ) is not marked. Therefore the wall M (n-1) 3 can no longer be defined as before. This time, we define a wall M


S Youssef Ai

S Nicaise
email: sylvie.nicaise@inrae.fr

N Benahmed
email: nadia.benahmed@inrae.fr

P Philippe
email: pierre.philippe@inrae.fr

A Bennabi
email: abennabi@estp-paris.euapoupardin@estp-paris.eu

Soil surface erosion: influence of soil properties on the erosion resistance

Keywords: Internal erosion, parametric study, hole erosion test, cohesive soils, erosion characteristics

L'érosion de conduit se produit lorsqu'il y a départ de particules de sol le long d'un chemin d'écoulement préférentiel existant dans un ouvrage hydraulique. Pour évaluer la vulnérabilité des sols vis-à-vis de ce phénomène, et plus généralement à l'érosion de surface, plusieurs dispositifs expérimentaux de laboratoire existent, parmi lesquels le Hole Erosion Test. Nous étudions ici l'impact de plusieurs paramètres constitutifs du sol sur sa résistance à l'érosion, à l'aide d'une campagne expérimentale réalisée sur des sols de teneur en argile, forme des grains grossiers et taux de compactage variables. En accord avec des résultats antérieurs, il est montré que la présence d'argile confère une meilleure résistance à l'érosion ; de même, plus l'échantillon est dense, plus il est résistant. Concernant la forme des grains, les essais réalisés montrent deux résultats : d'une part, et contre intuitivement, les échantillons avec des particules sub-arrondies à arrondies présentent une meilleure résistance vis-àvis de l'initiation de l'érosion, cependant, une fois l'érosion enclenchée, ces échantillons s'érodent plus rapidement que ceux ayant des particules de forme sub-angulaire à angulaire.

Introduction

The erosion of hydraulic earthworks by piping, among different kinds of so-called internal erosion, is widely considered as one of the main mechanism causing failures and accidents of embankment dams as described in Foster et al. 2000 [START_REF] Foster | The statistics of embankment dam failures and Figure 4 : Effect of clay mass fraction on the erosion characteristics, critical shear stress 𝜏 𝑐 and erosion coefficient 𝐾 𝑒𝑟 , for 7 various KA-HN mixtures. accidents[END_REF]. This topic has been broadly investigated in last decades and, owing that different soils show different erodibility, previous researches have attempted to investigate the effect of several parameters related to soil properties on the resistance of soil to surface erosion. However, data on this are not abundant, making difficult to draw reliable relationships between erosion characteristics and soil properties.

Assessing the resistance of soil to internal erosion lie in investigating its erodibility, which is quantified in terms of erosion rate and critical shear stress. The erodibility parameters are expressed by the following equation:

𝑚̇= 𝑘 𝑒𝑟 (𝜏 -𝜏 𝑐 )
Where 𝑚̇ is the mass rate of erosion in kg/m/s, 𝑘 𝑒𝑟 is the erosion kinetics coefficient in s/m, τ is the applied shear stress to soil surface in Pa, and 𝜏 𝑐 is the critical shear stress for initiation of erosion in Pa.

The present study focuses specifically on materials consisting of a mixture of coarse grains and fine particles, to investigate the effect of nature and state of soil on (i) the initiation conditions of erosion process, (ii) the erosion kinetics under stationary or near-stationary conditions. To this end, HET (Hole Erosion Test) device is used to determine the erosion parameters, 𝜏 𝑐 and 𝑘 𝑒𝑟 . A particular emphasis was put in studying the influence of different material properties on the process, such as coarse particles shape, compaction degree, and fines content on the erodibility of the tested soils.

HET apparatus and studied materials

Nowadays, Hole Erosion Test (HET) device is being used increasingly for studying the resistance of cohesive soils to surface erosion, in a configuration reminiscent, on a smaller scale, to concentrated leak erosion, one specific type of internal erosion. During this test, erosion takes place along the wall of a pipe, which was initially drilled with a cylindrical shape through the soil sample length. This pipe is assumed to remain cylindrical during the test with radial expansion [START_REF] Wan | Laboratory tests on the rate of piping erosion of soils in embankment dams[END_REF]. The apparatus used in the present study is shown in Figure1. It is composed by three chambers: upstream chamber, downstream chamber, and central chamber, which holds the soil sample. Flowmeter is installed at the inlet of the water flow to record the flowrate and pressure transducers are mounted on both sides of the soil sample to measure the differential pressure for calculating the hydraulic gradient. The flow rate is controlled by an outflow micrometric vane at the upstream side. Eroded mass is estimated using a turbidity meter at the downstream side, after specific calibration.

Experimental program and protocol

The influence of particle shape is investigated throughout 15 compacted soil samples (five samples per mixture) of 150 mm height each. All samples are composed, in mass, of 50% of fine cohesive material (KA) and 50% of coarse material (HN, TEN, or GB). Considering the influence of fines mass fraction on the resistance of soil to erosion, 7 soil samples are tested with different fine contents, namely 20%, 30%, 40%, 50%, 70%, 90%, and 100%. Due to the effect of sample preparation on the structure of sandclay mixtures [START_REF] Yin | Influence of sample preparation on the multi scale structure of sand-clay mixtures[END_REF], one mixing procedure of the different components is chosen: firstly, dried granular material is mixed with fine material, then water is sprayed over the mixture with an amount corresponding to the optimum water content determined by Proctor test. After 72 hours of storage in hermetic bags for moisture content homogenisation, samples are compacted dynamically in 5 layers to the desired bulk density. As regards both particle shape and density effects, the soil of each mixture is compacted at 5 different degrees of about 88%, 92%, 95%, 98%, and 100% of Standard Proctor maximum dry density, corresponding to dry densities varying from 1.51 t/m 3 to 1.72 t/m 3 . Furthermore, as far as the effect of fine percentage is concerned, all soil samples are identically compacted at 95% of maximum dry density. Each hole erosion test is performed immediately after sample compaction. Soil sample compaction and testing procedure are described in more detail in [START_REF] Benahmed | Investigating concentrated leak erosion behaviour of cohesive soils by performing hole erosion tests[END_REF].

In order to determine erosion law's parameters, namely 𝜏 𝑐 and 𝐾 𝑒𝑟 , we followed two of the main methods of interpretations currently used in the literature [START_REF] Bonelli Stephane | Erosion of Geomaterials[END_REF]- [START_REF] Lachouette | One-dimensional modeling of piping flow erosion[END_REF]. The application of these methods requires some input data such as soil density, water content, the evolution over time of the outlet turbidity, the total pressure drop, the average flowrate, and the final diameter of the eroded hole. One should point out that the latter is the most critical one, which is determined post-mortem at the end of the test.

Experimental results

Results are obtained for the variation of the erosion parameters, critical shear stress 𝜏 𝑐 and erosion coefficient 𝐾 𝑒𝑟 , with respect to the coarse grain shape, compaction degree, and the fines mass fraction of the soil sample. These parameters are determined considering three regression models, which use different input data of the experimentations. Each point represents a weighted average of the three values determined throughout these methods of interpretation, the weight of each being based on the corresponding R 2 value. Error bars indicate the standard deviation for each determined critical shear stress and erosion coefficient by means of the three interpretation models.

Influence of coarse grain shape

The experimental results obtained for the 15 tests focusing on the effect of coarse grain shape on the resistance of soil to surface erosion are presented in Table 1 and plotted in Figure 3. As can be observed, the mixture with rounded particles (KA50GB50) shows the highest critical shear stress 𝜏 𝑐 , whereas the mixture with angular particles (KA50HN50) presents the lowest one. The difference in behavior of these two types of materials become more pronounced as far as the compaction degree increases. Moreover, for compaction degrees less than 95%, the critical shear stresses of KA50TEN50 are quite close from those of KA50HN50, while for compaction degrees above 95%, there is a substantial divergence between the two mixtures to such an extent that, at 100% compaction degree, KA50TEN50 reaches the KA50GB50 value.

Concerning the erosion coefficient 𝐾 𝑒𝑟 , and considering global tendencies, the rounded particles mixture displays also the highest values compared to the two other mixtures, and, similarly to 𝜏 𝑐 , sub-angular mixture show the lowest erosion coefficient. Besides, as far as the degree of compaction increases, the erosion coefficient of the three mixtures tends to converge. Accordingly, contrary to 𝜏 𝑐 , the variation of 𝐾 𝑒𝑟 for the different mixtures demonstrates that the presence of glass beads seems to lead to the highest soil's erodibility, once erosion is initiated. The effect of compaction state on the resistance to erosion of the tested samples can also be discussed based on the results presented in Figure 3. In the whole range of compaction degrees and especially above 92%, the mixtures made of glass beads and Loire sand, namely KA50GB50 and KA50TEN50 respectively, reveal a significant increase of the critical shear stress 𝜏 𝑐 with the degree of compaction, particularly for Loire sand mixture which varies from 8 Pa to 35 Pa at maximum density. For the other mixture (KA50HN50), the increase of 𝜏 𝑐 remains much more limited below 95% and then becomes more marked beyond. The dependence of the erosion coefficient 𝐾 𝑒𝑟 on the compaction degree is also shown on Figure 3. Similarly, to 𝜏 𝑐 , there is a change in behavior around 95%, where the value of 𝐾 𝑒𝑟 decreases as the compaction degree increases, with the coefficient obtained for the glass beads mixture being systematically much higher than the one of the two other mixtures. Above 95% compaction degree, there is almost a plateau, or even a very slight decrease of ker, with values fairly comparable for the three mixtures.

Influence of compaction degree

Considering the two parameters 𝜏 𝑐 and 𝐾 𝑒𝑟 , this analysis leads to the conclusion that an increase in soil's density leads to a higher resistance to surface erosion. Nevertheless, there is a quantitative discrepancy when considering the shape of the coarse grains: regardless of their nature, the mixtures with sands show very similar values (except for 𝜏 𝑐 at the highest compaction degrees) while the glass beads mixture has significantly higher erosion parameters, both the critical stress and the erosion coefficient.

Influence of fines mass fraction

The experimental results related to the effect of fines fraction on the erosion parameters of samples prepared with mixtures of Armoricaine kaolinite and Hostun sand are given in Table 2 and Figure 4.

From the variation of both the critical shear stress 𝜏 𝑐 and the erosion coefficient 𝐾 𝑒𝑟 with respect to the percentage of fine material in the sand soil mixture, it can be seen that, for fines content ≤ 40-50%, there is a significant increase of the critical shear stress 𝜏 𝑐 , whereas the erosion coefficient 𝐾 𝑒𝑟 shows a steep decrease. For mixtures composed of more than 50% of fines, the erosion coefficient remains very small, on a kind of plateau about 1 s.m -1 with a minimum value of 0.74 at 70% of fines content. In the same time, 𝜏 𝑐 keeps increasing, but less sharply, until reaching almost 15 Pa at 100% of fines where the soil sample is essentially constituted of pure Armoricaine kaolinite. Overall, the trend is therefore towards a more resistant material as the fine fraction is increased. However, as highlighted here and clearly observed in Figure 4, one should point out that there are two distinct regimes of the erosion parameters evolution with a demarcation point around 50% of fines content.

Discussion & Conclusions

Based on HET test, the present study investigated the impact of the coarse particles shape, the soil's compaction degree, as well as the fines content on the resistance to surface erosion of samples constituted by mixtures of fines and coarse grains. The main outcomes are as follow:

1. The increase in soil's density leads, generally, to a decrease in soil's erodibility, which is in agreement with previous results [START_REF] Benahmed | Investigating concentrated leak erosion behaviour of cohesive soils by performing hole erosion tests[END_REF], [START_REF] Choo | Laboratory and theoretical evaluation of impact of packing density, particle shape, and uniformity coefficient on erodibility of coarse-grained soil particles[END_REF]. The response to the increase of soil's density seems to depend on the shape of coarse material. Indeed, it appears that the more rounded the coarse particles in the mixture the sharper the response for the erosion coefficient 𝐾 𝑒𝑟 . 2. Concerning particles shape, the experimental results point out two counterintuitive conclusions since, on the one hand, the highest 𝜏 𝑐 value is found for the more rounded coarse particles. This is in contradiction with the findings of Choo et al. 2020, who pointed out that the increase in particles angularity leads to an increase of resistance to erosion [START_REF] Choo | Laboratory and theoretical evaluation of impact of packing density, particle shape, and uniformity coefficient on erodibility of coarse-grained soil particles[END_REF]. However, Choo's study concerned purely coarse material erosion and cannot be strictly compared to our study on cohesive soils. This fact might probably be attributed to the manner that the coarse particles were coated by fines. On the other hand, it was showed that 𝐾 𝑒𝑟 is higher for rounded particles mixtures, which defines a more erodible soil in comparison to angular particles mixtures. Considering a non-cohesive soil, Guo et al. 2018 [START_REF] Guo | Influence of particle shape on the erodibility of noncohesive soil: Insights from coupled CFD-DEM simulations[END_REF] gave a possible explanation to this observation as they demonstrated numerically that the ellipsoidal particles show more resistance to erosion in comparison to spherical particles due to the lack of rotation resistance for spheres. 3. The increase in fines content results in less erodible soils, which is fully consistent with previous studies [START_REF] Haghighi | Improvement of Hole Erosion Test and Results on Reference Soils[END_REF]- [START_REF] Panagiotopoulos | The influence of clay on the threshold of movement of fine sandy beds[END_REF]. A change in behavior is observed around 40-50% of fines content, in line with a previous study by Panagiotopoulos et al. 1997 [13] which revealed that the increase in mud content results in two different regimes of increase for the threshold shear stress, with 30% mud content being the frontier between the two regimes.

Figure 1 :

 1 Figure 1 : Photograph of the hole erosion test apparatus.

Figure 2 :

 2 Figure 2 : Particle Size Distribution curves for the tested materials.

Figure 3 :

 3 Figure 3 : Effect of coarse particles shape on the erosion parameters, critical shear stress 𝜏 𝑐 and erosion coefficient 𝐾 𝑒𝑟 for three types of mixtures.

Table 1 :

 1 Erosion parameters of soil samples obtained from HET tests.

	Sample n°	Test ID	𝑾(%)	Compaction degree (%)	𝝆 𝒅 (t/𝒎 𝟑 ) 𝝉 𝒄 (Pa)	Ker (s/m)
		Mixture n°1 : 50% Armoricaine Kaolinite (KA) + 50% Glass Beads (GB0.4/0.8)
	1	KA50GB50-88%PROCTOR-15.6%	15.6	88		14	11.3
	2	KA50GB50-92%PROCTOR-15.6%	15.6	92	1.58	17	7.3
	3	KA50GB50-95%PROCTOR-15.6%	15.6	95	1.63	22	2.2
	4	KA50GB50-98%PROCTOR-15.6%	15.6	98	1.68	29	2.8
	5	KA50GB50-100%PROCTOR-15.6%	15.6	100	1.72	34	2.4
		Mixture n°2 : 50% Armoricaine Kaolinite (KA) + 50% Loire sand (TEN0.5)
	6	KA50TEN50-88%PROCTOR-16.6%	16.6	88	1.51	8	4.5
	7	KA50TEN50-92%PROCTOR-16.6%	16.6	92	1.58	9	2.3
	8	KA50TEN50-95%PROCTOR-16.6%	16.6	95	1.63	13	2.9
	9	KA50TEN50-98%PROCTOR-16.6%	16.6	98	1.68	27	0.9
	10	KA50TEN50-100%PROCTOR-16.6%	16.6	100	1.72	35	1.1
		Mixture n°3 : 50% Armoricaine Kaolinite (KA) + 50% Hostun sand (HN0.4/0.8)
	11	KA50HN50-88%PROCTOR-15.6%	15.6	88	1.50	7	5.5
	12	KA50HN50-92%PROCTOR-15.6%	15.6	92	1.57	7	3.1
	13	KA50HN50-95%PROCTOR-15.6%	15.6	95	1.62	9	2.3
	14	KA50HN50-98%PROCTOR-15.6%	15.6	98	1.67	17	1.2
	15	KA50HN50-100%PROCTOR-15.6%	15.6	100	1.71	19	1.4

Table 2 :

 2 Erosion parameters of soil samples from HET tests.

	Sample n°	Test ID	% Clay 𝑾 (%)	Compaction degree (%)	𝝆𝒅 (t/𝒎 𝟑 ) 𝝉 𝒄 (Pa)	Ker (𝒔/𝒎)
	1	KA20HN80-95%PROCTOR-11.9%	20	11.9	95	1.82	3	26.5
	2	KA30HN70-95%PROCTOR-12.5%	30	12.5	95	1.80	7	9.7
	3	KA40HN60-95%PROCTOR-14%	40	14	95	1.73	10	8.4
	4	KA50HN50-95%PROCTOR-15.6%	50	15.6	95	1.62	9	2.3
	5	KA70HN30-95%PROCTOR-22.6%	70	22.6	95	1.47	13	0.7
	6	KA90HN10-95%PROCTOR-28.2%	90	28.2	95	1.32	13	0.6
	7	KA100HN0-95%PROCTOR-32%	100	32	95	1.23	15	0.6