Marie-Pierre Ledru
email: marie-pierre.ledru@ird.fr

Vivian Jeske-Pieruschka

Laurent Bremond

Anne-Lise Develle

Pierre Sabatier

Eduardo Sávio Passos

Rodrigues Martins

Manuel Rodrigues

Diógenes Freitas Filho

Diógenes Passos Fontenele

Fabien Arnaud

Charly Favier

Eduardo Sávio

Passos Rodrigues Martins

Francisco Rony Gomes Barroso

Francisca Soares

When archives are missing, deciphering the effects of public policies and climate variability on the Brazilian semi-arid region using sediment core studies

Keywords: Caatinga, lake reservoir, sediment, drought, land degradation, human impact

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

Introduction

The brazilian semi-arid (BSA) is the world's most densely populated semi-arid region with more than 54 million inhabitants. It also hosts 1 500 plant, 178 mammal, 591 bird, 177 reptile, 79 amphibian, 241 fish and 221 bee species (Ministerio do Meio Ambiente, 2007) that mostly live in the dry deciduous forest biome known as Caatinga. The combination of rainfall variability and intensive land use will make this region one of the world's most vulnerable to climate change in the coming century (IPCC Climate Change Synthesis Report, 2014). The rainfall regime has been highly variable since the 16 th century with successive periods of high precipitation and droughts. For example, 48 drought episodes have been reported of which 19 lasted two years or more, and in 1958, millions of people were obliged to flee starvation [START_REF] Marengo | Regional differences in aridity/drought conditions over Northeast Brazil: present state and future projections[END_REF]. Climate variability, deforestation, the creation of pastureland, irrigation, extensive farming activities, mining since the 17 th century [START_REF] Neto | Os primórdios da organização do espaço territorial e da vila cearensealgumas notas[END_REF] and, more recently, overpopulation, have encouraged desertification (Centro de Gestão e Etudos Estratégicos, 2016). Brazilian policies in the 20 th century atttempted to halt this process through several land use initiatives and recovery of degraded areas (Centro de Gestão e Etudos Estratégicos, 2017) although such projects were often limited to a single region, one such example being « The Caatinga protected areas program » (Portuguese acronym UCCA) developed by the Environmental Municipal Agency for a single municipality. About 447 reservoirs were built with a total capacity of 30 261 hm 3 although today they only function at 18% of their capacity (https://olhonagua.insa.gov.br).

Since the beginning of the 21 st century, in addition to natural variability and anthropogenic land use change, ongoing climate change is again challenging the BSA. Both global and regional climate change scenarios suggest rainfall deficits and increased aridity will prevail in the region by the second half of the 21 st century [START_REF] Marengo | Regional differences in aridity/drought conditions over Northeast Brazil: present state and future projections[END_REF]. Aware of these risks, the Brazilian government set up the "Recovery of degraded areas and reduction of climate vulnerability" (Portuguese acronym URAD) initiative to address the control of the main drivers of land degradation in the Caatinga biome.

Land cover is a common proxy for the long-term impacts of climate change and land use [START_REF] Lebel | Au Sahel pas de retour à la normale après la "grande sécheresse[END_REF]. For an integrated assessment of the driving forces likely to affect the landscape, data are collected in field surveys [START_REF] Ferrenberg | Climate change and physical disturbance cause similar community shifts in biological soil crusts[END_REF], from agricultural yields, and from aerial photographs [START_REF] Barbosa | Influence of rainfall variability on the vegetation dynamics over Northeastern Brazil[END_REF], as the latter provide an accurate view of the effects of forcing on a landscape and can distinguish anthropogenic from natural forcing in one or several specific time windows (IPCC, 2018). However, in the BSA, the scarcity and the irregularity of such data make it impossible to analyze these effects on the landscape, thereby preventing the production of the crucial information required to take sustainable resource management decisions. Moreover, field surveys do not allow continuous long-term observations (at the scale of several decades) that would be necessary to account for the impacts of current and future climate change when evaluating the future ecosystem services provided by this semi-arid biome. When such archives are absent, sediment cores from lakes (e.g. [START_REF] Levine | The eutrophication of Lake Champlain's northeastern arm: Insights from paleolimnological analyses[END_REF] or reservoirs [START_REF] Cardoso-Silva | Temporal and spatial accumulation of heavy metals in the sediments at Paiva Castro Reservoir (São Paulo, Brazil)[END_REF] are commonly used to reconstruct the chronology of pesticides [START_REF] Sabatier | Long-term relationships among pesticide applications, mobility, and soil erosion in a vineyard watershed[END_REF], organic polluants [START_REF] Zhang | Apparent relationships between anthropogenic factors and climate change indicators and POPs deposition in a lacustrine system[END_REF], trophic history [START_REF] Levine | The eutrophication of Lake Champlain's northeastern arm: Insights from paleolimnological analyses[END_REF] to estimate degradation and propose some alternatives. However these studies do not inform about total landscape and environmental histories. Here to evaluate the combined effects of public policies, human activities and natural climate variability on the total environment and landscape we aim testing a new approach based on multiproxy analyses geochemistry and biological proxies crossed with the instrumental data and the launching of the main public policies and, evaluate how well the sediment cores are able to reconstruct an integrated history of the landscape in the BSA (Appendix A).

For our analysis, we proposed an approach based on indirect reconstructions and produced continuous data from archives retrieved from reservoir sediments, using paleo-limnological techniques. Multiproxy analyses (sedimentological, geochemical, pollen, charcoal) associated with a chronology based on short-lived radionuclides were performed to characterize the signature of both climate variability and land use in the landscape continuously over time.

The aim of the paper is twofold. Firstly, to test a new approach for historical reconstruction of a landscape in the absence of archives. Secondly, to provide reliable data to understand the changes in landscape resources in recent decades in order to enable policy makers to incorporate these factors to insure the sustainable development of the region.

Study area 1.1 Acaraú-Mirim and Araras dams

Our study focused on two dam reservoirs, Acaraú-Mirim and Araras, located in the hydrographic basin of Acaraú, which covers 10% of the total surface area of the State of Ceará with 15 dams along the Acaraú River alone (Fig. 1). The Acaraú River is 320 km long, its headwaters are in the Serra da Mata (700 m asl) near the city of Monsenhor Tabosa and the river flows into the Atlantic Ocean. The Acaraú-Mirim reservoir (03°30.268´S; 40°16.766´W) (ACA) is located in the district of Ipaguaçu, 84 m asl, at the foot of the Serra de Meruoca that feeds the dam; the dam spillway feeds the Acaraú River just below the dam. The montane is covered with a dense seasonal deciduous (Caatinga) and semi-deciduous forest. Annual precipitation at Meruoca (990 m asl) is 1 400 mm. The dam has a capacity of 52 000 000 hm 3 and is used to supply the city of Massapê and farms in the vicinity with water for irrigation and fishery activities [START_REF] Gurgel | Fisheries in semi-arid Northeast Brazil with special reference to the role of tilapias[END_REF]. The Araras reservoir (04°12.559´S; 40°27.254´W) (ARA) (183 m asl) is located in the seasonal tropical deciduous forest landscape in the district of Varjota and provides water for ~141,000 people. Mean annual precipitation is 880 mm. The dam has a capacity of 891 000 000 hm 3 and is intended to supply drinkable water to three cities and their surroundings, as well as water for intensive irrigation (Araras Norte Project), local irrigation and fishery. These reservoirs normally capture a large volume of water during the rainy season, which is irregular in intensity and duration because of the complex climatic pattern [START_REF] Marengo | Regional differences in aridity/drought conditions over Northeast Brazil: present state and future projections[END_REF].

Climate

Between 1990 and 2010, mean annual rainfall in the Acaraú basin was 913 mm. For 45% of the period, annual rainfall was above the mean and for 15% of the period, annual rainfall was below the mean. Rainfall is concentrated between March and May. Thirteen drought events have been recorded in the BSA region in the last 60 years : 1958: , 1966: , 1970: , 1976: , 1979[START_REF] Andrade-Lima | The Caatinga dominium[END_REF]: , 1982: -1983: , 1992[START_REF] Davison | Iron and manganese in lakes[END_REF][START_REF] Araújo Filho | Desenvolvimento Sustentado da Caatinga[END_REF]: -1998[START_REF] Heiri | Loss on ignition as a method for estimating organic and carbonate content in sediments: Reproducibility and comparability of results[END_REF]: -2002: , 2005: , 2007: , 2010: , and 2012[START_REF] Boyin | Extended Reconstructed Sea Surface Temperature (ERSST), Version 4[END_REF][START_REF] Alexander | Projected sea surface temperatures over the 21 st century: Changes in the mean, variability and extremes for large marine ecosystem regions of Northern Oceans[END_REF][START_REF] Marengo | Regional differences in aridity/drought conditions over Northeast Brazil: present state and future projections[END_REF]. The February-May rainy season in The climate of the mountain area near ACA is moister than the inland climate at ARA due to the relief and the proximity of the ocean. Consequently, differences in rainfall anomalies are observed between our two sites. ACA with 1 400 mm rainfall per year is less sensitive to drought than the whole hydrographic basin comprising Acaraú and ARA.

At ACA, the driest years were 1983were , 1992were -93, 1997were -98, 2005were and 2012were -2015were . High moisture rates were recorded in 1985were , 1994were and 2009. . At ARA, monitoring started in 1987. The driest years recorded were 1992-93, 1997-98, 2012-2015 and the wettest years were 1989, 1995 and 2009. Considering the Acaraú hydrographic basin as a whole, the driest years were 1979-83, 1992-93, 1998, 2012-2016. Based on these observations, the year 2005 was considered « dry » only at ACA, near the coast. Correlations with the Pacific Decadal Oscillation (PDO 1900(PDO -2013) ) anomaly show that it was negative during the driest years : 1992, 1997-2003, and 2011-2013, while no correlation was observed with the AMO anomaly. A negative anomaly in tropical North Atlantic SST was observed in 1970-75, 1982-84, 1992-93 and in 1997-98 a sharp decrease was observed but the anomaly remained positive. Another correlation between SST and rainfall is observed when considering the SST gradients between tropical North and South Atlantic, the interhemispheric SST dipole, with the driest /wettest years when the S/N Atlantic SST was colder/warmer. Indeed, Northeastern Brazil is one of the key regions that are severely affected by the meridional mode of tropical Atlantic SST variability [START_REF] Alexander | Projected sea surface temperatures over the 21 st century: Changes in the mean, variability and extremes for large marine ecosystem regions of Northern Oceans[END_REF]. During at least the last 50 years, the changes in SST have been closely correlated with rainfall variability.

Vegetation

The Caatinga dominium is one of the three arid and semi-arid cores of the South American continent (Andrade- [START_REF] Andrade-Lima | The Caatinga dominium[END_REF]. The Caatinga dominium is a mosaic of tropical thorn woodland and very dry thorn forest at altitudes below 500-600 m asl and dry forest at higher altitudes (sensu the Holdridge life zone ecology). The phytogeographic classification system of Brazil is called savanna-stepic (IBGE, 2010) and represents the dryland biome in Brazil. It is located in northeastern semi-arid region, which occupies 11% of the country with an area of about 1.3 million km 2 . Delimitation was based on rainfall, mean anual rainfall below 800 mm, aridity index up to 0.5, and risk of drought > 60% relative to the 1970-90 climatology [START_REF] Joly | Evolution of the Brazilian phytogeography classification systems: implications for biodiversity conservation[END_REF]. Despite the semi-arid physiognomy, the Caatinga hosts high biomass and high biodiversity.

Four main types of tropical vegetation [START_REF] Holdridge | Tropical Science Center[END_REF] are observed in the Acaraú River Basin: very dry thorn forest (Caatinga) and three types of seasonal dry forests on the slopes and tops of mountains (deciduous, semi-deciduous with babaçu and evergreen forest). The dominant ecosystem is the very dry thorn forest or Caatinga associated with stony impermeable crystalline soils adapted to hydric deficiency. At Acaraú, the Caatinga has been severely impacted by human activities and hotspots of ongoing desertification.

Intensive cottton crops, introduction of grazing in the 16 th century in the Acaraú Basin and wood cutting have greatly modified the Caatinga. The degradation of the arboreal Caatinga led to expansion of the shrubs and grasses observed today. The limiting factors for agriculture, i.e. the climate, shallow soil, stony surface, hydric deficiency) obliged the inhabitants to practice itinerant agriculture, in which the land is abandoned after 2 or 3 years, leaving behind secondary vegetation of no economic value (Araújo Filho, 1997) (Appendix D).

Ecological indicators

The following species were assigned to the drought tolerant group Mimosa ceasalpiniifolia, M tenuiflora, Piptadenia stipulacea, Anadenanthera colubrina, Poincianella bracteosa, Combretum, Myracrodruon urundeuva, Spondias, Zizyphus joazeiro. Mimosa caesalpiniifolia is the dominant arboreal taxon in the ecosystem despite being very sensitive to degradation in areas susceptible to desertification. These trees occur naturally in the Caatinga and are one of the main sources of stakes for fences in Northeastern Brazil. Considering the climatic conditions of the Brazilian semi-arid region, Mimosa caesalpiniifolia is assumed to be a fastgrowing plant species, and is the wood used as fuelwood and for charcoal production [START_REF] Albuquerque | How ethnobotany can aid biodiversity conservation: reflections on investigations in the semi-arid region of NE Brazil[END_REF][START_REF] Milliken | Impact of management regime and frequency on the survival and productivity of four native tree species used for fuelwood and charcoal in the caatinga of northeast Brazil[END_REF]. Combretum leprosum is a good pollen producer although not abundant in the vegetation. It grows easily in degraded areas.

Alternanthera, Poaceae, Cyperaceae, Mitracarpus and Borreria are common in the herbaceous strata of the Caatinga, including on the margins of the reservoirs.

Material and methods

Coring

Two sediment cores were collected in two reservoirs, Acaraú Mirim (ACA) and Araras (ARA) with a gravity corer UWITEC pushed into the sediment in a stable area of the reservoir where the height of the water column was between 6 and 7 m, to guarantee permanent water. The ACA 15-2 core went down to a depth of 97 cm. The ARA 15 core went down to a depth of 60 cm. Cores were split in half using a cutter and a thin metal wire. One half core was kept intact and sent to the EDYTEM laboratory in Chambéry, France for XRF analyses the other half, used for proxy analyses, was cut into 1-cm thick sections, that were immediately placed in plastic bags for storage and kept refrigerated in the laboratory until processing and analysis. Proxy analyses were performed on alternate 1-cm samples.

The base of core ACA 15-2 is composed of compact organic sandy clay up to 80 cm, a layer of fine dark clay between 80 and 77 cm, and fine brown organic clay to the top of the core.

The ARA 15 core is composed of fine organic clay with loose wet sediment in the top 10 cm.

To define an age depth relation for each core, 13 samples from ACA and 11 samples from ARA were sent to Laboratoire Souterrain de Modane France for short-lived radionuclide measurements.

Dating

210 Pb, 226 Ra, and 137 Cs activities were determined using planar Broad-Energy Germanium (BEGe) detectors placed at the Laboratoire Souterrain de Modane following [START_REF] Reyss | low background well-type detectors for measurements of environmental radioactivity[END_REF] on 1-to 3-cm-thick samples with 12 samples on ARA15 sediment core and 10 samples on ACA15 sediment core. The levels of 226 Ra activity were determined using its short-lived daughters 210 Pb (295-and352-keV peaks) and214 Bi (609-keV peak), assuming secular equilibrium with 226 Ra. 210 Pb (22.3 y), 241 Am (432.2 y), 7 Be (53.4 d), and137 Cs (30.2 y) activity levels were directly measured by their gamma emissions at 46.5, 60, 477, and 662 keV, respectively. Excess 210 Pb activity was calculated by subtracting 226 Ra-supported activity from total 210 Pb activity.

Sedimentological and geochemical analyses

The grain size distributions of core ACA15 were determined using a Malvern Mastersizer S (EDYTEM lab) at a 1 cm continuous interval. After inserting the bulk sediment into the granulometer, ultrasound was applied to minimize particle flocculation. Core ACA15 was also sampled at 1-cm steps and dried at 60 °C during 4 d to obtain its dry bulk density, and then the LOI of each sample was measured using the protocol of [START_REF] Goldberg | Geochronology with 210Pb[END_REF] and [START_REF] Heiri | Loss on ignition as a method for estimating organic and carbonate content in sediments: Reproducibility and comparability of results[END_REF]. The LOI at 550 °C and 950 °C corresponds to the organic and carbonate components of the sediment, respectively.

The relative concentrations of major elements were analyzed on the surface of the sediment core by X-ray fluorescence (XRF) at high resolution (5 mm sampling step) with an Avaatech Core Scanner (EDYTEM Laboratory, CNRS-University Savoie Mont Blanc). The X-ray beam was generated with a rhodium anode and a 125-μm beryllium window, which allows a voltage range of 7-50 kV and a current range of 0-2 mA. Element intensities were expressed in counts per second (cps). Geochemical data were obtained with different settings depending on the elements analyzed. They were adjusted to 10 kV and 1 mA for 20s to detect Si, Ca, Al, Fe, Ti, K, Mn, and S. For heavier elements (i.e. Sr, Rb, Zr, Br, and Pb), measurements were performed at 30 kV and 0.75 mA for 30 s [START_REF] Jansen | CORTEX, a shipboard XRF-scanner for element analyses in split sediment cores[END_REF]. Identification of relationships between element depth-series, the data on each core was subjected to principal component analyses (PCA) using vegan package (Oksanen et al., 2018) on the R platform (R Core Team, 2018).

Pollen analyses

Pollen analyses were performed at 1 cm intervals on 76 samples from ACA and 50 samples fom ARA. Lycopodium spikes were used to estimate pollen concentrations [START_REF] Stockmarr | Tablets with spores used in absolute pollen analysis[END_REF]).

The samples treatment followed standard protocol [START_REF] Faegri | Textbook of Pollen Analysis[END_REF][START_REF] Kümmel | Handbook of paleontological techniques[END_REF]. A minimum of 300 terrestrial grains were counted per sample. Pollen grains and spores were identified using the reference pollen collection of about 130 Caatinga taxa held at the Prisco Bezerra Herbarium of the Federal University of Ceará, and several pollen atlases [START_REF] Miranda | Fundamentos de Palinologia: principais tipos de polen do litoral cearense[END_REF][START_REF] Oliveira | Prospecção palinológica em méis da Bahia (Print Mídia, Feira de Santana[END_REF][START_REF] Radaeski | Pólen nas angiospermas: diversidade e evolução[END_REF][START_REF] Salgado-Labouriau | Contribuição à palinologia dos Cerrados[END_REF][START_REF] Silva | Flora polínica das caatingas: Estação Biológica de Canudos [Canudos, Bahia, Brasil[END_REF]RCPol).

Burning biomass

The two cores (ACA 15 and ARA 15) were sampled continuously in each 1-cm layer for the charcoal analyses and prepared using the following steps: (1) Samples were soaked in 5% KOH solution; (2) particles were deflected and bleached in a 10% NaOCl solution; (3) 160μm mesh was used to separate microscopic coal particles; and (4) microcharcoals were analyzed under a stereomicroscope and using the image analysis software WindSeedle [START_REF] Umbanhowar | Experimental production and analysis of microscopic charcoal from wood, leaves and grasses[END_REF].

Sea Surface Temperature

Sea surface temperature (SST) data were taken from the NOAA's Extended Reconstructed Sea Surface Temperature, version 4 (ERSSTv4) dataset [START_REF] Boyin | Extended Reconstructed Sea Surface Temperature (ERSST), Version 4[END_REF], available on 2º x 2º global grids, for the years 1950-2017. In order to examine the role of SST, we computed an index of area-averaged SST anomalies within the northern Tropical Atlantic (1°-20° N, 48º-20º W) for the December-February season. Anomalies are relative to the 1981-2010 monthly average.

Remote sensing

We selected the LANDSAT images with the lowest cloud cover index between the years 1973 and 2016. The images were analyzed and interpreted using the functions available in a Geographic Information System (GIS) software for vector editing, cartographic projection conversion and digital image processing. All images were redesigned for the UTM projection, SIRGAS 2000 datum, zone 24 of the southern hemisphere.

Results

Age model

A chronological framework was established using short-lived radionuclides measurements.

The 210 Pb excess profiles were plotted on a logarithmic scale (Fig. 2a) and we applied the Constant Flux Constant Sedimentation model (CFCS) [START_REF] Goldberg | Geochronology with 210Pb[END_REF] for our age model.

For ACA 15 we obtained a mean accumulation rate of 13.8 +/-2.25 mm.yr -1 and the base of lake sediment dated to 1963 +/-8.5 yr AD. For ARA15, the 210 Pb excess profile (Fig. 2b) showed a regular decrease interrupted by drops in 210 Pbex activity. Following [START_REF] Arnaud | Flood and earthquake disturbance of 210Pb geochronology (Lake Anterne, NW Alps)[END_REF], the low values of 210 Pbex were disregarded for the construction of the synthetic sedimentary record, because the values are related to sedimentary facies that are considered to be instantaneous deposits (Fig. 2b) illustrated by an increase in Si content, corresponding to facies 2. Plotted on a logarithmic scale, the 210 Pbex activity vs corrected depth revealed a linear trend. Applying the CFCS model, we obtained a mean accumulation rate of 7.17 +/-1.45 mm.yr -1 . Ages were then calculated using the CFCS model applied to the original sediment sequence to provide a continuous age-depth relationship and an age interval of the four instantaneous deposits corresponding to facies 2 (1980±7 yr AD, 1976±8 yr AD, 1972±9 yr AD, 1964.5±10.5 yr AD). The base of the lake sediment (facies 1 and 2) dated to 1959.5 +/-11.5 yr AD in good agreement with the dam construction in 1958 AD. Values for 137 Cs activities were under the limit of detection in both cores.

Depositional environment

Color, grain size, loss on ignition (LOI), and sedimentary structure of the sediments in the reservoir (Fig. 3, 4, 5) allowed to identify three different sedimentary facies in both cores. A brown to green fine-grained sediment (facies 1) in the upper part of both cores (0-74.5 cm in ACA15; and 0-45.5 cm in ARA15), a grey fine-grained sediment (facies 2) interbedded at different depths in facies 1 and a brown coarse sediment (facies 3) at the base of each core (95.5-74.5 cm in ACA15; 60-45.5 cm in ARA15). At ARA, the top 10 cm of the record was a mixture of sediments/water interface and was considered as a single sample.

PCA variables and individual projections highlighted the correlation between different elements in the two cores (Fig. 3). Dimensions 1 and 2 (hereafter Dim1 and Dim2) represent 71% and 78% of the total variance for ARA15 and ACA15, respectively (Fig. 4). In both lakes, the ordination of variables among the first two axes are quite comparable and identify the same main factors of variations, probably linked to the common geological context and origin of the lakes. Three main chemical end-members or units were identified (Fig. 3). The first one was positively correlated with Dim1 and yielded high positive loadings for major terrigenous elements such as K, Zr, Sr and to a lesser extent Ca in both cores, plus Rb in ARA15, and plus Si in ACA15. The second pole with negative loadings on Dim1 made it possible to differentiate Br, S, Mn probably linked to in-lake processes: organic matter (Br) and changes in oxic/anoxic conditions (Mn) in the lake system. The third end-member with a positive loading on Dim2 included other terrigenous elements: Ti and Al in both cores, plus Si and Pb in ARA15 and Rb in ACA15. We were able to separate the two terrigenous endmembers thanks to their grain size, the size corresponding to K, Zr, Sr was related to the sand fraction of the sediment and the one enriched in Ti, Al was related to the fine-grained fraction (< 4 µm) of the sediment. Fe was located between organic and the fine terrigenous endmember, probably linked to its redox sensitive behavior. Adding the lithological facies 1, 2 and 3 allowed us to map the geochemical distribution of the data (Fig. 3,7,9). Mapping the units in the PCA revealed a clear link between facies 1 and the second end-member related to organic rich sediment (LOI550°C around 15%) linked to long-term sedimentation in the lake system (Fig. 5). Facies 2 was positively correlated with Dim2 and hence with the fine-grained terrigenous end-member and probably with the high terrigenous input from the watershed due to erosion during high rainfall events/periods. Facies 3 was positively correlated with Dim1 linked to the coarse sediment at the base of each sediment sequence, probably corresponding to soil sediment that existed before the dam was constructed at each site.

Reconstruction of the tree cover

The composition of the vegetation was reconstructed from pollen analyses (Fig. 7, 8, 9, 10; Apendix). Among the tree taxa of the Caatinga Mimosa caesalpiniifolia, M tenuiflora, Piptadenia stipulaceae, Anadenanthera colubrina, Poincianella bracteosa, Combretum leprosum, Myracrodruon urundeuva, Spondias sp., Zizyphus joazeiro were assigned to the drought tolerant group [START_REF] Lima | Lignin composition is related to xylem embolism resistance and leaf life span in trees in a tropical semiarid climate[END_REF]. M caesalpiniifolia, one of the dominant arboreal taxa in the Caatinga is very sensitive to degradation [START_REF] Milliken | Impact of management regime and frequency on the survival and productivity of four native tree species used for fuelwood and charcoal in the caatinga of northeast Brazil[END_REF]. C leprosum is a good pollen producer that grows easily although not abundantly in degraded areas [START_REF] Milliken | Impact of management regime and frequency on the survival and productivity of four native tree species used for fuelwood and charcoal in the caatinga of northeast Brazil[END_REF]. Alternanthera, Poaceae, Cyperaceae, Mitracarpus and Borreria are common in the herbaceous strata of the Caatinga (Costa et al., 2017), including on the margins of the reservoirs [START_REF] Tabosa | Live fast and die young: The aquatic macrophyte dynamics in a temporary pool in the Brazilian semiarid region[END_REF].

For ACA15, five pollen zones and 81 different pollen and spore types provide evidence for several changes in the landscape during the last 55 years (Fig. 6, 7, Appendices B, C).

The basal pollen zone (ACA-I, 75 -64.5 cm; 11 samples; 1961 to 1969 AD) comprises taxa of the Caatinga with abundant arboreal pollen (27-48%). Fabaceae pollen, mainly Mimosa caesalpiniifolia (21-44%) and Mimosa tenuiflora (0.4-4%), was well represented.

Drought tolerant tree pollen taxa ranged from 1% to 6%. This group was represented by Mimosa tenuiflora, Piptadenia stipulacea, Anadenanthera colubrina, Poincianella bracteosa, Combretum leprosum, Myracrodruon urundeuva, Spondias and Zizyphus joazeiro pollen grains. The herb taxa of the Caatinga were also well represented. A single pollen grain of Zea mays was counted at a depth of 67 cm. Scattered pollen grains belonging to taxa of the seasonal semi-deciduous forest -Serra da Meruoca -were observed in this zone. Aquatic taxa (≤2%) and fern spores (<2%) were scarce. The number of carbonized particles was low and infrequent in this period.

Zone ACA-II (64.5 -48.5 cm; 16 samples; 1969 to 1981 AD) was characterized by higher values of tree pollen taxa (44-62%). Although Mimosa caesalpiniifolia (31-51%) was the most abundant taxon, greater numbers of Mimosa tenuiflora (2-10%), Piptadenia stipulacea (0.3-3%) and Senna (0.3-3%) accounted for the increase in arboreal pollen. Drought tolerant tree pollen taxa reached higher values (5-13%), mostly due to increased percentages of Mimosa tenuiflora (2-10%) and Piptadenia stipulacea (0.3-3%) pollen. While pollen of Anadenanthera colubrina appeared in more samples than previously, percentages of Combretum leprosum pollen were still very low (0-1%). Pollen from Alternanthera (15-25%) and Poaceae (4-15%) were still well represented, whereas pollen from Borreria (2-18%) and Mitracarpus (1-4%) declined and both Asteraceae (0-2%) and Urticaceae (0-2%) pollen were still only present at lower proportions. Abundances of Cyperaceae pollen remained stable (2-5%) but both Cleome (0-4%) and the aquatic plant Echinodorus (0-2%) increased slightly in this zone. The amount of charcoal particles was still low in this zone.

The pollen assemblages of the ACA-III zone (48.5 -37.5 cm;11 samples;[START_REF] Andrade-Lima | The Caatinga dominium[END_REF]to 1988 AD) AD) were characterized by a significant decrease (from 42% to 16%) in arboreal pollen. NAP taxa of the Caatinga were well represented. Pollen of Zea mays (<1%) appeared more frequently.

There was a marked decrease in arboreal pollen, mostly due to the decline in Mimosa caesalpiniifolia (7-34%). Drought tolerant tree pollen taxa decreased (1-8%), mainly due to the low percentages of pollen of Mimosa tenuiflora (0.2-5%) and the low proportions of Piptadenia stipulacea (0-2%) and Myracrodruon urundeuva (0-1%). However, while the proportions of Combretum leprosum pollen remained low at the base of the zone (0.3%) they increased to 3% toward the top. Dam bank taxa were more abundant (4-12%) than in zone ACA-II, mainly due to higher proportions of Cyperaceae (1-8%) but also Cleome (2-4%) pollen. Aquatic taxa (<2%) were slightly less frequent than in the previous zone. Charcoal particles were very scarce in this period. ACA-IV (37.5 -25.5 cm;12 samples;1988to 1998 AD) AD) showed an increasing trend in arboreal pollen abundances (30-61%), whereas the opposite was observed for non-arboreal pollen percentages (35-63%). Arboreal pollen was continuously represented by Mimosa caesalpiniifolia (26-54%) in which there was a marked increase compared to the previous zone. Among non-arboreal pollen, there was a significant decrease in Alternanthera (12-25%), Mitracarpus (0-5%), Asteraceae (<2%), Urticaceae (0-2%) and in Acalypha (0-2%).

Pollen of Poaceae (2-9%), Zea mays (<1%) and Borreria (4-23%) decreased toward the top of the zone. The increase in the pollen curve for drought tolerant tree taxa (2-11%) was mainly due to higher pollen counts of Mimosa tenuiflora (1-9%), but Piptadenia stipulacea (<2%), Combretum leprosum (<2%) and Myracrodruon urundeuva (0-1%) remained low. Pollen of Poincianella bracteosa appeared for the first time at 35 cm. Abundances of plant species growing on the dam banks (2-8%) decreased in comparison to the previous zone. No charcoal particles were counted in this zone.

The uppermost pollen zone (ACA-V, 25.5 -0 cm; 26 samples; 1998 to 2015 AD) was characterized by stable arboreal (37-63%) and non-arboreal (32-56%) pollen taxa. Mimosa caesalpiniifolia (28-52%) dominated the arboreal pollen group, while Alternanthera (11-21%), Poaceae (3-11%), Borreria (10-24%) and Mitracarpus (1-8%) accounted for the highest values in the non-arboreal pollen group. Pollen of Zea mays (<2%) was still present.

Proportions of pollen of drought tolerant tree pollen taxa increased from 3% to 14%.

Abundances of Mimosa tenuiflora (1-11%) and Myracrodruon urundeuva (0-2%) pollen remained stable, while percentages of Combretum leprosum pollen continued to increase (0.3-4%). Pollen of Piptadenia stipulacea (0-2%), Anadenanthera colubrina (<1%) and Senna (0-3%) pollen increased slightly. Pollen of Cyperaceae (2-10%) and Cleome (1-6%) reached their highest values in the whole record. The percentage of moist semi-deciduous forest increased to 2% due to higher proportions of pollen of Cecropia (<1%), Myrtaceae (<1%), Anacardium occidentale (<1%) and Attalea speciosa (0-1%). Pollen of Sapindaceae (not shown in the diagram) and Copaifera appeared for the first time in this zone. The number of carbonized particles increased markedly in this period, but gradually decreased again toward the top.

For ARA15, constrained cluster analysis divided the pollen diagrams into four main zones (ARA-I to IV) and 74 pollen taxa and four spore types were identified (Fig. 8,9; Appendices B, C).

The top 10 cm of the record was a mixture of sediment/water interface and was considered as a single sample. Based on the constrained cluster analysis as well as the 74 pollens and four different spore types identified, the pollen diagrams were divided into four main zones (ARA-I to IV).

Zone ARA-I (60 -45.5 cm;15 samples;before 1960before to 1966 AD) AD) was characterized by high percentages of NAP (62-84%) and low percentages of AP (9-31%). Pollen from drought tolerant tree taxa ranged from 3% to 13%. This group was represented by Mimosa tenuiflora, Piptadenia stipulacea, Anadenanthera colubrina, Poincianella bracteosa, Combretum, Myracrodruon urundeuva, Spondias and Zizyphus joazeiro pollen grains. Low abundances of AP were mostly due to lower percentages of Mimosa caesalpiniifolia pollen (1-16%) and Mimosa tenuiflora pollen (1-8%). Pollen of Croton reached higher values, between 1% and 7%. The NAP were dominated by pollen of Alternanthera (33-54%), Poaceae (5-12%), Borreria (1-4%), Mitracarpus (4-10%), Asteraceae (4-16%) and other herbaceous plants (Froelichia/Gomphrena, Euphorbia and Microtea) not shown on the diagram. Pollen grains of the moist seasonal semi-deciduous forest occurred sporadically. This zone contained high pollen abundances of sedges (Cyperaceae 7-20%) as well as herbs such as Cleome (0-2%) and Cuphea (0-2%) that grow on the banks of the dam. The aquatic group was well represented by pollen of Echinodorus (6-19%) and spores of Salvinia (0-2%). Carbonized particles were only present in the uppermost part of the sample. The seasonal semi-deciduous forest group increased up to 14% in this zone, because of higher proportions of Myrtaceae (0-11%) pollen.

Zone ARA-II (45.5 -39.5 cm; five samples; 1966 to 1972 AD) was marked by an increase in AP from 40% to 53%, mostly due to an increase in Mimosa caesalpiniifolia with values up to 33%. The higher representation of pollen from drought tolerant tree pollen taxa, which varied between 9% and 25%, was due to the high pollen frequencies of Mimosa tenuiflora (5-17%), Piptadenia stipulacea (1-5%) and Myracrodruon urundeuva (0.3-3%) while Combretum pollen frequencies remained stable (0.3-2%). Abundances of Croton (0-2%) pollen decreased in this zone. Herbs, mostly represented by pollen of Alternanthera (5-21%), Poaceae (3-15%), Borreria (1-2%), Mitracarpus (2-12%), Asteraceae (2-10%) and Urticaceae (1-2%) were less frequent. The moist semi-deciduous forest taxa increased to 14% in this zone, because of higher proportions of Myrtaceae (0-11%) pollen. The abundance of herbs growing on the dam banks decreased (4-10%), due to less frequent Cyperaceae pollen (2-8%). The proportion of aquatics (0.6-6%) was lower than in the previous zone due to the decline of Echinodorus pollen (0.6-4%) and the absence of Salvinia spores. Charcoal particles were observed at depths of 42 and 43 cm.

Zone ARA-III (39.5 -18.5 cm;21 samples;1972to 1989 AD) AD) -AP increased continuously from 30% to 82% but then dropped to 35% at the top of the zone. The most abundant arboreal pollens were Mimosa caesalpiniifolia (10-71%), Mimosa tenuiflora (5-40%), Piptadenia stipulacea (2-7%), Anadenanthera colubrina (0-2%), Combretum (0-3%) and Myracrodruon urundeuva (0-9%). NAP had a minimum of 13% in the lower part of the zone but increased to 64% toward the top. Among herbaceous taxa, Alternanthera (5-32%) and Poaceae (2-22%) pollen were abundant. Pollen of Borreria (0-3%), Mitracarpus (0.3-5%), Asteraceae (0.3-5%), Urticaceae (0-3%), Acalypha (0-4%) and Zea mays (<1%) were also well represented.

Pollen from drought tolerant tree taxa showed an increasing trend from 11% to 46%. In this zone, taxa from the moist semi-deciduous forest group were represented by lower pollen frequencies, between 0% and 4%, mainly because of very low values of Myrtaceae (<1%).

Percentage values of Attalea speciosa pollen (<1%) and Copaifera (<1%) remained low, whereas Arecaceae (0-3%) pollen grains became more frequent. Taxa that grow on the dam banks also reached higher values, mostly due to increasing percentages of Cyperaceae (2-10%) and Cleome pollen (0-4%). Aquatics taxa were represented by increasing Echinodorus (0.3-5%) pollen and Salvinia (<1%) spores. The amount of microscopic charred particles increased during this period.

Zone ARA-IV (18.5 -0 cm; 9 samples; 1989 to 2002 AD) -AP predominated with frequencies ranging from 53% to 71%, whereas NAP decreased, ranging from 26% to 38%.

Percentages of Mimosa caesalpiniifolia pollen remained high (17-46%). The increase in pollen from drought tolerant tree taxa from 20% to 38% was mainly due to higher values of Mimosa tenuiflora (10-23%) and Piptadenia stipulacea pollen (1-13%). Pollen abundances of Anadenanthera colubrina (0-3%) and Myracrodruon urundeuva (1-4%) increased slightly, while percentages of Combretum pollen (up to 2%) remained low. Poincianella bracteosa pollen appeared for the first time at a depth of 15 cm. Herbs were characterized by decreasing values of Alternanthera (1-8%) and Poaceae (5-12%), and percentages of pollen from other herbaceous taxa continued to be low. Pollen of Borreria (1-4%) and Mitracarpus (2-6%) increased in the upper part of this zone. Proportions of seasonal semi-deciduous forest taxa decreased from 2% to 0%. Pollen frequencies of Cyperaceae increased to 17% while percentages of Cleome (<2%) decreased in this zone. Aquatic taxa remained stable at (2-6%).

No charcoal particles were counted in this zone.

Discussion

Effects of climate variability

Interannual rainfall variability does not reflect unequivocally in any of both records. Changes in the geochemical proxies Si and Zr identified in facies 2, revealed four terrigenous events in both records although the timing was different, showing an increase in erosion probably linked to either extreme precipitation events or to agricultural activities [START_REF] Hounsou-Gbo | Tropical Atlantic contributions to strong rainfall variability along the Northeast Brazilian coast[END_REF][START_REF] Wilhelm | 1.4 kyrs of flash flood events in the Southern European Alps : implications for extreme precipitation patterns and forcing over the north western Mediterranean area[END_REF]. At ACA15, these events match periods of increased precipitation in the years 1982-83, 1987-88, 1992-96 and 2006-2012, with associated age model uncertainties, that can be clearly observed in the sediment sequence but less in the vegetation composition (Fig. 4). Between AD 1969 and 1981 (pollen zone ACA-II; Fig. 6) there was a significant increase in the tree taxon M tenuiflora although this zone includes two climatic intervals: the interval 1969 to 1974 with above normal rainfall, with six months of more than 50 mm/month, and the interval 1975 to 1980 with below average rainfall (Appendix B).

Between AD 1997 and 2015, the arboreal (37-63%) taxa remained stable even though this zone includes the driest phase of the entire study with only two years with consecutive aboveaverage rainfall (1999 -2000 and 2008 -2009), of which only 2009 had significantly more rain. Over the last 15 years, above average rainfall was only recorded in three years (2008, 2009 and 2011). In the years 2012, 2014 and 2015, the recorded rainfall values resembled those during the severe droughts in 1966 and 1983. But despite several consecutive years of below normal rainfall, rains were well distributed over five to six months with ≥ 50 mm/month (Appendix B) unlike in 2012 and 2015 when the little rain that fell was distributed over three months. In 2014, even with four rainy months, the standard deviation of the annual mean precipitation was the lowest since 1963.

In ARA15, the four observed terrigenous events in 1964.5±10.5 yr AD, 1972±9 yr AD, 1976±8 yr AD, 1980±7 yr AD (Fig. 4) do not match the major reservoir overflows that were observed twelve times between 1978 and2011 (Fig. 10, 11). Between 1960 and 1980 the tree cover progressively increased either due to the absence of major drought episodes during this interval or to an increase in local moisture due to the construction of the dam (Fig. 3, 10, 11).

We observed changes in the vegetation with a sharp decrease in M caesalpiniifolia during the 1982 drought followed by re-expansion during the rainy interval in 1984-85 (Fig. 11).

Between AD 1990 and 2002, the Caatinga tree taxa increased and showed the highest values of the whole record. No major changes in vegetation were observed during the drought in 1992-1993, whereas a sharp decline in Piptadenia stipulacea was observed during the drought in 1998, a decrease in M caesalpiniifolia was observed during the drought in [START_REF] Heiri | Loss on ignition as a method for estimating organic and carbonate content in sediments: Reproducibility and comparability of results[END_REF]-2002, whereas M tenuiflora continued to expand progressively (Fig. 8, 9, 11), whereas M tenuiflora continued to expand progressively (Fig. 8, 9, 11). These results show that at the decadal scale Caatinga tree taxa have generally always been well represented irrespective of the mean annual climatic conditions and highlight the resilience of this ecosystem, which is well adapted to the alternance of drought and high moisture episodes. Therefore, rainfall variability was not the primary driver of vegetation and land cover and human activies have to be taken into consideration.

At a multidecadal scale, the progressive increase in tree pollen taxa C leprosum at ACA and M tenuiflora at ARA (Fig. 7,9) followed the same pattern as tropical Atlantic SST (Fig. 11).

Given the close link between SST and climate in the Ceará [START_REF] Alves | Relationship between ocean climatic variability and rain-fed agriculture in northeast Brazil[END_REF][START_REF] Hounsou-Gbo | Tropical Atlantic contributions to strong rainfall variability along the Northeast Brazilian coast[END_REF], we can infer that these plant taxa responded to regional changes in temperature and/or rainfall. Thus, the internal floristic composition of the Caatinga is likely to evolve progressively as a function of changes in the SST.

Human impacts

Bromine (Br) in the sediment was associated with organic matter and used as a proxy for organic content [START_REF] Bajard | Erosion record in Lake La Thuile sediments evidences montane landscape dynamics through the Holocene[END_REF], confirmed by the correlation between LOI550 and Br contents (Fig. 5) with a similar increase in the upper part of ACA15 around 1987 and since 2003. XRF concentration of manganese (Mn) is an indicator for redox processes and rapid changes in oxic-anoxic conditions at the water-sediment interface [START_REF] Davison | Iron and manganese in lakes[END_REF][START_REF] Deflandre | Early diagenetic processes in coastal marine sediments disturbed by a catastrophic sedimentation event[END_REF][START_REF] Elbaz-Poulichet | Sedimentary record of V, U, Mo and Mn in the Pierre-Blanche lagoon (Southern France) -Evidence for a major anoxia event during the Roman period[END_REF]. Abrupt changes can be caused by flood events that may favor re-oxygenation of bottom water [START_REF] Wilhelm | 1.4 kyrs of flash flood events in the Southern European Alps : implications for extreme precipitation patterns and forcing over the north western Mediterranean area[END_REF]. In organic rich lake environments with hypoxic conditions caused by the decomposition of organic matter, an increase in Mn may be observed during abrupt re-oxygenation of the sediment/water interface linked to mass movement of water [START_REF] Wirth | Combining sedimentological, trace metal (Mn, Mo) and molecular evidence for reconstructing past water-column redox conditions: The example of meromictic Lake Cadagno (Swiss Alps)[END_REF]. In ACA15, we observed both indicators (organic matter and oxygenation; Fig. 11) since 2003 suggesting that this lake system experienced a change in trophic status with an increase in in-lake productivity and more hypoxic conditions. Consequently two distinct lacustrine phases were observed: 1) a first one between 1961 and 2003 with fine sedimentation and flooding evidenced by terrigenous material and occasional precipitation of Mn caused by re-oxygenation due to flood events [START_REF] Wilhelm | 1.4 kyrs of flash flood events in the Southern European Alps : implications for extreme precipitation patterns and forcing over the north western Mediterranean area[END_REF]; 2) a second phase between 2003 and 2015 with a change in trophic status possibly related to the development of fisheries [START_REF] Gurgel | Fisheries in semi-arid Northeast Brazil with special reference to the role of tilapias[END_REF] and to the increase in population around the lake (20% increase between 1999 and 2003) (Fig.

Effects of public policies

The date of occurrence of the different events discussed below integrate the error interval defined by the age model when comparing with the launching of the public policies (Fig. 2).

In ARA, before 1958, when the reservoir was not yet active, low tree pollen frequencies (less than 5%) characterized a dry open landscape with abundant dry herbaceous taxa (mainly Amaranthaceae Alternanthera sp. and Rubiaceae Mitracarpus sp.) (Fig. 8,9). This interval reflects the nature of the landscape before and during the construction of the dam that started in 1951 and was inaugurated in 1958. The decrease in Zr content, a geochemical indicator associated with coarse element represented by medium sand (facies 3 and Fig. 4, 11) was evidence for the installation of the lacustrine system due to the construction of the dam. After the dam began functioning, the main changes in the vegetation between AD 1963 and[START_REF] Stockmarr | Tablets with spores used in absolute pollen analysis[END_REF] were reforestation with the expansion of the Caatinga tree species mainly M caesalpiniifolia, M tenuiflora and Myracrodruon sp. and an abrupt decrease in the herb taxa. Two more humid periods in 1987 and 1994 were visible through geochemical proxies (Fig. 11). In 1986, the irrigation law for northeastern Brazil was passed (Lei 10 Septembre 1986 Programa de Irrigação do Nordeste -PROINE) resulting in a significant extension of the cultivated area with, in particular in the hydrolic basin of Varjota, the Araras Norte Project that included 1 346 irrigated ha destined for fruit production [START_REF] Ximenes | O projeto Araras Norte em meio a seca no sertão revelando as fragilidades dos perímetros irrigados implantados no semiarido nordestino[END_REF]. The highest concentration of Br recorded in 1987 related to an increase in organic matter in the lake likely due to the erosion of soil litter. The implementation of the irrigation plan was associated with a sharp prolongation of the local deforestation as evidenced by the decrease in M caesalpiniifolia and M tenuiflora, more frequent microscopic charred particles associated with increased burning (Fig. 9) and the expansion of cultivated crops (Fig. 9). Remote sensing reconstructions for the year 1999 show an increase in natural vegetation cover related to high moisture rates while in 2016, after four years of very low rainfall, a significant reduction in vegetation cover and an increase in bare soil emphasize the limits of the above mentionned irrigation plan (Fig. 12 for the use of land located upstream and downstream of the dam. The purpose of the land concession agreement for agriculture and livestock was to settle farming families in the vicinity of the public dams. Confirming fossil data, the satellite image from 1985 clearly shows the impact of the law (Fig. 12, Table 1). Reforestation began in 1988, but the years 1990 to 1993 showed negative standard deviation of the annual mean precipitation, culminating in another severe drought in the year 1993. Here again, analyses of a satellite image acquired in 1991, which shows the return of a natural vegetation cover and a reduction in the extent of bare soil, support the significant rise in tree pollen taxa (30-61%) mainly M caesalpiniifolia between 1990 and 1997 (Fig. 12, Table 1). In 1989, the return of vegetation cover was favored by the abandonment of cultivated land due to another law (National Institute for Colonization and Agrarian Reform INCRA Law 97.844 19/07/1989).

Conclusions

Our study demonstrates the possibility of reconstructing the disturbances caused by human activities and climatic changes in a hydrographic basin using sedimentary archives when historical data are inexistent, rare or damaged. Despite similar climate variability, major differences were observed between the two sites due to specific local human occupation and development policies. The ability to monitor these trends at a variety of scales provides crucial information for sustainable resource management decisions that are essential for the maintenance and improvement of human wellbeing in the BSA [START_REF] Wu | Landscape sustainability science: Ecosystem services and human wellbeing in changing landscapes[END_REF]. Three levels of changes were observed depending on the duration of their impact on the landscape, annual, decadal or pluridecadal. In the Acaraú basin, the high climatic variability in recent decades did not seem to affect the composition of the Caatinga vegetation due to its extremely high resilience. As models have difficulty predicting the response of this natural rainfall regime to increasing land and sea temperatures [START_REF] Marengo | Regional differences in aridity/drought conditions over Northeast Brazil: present state and future projections[END_REF], accurate monitoring is needed before appropriate solutions for both people and the environment can be proposed. We emphasize that public policies play a significant role in changes in the landscape and in resources in the long term, whereas interannual climate variability and gradual climate change have not affected the vegetation cover of the Caatinga in the last 60 years. It is now crucial to define new forms of management that account for the physical, societal and sustainability specificities of the BSA and, to promote sustainable land use when scarce water resources and increasing temperatures threaten all the living organisms in these extreme conditions. For ARA, Si contents illustrate high terrigenous content related to high rainfall events that were considered as instantaneous events and removed for age modelling (grey bands). data in the ARA15 (left panels) and ACA15 cores (righ panels). Upper panels display correlation circles for each PCA in the plan 1-2, with interpretations of each group of elements as representative of inlake process/organic sediment; generic detrital sediment; fineand coarse-grained detrital sediment. Lower panels reprensent the loadings of each sample in the plan 1-2, colored according to the sedimentological facies they belong to: black dots, facies 1, brown to green organic rich fine-grained sediment ; red dots, facies 2, grey finegrained sediment ; green dots, facies 3, brown coarse sediment. 3). A) Zr/Ti represents sedimentation before the dam was constructed and Si/Ti (B) identifies erosive events after the dam was finished. Fe/Mn (C) has been added to highlight redox processes and Br/Ti (D) was used as a proxy of organic matter content. LOI550 (E) has been added to show the covariation between LOI550 and bromine. Finally, these all proxies are compared with annual precipitation rates in the Acaraú basin, Ipaguassu station close to ACA (dashed line) and Varjota station (full line) near ARA (F).

  BSA in 2012 was the driest between 1961 and 2012, and the two driest years since 1961 were 1982 and 2012. During the dry/wet years, inter hemispheric gradient is steep and the ITCZ stays in a northern/southern position while the SST of the equatorial Pacific is anomalously warm/cold. Model predictions show a scenario with a warmer SST in the northern Atlantic, less rainfall and warmer temperatures in the BSA.

  10, 11). Many reservoirs created in the BSA in the second half of the 20 th century currently suffer from eutrophication and no longer provide drinkable water (CGEE 2017). Today, the region around ACA is a vast area of intensive agriculture. Pollen of Arecaceae, Cecropia sp. (a pioneer indicator of deforestation of the semi-deciduous forest) and Anacardium occidentale (cashew) appeared often in the last 20 years, while pollen of Sapindaceae, Copaifera sp. and Attalea sp. (Babaçu a species of palm that produces oil) were found for the first time after 1963. Cashew and native babaçu palm are cultivated in the Serra da Meruoca for commercial purposes (Coradin et al. 2018). Despite the creation of the Serra da Meruoca Environmental Protection Area in 2008 (Law N° 11891), some families continued to practice subsistence farming there. At ARA local deforestation are shown by the increase in charcoal particles with a peak in 1983 and of cultivated taxa such as Zea mays (Fig. 8, 9).

  ).AtACA, between 1980 and 1988, the vegetation cover became more open with a decline in trees and shrubs, and more particularly in M caesalpiniifolia despite years with abundant rainfall after the severe drought of 1982-83 (ACA-III in Appendix B; Fig. 10, 11). Human activities were triggered by the award of public contracts (Portaria nº 3 -DGO 22/02/1985)

Fig 1 .

 1 Fig 1. Location of the study area. A) Map of South America showing the four main Brazilian biomes. B) Map of the State of Ceara in Brazil with the hydrographic basin of the River Acaraú (in grey). C) Representation of the landscape of the Acaraú basin (data from http://funceme.org) with the location of the two sediment cores ACA and ARA.

Fig 2 .

 2 Fig 2. Presentation of the age model. 210 Pbex activities for (A) the ACA and (b) ARA cores.

Fig. 3 .

 3 Fig. 3. Definition of the facies. Principal component analyses (PCAs) of geochemical (XRF)

Fig. 4 .

 4 Fig. 4. Downcore geochemical variations. Results for ARA15 (red curve) and ACA15 (blue curve) are plotted against time. The elements plotted here represent the processes highlighted by principle component analysis (Fig. 3). A) Zr/Ti represents sedimentation before the dam

Fig. 5 .

 5 Fig. 5. Graph showing the correlation between Br and LOI in core ACA 15.

Fig. 6 .

 6 Fig. 6. Pollen diagram of the ACA 15-2 core. Groups of dominant taxa are classified as woody plants, herbs, plants growing on the lake margin, seasonal semi-deciduous forest (Serra da Meruoca) and aquatic plants. Pollen zones were defined based on the results of the cluster analysis. Hollow curves correspond to an exaggeration x10.

Fig. 7 .

 7 Fig. 7. Summary pollen diagram of the ACA 15-2 core. Representation of the frequencies of ecological pollen groups together with the number of charcoal particles in each sample, as well as the cluster analysis dendrogram. The three lithological units were defined based on the geochemical analyses (Fig. 3, 4).

Fig. 8 .

 8 Fig. 8. Pollen diagram of the ARA 15 core. Groups of dominant taxa are classified as woody plants, herbs, seasonal semi-deciduous forest (SSDF) -Serra Ibiapaba plants growing on the lake margin, and aquatic plants. Pollen zones were defined based on the results of the cluster analysis. Hollow curves correspond to an exaggeration x10.

Fig. 9 .

 9 Fig. 9. Summary pollen percentage diagram of the ARA 15 core. Representation of the ecological groups, the number of charcoal particles and the cluster analysis dendrogram. The three lithological units were defined based on the geochemical analyses (Fig. 3, 4)..

Figure

  Figure 2

  Figure 2

Figure 10

 10 Figure 10

Acknowledgements

We are grateful to COGERH at Massapê, DENOCS and the fishermen of Varjota for letting us use their boats, Maria Iracema Loiola for information about plants of the Caatinga and for permission to collect flowers at the herbarium for our reference pollen collection, Ligia Queiroz Matias for information about plant species in the semi-arid wetlands, Vaneicia Gomes for preparing some pollen slides, Sandrine Canal for welcoming and supervising the students in the pollen laboratory at ISEM.

Fundings

This research was funded by a "sciences without boarder" (portugues acronym CSF) project « Evolution of biodiversity loss of areas in desertification processes » (CNPq 400890-2014-3), INCT IN-TREE: Interdisciplinary and transdisciplinary studies in ecology and evolution (CNPq CAPES FAPESBA) and covered a visiting foreign professor fellowship for MPL at the Federal University of Ceara (UFC) in the frame of the CSF project (2014-2018) and a post doctoral grant (CNPq/FUNCAP) for VJP at the UFC (2015-2017) and for FRGB (PNPD/CAPES Program) at the PPGERN-UFC (2017-2019). This research is also part of the Labex-CEBA.

Data availability

The sequences reported in this paper will be deposited at the french national data repository Cyber carothèque https://cybercarotheque.fr/ with a link to PANGEA data base.

Author contribution

MPL and FSA defined the research objectives VJP performed the pollen analyses LB collected the samples PS, ALD performed the geochemical analyses, FRGB performed the charcoal analyses, ESPRM provided the climatic data, MRFF, DPF performed the image analyses, CF built the synthetic figures All the authors wrote and commented on the manuscript. All the authors read and approved the manuscript for submission.

List of Table Table 1. Extent of natural vegetation, land use and reservoir from analysis of the satellite images acquired in the vicinity of ACA and ARA cores and presented in Fig. 12.

List of Figures

Fig. 10. Main results. In the Acaraú basin A) mean annual rainfall (http://funceme.org), B) demography (from the Brazilian institute for geography and statistics IBGE https://www.ibge.gov.br), C) in-lake productivity inferred from principal component analysis axis 1 scores (Fig. 3), D) arboreal pollen frequency in core ACA 15, E) launching of public policies in the region of ACA, and in the Araras basin with F) mean annual rainfall (http://funceme.org) G) demography (from the Brazilian institute for geography and statistics IBGE https://www.ibge.gov.br), H) in-lake productivity inferred from principal component analysis axis 1 scores (Fig. 3), I) arboreal pollen frequency of core ARA 15, J) launching of public policies in the region of ARA. Orange bars identify the drought years.

Table

Table S1


Myriam Tami

Marianne Clausel

Emilie Devijver

Adrien Dulac

Eric Gaussier

Stefan Janaqi

Meriam Chebre

Uncertain Trees: Dealing with Uncertain Inputs in Regression Trees

. To the best of our knowledge, no extension has been provided yet for dealing with uncertainties in the input variables, even though such uncertainties are common in practical situations. We propose here such an extension by showing how standard regression trees optimizing a quadratic loss can be adapted and learned while taking into account the uncertainties in the inputs. By doing so, one no longer assumes that an observation lies into a single region of the regression tree, but rather that it belongs to each region with a certain probability. Experiments conducted on several data sets illustrate the good behavior of the proposed extension.

Introduction

Most measures carried out in the real world, e.g., by sensors embedded in different machines or by analyses of samples, are uncertain if not erroneous in some cases. This uncertainty may be due to the generating process of the samples being measured or from the intrinsic limitations of any measurement process. Considering such measures, that constitute many of the data sets used in data science applications in both industry and academy, as certain is thus in the best case a naive position. Our experiments illustrate this point inasmuch as the method we propose here to handle uncertainty outperforms standard approaches on several, real benchmark data sets.

However, if several methods have been developed to obtain uncertainty estimates from data sets, very few studies have been devoted to designing data science methods that can deal with such uncertainties. We * Institute of Engineering Univ. Grenoble Alpes address this problem here in the context of regression trees, a popular machine learning method at the basis of widely used ensemble methods as Random Forests.

In this context, recent studies have focused on providing conditional quantiles, as opposed to conditional means, so as to better represent the output variable and avoid the uncertainty inherent to point estimates. The work on quantile regression forests developed by Meinshausen [START_REF] Meinshausen | Quantile regression forests[END_REF] is a good illustration of this. We take here a different approach and directly model the uncertainty of the input variables in the regression trees we consider. This leads to a regression tree in which observations no longer belong to a single leaf. Instead, each observation has a non-null probability of being assigned to any leaf of the tree. The construction process associated to such trees is slightly more complex than the ones of standard trees (it involves in particular the inversion of a K × K matrix, where K is the number of leaves of the tree), but the improvements obtained on the prediction fully justify this additional complexity, as shown in Section 4 on both benchmark and modified (with an additional noise) data sets.

The idea of including information on the uncertainty of the input data in a regression method is not new and is related to uncertainty propagation and sensitivity analysis. Several authors have indeed proposed to integrate uncertainties on the inputs in different regression methods, as multivariate linear models [START_REF] Reis | Integration of data uncertainty in linear regression and process optimization[END_REF] or neural networks [START_REF] Gal | Dropout as a Bayesian approximation: Representing model uncertainty in deep learning In international conference on machine learning[END_REF]. To the best of our knowledge, our approach is the first one to address this problem in the context of regression trees and ensemble methods based on such trees, as Random Forests.

The remainder of the paper is organized as follows: Section 2 presents the general model we rely on and its main properties; Section 3 then describes the algorithms for constructing the regression tree and the associated prediction rule, while Section 4 presents the experiments conducted and the numerical results obtained. Finally, Section 5 positions our work wrt previous studies while Section 6 concludes the paper.

Regression trees with uncertain input

Let Y be an output random variable and X = (X1 , • • • , X p ) a p-dimensional input random vector. A classical question is to find some relationship between X and Y , estimating the so called link function f involved in the model Y = f (X) + ǫ.

Tree-based ensemble methods, as Random Forests or Gradient Boosted Trees, are popular machine learning methods, developed to address the above regression problem [START_REF] Hastie | The Elements of Statistical Learning[END_REF]. In these methods, a set of regressors are constructed and aggregated in a convenient way. The building blocks are decision regression trees [START_REF] Breiman | Classification and regression trees[END_REF], which are defined from a partition of the space X of input variables into K regions (R k ) 1≤k≤K obtained by dyadic splits minimizing a risk function. A weight γ k is associated to each region leading to a piece-wise predictor, for a new input x ∈ X , of the form:

T (x; Θ) = K k=1 γ k 1 {x∈R k } (2.1)
where Θ = (R k , γ k ) 1≤k≤K is the set of parameters, learned from a training data set, defining the tree. Both categorical and quantitative inputs can in theory be considered. For the sake of simplicity, we however focus in this study on quantitative inputs, thus considering that X = R p .

To deal with uncertainty in the inputs, we introduce an auxiliary latent random vector U representing the true value of the data and consider the general regression function that relates U |X to Y through Y = f (U |X) + ǫ. The standard regression model is obtained from this general model by considering that U j |X j is a Dirac at X j = x j (or equivalently is Gaussian with mean x j and variance 0). We further assume here that the variables are independent of each other and that the true measure given the observation, on each variable, is Gaussian, leading to the following complete model:

(2.2) Y = f (U |X) + ǫ, U j |X j = x j ∼ N x j , σ 2 Uj for 1 ≤ j ≤ p.
The Gaussian distribution is widespread, can be used to approximate several other distributions and is easy to manipulate, hence its use here. Other distributions can nevertheless be considered for the measurement error, but this is beyond the scope of this study. In the remainder, we will denote σ 2 U = (σ 2 U1 , . . . , σ 2 Up ) the vector of variances of the Gaussian distributions. In practical situations, these variances may be given (for example when the data corresponds to measurements from machines for which the uncertainty is known) or may be directly learned from the data.

When one is dealing with uncertain inputs, and we want to stress again that this is the general situation, one can no longer assume that observations are hard assigned to regions. Instead, each observation has a probability of being associated to each region, leading to the following prediction rule that directly generalizes Eq. (2.1):

T un (x; Θ) = K k=1 γ k P(U ∈ R k |X = x, σ 2 U ). (2.3)
Note that the set of parameters Θ now includes the variances σ 2 U . In addition, because of the independence assumption at the basis of the model retained, one has, for 1 ≤ k ≤ K:

P(U ∈ R k |X = x, σ 2 U ) = p j=1 P(U j ∈ R j k |X j = x j , σ 2 Uj ),
where the interval R j k corresponds to the region R k projected onto the jth variable. We now turn to the procedure for learning the parameters of the model.

Estimation procedure

The estimation procedure of the parameters is based, in this study, on the minimization of the empirical quadratic risk, which is the standard risk considered in regression trees. More precisely, for the learning set defined by:

L n = x 1 i , • • • , x p i , y i 1≤i≤n , (2.4)
with (x i , y i ) 1≤i≤n the observed sample, we define the empirical risk:

R emp (T un (.; Θ), L n ) = 1 n n i=1 (y i -T un (x i ; Θ)) 2 (2.5)
where T un has been introduced in Eq. (2.3). This criterion has to be minimized on the training set wrt the parameters denoted by Θ = (R k , γ k ) 1≤k≤K , σ 2 U : regions of the tree, associated weights, and the variances of the vector U |X. To do so, we introduce the matrix P ∈ M n,K (R) 1 defined, for 1 ≤ i ≤ n and for 1 ≤ k ≤ K, by:

P i,k = P U i ∈ R k |X i = x i , σ 2 U = p j=1 P U j i ∈ R j k |X j i = x j i , σ 2 Uj .
Estimating γ -It is easy to see that when fixing the regions (R k ) 1≤k≤K and the vector of variances σ2 U , minimizing Eq. (2.5) with respect to γ corresponds to a weighted average of y 1 , . . . , y n , in a way similar to the linear regression model if P T P is not singular:

γ =argmin γ∈R K R emp (T un (.; Θ), L n ) (2.6) = P T P -1 P T y,
where y denotes the vector of n univariate outputs y 1 , . . . , y n . Indeed, by definition:

R emp (T un (.; Θ), L n ) = 1 n n i=1 y i - K k=1 γ k P i,k 2 .
Differentiating wrt γ k , for all 1 ≤ k ≤ K, one gets:

∂ ∂γ k R emp (T un (.; Θ), L n ) = 0 ⇔ 2 n i=1 P i,k • y i - K k ′ =1 P i,k ′ γ k ′ = 0 ⇔ n i=1 P i,k • y i = n i=1 K k ′ =1 P i,k P i,k ′ γ k ′ .
So that, if P T P is not singular:

P T y = P T P γ ⇔ γ = P T P -1 P T y,
which is a minimum. In Section 2.3, we derive assumptions under which P T P is indeed not singular. In practice, one can always use the pseudo-inverse of P T P , which we will do in our experiments. Lastly, note that γ depends on the regions through the matrix P .

Estimating (R k ) 1≤k≤K -The regions are constructed in a way similar to the construction process of standard trees with the difference that here, during the construction process of the regression tree, the number of regions K is not fixed and increases step by step, which implies that the size of the matrix P is also changing during the iterative process. Let us assume that K regions have already been identified, meaning that the current tree has K leaves, each leaf corresponding to a region (i.e., hyper-rectangle). As in standard regression trees, each current region R k , 1 ≤ k ≤ K, can be decomposed into two sub-regions wrt a variable X j and a splitting point s j k , for 1 ≤ j ≤ p:

R j k,L = {x ∈ R k | X j < s j k } R j k,R = {x ∈ R k | X j > s j k }.
This decomposition adds a new region for which the associated elements, P ∈ M n,K+1 (R) and γ ∈ R K+1 can readily be computed. P has now K + 1 regions corresponding to the current regions R k ′ (with k ′ = k) and the two new regions R j k,L and R j k,R . For each current region R k , one is looking for the best split, i.e., the best variable X j and the best splitting point s j k that minimize:

n i=1 y i - K+1 l=1 γ l P i,l 2 .
The sum includes all possible observations as each observation has a non null probability to belong to any region. Using Eq. (2.6), the above problem can be reformulated, for each current region R k , as:

argmin 1≤j≤p,s∈S j k    n i=1 y i - K+1 l=1 P T P -1 P T y l P i,l 2    , (2.7)
where S j k denotes the set of splitting points corresponding to the middle points of the j th coordinates of the observations in R k sorted according to X j .

The decomposition corresponding to the best split is then used to build the child nodes of R k , which is replaced, in the set of current regions, by its two children. In this process, that is repeated till a stopping criterion is met 2 , the number of regions, the matrix P and the weights γ are gradually updated. Section 3 provides the algorithm corresponding to this construction. Estimating σ 2 U -Lastly, the vector of variances, σ 2 U , can either be set according to some high level principles, or can be learned through a grid search on a validation set guaranteeing. The latter is more demanding wrt computational resources, but is likely to lead to better results. However, in our experiments, we use the former strategy, with the aim of showing that our approach is robust in the sense that it yields good results even when σ U is set a priori.

Final prediction

For a new observation x ∈ R p , one first computes its distribution over all the obtained regions:

∀k, 1 ≤ k ≤ K, P x k = P U ∈ R k |X = x, σ 2 U .
The prediction is then a direct application of (2.3):

T un x; Θ = K k=1 γk P x k , (2.8)
where γ, (R k ) 1≤k≤K and σ 2 U are estimated as described above.

2.3 A sufficient condition on the invertibility of P T P The formula used for the construction of the tree relies on the inverse on the matrix P T P ∈ M K,K (R). Even if numerically the Moore-Penrose pseudo-inverse might be used to approximate this inverse, we provide here a sufficient condition on the invertibility of P . Without loss of generality, the regions involved in the definition of the regression tree are of the form

R k = p j=1 [a j k , b j k ] for 1 ≤ k ≤ K, where (a j k , b j k ) ∈ (-∞, +∞) 2 .
As usual, we denote q α the α quantile of the standard Gaussian distribution. The general invertibility result is stated in Theorem 2.1.

Theorem 2.1. With the model defined in (2.2), if the following assumption is satisfied:

∀j, 1 ≤ j ≤ p, σ Uj < min 1≤k≤K (b j k -a j k ) 2q 1+0.5 1 p 2
, then the matrix P T P is invertible.

Roughly speaking, the matrix P T P is invertible provided that the standard deviation σ Uj is sufficiently small for all 1 ≤ j ≤ p. The smaller the regions and/or the larger the number p of input variables, the lower the uncertainty on the measurement has to be to ensure the invertibility of the matrix P T P .

To prove Theorem 2.1, we prove that P is of full rank. To do so, we first prove the following result:

Proposition 2.1. Let us fix 1 ≤ k ≤ K and consider 1 ≤ i ≤ n such that a j k < x j i < b j k . Assume that (2.9) ∀1 ≤ j ≤ p, σ Uj < b j k -a j k 2q 1+0.5 1/p 2 .
Then, P i,k > 0.5.

Proof. Observe that a sufficient condition to have

P i,k > 0.5 is ∀j, 1 σ Uj √ 2π R j k e - (u-x j i ) 2 2σ 2 U j du > 0.5 1/p .
We now search a sufficient condition to have the inequality just above. To get this inequality the following condition is sufficient: for all 1 ≤ j ≤ p,

P U j i -x j i σ Uj < b j k -x j i σ Uj |X i = x i -P U j i -x j i σ Uj < a j k -x j i σ Uj |X i = x i > 0.5 1/p . (2.10)
Note that this last condition is satisfied if we have, for all 1 ≤ j ≤ p,

   P U j i -x j i σU j < b j k -x j i σU j |X i = x i > 1+0.5 1/p 2 , P U j i -x j i σU j < a j k -x j i σU j |X i = x i < 1-0.5 1/p 2 ,
which is equivalent to

a j k -x j i σ Uj < q 1-0.5 1/p 2 and b j k -x j i σ Uj > q 1+0.5 1/p 2 . Note that q 1-0.5 1/p 2 = -q 1+0.5 1/p 2 < 0. A sufficient condition is then, for all 1 ≤ j ≤ p, σ Uj < min 1≤k≤K    min x j i -a j k , b j k -x j i q 1+0.5 1/p 2    .
Since a j k < x j i < b j k , a condition even more conservative is the following:

∀j, σ Uj < b j k -a j k 2q 1+0.5 1/p 2 .
This concludes the proof.

We assume in the following that Assumption (2.9) is satisfied. Then, the set I k = {1 ≤ i ≤ n, P i,k > 0.5} is non-empty. Let us consider, for all 1 ≤ k ≤ K, i k a representative of this set and let us introduce the matrix Q ∈ M K,K (R) defined, for all 1 ≤ l, k ≤ K, by:

Q l,k = P i l ,k .
Proposition 2.2. Assume that Assumption (2.9) holds. Then the K × K matrix Q is invertible.

Proof. We first show that Q is a strictly dominant diagonal matrix, i.e.:

∀k, Q kk > k ′ =k Q kk ′ .
Indeed, by Proposition 2.1, we know that P i k ,k > 0.5 and k is the only region where it is true:

l =k Q l,k ≤ l =k P i l ,k = 1 -P i k ,k < 0.5 < P i k ,k = Q k,k .
According to Hadamard's Lemma, we know that a strictly dominant diagonal matrix is invertible, which concludes the proof.

As Q is invertible, P is of full rank, leading to the fact that P T P is invertible which concludes the proof of Theorem 2.1.

Extension to Random Forests

It is straightforward to use the uncertain regression trees introduced above in Random Forests, leading to Uncertain Random Forests. Indeed, each tree of the forest is now an uncertain regression tree that can be constructed as outlined before. Assuming a forest of τ uncertain trees and denoting γt , 1 ≤ t ≤ τ , the weights estimated for each tree and P x,t kt the probability distribution of a new observation x ∈ R p over the K t regions (1 ≤ k t ≤ K t ) of the t th tree (1 ≤ t ≤ τ ) of the forest, the prediction rule of the uncertain random forest takes the form:

RF un x; Θt, 1≤t≤τ = 1 τ τ t=1 Kt kt=1 γt kt P x,t kt . (2.11)
As one can note, the above prediction rule is a direct extension of Eq. (2.8).

3 Associated algorithms Algorithm 3.1 describes the construction of uncertain regression trees. This construction parallels the one of standard regression trees except that we consider a matrix P encoding the probability distribution of training data points over regions, which is built dynamically (as the regions and the weights) by adding a new column when a given region is split into two sub-regions.

Each time a region R, corresponding to a current leaf of the tree being constructed, is considered (through the pop function in Algorithm 3.1), one identifies the best split (j ⋆ , s ⋆ ) that minimizes the empirical risk in Eq. (2.7) among all possible splitting points in S j (R) of each variable j. The set S j (R) is defined by S j (R) = {(x j i l+1 + x j i l )/2 | 1 ≤ l ≤ r -1}, with (x i1 , . . . , x ir ) corresponding to the r observations of the learning set L n belonging to the region R, sorted such that x j i1 ≤ x j i2 ≤ • • • ≤ x j ir . The k-th column of P corresponding to the current region is finally replaced by the probability distribution of the training data points over its left sub-region (denoted R L ) and an additional column is added to P for the right sub-region (denoted R R ). The weights γ are easily deduced from P at each step through (2.6).

The algorithm finally outputs the set of regions R(S) corresponding to the leaves of the tree and the weights γ.

Algorithm 3.1. Uncertain regression trees

Input: L n Initialization: K = 1, S ← (1, X ), P = 1 n while stopping criterion not met (k, R) = S.pop() for j = 1 to p Construct S j (R)

for s ∈ S j (R) R L = {x ∈ R, x j ≤ s} R R = {x ∈ R, x j ≥ s} Q L = (P(U i ∈ R L |X i = x i )) 1≤i≤n Q R = (P(U i ∈ R R |X i = x i )) 1≤i≤n P [j,s] ← Merge {P [, 1 : k -1], Q L , P [, k + 1 : K], Q R } end for end for (j ⋆ , s ⋆ ) = argmin 1≤j≤p,s∈Sj (R)
risk(P [j,s] ) #defined in (2.7)

P = P [j ⋆ ,s ⋆ ]
Update γ acc. to Eq. 2.6 K = K+1 S.append((k,R L ), (K, R R )) end while Output: (R(S), γ) Note that in this version the tree is constructed in a depth-first manner, and that the usual stopping criteria for regression trees can be used (as imposing a minimum number of observations in each leaf or a maximal depth for the tree). Algorithm 3.2 provides the pseudo-code for the prediction rule given in (2.8) for a new covariate x pred .

Algorithm 3.2. Prediction

Input: (R k ) 1≤k≤K , γ, x pred , σ 2 U ŷpred = 0 for k = 1 to K ŷpred + = P(U ∈ R k |X = x pred , σ 2 U )γ k end for Output: ŷpred 4 Experimental validation
The goal of our experiments is to assess how uncertain regression trees:

• compare to standard regression trees,

• behave in Uncertain Random Forests.

By standard regression trees we mean here regression trees based on the quadratic risk and the prediction rule given in Eq. (2.1). In addition, we consider a trade-off between standard and uncertain regression trees based on standard trees (and thus avoiding the complex construction process outlined before) but using the prediction rule given in (2.8). The matrix P and the weights γ are thus computed only once, when the standard trees have been built. The rationale for this method is to rely on the strengths of the two approaches: a simple construction tree process and a robust prediction rule. As one can conjecture, this method yields results in between the two approaches.

We make use in our experiments of four benchmark data sets commonly used in regression tasks and described in Section 4.1. We first use these data sets without any modification, to illustrate the fact that real data sets contain uncertain inputs. The uncertain regression trees proposed here indeed outperform standard trees and Random Forests on these data sets, as shown in Section 4.2. We then modify these data sets by adding a uniform noise bounded by a quantity proportional to the empirical variance of the data. This additional perturbation aims at assessing the robustness of the different methods (standard and uncertain trees) in situations where the input data is highly uncertain. Once again, the results obtained show the good behavior of the uncertain trees and Random Forests (Section 4.3). In all our experiments, the results are evaluated with the Root Mean Squared Errors (RMSE), which is a standard evaluation measure in regressions tasks. To ensure that all the available data is used for both training and testing, we further rely on 5-fold cross-validation and report the mean RMSE and its standard deviation over the 5 folds.

The stopping criterion for the trees, both standard and uncertain, is based on the fact that all leaves should contain at least ten percent of the training data. For Random Forests, three features are randomly selected for constructing a tree, which roughly corresponds to one third of the features on the data sets considered, a standard proportion in Random Forests.

In this study, the vector of variances is fixed according to some high level principle. In particular, when no additional noise is introduced on the data, the variance for a particular variable U j |X j is set to the empirical variance of X j (in other words, we assume that the variance of the true values is of the same order as the variance of the observed values). When some noise is added to the data, the variance of U j |X j is set to one half of the variance of the observed, noisy data (in this case, the variance of the true values should be lower than the empirical variance of the observed values; we arbitrarily chose one half in this study).

Lastly, our algorithms have been implemented from scikit-learn [START_REF] Varoquaux | Scikit-learn: Machine Learning Without Learning the Machinery[END_REF], using Cython [START_REF] Behnel | Cython: the best of both worlds[END_REF].

Data Sets

Experiments are conducted on 4 publicly available data sets, popular in the regression tree literature. The main characteristics of these data sets are summarized in Table 1. As we focus in this paper on quantitative variables, we simply deleted the categorical variables from the data sets. Several applications are considered, among which environment (concentration in ozone over a day for the data set Ozone introduced in [START_REF] Cornillon | R for Statistics[END_REF]), health (data set Diabete, introduced in [START_REF] Varoquaux | Scikit-learn: Machine Learning Without Learning the Machinery[END_REF] and used in [START_REF] Efron | Least angle regression[END_REF]), economy (data set BigMac about price of sandwiches, available in R package alr3 and used in [START_REF] Meinshausen | Quantile regression forests[END_REF]) or biology (data set Abalone introduced in [START_REF] Nash | Ford The Population Biology of Abalone ( Haliotis species) in Tasmania. I. Blacklip Abalone ( H. rubra ) from the North Coast and Islands of Bass Strait, Sea Fisheries Division[END_REF] and used in [START_REF] Meinshausen | Quantile regression forests[END_REF]

Results on benchmark data sets

As mentioned in the introduction, data sets are by nature uncertain, so we illustrate our method on the four data sets introduced in Section 4.1. The noise reflects the uncertainty, so we propose to use the empirical standard deviation vector as the input parameter σ U . Table 2 displays the results obtained for each data set in terms of mean and standard deviation of RMSE using 5-fold cross validation. Standard trees, standard Random Forests (with τ = 100 trees) as well as uncertain trees and standard trees with uncertain prediction are compared. For all data sets, the best performances are achieved by uncertain trees. As expected, the standard Random Forest are performing better than considering only one tree. However, as one can see in Table 2, uncertain trees yield better results than standard Random Forests. Performances of standard trees with uncertain prediction are better than the ones of standard trees, however not always better than the ones of standard Random Forests.

Finally, one can note that the standard deviation of uncertain trees is smaller than the standard deviation of other methods, meaning that uncertain trees yield more stable results.

Results on artificial uncertain data sets

To assess the robustness of uncertain trees to and uncertain Random Forests to uncertainty in the input variable, we add artificial noise, which could correspond to some measure error, to the input observations. By doing so, one can consider that the variability of the data is coming from two sources, on the one hand the variance in the latent variables (related to the variance of U with the notations of Section 2) and on the other hand the variance of the uncertainty (related to the variance of X|U), both leading to the variance of the observations. Then, we assume here that σ U is smaller than the standard deviation of the observations σ X , which can be estimated through the data set. the main part of the variance is due to the uncertainty.

To construct artificial uncertain data sets satisfying this condition, we consider in this section BigMac, Ozone, Diabete and Abalone data sets introduced in Section 4.1. A noise is added to each observation using the following process. For each observation from an input variable X j , 1 ≤ j ≤ p, we add a noise generated from the product of a Rademacher variable and a uniform variable coming from the interval [σ Xj /10, σXj /4].

We consider here standard trees, standard Random Forests with τ = 500 trees, standard trees with uncertain prediction, uncertain trees and uncertain Random Forests with τ = 15 trees. The results obtained are displayed in Table 3.

As one can note, uncertain trees outperforms here again standard trees. RMSE scores are a bit higher than in Table 2, which is not surprising given the noise added to the data sets.

Similarly, uncertain Random Forests outperform all the other methods, including standard Random Forests even though the number of trees is one order of magnitude less.

Discussion

Regression trees have been introduced in the 1980s through the popular CART algorithm [START_REF] Breiman | Classification and regression trees[END_REF], notably allowing one to deal with both categorical and numerical input variables. They constitute the basic building block of state-of-the-art ensemble methods based on the combination of random trees, notably Random Forests introduced by Breiman in [START_REF] Breiman | Random forests[END_REF] to circumvent the instability of Regression trees [START_REF] Gey | Boosting and instability for regression trees[END_REF]. Since Random Forests are particularly well suited for Big Data analysis (see [START_REF] Genuer | Random Forests for Big Data[END_REF]), many applications have been addressed in various domains with Random Forests, for example in ecology or genomics. [START_REF] Verikas | Mining data with random forests : A survey and results of new tests[END_REF] provides a review of the use of Random Forests for data mining purposes. An interesting feature of Random Forests is the possibility to quantify the importance of input variables in the model obtained (see [START_REF] Genuer | Variable selection using random forests[END_REF] for more details on that point). Another interesting feature, which was empirically established, is their robustness to noise. They are thus, to a certain extent, robust to uncertain inputs (even though no mechanism was specifically designed for that). This said, explicitly modelling the uncertainty as done in the uncertain regression trees proposed here allows one to outperform the standard version of Random Forests, as illustrated in our experiments.

Several adaptations of ensemble methods for quantile regression have been proposed, as quantile Random Forests or quantile boosting trees [START_REF] Fenske | Identifying risk factors for severe childhood malnutrition by boosting additive quantile regression[END_REF][START_REF] Kriegler | Boosting the quantile distribution: A cost-sensitive statistical learning procedure[END_REF][START_REF]Small area estimation of the homeless in Los Angeles: An application of cost-sensitive stochastic gradient boosting[END_REF][START_REF] Meinshausen | Quantile regression forests[END_REF][START_REF] Qboost | Predicting quantiles with boosting for regression and binary classification[END_REF]. These studies however focus on the uncertainty in the output variable (by producing conditional quantiles rather than a conditional mean) and not on the uncertainties in the inputs, as done here. It is of course possible to combine both approaches, which we intend to do in the future.

Lastly, the idea of including information on the uncertainty of the input data in a regression method is not new and is related to uncertainty propagation and sensitivity analysis. Several authors have indeed proposed to integrate uncertainties on the inputs in different regression methods, as multivariate linear models [START_REF] Reis | Integration of data uncertainty in linear regression and process optimization[END_REF] or neural networks [START_REF] Gal | Dropout as a Bayesian approximation: Representing model uncertainty in deep learning In international conference on machine learning[END_REF]. In each case, the methods have been improved, showing the benefits of explicitly modelling uncertainties in the input data. To the best of our knowledge, our approach is the first one to address this problem in the context of regression trees (and ensemble methods based on such trees). Our conclusion on the benefits of this approach parallels the ones of the above-mentioned studies.

Conclusion

We have introduced in this study uncertain regression trees with can deal with uncertain inputs. In such trees, observations no longer belong to a single region, but rather have a non-null probability to be assigned to any region of the tree. This extra flexibility leads nevertheless to a construction process that is more complex than the one underlying standard regression trees and that necessitates the inversion of a square matrix (for which we have theoretically provided a sufficient condition). In practice, we rely on the pseudo-inverse of ]. For each uncertain tree-based method, σ U is to half of the empirical standard deviations of the observed (modified) input variables. this matrix. The experimental results fully justify the approach we have proposed and show that uncertain regression trees improve the results of standard regression trees and Random Forests on several benchmark data sets. A similar conclusion is drawn on artificial uncertain data sets obtained from the standard ones by introducing additional uncertainty in the form of a uniform noise.

The methodology developed in this study can also be adapted to the case where some input data are categorical. We plan to work on such an adaptation in the future, by considering different types of uncertainties.

Furthermore, as mentioned before, the vector of variances for modelling uncertainties can easily be learned by grid search on validation sets (typically using cross-validation). One can expect by doing so that the results would further improve. We have set this vector in our experiments to values that we believe are reasonable, so as to show that our approach is robust in the sense that it yields good results even when σ U is set a priori. We nevertheless plan to run additional experiments to learn this vector. As a grid search can be easily parallelized, this learning should not impact the running time of the algorithm.

Lastly, exploring the extension of our method to Random Forests or boosting trees, as in [START_REF] Freund | A short introduction to boosting[END_REF][START_REF] Ridgeway | Generalized Boosted Models: A guide to the gbm package[END_REF], is another promising research direction. We also intend to explore the use of other empirical loss functions, as the quantile loss used in the definition of quantile Random Forests or quantile boosting trees [START_REF] Fenske | Identifying risk factors for severe childhood malnutrition by boosting additive quantile regression[END_REF][START_REF] Kriegler | Boosting the quantile distribution: A cost-sensitive statistical learning procedure[END_REF][START_REF]Small area estimation of the homeless in Los Angeles: An application of cost-sensitive stochastic gradient boosting[END_REF][START_REF] Meinshausen | Quantile regression forests[END_REF].

Table 1 :

 1 among others). Characteristics of data sets used in our experiments, ordered by sample size n.

	Data set sample size n number of variables p
	BigMac	69	9
	Ozone	112	9
	Diabete	442	10
	Abalone	500	7

Table 2 :

 2 Basically it means that Average RMSE based on 5-fold cross-validation for the 4 benchmark data sets. Standard deviations are given in parentheses. For each uncertain tree based method, σ U are set to the empirical standard deviations of the observed input variables.

	Data sets	BigMac	Ozone	Diabete	Abalone
	Standard tree	25.09 (12.3) 17.82 (4.3) 60.29 (4.3) 2.70 (0.3)
	Standard RF, τ = 100	19.78 (11.0) 15.86 (4.1) 58.18 (4.3) 2.65 (0.4)
	Standard tree with uncertain prediction 21.49 (8.7) 16.79 (4.1) 57.05 (3.7) 2.41 (0.2)
	Uncertain tree	18.74 (8.9) 15.39 (3.7) 56.56 (3.3) 2.33 (0.3)

Table 3 :

 3 (8.7) 19.37 (4.1) 60.47 (2.81) 2.54 (0.17) Standard tree with uncertain prediction 21.68 (9.8) 17.35 (4.9) 57.92 (3.43) 2.40 (0.15) Uncertain tree 19.28 (13.4) 16.87 (6.0) 58.56 (3.45) 2.34 (0.20) Standard RF, τ = 500 19.25 (7.8) 15.72 (3.0) 59.55 (4.39) 2.64 (0.18) Uncertain RF, τ = 15 18.06 (9.3) 15.49 (3.7) 55.66 (4.31) 1.98 (0.12) Average RMSE based on 5-fold cross-validation for the 4 modified data sets. Standard deviations are given in parentheses. On each data set, each observation is modified by a noise generated from the product of a Rademacher variable and a uniform variable coming from the interval [

	Uncertain data sets	BigMac	Ozone	Diabete	Abalone
	Standard tree	22.28 σX j 10 ,	σX j 4

P depends on the number of regions considered. It is thus a dynamic matrix that evolves during the construction of the tree. For clarity sake, we do not explicit this dependence in our notation but the reader has to keep this in mind.

Any standard stopping criterion can be used here.