Nikolay Sirakov

Bénédicte Fontez
email: benedicte.fontez@supagro.fr

Thérèse Libourel Rouge

Alessio Dos Santos

Danielle Mitja

Patrice Loisel

Thérèse Libourel

A stage-structured hierarchical Bayes model for the babassu palm tree population dynamics -Estimated from anthropogenic open area data sets

Keywords: mathematical modelling, Bayesian statistics, population dynamics, babassu palm tree

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The primary threat to primary forests, long before climate change, is human action, [START_REF] Jha | World's forests can adapt to climate change, study says[END_REF]. The majority of plant species that grow in tropical rainforests do not resist land use change at the time of cultivation, [START_REF] Fearnside | Deforestation in brazilian amazonia: History, rates, and consequences[END_REF]; [START_REF] Barlow | Anthropogenic disturbance in tropical forests can double biodiversity loss from deforestation[END_REF]. Some other tree and palm species are resilient, [START_REF] Anthelme | How anthropogenic disturbances affect the resilience of a keystone palm tree in the threatened andean cloud forest?[END_REF]; [START_REF] Norden | Resilience of tropical rain forests: tree community reassembly in secondary forests[END_REF]. The babassu palm tree (Attalea speciosa Mart. ex Spreng.) which can reach 30 m height, is endemic to the South American Amazonian forest. It is sometimes present in agrosystems installed following clearing of the primary forest, while many other woody species have disappeared, [START_REF] Mitja | Satellite images combined with field data reveal negative changes in the distribution of babau palm after clearing off amazonian forest[END_REF]; dos Santos and [START_REF] Dos Santos | Wooded Cattle Pasture in the Benfica Seetling Project in Itupiranga, Para, Brazil[END_REF]. It persists years later in secondary forests resulting from the abandonment of crops and pastures, and it can even become invasive and hinder the development of cultivated species. While a portion of the population wants to eliminate it, another part wants to maintain it because of its multiple uses, [START_REF] Araújo | Diversity of use and local knowledge of palms (Arecaceae) in eastern Amazonia[END_REF]. Some plot-scale work has demonstrated higher densities of adult babassu in pastures than in the original forest, [START_REF] Anderson | The biology of Orbignya Martiana (Palmae), a tropical dry forest dominant in Brazil[END_REF]; [START_REF] Campos | Population structure and fruit availability of the babassu palm (attalea speciosa mart. ex spreng) in human-dominated landscapes of the northeast region of brazil[END_REF], while other work has shown the opposite, [START_REF] Barot | Reproductive plasticity in an Amazonian palm[END_REF]. A territory-wide study of the PA-Benfica smallholder farming community revealed this variability in the density of babassu cauliflower palms in primary forests as well as in pastures and secondary forests. In these environments, the few plots with very high densities are only observed in pastures or secondary forests, [START_REF] Mitja | Satellite images combined with field data reveal negative changes in the distribution of babau palm after clearing off amazonian forest[END_REF]. In contrast, because the number of plots without cauliflower babassu is higher in pastures, the average babassu density in the study area is lower in pastures than in secondary or primary forests, [START_REF] Mitja | Satellite images combined with field data reveal negative changes in the distribution of babau palm after clearing off amazonian forest[END_REF]. Thus, babassu has the potential to invade agrosystems even if this does not occur systematically. In fact, there is a high variability of the palm density in pastures potentially that is related to the initial density in the forest, depending on the cultivation practices of farmers (who can conserve or eliminate palm trees), [START_REF] Mitja | Satellite images combined with field data reveal negative changes in the distribution of babau palm after clearing off amazonian forest[END_REF]. In pastures, babassu are influenced by factors such as sun exposure, [START_REF] Mitja | Establishement of Babassu in Pastures in Pará, Brazil[END_REF]; [START_REF] Anderson | The Subsidy from Nature: Palm Forests, Peasantry, and Development on an Amazon Frontier[END_REF], burning, chemical or manual cleaning of plots, competition with forage grasses, burial of fruit in the soil and clear-cutting. These factors sometimes act in a contradictory way on fruits, seedlings and young individuals and adults and can counteract the high resilience of this species at certain stages, [START_REF] Mitja | Satellite images combined with field data reveal negative changes in the distribution of babau palm after clearing off amazonian forest[END_REF]. In our study context, the studied pastures did not undergo significant anthropogenic actions such as chemical cleaning or clear-cutting.

To explain this ability of babassu to remain on pasture after forest clearing, the possibility of natural biological diversity within babassus (diversity of genetic, epigenetic origins or more generally due to genetic interaction with the environment) cannot be excluded, [START_REF] Bräutigam | Epigenetic regulation of adaptive responses of forest tree species to the environment[END_REF]; [START_REF] Richards | Ecological plant epigenetics: Evidence from model and non-model species, and the way forward[END_REF], although we do not currently have data to support these assumptions.

Previous works, focusing on all babassu individuals from seedlings to adults, have examined the life cycle divided into ontogenic stages: three stages for [START_REF] Campos | Population structure and fruit availability of the babassu palm (attalea speciosa mart. ex spreng) in human-dominated landscapes of the northeast region of brazil[END_REF]; four stages for [START_REF] Anderson | The biology of Orbignya Martiana (Palmae), a tropical dry forest dominant in Brazil[END_REF]; five stages for [START_REF] Barot | Reproductive plasticity in an Amazonian palm[END_REF]; six stages for dos [START_REF] Dos Santos | Estimating babassu palm density using automatic palm tree detection with very high spatial resolution satellite images[END_REF]. In the latter case, which concerns our field of study, the life cycle was divided into a seedling stage, four stages of young individuals and an adult stage. Stages 1, 2 and 3 are characterized by a rosette configuration and a buried terminal bud. Stage 4 corresponds to the presence of the terminal bud on the soil surface, protected by the leaf sheaths. Stage 5 (last young stage) corresponds to individuals with a terminal bud above ground that is very often still protected by the leaf sheath. Stage 6 is the adult stage, most often ascertained by a free stipe. Considering the stages, some authors have studied the population structure using a single "one shot" measure, [START_REF] Barot | Reproductive plasticity in an Amazonian palm[END_REF]; [START_REF] Campos | Population structure and fruit availability of the babassu palm (attalea speciosa mart. ex spreng) in human-dominated landscapes of the northeast region of brazil[END_REF]. [START_REF] Barot | Reproductive plasticity in an Amazonian palm[END_REF] obtained an inverted J-shaped curve for the primary tropical forest, characteristics of natural unharvested populations, [START_REF] Hall | Methods to assess the impact of extraction of nontimber tropical forest products on plant populations[END_REF], while in pastures this structure was found to be different with a lower seedling density. These latter results have not been verified by studies by [START_REF] Campos | Population structure and fruit availability of the babassu palm (attalea speciosa mart. ex spreng) in human-dominated landscapes of the northeast region of brazil[END_REF] and [START_REF] Anderson | The biology of Orbignya Martiana (Palmae), a tropical dry forest dominant in Brazil[END_REF], who showed an "inverted J" curve not only for primary forest but also for pastures. However, the accepted interpretation that species with inverted J population structures are in equilibrium, while those with a small number of young stages would correspond to declining populations, has been questioned. Indeed, the survival of the plant population would depend more on the growth rate of individuals in the various stages, [START_REF] Condit | Predicting population trends from size distributions: A direct test in a tropical tree community[END_REF], which, if rapid, can maintain a sufficient density of adults. To determine whether the babassu will persist in this newly set-up farmland in the long term, it will be necessary to assess the growth rate and thus monitor the evolution of the population over time, [START_REF] Johnson | The hazards of interpretation of static age structures as shown by stand reconstructions in a pinus contortapicea engelmannii forest[END_REF] by studying the population dynamics through periodic measurements over the long term.

Recently, individual-based models have received great interest and have been used to predict population dynamics and the effects of management strategies [START_REF] Buckley | Demography and management of the invasive plant species Hypericum perforatum. II. Construction and use of an individual-based model to predict population dynamics and the effects of management strategies[END_REF]. Following the ideas of Uchmaski and Grimm (1996) and [START_REF] Lomnicki | Individual-based models and the individual-based approach to population ecology[END_REF], the interest in individual-based models is to take the individual interactions (competition) and the resource dynamics explicitly into account. In our case, pastureland is a stable and non limited resource environment where competition between trees is negligible. An individual-based model could be interesting to ascertain the time spent at each stage of development and to define trajectories, but as the data are from a three-year experiment, the results may not be sufficiently extensive to be reliable. Because of the very low transition rates often observed between some stages, a longer observation period would be necessary to be of use. The population dynamics models appear to offer a simpler and better adapted framework for our study Zuidema and Zagt (2000). Two components must be taken into account when modelling: the ecological process and the observation process. In our context, the ecological process describes variation in abundance over transects and time (two periods of measurement).

Demographic models represent a standard approach for abundance data that provides insight into population growth [START_REF] Caswell | Analysis of life table response experiments I. Decomposition of effects on population growth rate[END_REF]. Deterministic matrix models have been used for trees in the following articles: [START_REF] Bullock | Demography of an Undergrowth Palm in Littoral Cameroon[END_REF][START_REF] Pinero | A Population Model of Astrocaryum Mexicanum and a Sensitivity Analysis of its Finite Rate of Increase[END_REF], [START_REF] Enright | Population dynamics of the nikau palm, Rhopalostylis sapida (Wendl. et Drude), in a temperate forest remnant near Auckland, New Zealand[END_REF], [START_REF] Pinard | Impacts of Stem Harvesting on Populations of Iriartea deltoidea (Palmae) in an Extractive Reserve in Acre, Brazil[END_REF], [START_REF] Olmsted | Sustainable Harvesting of Tropical Trees: Demography and Matrix Models of Two Palm Species in Mexico[END_REF], [START_REF] Barot | Demography of a Savanna Palm Tree in Ivory Coast (Lamto): Population Persistence and Life-History[END_REF], [START_REF] Holm | Population Dynamics of the Dioecious Amazonian Palm Mauritia flexuosa: Simulation Analysis of Sustainable Harvesting[END_REF], and also [START_REF] Guarín | Modeling the stipe growth of the Oenocarpus bataua palm in the Central Cordillera of the Andes, Colombia[END_REF]. Demographic models include a deterministic part but stochastic elements can be added to accommodate any sources of variability (transects, time, individuals) that affect the observations. To take into account the geographical structure of our data, palm trees are observed on different transects; thus defining the hierarchical model. This approach was advocated by [START_REF] Clark | Uncertainty and Variability in Demography and Population Growth: A Hierarchical Approach[END_REF] who insists that variability associated with space, time, and among individuals who are not accommodated by demographic models can make parameter estimates and growth rate predictions unrealistic. Stochastic models can be defined in different ways. In the method defined by Tuljapurkar and Caswell (1997), the stochastic matrix involves elasticity and stochasticity of the vital rate. In our case, we did not focus on one parameter only but on a global approach that allows one to deal with broader issues and topics from different fields (economy, social impact, biology, technical aspect of the management and so on).

IPM (integrated population models) is another possible approach, which is also known as hierarchical or state-space models [START_REF] Eacker | Assessing the importance of demographic parameters for population dynamics using Bayesian integrated population modeling[END_REF], [START_REF] Schaub | Integrated population models: a novel analysis framework for deeper insights into population dynamics[END_REF]. More generally than IPM, the Bayesian hierarchical model is well adapted and sufficiently flexible to incorporate hierarchical data and expert knowledge, [START_REF] Hobbs | Bayesian Models -A Statistical Primer for Ecologists[END_REF]. It can be adjusted over successive years of observation and updated over time, posterior estimation is used from the first year as prior information for the second year estimation, and so on, as suggested by [START_REF] Parent | Introduction to Hierarchical Bayesian Modeling for Ecological Data[END_REF].

To understand how babassu, which is a forest species, adapts to grazing and under certain conditions becomes invasive, this study proposes to develop a population dynamics model that does not currently exist for this species. This model will be a hierarchical stochastic matrix population model in the Bayesian framework. Beyond the knowledge produced, this model can be used in future species management projects to make forecasts. This article covers the following topics:

1. "Materials" provides a description of the biology of the babassu, as well as a description of the study site and of the sampling.

"

Model" provides a comprehensive description of all aspects of the modelling (matrix model, hierarchy and Bayesian framework).

3. "Results-Discussion" provides comments on the model validation and on population dynamics and the existence of a bottleneck. These sections are followed by a global conclusion.

Materials

Biology of babassu in pastures

The babassu life cycle can be organized in biological stages of development.

Six stages have been defined varying from the seedling to the adult palm tree (Fig. 1, see [START_REF] Mitja | Viability of the babassu palm eco-sociosystem in Brazil: the challenges of coviability[END_REF] for more details). In addition, the reproductive stage and fructification of the babassu has been studied by [START_REF] Barot | Reproductive plasticity in an Amazonian palm[END_REF] and [START_REF] Anderson | The Subsidy from Nature: Palm Forests, Peasantry, and Development on an Amazon Frontier[END_REF] providing knowledge for understanding the reproductive process and recruitment in our modelling.

Study site

Our data was collected over three years of field monitoring (from 2013 to 2015). This monitoring was carried out on the Benfica Agriculture Installation Project (Projeto de Assentamento) site which is located in the Brazilian State of Pará (Fig. 2). This territory is located within the Amazonian rain forest on the territory of the Itupiranga municipality and covers 10,026 ha. The shape of Benfica is irregular and lies within a quadrilateral space formed by the straight line connection at the coordinates 05 • 12'20" south latitude, 05 • 20'40" south latitude and 49 • 48'00" west longitude, 49 • 56'40" west longitude.

The particularity of this site is that it includes man-made spaces that, since 1989, have been progressively colonized by one thousand people, Sampaio (2008).

Field sampling

The fieldwork protocol was developed in early 2012. Seven surface areas of pastureland were designated for this research project. Within each piece (2017). We obtained the total number of seventeen transects from which we collected data regarding growth, sprouting and fruiting processes. In addition, we also used seven supplementary transects to obtain more data about the growth process, especially for stages 5 and 6 (Table 1). The counting and marking of all individuals at all six stages of development was performed within these transects. In 2013, 2,987 palm trees of all stages were identified. In 2014 and 2015, we surveyed the same palm trees, to which we respectively , added 471 and 482 newly recruited individuals of stage 1 (Table 2).

Model

The modelling was constructed based on three distinct experimental data sets (see flowchart Fig. 3) that were observed during two observation periods We decided to infer three independent models based on three different observational units for the following reasons:

• The transition-survival model was our main concern for answering demographic questions such as the existence of a bottleneck or a reverse J-shape curve. The survey was constructed for that purpose. The sample land was chosen to be representative of the different ages of pasture in Para. The 24 transects were chosen in order to sample a sufficient number of babassu • The fruiting and sprouting processes were less precisely observed and only on the 17 transects. We wanted to avoid error propagation of the measurement or modelling approximations to the transition-survival modelling.

Therefore, we decided to construct independent fruiting and sprouting models. Moreover, estimation of these models on the 17 transects might not be adapted to the 7 specific transects of older palm trees. However, even if the three models were independently estimated, they could be combined to simulate a life cycle.

• The survey did not focus on the seed dispersal (even if it exists, see [START_REF] Mitja | Establishement of Babassu in Pastures in Pará, Brazil[END_REF] for more detailed) or the germination. No information was collected on the germination of fruits or seedlings. No information was available for outside the transects, (linear paths). Without this information, it was not possible to fit and use seed dispersal models such as the one proposed by [START_REF] Clark | Seed dispersal near and far: Patterns across temperate and tropical forests[END_REF]. Localization of adult palm trees around or adjacent to the transect would be necessary to achieve such modelling, and the new recruits were considered to be of unknown origin.

Brief description of the three independent models

• A demographic model with two levels (transect -studied area in Benfica) was defined. The observation were performed at an individual level, but the palm trees were considered to be independent, and the observations were pooled at the transect level for each stage. Demographic stochasticity was modelled with hierarchical Dirichlet (prior) and multinomial (likelihood) laws.

• The annual number of fruits from a palm tree was modelled with a log-Normal distribution to restrict them to positive values and to maintain flexibility in the meantime. A non-conjugated hierarchical model is used to model the annual mean number of fruits per palm tree across the studied area level (Normal prior) in addition to the observation/individual level (logNormal likelihood). This model is conditional on the probability that a palm tree can have fruits. This probability is modelled with a classical conjugate Beta-Binomial distribution.

• We used a classical conjugate Beta-Binomial distribution to model the probability of sprouting across the studied area. Recall that the sprouting process was observed at the transect level, from a number of pooled observed fruits (including contribution from Babassu palm trees besides/along and outside the transect).

Hierarchy

The ecological process was observed on transects belonging to the designated area of land. Only two levels were considered for the transition survival process: the transect and the whole studied area in Benfica. Moreover, our estimation concerning the probability of transition between stages (see Fig. 4) did not exhibit any structure at the level of the farmer's piece of land (each piece of land is separated by a red line). No spatial structure was observed, except for the transition from stage 2 to stage 3 where a link with spatial coordinates might be possible (see appendix and discussion) but it was difficult to estimate or ensure with our dataset. Therefore, we preferred to employ a multilevel model using the level of the transects and of the whole studied area in Benfica.

Observation process and definition of the observed variables

Survival-transition process

Observations were performed on 17 (+7 for stages 5 and 6) transects (see Table 1). For transect l, we observed the total number of plants in stage i at year t -1, n i,t-1,l . We also observed the number of survivors at year t in stage i, n i,t,l . This number was the sum of the observed number of survivors that remained in this stage, denoted as s i→i,t,l and of the observed number of survivors from stage i -1 who had moved to stage i, denoted as s i-1→i,t,l , n i,t,l = s i→i,t,l + s i-1→i,t,l . Moreover, in a given transect l at time t, the number of deaths among those in stage i was the total number of plants in stage i from the previous year minus the number of survivors from the previous year that remained in stage i minus the survivors from the previous year that changed stages, dead i,t,l = n i,t-1,l -s i→i,t,l -s i→i+1,t,l . From these observations, we could infer proportions for each transect l. We assumed that the expected proportions remained constant over time and denoted π i→i,l as the expected proportion of plants that remained at the same stage i, π i→i+1,l as the expected proportion moving from stage i to the next stage i + 1 and 1 -π i→i,l -π i→i+1,l as the expected mortality proportion in stage i. These three parameters allowed us to define a categorical distribution for each transect l:

π i,l = (π i→i,l , π i→i+1,l , 1 -π i→i,l -π i→i+1,l )
We also defined a categorical distribution at the level of the whole studied area:

π i = (π i→i , π i→i+1 , 1 -π i→i -π i→i+1 )
No observations performed at this level, therefore π i→i and π i→i+1 are latent variables.

Fruiting and sprouting processes.

Fruiting First, the fruit production per palm tree was observed on the 17 transects for each palm tree in stage 6. The reproduction process begins with the definition of the sexual state (no inflorescence, male or female inflorescence or both inflorescences). The babassu palm tree has a complex reproductive process that is not well understood and has not been well measured. The female and male inflorescences change over undefined periods of time.

Therefore the reproductive process has been simplified and resumed with the observed reproductive state variable u j,t,l , which indicates the presence of fruits at year t on palm tree j from transect l. This variable was coded as 1 when fruits were present and 0 if they were not. We assumed that the expected proportion φ of palm trees of stage 6 in a reproductive state 1 was constant over time and was common to all transects.

When a palm tree j of year t -1 was in a reproductive state 1, its number of fruits, denoted F ruits j,t-1,l , was the sum of the observed number of fruits on the palm tree and on the ground immediately under the palm tree. This observation was performed directly in the field but also validated later by using photography of the palm tree. We approximated fruit observations using a normal distribution with a mean parameter µ and a variance parameter σ 2 . Sprouting Next, the seedlings were not observed outside the transect, but only the number of recruited seedlings in stage one per transect l and per year t was measured in the field. We denoted this variable r t,l , but the origin of these new palm trees was unknown. They could represent fruit from an adult palm tree (stage 6) adjacent to the transect. Therefore, the total number of fruits on transect l, year t -1 was counted (regardless of their origin), and we denoted this variable T otalF t-1,l . We assumed that the recruitment occurred only in the first stage and that the expected proportion of fruit providing a new palm tree in one year was constant over time and was common to all transects. We denoted this parameter λ.

Remark: All variables were directly measured/observed on the ground by expert biologists. We assumed that there were no errors and our observation process was identity.

Lefkovitch Matrices for the ecological process

Six growth stages have been defined for babassu [START_REF] Mitja | Viability of the babassu palm eco-sociosystem in Brazil: the challenges of coviability[END_REF]. Vector n t,l = (n 1,t,l , . . . , n 6,t,l ) is the state vector, where n i,t,l is the number of individuals in stage i at year t in transect l. The ecological process is modelled with a Markovian process:

E      n 1,t,l . . . n 6,t,l n 1,t-1,l . . . n 6,t-1,l      =               π 1→1,l 0 0 0 0 0 π 1→2,l π 22,l 0 0 0 0 0 π 2→3,l π 3→3,l 0 0 0 0 0 π 3→4,l π 4→4,l 0 0 0 0 0 π 4→5,l π 5→5,l 0 0 0 0 0 π 5→6,l π 6→6,l                    n 1,t-1,l . . . n 6,t-1,l      + E      r t,l . . . 0      (1)
The projection of the population was modelled at the transect level with a Survival-Transition matrix following [START_REF] Caswell | Matrix population models : construction, analysis, and interpretation[END_REF].

The model must take into account the demographic stochasticity and hierarchy. Integration of different kinds/sources of variability was performed in a Bayesian framework. In this framework, all parameters except σ 2 were considered as random variables (see the DAG-Directed Acyclic Graph, Fig. 5, for details). This parameter was fixed in the observation in an empirical Bayes approach to ensure MCMC converged because we omitted the observation process error (lack of measurements allowing estimation of the observation error apart from the variability in the number of fruits) and also due to the very small sample size. Three independent Bayesian sub models were defined:

Survival-Transition for the six stages

According to [START_REF] Buckland | Embedding Population Dynamics Models in Inference[END_REF], demographic stochasticity is modelled with multinomial laws (extension of the binomial law to a context with more than two possibilities). Fruiting Sprouting Survival -Transition, i = 1, . . . , 6

     s i→i,t,l s i→i+1,t,l dead i,t,l      | n i,t-1,l π i,l =      π i→i,l π i→i+1,l 1 -π i→i,l -π i→i+1,l      π i =      π i→i π i→i+1 1 -π i→i -π i→i+1      η i r t,l | T otalF t-1,l F ruits j,t-1,l u j,t-1,l φ µ λ
• Evolution from year t -1 to year t of the six stages of transect l:

(s i→i,t,l , s i→i+1,t,l , dead i,t,l ) | π i,l , n i,t-1 ∼ M ult(n i,t-1,l , π i,l ), i = 1, . . . , 5 (s 6→6,t,l , dead 6,t,l ) | π 6,l , n 6,t-1,l ∼ M ult(n 6,t-1,l , π 6,l ) (2)
To add to the hierarchy for each stage of the preceding model, we adapted a hierarchical Dirichlet distribution from [START_REF] Cowans | Information Retrieval Using Hierarchical Dirichlet Processes[END_REF]. or information retrieval from documents [START_REF] Haffari | Hierarchical Dirichlet Trees for Information Retrieval[END_REF].

• Hierarchical Dirichlet distribution to add to the hierarchy for each stage i

π i,l | π i , η i ∼ Dirichlet(η i × π i ), i = 1, . . . , 6
π i ∼ Dirichlet(W eights i ), i = 1, . . . , 6 (3)
η i ∼ U nif orm(u 1 , u 2 ), i = 1, . . . , 6
The shared mean parameter by stage i induces sharing of information across transects. W eights i of the Dirichlet hyperprior were calibrated with both observed proportions (as in an empirical Bayes approach) and uniform distribution (following [START_REF] Cowans | Information Retrieval Using Hierarchical Dirichlet Processes[END_REF] idea) to test the sensitivity of the model to this hyperprior. The latent variable η i could be considered to be a size effect which indicated the level of confidence in the realization of the latent vector π i .

Fruiting and number of fruits.

The annual number of fruits of a palm tree in stage 6 is conditioned by its reproductive status u j,t-1,l :

   F ruits j,t-1,l | u j,t-1,l = 1 ∼ logN ormal(µ, σ 2 ) F ruits j,t-1,l | u j,t-1,l = 0 ∼ δ 0 (4)
where δ 0 indicates the Dirac mass at 0. The logNormal distribution was preferred to the usual Poisson or negative binomial distributions to increase the flexibility of the modelling to address over dispersion and add explicitly to the population latent variable µ (the mean annual number of fruits inferred at the level of the studied area). The logNormal distribution provided the best fit our data (see supplementary materials, Fig. 10). Next, we assumed that the palm trees in stage 6 shared the same parameter of fruiting φ.

u j,t-1,l ∼ Bernoulli(φ) (5)
To add to the hierarchy allowing the sharing of information across palm trees (and transects), we defined:

φ ∼ Beta(a, b) µ ∼ N (m, V ) (6)
Sprouting (or recruitment in the first stage).

Recruitment in the first stage r t,l was stochastic and was defined Conditional on the total number of fruits present on a transect l, year t -1 :

r t,l | T otalF l,t-1 ∼ Binomial(T otalF l,t-1 , λ) (7)
To add to the hierarchy allowing the sharing of information across palm trees (and transects), we defined:

λ ∼ Beta(c, d) (8)

Hyperparameters

Two sets of hyperparameters were used, with one per year of estimation.

For the first year, the prior was less informative as no expert knowledge was available. For the second year, we used the marginal posterior mean value and variance of the variables (φ, λ, π i ) to calibrate the hyperparameters. For both years the hyperparameters m, V and σ 2 were calibrated from observed field data.

The values are presented in Table 3.

Weights 1 2 3 4 5 6

Non info (1/3,1/3,1/3) . . . . . . (1,200)

(1/3,1/3,1/3) (1/2,
Table 3: Hyperparameter values for non-informative and informative priors (informative:

hyperparameters were calibrated with the results of the first year posterior estimation)

Bayesian inference and model validation

In Bayesian frameworks, estimation of the variables of interest is performed for the posterior distribution. In our case, the variables of interest were for transition-survival θ 1 = ((π 1,l , . . . , π 6,l ) l=1,...,24 , π 1 , . . . , π 6 , η 1 , . . . , η 6 ), recruitment in the first stage θ 2 = λ and fruiting θ 3 = (φ, µ). Whether the prior distribution of each θ p is conditional/depends on the hyperparameter values. The vectors of those hyperparameters are param 1 = (W eights 1 , . . . , W eights 6 , u 1 , u 2 , n t-1 ), param 2 = (c, d, T otalF t-1,1 , . . . , T otalF t-1,17 ) and param 3 = (a, b, m, V, σ 2 ).

Each vector is considered fixed, and its values are given in Table 3. Summaries for (n t-1,l ) l=1,...,24 and the histogram for (T otalF t-1,l ) l=1,...,17 , are given in Table 2 and Supplementary Material, Fig. 10, respectively. Observations or data respectively are the realization of the following variables:

data 1 = (s 1→1,t,l , s 1→2,t,l , . . . , s 6→6,t,l ) l=1,...,24 , data 2 = (r t,l ) l=1,...,17 and data 3 = (u j,t-1,l , F ruits t-1,l ) j∈transect l ,l=1,...,17 .

Based on those definitions, from Bayes' rule, the posterior distribution is proportional to

p(θ p | data p , param p ) ∝ p(data p | θ p , param p ) p(θ p | param p ) (9)
To simplify the notation, we will omit the dependence on param in the following equations. Of course, we are clearly aware that posterior estimation is conditional on param and therefore sensitive to the values of the hyperparameters, and this point will be discussed and studied later.

Using the DAG (Fig. 5) and conditional dependencies, Equation ( 9) can be decomposed into a product, depending on the type of observation: Transition-Survival, Fruiting and Sprouting:

p(θ 1 | data 1 ) ∝ p(Survival data | (π 1,l , .
. . , π 6,l ) l=1,...,24 , π 1 , . . . , π 6 , η 1 , . . . , η 6 )× p((π 1,l , . . . , π 6,l ) l=1,...,24 , π 1 , . . . , π 6 , n 1 , . . . , n 6 ) =

Π i Π l p(s ii,t,l , s i i+1,t,l , dead i,t,l | π i,l ) p(π i,l | π i , η i ) p(π i ) p(η i ) p(θ 2 | data 2 ) ∝ p(Sprouting data | λ) × p(λ) = Π l p(r t,l | λ) p(λ) p(θ 3 | data 3 ) ∝ p(Fruiting data | µ, φ) × p(µ, φ) = Π l Π j∈transect l p(F ruits j,t-1,l | µ, u j,t-1,l ) p(u j,t-1 | φ) p(φ) p(µ)
Sensitivity to Informative versus Non-Informative priors in a learning approach Estimation based on the first year of observation was performed according to [START_REF] Helser | Bayesian hierarchical modeling of Pacific geoduck growth increment data and climate indices[END_REF] since no study to date has provided prior information on babassu growth or reproduction processes in the babassu environment. We endorsed the approach of [START_REF] Smith | The hierarchical Bayesian approach to population pharmacokinetic modelling[END_REF], who suggest that priors should be non-informative to ensure that the likelihood dominates the prior.

The hyperparameters characterizing the distribution of π i , η i , φ, λ and µ were chosen to provide a prior that was as flat and non-informative as possible. Then, in a learning approach, the value of the hyperparameters W eights i was set to the posterior estimates of π i . We compared the estimation based on the second year of observation obtained with each set of hyperparameters. This comparison allowed us to study the sensitivity of the model and its stability over the years.

Model Validation

The predictive posterior distribution can be used to check whether the model is consistent with the observed data [START_REF] Gelman | Posterior predictive assessment of model fitness via realized discrepancies[END_REF] of the omitted transects. The total number per stage of the omitted transects was considered to be consistent with a prediction at the level of the studied area.

Computation with R and JAGS

JAGS which uses the dialect of BUGS Lunn et al. (2000), [START_REF] Plummer | JAGS: A program for analysis of bayesian graphical models using gibbs sampling[END_REF], [START_REF] Lunn | The BUGS book: A practical introduction to Bayesian analysis[END_REF] was used to fit the model to the data with the package rjags.

The estimates of parameters were evaluated based on 200,000 samples after 50,000 burn-in samples, thinning to one draw for each sample, by MCMC simulation of the joint posterior distribution. The MCMC inference samples analysis was processed with the coda package [START_REF] Plummer | CODA: convergence diagnosis and output analysis for MCMC[END_REF]. Convergence diagnostics consisted of visuals (traces and autocorrelations) and statistical calculations using the test of [START_REF] Gelman | Bayesian Data Analysis[END_REF] based on [START_REF] Gelman | Evaluating and using statistical methods in the social sciences -A discussion of "A critique of the Bayesian information criterion for model selection[END_REF]. For all variables in the model, the Gelman criteria of between-chain variance over within-chain variance converged into one. More information on the convergence diagnosis and on the script of the model are available in the supplementary materials.

Results and Discussion

Model validation

We focused the validation analysis on the Survival-Transition matrix at the level of the studied area. The validation concerned the second year of observation 2014-2015 with the informative prior.

First, we computed the posterior predictive p-value for estimation of the survival-transition probabilities obtained for the second year of observation (2014)(2015) and with the informative prior (see Table 4). When the ppp-value was much smaller than 0.5 (e.g., 0 or 0.1), then the probability simulated with the model was less variable than the observed probability, and the model failed to reproduce the variability of the data. Thus, the model underestimated the data variability. We imputed this phenomenon as a higher variability observed in the data during 2014-2015 than during 2013-2014, which may indicate a year effect. A longer time series of observation may well be needed to define an informative prior, as suggested by [START_REF] Bierzychudek | Looking Backwards: Assessing the Projections of a Transition Matrix Model[END_REF].

ppp-values

Stage Next, we checked the adequacy of the model with the cross-validation approach (see Table 5) for the second year of observation 2014-2015 and with the informative prior. The young stages (namely one and two) showed a poor prediction, as the observed values were not in the 95% credible interval of the prediction for all variables. Again, this mismatch for the young stages was a result of the extreme variability observed in the field for the dynamics of the young during the two periods of observation (2013-2014 and 2014-2015).

i 1 -π i→i -π i→i+1 π i→i π i→i+1 1 0 0.
dead i,
3 0 0 [0,2] 83 82 [79,83] 0 1 [0,3] 4 0 0 [0,1] 3 4 [2,5] 2 1 [0,3] 5 0 0 [0,1] 8 9 [6,10] 2 1 [0,4] 6 1 1 [0,3] 58 58 [56,58] - -
Table 5: Cross-validation for the year 2015 and for each stage i = 1, . . . , 6, the observed (obs.) and predicted (pred.) numbers of palm trees that died (dead i,2015 ), survived in the same stage (s i→i,2015 ) and survived and evolved toward the upper stage (s i→i+1,2015 ). The credible interval at a 95% level was computed for each prediction.

Eventually, the learning approach of using the posterior of the first period 2013-2014 to calibrate the prior of the second period 2014-2015 was not satisfactory for the young stage (namely 1 and 2).

Survival and transition between states: identification of a bottleneck

Table 6 summarizes the estimates of the survival-transition probabilities for both periods (2013-2014 and 2014-2015). The dynamic can be divided into two parts. Before compared with after stage 3, mortality was significantly elevated (5% to 27%). The first stage exhibited a high year variability, as underlined in the model validation section. This competition disappeared at a later stage of palm tree growth. The death probability for the first stage was estimated at 12% for the first period and 27% for the second period. These values were both lower than the estimate of 38% published by [START_REF] Anderson | The biology of Orbignya Martiana (Palmae), a tropical dry forest dominant in Brazil[END_REF] and [START_REF] Anderson | Taxonomy of the babassu complex (Orbignya spp. palmae)[END_REF] concerning another experiment in Brazil. After stage 3, the estimates were not sensitive to the prior or the year effect. The mortality rates were lower and difficult to estimate with small sample sizes. Except for the transition between stage 3 and stage 4 (lower than 1%), the probability of transitions were significant (10% to 28%). Hence, these results showed a clear bottleneck at stage 3. Due to its extremely low probability, the transition between stage 3 and 4 had a large impact on the "needed" time for a babassu to reach the adult stage 6. Thus, a supplementary sampling effort should be performed for stages 3 and 4 to better estimate and understand this transition.

Stage 3 is critical in the biology of the babassu because the palm terminal bud rises from the ground during this stage.

Link between the bottleneck and population structure

Finally, we observed a bottleneck at the end of stage 3 in pastures. This bottleneck was observed in a primary forest by [START_REF] Barot | Reproductive plasticity in an Amazonian palm[END_REF]. Moreover, a comparison of our bar plot of densities by stage (see Fig. 6) revealed a shape similar to [START_REF] Barot | Reproductive plasticity in an Amazonian palm[END_REF] for stages 1 to 3 (corresponding to ELS, SLS and Juv1, respectively, in their article); however in the next stage situated after the bottleneck, the densities are collapsed toward zero in our case. We expected, a priori, high densities in all young stages for this palm specie,s which is considered to be an invasive plant species in pastures, [START_REF] Barot | Reproductive plasticity in an Amazonian palm[END_REF]. A pasture is an open field/environment with less competition for light or with other species. In open areas, the number of germinations per fruit is lower but growth in stage 1 is more rapid, [START_REF] Mitja | Establishement of Babassu in Pastures in Pará, Brazil[END_REF].

Fig. 6 shows an increase in density in stage 3 over the time. In addition, the mortality rates for stages 3 to 6 were very low during this period. Therefore, we may assume that the observed bottleneck was due to the low transition rate from stage 3 to stage 4. Multiple causes could have contributed or explained this low transition.

First, competition for light and nutrients with fodder plants or introduced grasses might not be negligible in the younger stages. As stated in the article by [START_REF] Sankaran | Tree-grass coexistence in savannas revisited-insights from an examination of assumptions and mechanisms invoked in existing models[END_REF], "tree grass coexistence in savannas", competitionbased models or demographic-bottleneck models, or a mixture of the two, could explain the population structure of babassu in the pasture.

Second, the management of the pasture, with human intervention, has an impact on what is believed to be an appropriate density of babassu palm trees for cattle shelter. [START_REF] Runk | Productivity and sustainability of a vegetable ivory palm (phytelephas aequatorialis, arecaceae) under three management regimes in northwestern ecuador[END_REF] observed this kind of human intervention in the regulation of Ivory palm tree, where "the lack of juveniles to the sub adult stage appears to be a population bottleneck in managed stands, limiting the number of adults". However, during the period of our study, management was limited to an annual clearing to respect our research on the population dynamics and estimation of the natural mortality. Annual clearing of pasture plots by cutting leaves from stages 2 to 4 and sometimes even stage 5 can delay the transition between stages because the palm tree must reinvest in the production of new leaves, [START_REF] Mitja | Satellite images combined with field data reveal negative changes in the distribution of babau palm after clearing off amazonian forest[END_REF].

Third, the impact of cattle (for example, trampling, grazing) might provide an explanation for a demographic bottleneck model. The terminal meristem of the plant is below ground until stage 3, and the plant can regenerate even if the aerial part is damaged (grazed, etc.).These non-lethal actions (as we estimated a very low mortality in stage 3) might delay the transition and thus create a bottleneck for the transition from stage 3 to stage 4.

Fourth, none of these suggestions explain why this bottleneck or at least a reverse J-shaped curve (with a lower density in stage 4) was observed in primary rain forests. In addition to the competition for light with other trees (or babassus), the population structure could be due to different times of residence in the different stages. The prolonged juvenile phase of trees is obviously a central component of their life strategy. [START_REF] Brunner | Revisiting tree maturation and floral initiation in the poplar functional genomics era[END_REF] have suggested a different mechanism (mainly related to chromatin-based repression) to explain this prolongation. One hypothesis is that epigenetic mechanisms play a central role in the regulation of the tree phases. [START_REF] Bräutigam | Epigenetic regulation of adaptive responses of forest tree species to the environment[END_REF] have explored the epigenetic variation and adaptive capacity of plant species, including forest trees. Considered in this way, the prolonged juvenile phase might be a contributory strategy to the babassu's ability to colonize or persist in pastures.

Other biological information: fruiting and sprouting

The fruiting and sprouting parts of the model were simplifications of the complex reality. Nevertheless, rate estimation of fruit production (55%) and sprouting (12%) exhibited good statistical properties, and these estimations were enriched by the Bayesian learning period. We suggest that these processes (fruiting and sprouting) remained unchanged during the two periods of observation.

Fruiting

The mean estimate values ± the standard error of the fruiting rate (φ) were 0.52 ± 0.10 for 2013 and 0.54 ± 0.07 for 2014. The estimation for 2014 was performed with an informative prior using a Bayesian learning approach (the posterior of the first period was used as the prior for the second period).

Biologically speaking, these values indicate that a little more than fifty percent of the adult palm trees contribute annually to the fruiting process.

The mean estimate values ± the standard error of the number of fruits per palm tree (µ) respectively were 313 ± 50 for 2013 and 317 ± 48 for 2014.

Sprouting

The mean value estimates for the sprouting probability λ were approximately ≈ 0.12 for both periods, with standard error of 0.005 for the first period 2013-2014 and 0.004 for the second period 2014-2015. λ represents the probability that a fruit produced in year t becomes a seedling (stage 1) in the next year t + 1.

Bayesian Learning

Fruiting/sprouting

The estimation for the first period 2013-2014 was performed with a noninformative prior. A flat distribution was used as illustrated in blue, for the fruiting rate (φ ∼ Beta(a = b = 1.5)), Fig. 7a. The distribution obtained a posteriori (in red) was concentrated on the estimated mean value of 54%. This figure illustrated the Bayesian learning approach with the evolution between the flat distribution (non-informative prior) to the conditional distribution (knowledge of the information provided by the data). This informative distribution (in red, Fig. 7a) obtained during the first period of observation was then used as the prior (in blue, Fig. 7b) for estimation of the fruiting rate during the second period of observation, 2014-2015. Again, the distribution obtained a posteriori (in red) exhibited a concentration around the same mean estimated value. Bayesian learning helped to increase the precision of the estimation for the fruiting rate.

We observed similar figures for the sprouting rate and the number of fruits per palm tree. Concerning transition-survival probabilities, this learning approach was not so fruitful, because of important variability-fluctuation between the two periods of observation mainly during the early stages (1 and 2).

Conclusion and perspectives

To better understand the babassu population dynamics in pastures, we used a Bayesian hierarchical matrix model. This model with two levels of hierarchy (studied area in Benfica and transect) allowed a more reliable estimation of the population dynamics because information is shared between transects to improve the estimation at the level of the studied area and then retrospectively at the transect level. This is a significant improvement over direct empirical estimations using small sample sizes, with missing data and an unbalanced design, [START_REF] Fontez | Bayesian hierarchical model used to analyze regression between fish body size and scale size: application to rare fish species zingel asper[END_REF]. In our case, the observed stages in transects were unbalanced, as shown in Fig. 4. The model did not reveal unobserved empirical facts but provided a direct and reliable estimation at the level of the all the palms have an equal probability of passing from stage 3 to stage 4 each year; 2. only a portion of the palm trees are able to undergo this transition.

• In the first case, all the palm trees are impacted by the bottleneck and a rough estimation of the period in stage 3 with a geometrical model could be π 3→4 /(1 -π 3→3 ) 2 = 0.005/(1 -0.99) 2 ≈ 50 years. This bottleneck dramatically raised the age at which a palm tree becomes an adult. This result contradicts information from experiments (some palm trees became adults only after 10-12 years).

• In the second case, a fraction of individuals f can successfully transit and a fraction 1 -f will never transit. Therefore, π 3→4 is a compromise (a mean value) between those two kinds of palm trees. The very low would make it possible to determine for an individual whether the past history of the passages from stage 1 to stage 2 and stage 2 to stage 3 affected the probability of passage from stage 3 to stage 4.

In open environments (such as pastureland), the number of seedlings per fruit is lower for babassu, [START_REF] Mitja | Establishement of Babassu in Pastures in Pará, Brazil[END_REF]. In the present work, we highlighted a high and year-dependant mortality rate at stage 1 and a bottleneck at stage 3. These three points highlight a new perspective maybe in contradiction with the common idea that the babassu is an invasive species, [START_REF] Barot | Reproductive plasticity in an Amazonian palm[END_REF]. These three points in conjunction with harvesting might explain the negative evolution of density in pastures, observed in [START_REF] Mitja | Satellite images combined with field data reveal negative changes in the distribution of babau palm after clearing off amazonian forest[END_REF], and the need for reliable preventive scientific information to help manage babassu population dynamics in pastures.

Our model could be extended to a larger scale (region, country) by incorporating palm density information from high spatial resolution satellite images dos [START_REF] Dos Santos | Estimating babassu palm density using automatic palm tree detection with very high spatial resolution satellite images[END_REF]. However, to achieve the sustainable management of the species, a more comprehensive study of the mechanisms underlying the low transition rate between stage 3 to 4 is needed. Some factors, such as competition for light and nutrients with fodder plants as well as information on the genetic diversity (or a potential epigenetic mechanism) should first be adressed.

Author's contributions

Nikolay Sirakov drafted this article, and his PhD thesis is on this work. B.

Fig. 1 :

 1 Fig. 1: Stage definition from Mitja et al. (2019)

Fig. 2 :

 2 Fig. 2: PA Benfica in Brazil. Magnification of the localization of transects within the farmer's piece of land. The localization of some transects (green circle) are superimposed.

Fig

  Fig. 3: Flowchart of the whole process

Fig. 4 :

 4 Fig. 4: Estimation of the probabilities of transition from one stage to the next per transect (transects are grouped by farmer's pieces of land, with each separated by an underlying red line)

Fig. 5 :

 5 Fig. 5: Representation of the hierarchy and dependences between variables with a Directed Acyclic Graph constructed for stage i and transect l. Three processes, fruiting, recruitment and survival-transition are represented from left to right. Observed variables are represented in grey.

  Cowans proposed a hierarchical Dirichlet document model where each document (transect for us) is modelled using a multinomial distribution whose parameters are given a Dirichlet prior. The mean value of the Dirichlet prior is itself assumed to be random and given a Dirichlet hyperprior. This kind of modelling with Dirichlet priors and Dirichlet hyperpriors is used for topic modelsWallach et al. (2009)

Fig. 6 :

 6 Fig. 6: Mean and standard error of the transect densities by stage and by year of observation.

Fig. 7 :

 7 Fig. 7: Estimation of the fruiting rate. Prior and posterior for both periods.

Fontez

  and P. Loisel provided the modelling and co-drafted the article. A. dos Santos, provided the field data and proof-read the article. D. Mitja drafted the article and helped to define the babassu stages and interpret the results. T. Libourel contributed to the modelling and implementation. lis Foundation and by the French National Institute of Agronomical Research (INRA) and Montpellier SupAgro. We thank all our partners who participated in the Open Science research project (especially É. Delaître and L. Demagistri from IRD and our Brasilian partners, with special thanks to I. de Souza Miranda and D. Deurival da Costa Carvalho) and E. Bonnier for providing language support, proof-reading and reformulation.

Fig. 8 :Fig. 9 :

 89 Fig. 8: Example of convergence for π 1 2 (2014-2015).

Fig. 11 :

 11 Fig. 11: Estimation of the probability of transition from stage 2 to stage 3 according to the coordinates of the pieces of land.

Table 1 :

 1 Design of the field sampling

	Biological	Stages	Nb of
	process		transects
		1 to 4	17
	Growth
		5 and 6	17 + 7
	Sprouting new stage 1	17
	Fruiting	6	17

Table 2 :

 2 Empirical probabilities and sample size for babassu in pastures

	Stage		2013-2014 2014-2015
		death	0.08	0.21
		idle	0.66	0.66
	1
		transition	0.26	0.13
		Sample size	736	956
		death	0.04	0.03
		idle	0.78	0.87
	2
		transition	0.18	0.1
		Sample size	1737	1548
		death	0.006	0
		idle	0.988	0.98
	3
		transition	0.006	0.02
		Sample size	340	638
		death	0	0
		idle	0.79	0.76
	4
		transition	0.21	0.24
		Sample size	19	21
		death	0	0
		idle	0.83	0.87
	5
		transition	0.17	0.13
		Sample size	52	47
		death	0.01	0.01
	6	idle	0.99	0.99
		Sample size	103	111

(one in

2013-2014 and one in 2014-2015)

on 17 or 24 transects.

  measurement, y obs is the observed value and y rep is a replicate sampled from the predictive posterior p(y rep | θ) p(θ | data) dθ Steinbakk and Storvik (2009). Posterior p-values of approximately 0.5 indicate that the distributions of replicated and actual data are close, while values near zero or one indicate a discrepancy between them.Moreover, validation was possible by leaving out 25% of the transects for the second year of observation 2014-2015. The posterior population estimate π i was then used to predict the total expected number of survivors (s i→i , s i→i+1 , dead i )

Table 4 :

 4 Posterior predictive p-values of the survival-transition probabilities (π i , i = 1, . . . , 6) obtained for the second year of observation 2014-2015 and with the informative prior

	5	0.22

Table 6 :

 6 Sensitivity analysis of the prior. The given estimates are posterior mean values for the transition-survival probabilities.

	Probabilities	2014 non info. 2015 non info. 2015 info.
		Death	0.12	0.27	0.21
	π 1,t	Idle	0.6	0.52	0.57
		Transition	0.28	0.21	0.22
		Death	0.05	0.07	0.06
	π 2,t	Idle	0.79	0.83	0.81
		Transition	0.16	0.10	0.13
		Death	0.005	0.001	0.005
	π 3,t	Idle	0.993	0.994	0.989
		Transition	0.002	0.005	0.006
		Death	0.02	0.02	0.01
	π 4,t	Idle	0.77	0.75	0.77
		Transition	0.21	0.23	0.22
		Death	0.01	0.01	0
	π 5,t	Idle	0.85	0.87	0.85
		Transition	0.14	0.12	0.15
		Death	0.01	0.01	0.01
	π 6,t	Idle	0.99	0.99	0.99

Table 9: Posterior densities of the transition survival probabilities at the level of the whole studied area, respectively, over 2013-2014 on the left and over 2014-2015 on the right: death is shown in red, transition in green, and idle in blue. Stage 5 Stage 6

Acknowledgments

This PhD research work has been accomplished as a part of the Open Science research project (under the reference ID 1403-026) financed by the Agropo-

value estimated for π 3→4 could be due to this fraction 1 -f that never transited. According to [START_REF] Barot | Reproductive plasticity in an Amazonian palm[END_REF] fraction f that transits is probably composed of "juvenile with a high survival and growth rate to be quickly recruited into the adult stage and hence representative of a small percentage of the population". For this kind of juvenile, using a geometrical model, we deduced the expected time spent in the stages that are above soil and observable (4-6). Using estimates presented in Table 6 we obtained (3.68 to 4.16) for stage 4 and (6.22 to 7.1) for stage 5 for a total sum of these stages of 9.9 to 11.3 years. In the field, the people interviewed suggested 10 to 12 years. If we consider the stages that are above the soil and observable (4-6), the predictions are coherent with our model estimations.

The reality is undoubtedly between these two extreme cases. To arbitrate between these two hypotheses, follow-up sequences should be available in the future, which will soon be possible as new surveys are underway.

Concerning the second case, one possible explanation could be the high genetic variability of the babassu, which is a wild plant. Thus, important geneticenvironment interactions may be involved. A genetic difference for the f fraction of successful babassu is coherent with the observed increase in density in stage 3 over the years. At the population level, clearing of the forest and the use of fire to clear brush have resulted in a population bottleneck. Following [START_REF] Rapp | Epigenetics and plant evolution[END_REF], these anthropic actions may have created genomic stress leading to the induction of epigenetic phenomena. From an epigenetic perspective, a population bottleneck leads to a decrease in genetic diversity but an increase in phenotypic variability. In fact, we observed a great density variability between pieces of land, which could be due to phenotypic variability (see, for example, the probability of transition from stage 2 to 3 in the supplementary material).

Based on our study of very local and limited points of observation in Pará, it is impossible to comment further without additional research.

Another explanation could be the past history of each plant. New sequences Uchmaski, J., Grimm, V., 1996. Individual-based


Yassine Bentahar
email: bentaharyassine@gmail.com

Khalid Draoui

Charlotte Hurel

Omar Ajouyed

Slimane Khaïroun

Nicolas Marmier

Slimane Khairoun

Physico-chemical characterization and valorization of swelling and non-swelling Moroccan clays in basic dye removal from aqueous solutions

Keywords: Effect of temperature, Intra-particle diffusion, Kinetics of adsorption, Moroccan clays, Methylene blue

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Industries such as textile, paper, plastic, and leather are the main sources of dyes discharge in aquatic systems [START_REF] Doğan | Adsorption kinetics and mechanism of cationic methyl violet and methylene blue dyes onto sepiolite[END_REF]. Due to their important molecular weight and their complex structure, dyes resist to biodegradation. Hence, they accumulate in water stream and present a potential risk for aquatic organisms and humans. Dyes toxicity depends on their chemical form. For example, azo dyes are the most toxic since reduction products are aromatic amines which are carcinogenic compounds [START_REF] Chen | Toxicity assessment of aromatic amines to Pseudomonas luteola: Chemostat pulse technique and dose-response analysis[END_REF]. Indigo dyes also exhibit significant toxicity. Human exposure to this family of dyes causes skin irritation, gastric disorder, eczema and breathing difficulties [START_REF] Srivastava | Toxicological effects of malachite green[END_REF]. Basic dyes like methylene blue (MB) present also an important menace to human health [START_REF] Rafatullah | Adsorption of methylene blue on low-cost adsorbents: a review[END_REF]. In Morocco, textile industry involves 160 000 employees and represents the first industrial employer in the country [START_REF] Gebrati | Detoxification of textile industry effluents by photocatalytic treatment[END_REF]. Among the textile companies, 29 % are specialized in dying activities. The Moroccan scientific community reported numerous studies on the treatment of effluents from dying industry by using several technologies such as coagulation-flocculation, nanofiltration, and photocatalysis. [START_REF] Abouzaid | Traitement d'effluents de l'industrie textile par combinaison de procédés de coagulation-floculation-nanofiltration[END_REF] studied the combination of coagulation-flocculation-nanofiltration and found that this process exhibited a good ability to remove organic pollutants. Total discoloration of textile effluents was observed with rejection of 98 % of soluble organic matter and 99 % of salts.

Photocatalysis was applied by [START_REF] Gebrati | Detoxification of textile industry effluents by photocatalytic treatment[END_REF]. The results showed a removal rate of 53 % for chemical oxygen demand COD and 89 % for discoloration of industry effluents.

Adsorption process is the most efficient technology for dyes removal. Adsorbents can be of biomass origin such as grass waste, walnut shells, rice husk, algal waste or industrial origin such as red mud and fly ash ( Miyaha et al 2018;[START_REF] Hameed | Grass waste: A novel sorbent for the removal of basic dye from aqueous solution[END_REF][START_REF] Chowdhury | Rice husk ash as a low cost adsorbent for the removal of methylene blue and congo red in aqueous phases[END_REF][START_REF] Gupta | Removal of Rhodamine B, Fast Green, and Methylene Blue from Wastewater Using Red Mud, an Aluminum Industry Waste[END_REF][START_REF] Kumar | Modeling the mechanism involved during the sorption of methylene blue onto fly ash[END_REF][START_REF] Vilar | Methylene blue adsorption by algal biomass based materials: Biosorbents characterization and process behaviour[END_REF]. All of them exhibited a high capacity to remove basic dyes but adsorbents from biomass need preliminary treatments to enhance their adsorption efficiency, and adsorbents from industrial origin can generate a secondary source of pollution. Considering this, clays are considered as good candidates to remove cationic dyes since they occur naturally and abundantly, no supplementary treatment is needed, and they do not provide any secondary pollution. Clays present interesting physico-chemical properties such as low permeability, high porosity, high cation exchange capacity and large specific surface area, allowing efficient dyes adsorption [START_REF] Bentahar | Adsorptive properties of Moroccan clays for the removal of arsenic (V) from aqueous solution[END_REF][START_REF] Christidis | Acid activation and bleaching capacity of bentonites from the islands of Milos and Chios, Aegean, Greece[END_REF]. Dyes adsorption on clays is generally attributed to an ion exchange mechanism and linked to the surface specific area of the adsorbent (Li et al., 2011, Bagane and[START_REF] Bagane | Removal of a dye from textile effluents by adsorption[END_REF].

The permanent negative charge related to isomorphic substitution in the interfoliar space of swelling clays improved the adsorption capacity of cationic elements [START_REF] Azejjel | Natural and modified clays from Morocco as sorbents of ionizable herbicides in aqueous medium[END_REF] or dyes [START_REF] Hong | Adsorption thermodynamics of Methylene Blue onto bentonite[END_REF]. Regarding non-swelling clays, an increased adsorption of MB was reported when clay was treated with NaOH [START_REF] Ghosh | Adsorption of methylene blue on kaolinite[END_REF]. Whatever the kind of clay, the electrostatic interaction between adsorbate and adsorbent is an important parameter conditioning adsorption efficiency [START_REF] Al-Ghouti | The removal of dyes from textile wastewater: a study of the physical characteristics and adsorption mechanisms of diatomaceous earth[END_REF]. A nowadays challenge is to promote the use of naturally occurring solids, without any preliminary treatment for environmental applications. This is both for an economical point of view, since naturally occurring materials are less expensive than manufactured ones, and for an ecological point of view, since these materials will prevent secondary pollution effects. Recent papers investigated adsorption of dyes on natural minerals (illite or kaolinite) [START_REF] Mouni | Removal of methylene blue from aqueous solutions by adsorption on kaolin: kinetic and equilibrium studies[END_REF][START_REF] Omer | Adsorption thermodynamics of cationic dyes (methylene blue and crystal violet) to a natural clay mineral from aqueous solution between 293[END_REF][START_REF] Özmetina | Removal of Methylene Blue dye from aqueous solutions by illite clay[END_REF]. One difficulty with data provided in these papers is that studies were carried out on only one kind of solid, illite on one hand [START_REF] Mouni | Removal of methylene blue from aqueous solutions by adsorption on kaolin: kinetic and equilibrium studies[END_REF] and kaolinite on the other hand (Özmetina et al., 2018), some other authors investigated the effect of temperature on natural clay mainly composed of a mixture of illite and kaolinite [START_REF] Omer | Adsorption thermodynamics of cationic dyes (methylene blue and crystal violet) to a natural clay mineral from aqueous solution between 293[END_REF].

The objective of this work was to provide a complete adsorption study (kinetic effects, equilibrium effects, temperature effects) for two naturally occurring Moroccan clays (Rhassoul-clay and Red clay). Rhassoul-clay and Red clay have different physicochemical properties, and could be promising candidates to remove MB from water. Rhassoul-clay is mainly composed of stevensite, a swelling clay that belongs to the smectite family characterized by a large surface specific area and high cationic exchange capacity. Its texture presents both micro and mesopores. Furthermore, an important negatively charged reaction site was found on its surface related to the isomorphic substitution [START_REF] Bouna | Correlation between eletrokinetic mobility and ionic dyes adsorption of Moroccan stevensite[END_REF].

Therefore, a strong capacity to attract cationic dyes could be expected, as previously shown for methyl violet [START_REF] Elass | Removal of methyl violet from aqueous solution using a stevensite-rich clay from Morocco[END_REF]. For this kind of clay, the electrostatic attraction was the dominant mechanism controlling the adsorption processes. Other swelling clays like bentonite and montmorillonite also exhibited a good affinity toward cationic dyes [START_REF] Hong | Adsorption thermodynamics of Methylene Blue onto bentonite[END_REF][START_REF] Li | Mechanism of methylene blue removal from water by swelling clays[END_REF]. Red clay consists mainly of kaolinite, a non-swelling clay with a low CEC.

Previous study reported that kaolinite with low CEC can be valorized as adsorbent for MB in relatively low concentration [START_REF] Ghosh | Adsorption of methylene blue on kaolinite[END_REF]. Red clay is very abundant in North of Morocco and essentially used in ceramic industry. Its valorisation in water treatment should give an added value to this low-cost mineral resource and give more information about the influence of clay composition for dyes adsorption.

Studying adsorptive behavior of a swelling clay (Rhassoul-clay) and a non-swelling clay (Red clay) towards MB in identical experimental conditions and under various physico-chemical parameters (pH, contact time, temperature, MB initial concentration) could allow us to state about the importance of edge sites and exchange sites reactivity, and to evaluate which parameter is important regarding each kind of surface sites. Finally, considering the whole results one could be able to select which one of natural, low-cost and easy to provide clays could be used in Morocco for dyes removal from wastewater.

Materials and methods

Materials

The adsorbents are clayey materials, extracted from two different sites of Morocco. Red clay was collected on the Mediterranean coast of Morocco in Oued Laou, at 44 km from Tétouan while Rhassoul-clay was collected from Jbel Rhassoul located at the East side of the Middle Atlas Mountains. The samples were used in their natural state without chemical modification.

They were crushed and sieved to obtain particle below 0.1 mm. The mineralogical content of solids were determined by diffractometer model Phillips Xpert-pro using CuKᾳ radiation at 45 kV and a scanning speed of 0.039° s -1 . Chemical characterization of the samples was carried out using an X-ray fluorescence spectrometer (WDXRF, S4 Pioneer from BRUKER).

BET surface area was measured by a Multi-point Beckman Coulter surface analyzer SA 3100 with N2 gas as adsorbate. The morphology of the clays was investigated using Scanning Electron Microscopy (SEM) (Quanta FEG 250). The cation-exchange capacity (CEC) was determined according to hexamino-cobalt method [START_REF] Orsini | Utilisation du chlorure de cobaltihexammine pour la determination simultanee de la capacite d'echange et des bases échangeables des sols[END_REF]. Zeta potentials of solid particles were measured in clay/water system at (m/v = 1 g L -1 ) in pH of suspension and in function of different pH values using a Zetasizer Nano Z (Malvern Instruments).

Reagents

The MB was used in this study as adsorbate and was purchased from Fisher scientific. Its chemical formula is C16H18ClN3S. pH was adjusted using micro-volumes of nitric acid (Fluka, 69 % first grade) and sodium hydroxide (Fluka, purity 98%) solutions at 0.01 M.

Hexaamminecobalt was used to determine the cation exchange capacity. It was prepared at 5 × 10 -3 M from hexaamminecobalt chloride (Fluka, purity ≥ 99%). All reagents used in this study were prepared using a high purity deionized water .

Adsorption studies

The adsorption experiments were realized using batch method at room temperature.

Polypropylene tubes of 50 ml were used as reactors. The solid/liquid ratio used for all adsorption studies was 2.5 g L -1 and 1.25 g L -1 for Red clay and Rhassoul-clay, respectively. These ratios were determined according to recommendations made in [START_REF] Roy | Batch-type procedures for estimating soil adsorption of chemicals[END_REF] EPA procedure. When necessary, the pH of the suspensions was adjusted all along the duration of experiment by adding few drops of dilute HNO3 or NaOH in order to get a constant pH around 8 for Rhassoul-clay and 5 for Red-clay, corresponding respectively to the natural equilibrium pH of each clay. At the end of experiments, the equilibrium pH was measured and found to be constant for each clay. A tube rotator was used to shake the suspensions at 50 rpm. The centrifugation was made at the end of experiment using a Sigma 3-30 K centrifuge.

The supernatants obtained were filtered using 0.45 µm cellulose acetate filters. The concentrations of methylene blue were measured by a UV-vis spectrophotometer (Metertech) at λmax = 660 nm. The adsorbed quantity of methylene blue was calculated as follows:

= ( ) (1)
Where Ci (mg L -1 ) is the initial MB concentration, Cf (mg L -1 ) is the concentration of MB in final solution, qe (mg g -1 ) is the amount of MB adsorbed at equilibrium, m (g) is the mass of adsorbent and V (L) is the volume of solution.

Kinetic experiments

The initial concentration of MB was fixed at 40 mg L -1 and 130 mg L -1 for Red-clay and Rhassoul-clay, respectively. These conditions were imposed by the differences in the adsorption capacity. The experiments were carried out without adjusting the initial pH of the solutions. Experiments with increasing contact times, ranging from 10 min to 4 hours, were carried out.

Effect of pH

The initial MB concentration used is 40 mg L -1 and 130 mg L -1 for Red-clay and Rhassoulclay, respectively. The pH values varying from 3 to 12. The contact time was fixed at 60 min and corresponded to the time necessary to reach adsorption equilibrium for both adsorbents.

Effect of temperature

The effect of temperature was investigated at four temperatures 20, 25, 30, and 35 °C. The initial MB concentration used is 40 mg L -1 and 130 mg L -1 for Red-clay and Rhassoul-clay, respectively. The experiments were carried out on a shaking thermostated plate at 50 rpm.

The experiments were carried out without adjusting the initial pH of the solutions.

Adsorption isotherms

The adsorption isotherms were realized as a function of MB concentration. The experiments were carried at the natural equilibrium pH of each clay. The initial MB concentration varied from 20 mg L -1 to 100 mg L -1 for Red-clay and 100 mg L -1 to 200 mg L -1 for Rhassoul-clay.

Kinetics and isotherm model

Kinetic models

Adsorption kinetic is primordial to evaluate the equilibrium time necessary for the retention of pollutants from aqueous solution and to select the appropriate adsorption conditions.

In this study, the pseudo-first-order, pseudo-second-order and intraparticle diffusion models were used to investigate the mechanism of adsorption and rate controlling steps. The aim of these models is to formulate the general expression to describe the adsorption kinetics on solid surfaces for the solid-liquid adsorption system. These models were employed in many studies to describe the adsorption of dyes on solid sorbents [START_REF] Albadarin | Mechanisms of Alizarin Red S and Methylene blue biosorption onto olive stone by-product: isotherm study in single and binary systems[END_REF][START_REF] Borah | Adsorption of methylene blue and eosin yellow using porous carbon prepared from tea waste: Adsorption equilibrium, kinetics and thermodynamics study[END_REF][START_REF] Doğan | Adsorption kinetics and mechanism of cationic methyl violet and methylene blue dyes onto sepiolite[END_REF][START_REF] Eren | Adsorption performance and mechanism in binding of azo dye by raw bentonite[END_REF][START_REF] Sun | Performance, kinetics, and equilibrium of methylene blue adsorption on biochar derived from eucalyptus saw dust modified with citric, tartaric, and acetic acids[END_REF]Yan et al., 2015). [START_REF] Ho | Pseudo-second order model for sorption processes[END_REF] indicated that application of a sample kinetic model such as first or secondorder models to describe the adsorption on solid surface is often incorrect because of the heterogeneity of the solid and the diversity of sorption phenomena. On one hand, they have also reported that the pseudo first order kinetic model was widely applied and they showed that the rate limiting step was mainly diffusion when the process follows the pseudo-first order kinetic model. On the other hand, the pseudo-second order model is based on the assumption that the rate-limiting step may be chemical adsorption involving valence forces through electron sharing or exchanging between sorbent and sorbate.

Pseudo-first-order model

The pseudo-first-order kinetic model was proposed by [START_REF] Lagergren | Zur theorie der sogenannten adsorption gelöster stoffe[END_REF] and is expressed by [START_REF] Yalçin | The adsorption kinetics of cethyltrimethylammonium bromide (CTAB) onto powdered active carbon[END_REF] as follows:

log(q -q ) = logq -. t (2)
Where K1 is the rate constant (min -1 ), qe and qt represent the MB adsorbed per unit weight (mg g -1 ) at equilibrium and at any time t, respectively.

When the pseudo first order kinetic is applicable, the plot of log (qe-qt) against t should give a linear relationship, and consequently K1 and qe are obtained from the slope and intercept respectively.

Pseudo-second-order model

The pseudo-second-order kinetic is expressed as follows [START_REF] Tseng | Liquid-phase adsorption of dyes and phenols using pinewood-based activated carbons[END_REF]):

= ! + ! # (3)
Where k2 is the rate constant of pseudo-second order (g mg -1 min -1 ). If the pseudo secondorder kinetic is applicable, plot of (t/qt) against t should give a linear relationship. From this representation, the constants K2 and qe can be determined.

Intraparticle diffusion model

The intraparticle diffusion model assumes that intraparticle diffusion is the only ratecontrolling step. The mathematical expression for this model is [START_REF] Kismir | Adsorption characteristics of the hazardous dye Brilliant Green on Saklıkent mud[END_REF] :

= $ # .% + c (4)
Where k3 is the intraparticle diffusion rate constant (mg g -1 min -0.5 ) and c is a constant. The values of k3 and c can be evaluated from slope and intercept of the plot of qt versus t 0.5 . The value of c gives an indication of the thickness of the boundary layer. High intercept value means important effect of boundary layer [START_REF] Banerjee | Adsorption characteristics for the removal of a toxic dye, tartrazine from aqueous solutions by a low cost agricultural by-product[END_REF]. In some cases the plot qt versus square root time can show multi-linearity which indicates that several steps occur in the process.

Isotherms models

Several empirical models were employed to characterize interaction adsorbent/adsorbate.

Langmuir and Freundlich models are often useful for describing adsorption process and were consequently applied in this study.

The Langmuir model supposes that the adsorbent surface is ideally uniform and there are no interactions between adsorbate molecules on adjacent sites. While Freundlich model consider that the adsorbent surface is heterogeneous.

The expression of Langmuir isotherm model can be expressed as:

! = & ' ( + ! ( ( 5
)
Where qm is the maximum monolayer adsorption, KL is the equilibrium constant of Langmuir related to the energy of adsorption and Ce is the concentration of solute at equilibrium.

The expression of Freundlich isotherm model can be represented as:

q= ) * + , (6)
Where Kf is Freundlich constant gives indication of the adsorption affinity of adsorbents and 1/n is an another constant that represents the intensity of adsorption.

Results and discussion

Characteristics of Red-clay and Rhassoul-clay

The mineralogical composition of Red-clay and Rhassoul-clay was characterized in a previous work [START_REF] Bentahar | Adsorptive properties of Moroccan clays for the removal of arsenic (V) from aqueous solution[END_REF]. As shown in Figure 1, kaolinite, illite, and quartz are the major phases in Red-clay while in Rhassoul-clay a presence of stevensite, quartz, calcite and dolomite was observed. The SEM micrograph of the fine fraction particles of Rhassoul and red clay reveal roughness, irregular layer structure and large aggregates consisting of micrometre-sized platelets. The particles seem to be stuck together. The microstructure is characterized by a randomly oriented platelets and porous surface morphology (Fig. 2). The data given in table 1 shows the proportion of different elements existing in the Red-clay and Rhassoul-clay. These results were obtained by X-ray fluorescence analysis. Red-clay is mainly composed of Al2O3, SiO2, Fe2O3, and CaO, while Rhassoul-clay is mainly composed of SiO2, and MgO. The specific surface area (SSA) and cation exchange capacity (CEC) are

given in Table 2. The swelling propriety of Rhassoul-clay was clearly demonstrated by the high values of its CEC (83,3 meq.100 g -1 ) and its SSA (119 m 2 g -1 ). In aqueous solution the SSA of swelling clay can reach 800 m 2 g -1 (Yukselen and Kaya, 2008). For red clay the measured SSA and CEC are 22.4 m 2 g -1 and 35 meq.100g -1 , respectively. Zeta potential was measured at equilibrium pH for the two solids, as illustrated in the Table 2. A net negative charge was detected on the surface of both adsorbents. This parameter can promote the adsorption of MB by an electrostatic attraction between negatively charged surface and positively charged dye molecule.

Adsorption kinetics

Figure 3 shows the effect of contact time on the adsorption of MB onto Red-clay and Rhassoul-clay. No initial adjustment of the pH-value was performed. Naturally, Rhassoul-clay tends to equilibrate at pH 8 and Red-clay tends to equilibrate at pH 5. The adsorption of the two samples was initially rapid, and then became nearly constant after 60 min. These results are correlated with previous studies which have shown that adsorption equilibrium was reached rapidly, sometimes before 30 min of contact [START_REF] Almeida | Removal of methylene blue from colored effluents by adsorption on montmorillonite clay[END_REF][START_REF] Hajjaji | Removal of methylene blue from aqueous solution by fibrous clay minerals[END_REF]. This fast reaction is interesting when considering natural clays as potential dyes adsorbant. Compared to Red-clay, Rhassoul-clay exhibited a higher uptake of MB. To study the adsorption kinetic of MB onto both clays, three models were used: the pseudo first-order kinetic, the pseudo-second-order kinetic (Fig. 4), and the intraparticle diffusion (Fig. 5). The kinetics parameters and the correlation coefficients of kinetics models were presented in Table 3. For pseudo-first-order kinetic the correlations coefficients obtained were 0.97 and 0.71 for Red-clay and Rhassoul-clay, respectively. While the correlations coefficients of pseudo second-order kinetics were 0.99 for Red-clay and 1 for Rhassoul-clay. Furthermore the qm values calculated by pseudo-second-order (Table 3) are very consistent with experimental qm values which confirm that the pseudo-second-order is more suitable to model the adsorption kinetics of Red-clay and Rhassoul-clay than the pseudo-first-order model.

These results are coherent with previous studies on methylene blue adsorption kinetics on sepiolite and pyrophyllite [START_REF] Küncek | Adsorption of methylene blue onto sonicated sepiolite from aqueous solutions[END_REF]Zhang et al., 2015). Application of the pseudo second-order to ours two solids implies that the rate -limiting step may be a chemical adsorption involving valency forces through sharing or exchanging of electrons between adsorbent and adsorbate [START_REF] Ho | Pseudo-second order model for sorption processes[END_REF]. Intraparticle diffusion kinetic model is illustrated in figure 5. Presence of two straight lines was observed for Rhassoul-clay and Red clay, indicating that adsorption process of MB took place in two stages. Similar finding was observed by [START_REF] Gil | Removal of dyes from wastewaters by adsorption on pillared clays[END_REF] on pillared clays. The authors observed a decrease of the diffusion rate when the contact time was increased, and they concluded that the free path of the molecules in the pore was reduced, blocking some molecules during the diffusion process.

The first linear section was attributed to the surface adsorption into the micropores and the second linear section corresponded to adsorption on the external surface. The rate constant corresponding to the first linear section was higher for Rhassoul-clay than for Red clay, and the rate constant for the second linear section was lower for Rhassoul-clay (almost zero) than

for Red clay (Table 3)

Effect of pH

The effect of pH on MB adsorption onto Red-clay and Rhassoul-clay was investigated under varying pH value from 3 to 12. As presented in the Fig. 6, adsorption capacity of MB by the two adsorbents increased slightly with increasing pH. For Red-clay, a drop of adsorption is observed between pH 3 and 7, while for Rhassoul-clay a slight increase occurs between pH 3 to 5, then a plateau is reached between pH 5 to 7, and a second increase occurs from pH 7 to 12. It is well known that alkaline pH improves significantly dye adsorption [START_REF] Almeida | Removal of methylene blue from colored effluents by adsorption on montmorillonite clay[END_REF]. This behavior is due to deprotonation of the reactive edge surface sites, provoking an increase in the number of negative charge on the adsorbent surface and consequently improves the adsorption capacity of the cationic dye. In the case of swelling clays, this deprotonation of edge surface sites is combined to the permanent negative charge of the Rhassoul-clay.

To evaluate the impact of pH on surface solid, a measure of surface charge of adsorbents was realized as a function of pH, and without MB. The fig. 7 shows that surface charge of Red clay is highly dependent on the pH, and becomes more negative at high pH. This variation of charge is attributed to the deprotonation of silanols (SiOH) and aluminols (AlOH) groups on edge sites which are highly pH-dependent surface sites. Similar behavior was observed for other aluminosilicates clays [START_REF] Gücek | Adsorption and kinetic studies of cationic and anionic dyes on pyrophyllite from aqueous solutions[END_REF]Zhao et al., 2003). For Rhassoul clay, pH evolution does not influence significantly the zeta potential. This is due to the permanent negative charge arising from the interfoliar space. From these results we can conclude that adsorption of MB occurs at the surface of edge sites of Red-clay, whereas it occurs both at the surface of edge sites and in the interfoliar space of Rhassoul-clay.

Effect of temperature

The determination of thermodynamic parameters is important to understand the effect of temperature on adsorption. ΔRH° and ΔRS° of MB adsorption on Rhassoul-clay and Red-clay were determined consecutively to experimental data collected at different temperatures.

Knowing that Gibbs Energy depends on enthalpy, entropy and temperature:

ΔRG°= ΔRH° ̶ T ΔRS° (7)
And

) -= . ! ! ( 8
)
where Kd is the distribution coefficient. The adsorption enthalpy ΔRH° and the adsorption entropy ΔRS° can be calculated with the Van't Hoff equation:

ln () -) = 0 1 2° 3 ̶ 0 1 4°3 5 ( 9
)
Where R is the universal gas constant (8.314 J mol -1 K -1 ). The plot of ln(Kd) against 1/T represented in Fig. 8, is a straight line, hence ΔRH° and ΔRS° (supposed to be independent of T) can be calculated from slope and intercept respectively.

The effect of temperature of MB adsorption on Red clay and Rhassoul clay was examined at 20, 25, 30 and 35 °C without adjusting pH. Naturally, the pH equilibrated around 8 for Rhassoul-clay and around 5 for Red-clay. Temperature change from 20 to 35°C did not influence the equilibrium pH measured after each experiments. The results (Table 4) show that the temperature increase leads to an increased adsorption capacity of MB for both clays.

From 20°C to 35°C, the adsorption capacity (Q) increased of 11% for Red-clay and 15% for Rhassoul-clay. Considering this evolution, the thermal effect is benefic for both solids, and slightly more favorable for Rhassoul clay. This is explained by the swelling properties of this clay. Indeed, a temperature increase can enlarge the expandable layer, consequently allowing an easier diffusion of MB molecules [START_REF] Weng | Adsorption of a cationic dye (methylene blue) onto spent activated clay[END_REF]. Negative values of ΔRG° were observed for the two clays indicating the spontaneous nature of the adsorption. The values of ΔRG° were found to decrease with increase of temperature indicating that the adsorption processes of MB onto the Rhassoul-clay and Red-clay were more favorable at high temperature. Enthalpies values for Red-clay and Rhassoul-clay were positive and of the same order of magnitude, indicating (i) that adsorption process is endothermic [START_REF] Almeida | Removal of methylene blue from colored effluents by adsorption on montmorillonite clay[END_REF][START_REF] Doğan | Adsorption kinetics and mechanism of cationic methyl violet and methylene blue dyes onto sepiolite[END_REF], and (ii) adsorption is equivalent in terms of energy for both clays.

The values of ΔS° calculated by the intercept of the plot ln kd versus (1/T) were given as 132.35 J mol -1 for Red-clay and 168.01 J mol -1 for Rhassoul-clay. The positive value of ΔRS° reflects the affinity of these two adsorbents for the MB. Entropy values were of the same order of magnitude for both solids.

From these data, we can conclude that adsorption of MB on Red-clay and Rhassoul-clay is a spontaneous reaction, that is favored when temperature is increased. Enthalpy and entropy values were of the same order of magnitude for both solids, indicating that adsorption is thermodynamically equivalent at the surface of both systems. The main difference between the two solids resides in their adsorption capacity (Q).

Adsorption isotherm

Adsorption of MB on Rhassoul-clay and Red-clay was examined as a function of initial concentration of MB at naturally equilibrated pH (around 8 for Rhassoul-clay and 5 for Redclay). Initial MB concentration has a negligible effect on equilibrium pH. The results obtained were modeled using two empirical models: Langmuir and Freundlich models. The adsorptions isotherms are presented in Fig. 9. Rhassoul-clay and Red-clay exhibited a good correlation with Langmuir model, suggesting that adsorption is limited with monolayer coverage and the surface is energetically homogeneous. The constants and the correlation coefficients calculated by Langmuir and Freundlich are given in the table 5. From these results Rhassoulclay is the most effective adsorbent for methylene blue. This clay possesses the highest maximal adsorption capacity (166 mg g -1 ) compared to the adsorption capacity of MB on Red-clay (18.83 mg g -1 ). The ability of Red clay to remove MB is coherent with its CEC relatively low. Indeed, adsorption occurs only on the edge surface sites of the clay. The comparison of adsorption capacity of Rhassoul-clay and Red clay with others adsorbents reported in literature is presented in table 6. Adsorption capacity of Rhassoul clay was high when compared to the majority of adsorbents. While Red clay exhibited a comparable adsorption capacity with some natural adsorbents. [START_REF] Al-Futaisi | Aspects of cationic dye molecule adsorption to palygorskite[END_REF]. The adsorption properties of Rhassoul-clay is clearly linked to its swelling properties: its important cation exchange capacity resulting from the high negative charge at the clay surface due to the isomorphous substitution [START_REF] Chen | Adsorption of methylene blue on montmorillonite[END_REF][START_REF] Hong | Adsorption thermodynamics of Methylene Blue onto bentonite[END_REF]. Adsorption mainly occurs by ion exchange on the interfoliar space of the Rhassoul clay. These finding are in agreement with some previous studies [START_REF] Al-Futaisi | Aspects of cationic dye molecule adsorption to palygorskite[END_REF][START_REF] Li | Mechanism of methylene blue removal from water by swelling clays[END_REF]. However, some researchers highlighted a discrepancy between methylene blue adsorption and cation exchange capacity [START_REF] Hajjaji | Adsorption of methylene blue and zinc ions on raw and acid-activated bentonite from Morocco[END_REF][START_REF] Ruiz-Hitzky | Adsorption of monovalent organic cations on sepiolite:experimental results and model calculations[END_REF]. The nature of exchangeable cation affects drastically the distribution of methylene blue on the clays surface. Cations with an important hydration capacity such as Na + , Li + , Mg 2+ , Ca 2+

improve adsorption of cationic dye by ion exchange [START_REF] Cenens | Visible spectroscopy of methylene blue on hectorite, laponite B, and barasym in aqueous suspension[END_REF]Yukselen and Kaya, 2008). On the other hand, MB intermolecular iteractions can influence adsorption process [START_REF] Klika | Composition, structure, and luminescence of montmorillonites saturated with different aggregates of methylene blue[END_REF]. Indeed, MB has a tendency to form aggregates in the solution at high equilibrium concentrations. At low concentration (7 × 10 -6 M) the monomer species are predominant [START_REF] Hajjaji | Removal of methylene blue from aqueous solution by fibrous clay minerals[END_REF]. At higher concentrations, MB aggregates as dimers and trimers in the contact of clay particle surfaces [START_REF] Li | Mechanism of methylene blue removal from water by swelling clays[END_REF]. Formation of dimers and trimers of MB at the surface of Rhassoul-clay could explain that adsorption do not reach a plateau in figure 9b (as a function of initial concentration), contrary to figure 9a.

Regarding the experimental conditions used for adsorption isotherms on Rhassoul-clay, MB was introduced in concentrations ranging from 3.10 -4 M to 6.10 -4 M, which are significantly higher than the concentration recommended by [START_REF] Hajjaji | Removal of methylene blue from aqueous solution by fibrous clay minerals[END_REF]. Since the nature of clay and nature of exchangeable cation influence relatively the ratio of dimers and trimers species, this should be deeply investigated. The dimension of MB molecule and its orientation are the key parameters to understand the adsorption mechanism. It is assumed that MB molecule presents a rectangular shape with dimensions 17 Å×7.6Å×3.25Å (Hang and Brindly, 1970).

The orientation of this molecule was found to be arranged in tilted or perpendicular fashion [START_REF] Kaneko | Molecular orientation of methylene blue intercalated in layer-charge-controlled montmorillonites[END_REF][START_REF] Bujdàk | Molecular orientation of methylene blue cations adsorbed on clay surfaces[END_REF]. This orientation is appropriate to the basal plan of clay surface, particularly swelling clays.

Conclusion

. The main factor controlling the adsorption processes of methylene blue on Rhassoul clay and Red clay is adsorption capacity of the clay, pH, and initial concentration of methylene blue. The kinetic adsorption of MB on the two adsorbents is important and equivalent.

Adsorption was found to be spontaneaous and thermodynamically equivalent for Rhassoul

  . Considering Gil et al. (2011) interpretations, we can conclude that adsorption on Rhassoul-clay occurred mainly onto the microporous part of the clay, and occurred both on microporous and on the external surface of Red clay. Compared to pillared clays investigated by Gil et al. (2011) Rhassoul-clay exhibited an interesting diffusion rate, without any surface modification.

Fig. 3

 3 Fig. 1 X-ray patterns of Red-clay (a) and Rhassoul-clay (b)

Fig. 4 Fig. 5

 45 Fig. 4 (a) Pseudo-first-order kinetics and (b) Pseudo-second-order kinetics for MB adsorption onto Red clay and Rhassoul-clay.: Red clay (m/v= 2.5 g L -1 , [MB] = 40mg L -1 , pHeq= 5), Rhassoul-clay ( m/v= 1.25 g L -1 , [MB] = 130 mg L -1 pHeq = 8).

Fig. 6

 6 Fig. 6 Effect of pH on the adsorption of MB onto Red-clay (a) and Rhassoul-clay (b). 60 min of stirring. Red clay (m/v= 2.5 g L -1 , [MB] = 40mg L -1 ), Rhassoul-clay (m/v= 1.25 g L -1 , [MB] = 130 mg L -1 ).

Fig. 7

 7 Fig. 7 The variation of zeta potential of Red-clay and Rhassoul-clay in function of pH. (m/v = 1g L -1 , time 60 min, stirring speed 50 rpm).

Fig. 8

 8 Fig. 8 Plot of ln Kd versus T -1 : estimation of thermodynamic parameters for the adsorption of MB onto Red-clay and Rhassoul-clay. Red clay (m/v= 2.5 g L -1 , [MB] = 40 mg L -1 , pHeq≈ 5), Rhassoul-clay ( m/v= 1.25 g L -1 , [MB] = 130 mg L -1 pHeq ≈ 8).

Fig. 9

 9 Fig. 9 Adsorption isotherms of MB onto Red-clay (a) and Rhassoul-clay (b). 60 min of stirring. Red clay (m/v = 2.5 g L -1 ), Rhassoul-clay (m/v = 1.25 g L -1 ).

  clay and Red clay. Nevertheless, this point should be deeply investigated since two different kinds of surface sites are concerned (exchange sites for Rhassoul clay and edge surface sites for Red clay), justifying that we were expecting thermodynamic discrepancies. Adsorption was found to be more favorable in the pH range from 7 to 12 for both clays. Regarding the results obtained in this study, Rhassoul clay is the most efficient adsorbent. Red clay considered also as good adsorbents for MB, This clay exhibited an adsorption capacity better than some natural adsorbents recognized by their high adsorption potential like zeolite(see table6). Thus the both of clays can be used as a low cost adsorbent for the removal of cationic dye from wastewater. Nevertheless, adsorption reversibility should be deeply investigated to confirm the efficient of these clays are, and to confirm which of surface sites (ion exchange or edge surface) are sufficiently strong to prevent release of cationic dye according to physicochemical conditions. Yan, B.,Chen, Z., Cai, L., Chen, Z., Fu, J., Xu, Q., 2015. Fabrication of polyaniline hydrogel:

	Synthesis, characterization and adsorption of methylene blue. Applied Surface Science 356,
	39-47.
	Yukselen, Y., Kaya, A., 2008. Suitability of the methylene blue test for surface area, cation
	exchange capacity and swell potential determination of clayey soils. Engineering Geology
	102, 38-45.
	Zhang, J., Zhou, Y., Jiang, M., Li, J., Sheng, J., 2015. Removal of methylene blue from
	aqueous solution by adsorption on pyrophyllite. Journal of Molecular Liquids 209, 267-271.
	Zhao, S.M., Wang, D.Z., Hu, Y.H., Liu, B.D., Xu, J., 2003. The flotation behaviour of N-(3-
	aminopropyl)-dodecanamide on three aluminosilicates. Minerals Engineering 16, 1391-1395.

Table 4

 4 Thermodynamic parameters for adsorption of methylene blue on clays

	1

Table 5

 5 Freundlich and Langmuir isotherm constants determined for Red-clay and Rhassoul-clay

		Freundlich isotherm		Langmuir isotherm
	Adsorbents	1/n	Kf	R 2	Kl	Qmax (mg g -1 )	R 2
	Red-clay	0.09 13.5	0.92	2.3	18.83	0.99
	Rhassoul-clay	0.15 119.0	0.96	3.3	166.7	0.99

Acknowledgements

Authors would like to thank Mr. Gabriel Monge (Ecole des mines de Paris -Sophia Antipolis) for XRD measurements.