Adriana Reatto-Braga
email: adriana.reatto@univ-orleans.fr

Ary Bruand
email: ary.bruand@univ-orleans.fr

Fabrice Muller
email: fabrice.muller@univ-orleans.fr.

Euzebio Medrado Da Silva
email: euzebio@cpac.embrapa.br.

Osmar Abilio De Carvalho
email: osmarjr@unb.br.

Michel Brossard
email: michel.brossard@mpl.ird.fr.

Guy Richard
email: guy.richard@orleans.inra.fr

Adriana Reatto
email: reatto@cpac.embrapa.br

Eder De Souza Martins

Osmar Abílio De Carvalho

Development and origin of the microgranular structure in Latosols of the Brazilian Central Plateau: significance of texture, mineralogy, and biological activity

Keywords: Oxisol, Porosity, Termite, Soil structure, Brazil, Micromorphology, Microstructure

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

Introduction

Soil structure corresponds to the organization of the primary soil particles into secondary units that are called aggregates or peds [START_REF] Allaby | A dictionary of Earth Sciences[END_REF]. These secondary units are defined on the basis of their shape, size, and grade. They can be easily affected by the climate, biological activity, and soil management. In Brazil, Latosols occupy about 2,835,000 km 2 and 1/3 of the national territory [START_REF] Ker | Latossolos do Brasil: uma revisão[END_REF]. They are characterized by little differentiation of their horizons, low activity clay, a weak macrostructure, and a strong submillimetric granular structure [START_REF] Embrapa | Sistema Brasileiro de Classificação de Solos[END_REF]. In the Brazilian Soil Taxonomy [START_REF] Embrapa | Sistema Brasileiro de Classificação de Solos[END_REF], Latosols are defined by the presence of a diagnostic horizon Bw between the surface and 200 cm depth. This horizon is characterized by 1:1 clay minerals, variable amount of iron content, and SiO 2 and Al 2 O 3 removable by H 2 SO 4 , such as SiO 2 /Al 2 O 3 ≤ 2.2.

Most Latosols are Oxisols according to the Soil Taxonomy (Soil Survey Staff, 2006) or Ferralsols in the World Reference Base (IUSS Working Group WRB, 2006). The development of the microgranular structure was discussed in numerous studies based on the clay mineralogy, location in the landscape, parental material, and biological activity. [START_REF] Lima | Micromorfologia de "horizonte B" de latossolos do Sudeste e Sul do Brasil. In: Anais da III Reunião de Classificação, Correlação de Solos e Interpretação da Aptidão Agrícola[END_REF] stated that the development of the microgranular structure in the Bw horizon would increase with the hematite content in the clay fraction. Vidal-Torrado and Lepsch (1993) found a microgranular structure more developed on the upper slope than on the lower slope where the microgranular structure was weak or even absent in combination with a massive argillic horizon. [START_REF] Bennema | Micromorphology of some oxic and argillic horizons in south Brazil in relation to weathering sequences[END_REF] studied oxic and argillic horizons and showed an increase in the microgranular structure development with the clay content regardless clay mineralogy variation. The independency of the microgranular structure from the clay mineralogy was also observed by [START_REF] Gomes | Mineralogia, morfologia e análise microscópica de solos do Bioma Cerrado[END_REF]. On the other hand several authors showed that the microgranular structure was less developed in kaolinitic Latosols than in gibbsitic Latosols [START_REF] Ferreira | Mineralogia da fração argila e estrutura de latossolos da região sudeste do Brasil[END_REF][START_REF] Viana | Efeitos de ciclos de umedecimento e secagem na reorganização da estrutura microgranular de Latossolos[END_REF][START_REF] Schaefer | EDS/SEM study on microaggregates of Brazilian Latosols, in relation to P adsorption and clay fraction attributes[END_REF]. [START_REF] Schaefer | Brazilian latosols and their B horizon microstructure as long-term biotic constructs[END_REF] studied a set of Bw horizons with SiO 2 , Al 2 O 3 and Fe 2 O 3 removable by H 2 SO 4 such as 0.28 < SiO 2 /Al 2 O 3 < 1.17 and 80 < Fe 2 O 3 < 170 g kg -1 . They recorded a uniform microgranular structure with different degrees of roundness in spite of the great variability of parent material and mineralogy. The influence of the biological activity on microgranular structure was also discussed by [START_REF] Eschenbrenner | Contribution des termites à la micro-agrégation des sols tropicaux[END_REF]. [START_REF] Schaefer | Brazilian latosols and their B horizon microstructure as long-term biotic constructs[END_REF] and [START_REF] Schaefer | Pedogenesis on the uplands of the Diamantina Plateau, Minas Gerais, Brazil: a chemical and micropedological study[END_REF] said that the quartz grains in the centre of microaggregates were smaller than in the surroundings which indicates the biological origin including grain size segregation by termites and other soil animals. Based in this, [START_REF] Schaefer | Brazilian latosols and their B horizon microstructure as long-term biotic constructs[END_REF] concluded that the microgranular structure of Brazilian Latosols would result from long-term biological activity. [START_REF] Balbino | Comportement de la phase argileuse lors de la dessiccation dans dês Ferralsols microagrégés du Brésil: rôle de la microstructure et de la matière organique[END_REF][START_REF] Balbino | Changes in porosity and microaggregation in clayey Ferralsols of the Brazilian Cerrado on clearing for pasture[END_REF] and Volland-Tuduri (2005b) also studied the microgranular structure of this soil type and proposed that it might result from termite and ant activity. Thus the origin of the microgranular structure of Brazilian Latosols is still under discussion. The aim of this study is to analyze the characteristics of the microstructure development in Latosols that originated from different parental materials and are located on the two main geomorphic surfaces of the Central Plateau. The study also shows that the microstructure would be related to bioturbation by termites mainly.

Material and Methods

Site conditions

The Brazilian Central Plateau covers about 10% of the continuous Cerrados Biome and represents 176,000 km 2 [START_REF] Silva | Caracterização e distribuição das principais classes pedológicas do Planalto Central Brasileiro[END_REF]. Latosols occupy about 40% of its surface area (Fig. 1). The Central Plateau corresponds to two main geomorphic surfaces, the South American Surface (SAS) and Velhas Surface (VS) [START_REF] King | A geomorfologia do Brasil Central[END_REF][START_REF] Lepsch | Oxisol-Landscape relationships in Brazil[END_REF][START_REF] Motta | Relações solo-superfície geomórfica e evolução da paisagem em uma área do Planalto Central Brasileiro[END_REF][START_REF] Marques | Major element geochemistry and geomorphic relationships in Brazilian Cerrado soils[END_REF]. The SAS corresponds to a Brazilian landscape that begins as a vast peneplain produced by denudation between the lower Cretaceous and middle Tertiary. At this time, very humid climatic conditions were favorable for a deep weathering of rocks. The continent was uplifted afterwards and polycyclic stream incision carved valleys into the surface, transforming it into a dissected upland. [START_REF] Braum | Contribuição à geomorfologia do Brasil Central[END_REF] reported the tablelands (chapadas) 900 to 1200-m high of Central Brazil as corresponding to remnants of the SAS.

They appear as smoothly convex plane portions that terminates either as an abruptly escarpment or as a kind of a stepwise border on the sides of the plateau. The VS that initiates at the dissection of the SAS and is characterized by moderate slope, convex forms, and covers smaller extensions than the SAS. According to the [START_REF] Köppen | Grundriss der Klimakunde[END_REF] classification, the most representative climate of the Central Plateau is Megatermic or Humid Tropical (Aw) with the subtype savanna. It is characterized by a dry winter (medium temperature of the coldest month > 18 °C) and maximum rains in summer. The mean annual rainfall ranges from 1,500 to 2,000 mm, with the highest rainfall in January and the smallest in June, July and August (<50 mm month -1 ). The relative humidity of the air is about 75% between January and April and about 30% during the dry winter [START_REF] Assad | Veranicos na região dos cerrados brasileiros: freqüência e probabilidade de ocorrência[END_REF].

Soil selection and characterization

According to [START_REF] Reatto | Variabilidade mineralógica de latossolos da área da Embrapa Cerrados em relação aos do Bioma Cerrado[END_REF], ten Latosols (L) developed in different parent materials were selected along an approximately 350-km long regional toposequence across the SAS (L1 to L4) and VS (L5 to L10) (Fig. 2 and Table 1). Basic characteristics and hydraulic properties of these Latosols can be found in [START_REF] Reatto | Hydraulic properties of the diagnostic horizon of Latosols of a regional topossequence across the Brazilian Central Plateau[END_REF] and data on the clay mineralogy of their Bw horizons in [START_REF] Reatto | Variation of the kaolinite and gibbsite content at regional and local scale in the Latosols of the Brazilian Central Plateau[END_REF]. The Latosols L7 and L8 were studied by [START_REF] Volland-Tuduri | Direct analysis of microaggregates shrinkage for drying: Application to microaggregates from a Brazilian clayed Ferralsol[END_REF]2005a) and are similar to those studied by [START_REF] Balbino | Comportement de la phase argileuse lors de la dessiccation dans dês Ferralsols microagrégés du Brésil: rôle de la microstructure et de la matière organique[END_REF][START_REF] Balbino | Changes in porosity and microaggregation in clayey Ferralsols of the Brazilian Cerrado on clearing for pasture[END_REF][START_REF] Balbino | Change in the hydraulic properties of a Brazilian clay Ferralsol on clearing for pasture[END_REF]. The Latosols L5 and L6 were located on the upper VS, L7 and L8 on the intermediate VS, and L9 and L10 on the lower VS. Soils were described according to the field manual of [START_REF] Lemos | Manual de descrição e coleta de solo no campo[END_REF] and the Brazilian Soil Taxonomy [START_REF] Embrapa | Sistema Brasileiro de Classificação de Solos[END_REF]. About 2-m deep soil pits were dug and the top horizon (A), transitional horizons (AB and BA), and diagnostic horizon (Bw) were described. Disturbed samples were collected in every horizon. Basic soil characterization was performed on the air-dried <2-mm material according to the Brazilian standard procedures as described by [START_REF] Embrapa | Manual de métodos de análise de solo[END_REF] which in general follow the international procedures [START_REF] Klute | Methods of Soil Analysis, Part 1, Physical and Mineralogical Methods[END_REF][START_REF] Page | Methods of Soil Analysis, Part 2, Chemical and Microbiological Properties[END_REF]. The particle size distribution was determined using the pipette method after dispersion with NaOH 1N. The particle density was determined by using 95% hydrated alcohol with 20 g of air-dried soil material in a 50-ml pycnometer. The bulk density was calculated using the oven-dried mass of soil material contained in cylinders 100-cm 3 in volume (∅ = 5.1 cm, h = 5 cm). The soil pH was measured in distilled water and 1 N KCl using 1:1 mass soil to solution ratio. The cation exchange capacity (CEC) was computed as being the sum of the electric charges of Ca 2+ , Mg 2+ , and Al 3+ extracted with 1N KCl, of K + and Na + extracted with 0,05N HCl, and of H + and Al 3+ extracted with a tampon solution of Ca(CH 3 COO) 2 and CH 3 COOH at pH 7.0. The organic carbon content was analyzed by wet oxidation with 0.4N K 2 Cr 2 O 7 .

The SiO 2 , Al 2 O 3 , Fe 2 O 3 and TiO 2 content was estimated after dissolution in 1:1 H 2 SO 4 [START_REF] Vettori | As relações Ki e Kr na fração argila e na terra fina[END_REF][START_REF] Camargo | Classification, Characterization and Utilization of Oxisols, Part 2: Field Trip Background, Site and Pedon Descriptions, Analytical Data[END_REF][START_REF] Embrapa | Manual de métodos de análise de solo[END_REF][START_REF] Marques | Trade element distributions in Brazilian Cerrado soils at the landscape and micrometer scales[END_REF]. The SiO 2 ( 2 SiO S ) and

Al 2 O 3 ( 3 2 O Al S
) contents extracted with sulfuric acid were used to compute the Ki ratio and Kr ratio plus Fe 2 O 3 in according [START_REF] Resende | Uso das Relações Ki e Kr na estimativa da mineralogia para classificação dos Latossolos. In: Anais da 3. Reunião de Classificação, Correlação de Solos e Interpretação de Aptidão Agrícola[END_REF]) defined as follows:

) (

) ( 7 . 1 102 ) ( 60 ) ( 3 2 2 3 2 2 O Al SiO O Al SiO S S S S Ki × = = (1)
[ ]

) ( 6375 . 0 ) ( / ) ( 7 . 1 160 ) ( 102 ) ( 60 ) ( 3 2 3 2 2 3 2 3 2 2 O Fe O Al SiO O Fe O Al SiO S S S S S S Kr × + × = + = (2).
The kaolinite (Kt) and gibbsite (Gb) contents were determined as follows [START_REF] Resende | Mineralogia da argila de latossolos a partir do teor total de óxidos do ataque sulfúrico[END_REF][START_REF] Bruand | Analyse minéralogique quantitative d'un échantillon de sol: utilisation des donnés concernant la composition chimique de l'échantillon[END_REF]:

2 2 / SiO SiO Kt S Kt = (3) ( ) [ ] 3 2 3 2 3 2 / ) 1000 / ( O Al O Al O Al Gb Kt Kt S Gb × - = (4)
where Kt is the kaolinite content (in g kg -1 ) of the sample,

Gb

the Al 2 O 3 content of gibbsite and equaled to 654 in g kg -1 . Equations ( 3) and ( 4) assumed kaolinite and gibbsite without any substitution.

The goethite (Gt) and hematite (Hm) contents were computed by combining the two equations as follows:

Hm Gt S O Fe + × = 8989 . 0 3 2 (5) 33 . 8 / ) 50 . 3 ( ) /( - = + RI Gt Hm Hm (6) where 3 2 O Fe S
is the Fe 2 O 3 content (in g kg -1 ) of the sample recorded with H 2 SO 4 extraction, 0.8989 the specific proportion of Fe 2 O 3 in a goethite (goethite and hematite was assumed as being no Al-substituted), RI is the red index according [START_REF] Santana | Soil formation in a toposequence of oxisols from Patos de Minas region, Minas Gerais[END_REF][START_REF] Kampf | Óxidos de ferro em Latossolos do Brasil Sudeste e Sul[END_REF] and equaled to:

( )

V C M RI / + = (7)
with M a parameter related to the hue, C the chroma and V the value.

Undisturbed samples of the Bw horizon were also collected in triplicate by using cylinders of 100-cm 3 in volume. They were impregnated with a polyester resin that was diluted with styrene monomer and left 4 weeks to ensure complete polymerization. Vertical thin sections (4 x 6 cm 2 ) were produced following the method of [START_REF] Fitzpatrick | Micromorphology of Soils[END_REF] and examined under optical microscopy using a Zeiss polarizing microscope equipped with a digital camera. The pedological features of the Bw horizon were described using standard micromorphological techniques (Bullock et al., 1987;[START_REF] Stoops | Guidelines for analysis and description of soil and regolith thin sections[END_REF]. Estimation of the proportion of surface area occupied by each type of microstructure was performed visually by using the frequency charts proposed by [START_REF] Fitzpatrick | Micromorphology of Soils[END_REF]. The thin sections were also studied in scanning electron microscopy (SEM, Cambridge 90B) using the emission of backscattered electrons [START_REF] Bruand | Backscattered electron scanning images of soil porosity for analyzing soil compaction around roots[END_REF].

Results

Field observations

The Latosols studied showed a compound weak to moderate medium sub-angular blocky structure and strong fine to very-fine granular structure. Dense irregular sub-rounded volumes (DIRV) of soil with varying size and material similar to the surrounding soil were recorded in different proportions in the Latosols studied (Table 2). The DIRV were described as clay nodules in [START_REF] Reatto | Hydraulic properties of the diagnostic horizon of Latosols of a regional topossequence across the Brazilian Central Plateau[END_REF]. The biological macrocrovoids were represented by fine channels 1 to 2-mm in diameter resulting from root activity and channels varying in their size and cavities several centimeters in diameters both related to termite or ant activity (Table 2).

In L1, high content in DIRV < 10-mm in size was recorded between 100 and 200 cm depth.

Another characteristic of L1 was the high quartz grain content and iron nodules < 0.5 mm in size content in different depths. The biological macrovoids were concentrated within the upper 100-cm. In L2 and L3, the greatest DIRV content was detected between 30 and 70-cm depth that differed much in size (from 5 to 50-mm). Biological macrovoids were very few in these two profiles. In L4, any DIRV was evidenced. This profile showed also distinct mottles (7.5YR 5/6) below 130-cm depth and fragments of weathered plinthite (5YR 5/4) below 160cm depth (Bw 2 horizon). Numerous DIRV from 5 to 50-mm in size were recorded in L5 from the surface to 90-cm depth. Only a few biological macrovoids were observed in these two profiles. In L6, L7 and L8, numerous biological macrovoids were recorded in the whole profile and high DIRV content with sizes varying from 5 to 50-mm between 20 and 140-cm depth. Finally, in L9 and L10 biological macrovoids were very few and DIRV were recorded between 15 to 110-cm depth.

Physico-chemical properties

The clay content of the Bw horizons ranged from 300 to 780 g kg -1 (Table 3). A large range of bulk density (Db) was also recorded for the Bw horizons of these Latosols with 0.83 ≤ Db ≤ 1.21 g cm -3 . The pH w varied from 4.80 to 5.34 and the pH KCl from 4.00 to 6.20.

S

= 577 g kg -1 ), and small SiO 2 content, (average 2 SiO S = 165 g.kg -1 ) (Table 3). On the opposite, on the VS the Bw horizons showed small Al 2 O 3 content (average 3).

Observations in microscopy

The Bw horizons studied showed several types of microgranular structure according to their development and related inter-microaggregates voids. A strong microgranular structure with well developed microaggregates 50 to 200-µm in size that are combined in a loose assemblage with 40 to 50% of the surface area large compound packing voids (type A) was recorded in all the Bw horizons (Fig. 3a andb; Fig. 4a) except in L4. The microaggregates showed a double-spaced-porphyric and open-porphyric coarse/fine related distribution in the groundmass [START_REF] Stoops | Guidelines for analysis and description of soil and regolith thin sections[END_REF]. A microgranular structure developed as well as in the type A with similar internal characteristics but with microaggregates combined in a close assemblage (type B) with 20 to 35 % of the surface area corresponding to complex packing voids was also observed in all the Bw horizons but in smaller proportion than type A (Fig. 3c andd; Fig. 4b).

Besides types A and B, a very weak microgranular structure with a compact assemblage of microaggregates and 10 to 30 % of the surface area corresponding regular and star-shaped vughs, channels and planar voids (type C) was detected (Fig. 3e and f; Fig. 4c). Finally, a microstructure corresponding to the assemblage of partially clay coated quartz grains with small rounded microaggregates 10 to 30-µm in size (type D) with simple to compound packing voids corresponding to 30 to 40% of the surface area was also visualized (Fig. 3g and h; Fig. 4d).

Observations of the thin sections of the Bw horizons under optical microscopy showed around 95% of the surface area in L3 and L6 corresponded to type A (Table 4) with compound and complex packing voids in L3 and simple to compound packing voids in L6 [START_REF] Stoops | Guidelines for analysis and description of soil and regolith thin sections[END_REF]. In L2, the observations showed that type A was dominant (70 %) and associated with type C (30 %). The other Bw horizons studied showed a varying proportion of the types A, B and C (Table 4, Figure 5), except for L4 where the microstructure corresponded to type D alone.

Discussion

Mineralogical classification of the Latosols studied

Latosols are classified as kaolinitic Latosols or gibbsitic Latosols according to the value of the Ki ratio in the diagnostic Bw horizon. They can be also classified as sexquioxid Latosols or not according to the value of the Kr ratio [START_REF] Resende | Uso das Relações Ki e Kr na estimativa da mineralogia para classificação dos Latossolos. In: Anais da 3. Reunião de Classificação, Correlação de Solos e Interpretação de Aptidão Agrícola[END_REF]. The Ki ratio ranged from 0.3 to 1.4 and Kr ratio from 0.3 to 1.0 (Table 5). On the SAS, 0.3 ≤ Ki ≤ 0.7 and 0.3 ≤ Kr ≤ 0.6 when on the VS 0.7 ≤ Ki ≤ 1.4 and 0.6 ≤ Kr ≤ 1.0 for the Bw horizons studied. Hence, L1, L2, L3 and L4 located on the SAS and L10 on the VS were classified as gibbsitic-sexquioxidic Latosols, L5, L7, L8 and L9 on the VS as kaolinitic non-sesquioxidic Latosols and L6 on the VS as kaolinitic-sesquioxidic Latosol.

Mineralogy of the diagnostic Bw horizons

The results showed that the Kt content was smaller (196 ≤ Kt ≤ 412 in g kg -1 ) in most Bw horizons of the Latosols located on the SAS than in those of the Latosols located on the VS (399 ≤ Kt ≤ 645 in g kg -1 ) (Table 5). In contrast, the Gb content was higher (442 ≤ Gb ≤ 625 in g kg -1 ) for Latosols on the SAS in comparison with the Latosols on the VS (183 ≤ Gb ≤ 405 in g kg -1 ) (Table 5). We computed the ratio of the kaolinite content to both the kaolinite and gibbsite content (R Kt/Kt+Gb ). The data indicated that 0.2 ≤ R Kt/Kt+Gb ≤ 0.5 in the Latosols located on the SAS and 0.5 ≤ R Kt/Kt+Gb ≤ 0.8 for those located on the VS (Table 5). The Hm content ranged from 0 to 205 in g kg -1 and Gt content from 0 to 178 g kg -1 for the Bw horizons studied (Table 5). The averaged Hm and Gt content in the Bw horizon was respectively 87 and 106 in g kg -1 , on the SAS and 136 and 45 in g kg -1 , on the VS (Table 5).

We computed the ratio of the hematite content to both the hematite and goethite content (R Hm/Hm+Gt ). The results pointed that 0.0 ≤ R Hm/Hm+Gt ≤ 0.7 for the Latosols situated on the SAS, and 0.4 ≤ R Hm/Hm+Gt ≤ 1.0 for those situated on the VS (Table 5).

Types of microgranular structure, particle size distribution, and mineralogy

The analysis of the type of microgranular structure did not show any relationship with the composition of the Bw horizons studied, neither with their particle size distribution (Fig. 5a andb), nor with their mineralogy of the fine fraction (Fig. 5a andc). The Bw horizons of L1, L2, L3, L6 and L10 on the one hand, and those of L5, L7, L8 and L9 on the other hand showed kaolinite and gibbsite different contents (Fig. 5c) when the proportion of the types of microgranular structure varied independently of that distinction (Fig. 5a). These results are in line with many former studies that led to contradictory results when the development of the microgranular structure was discussed according to the clay content and clay mineralogy [START_REF] Bennema | Micromorphology of some oxic and argillic horizons in south Brazil in relation to weathering sequences[END_REF]Vidal-Torrado and Lepsch, 1999;[START_REF] Ferreira | Mineralogia da fração argila e estrutura de latossolos da região sudeste do Brasil[END_REF][START_REF] Gomes | Mineralogia, morfologia e análise microscópica de solos do Bioma Cerrado[END_REF][START_REF] Viana | Efeitos de ciclos de umedecimento e secagem na reorganização da estrutura microgranular de Latossolos[END_REF][START_REF] Schaefer | Brazilian latosols and their B horizon microstructure as long-term biotic constructs[END_REF][START_REF] Schaefer | EDS/SEM study on microaggregates of Brazilian Latosols, in relation to P adsorption and clay fraction attributes[END_REF]

Microgranular structure development and associated voids

The types of microgranular structure A, B and C were observed in different proportions in the Bw horizons studied (Table 4). They are characterized by a varying volume of voids between the microaggregates (Fig. 3 and Fig. 4a, b andc) except for L4 where the microgranular structure corresponds to type D with many voids between the clay-coated quartz grains and microaggregates (Fig. 3g and h and Fig. 4d). The volume of intermicroaggregate voids was calculated as follows:

ra v total v er v V V V int int + = (8)
where inter v V was the volume of inter-microaggregates voids (in cm 3 g -1 ), intra v V the volume of intra-microaggregates voids (in cm 3 g -1 ), and total v V , the total volume of voids (in cm 3 g -1 ). The latter was computed as follows:

) / 1 ( ) / 1 ( Dp Db V total v - = (9)
where Dp was the particle density and equal to the reciprocal of 2.65 g cm 3 . The value of ra v V int for each Bw horizon studied was computed by using the relationship recorded by [START_REF] Balbino | Changes in porosity and microaggregation in clayey Ferralsols of the Brazilian Cerrado on clearing for pasture[END_REF] for Latosols. The results showed that 0.156 ≤ ra v V int ≤ 0.234 cm 3 g -1 and 0.305 ≤ er v V int ≤ 0.585 cm 3 g -1 .

Comparison of er v

V int with the proportion of the different types of microstructure showed that both er v V int and the surface area proportion of type A were the greatest in L2, L3, L6 and L10 (Fig. 5a; Table 6). On the other hand, this comparison also showed that both er v V int and the surface area proportion of type A was the smallest in L1 and L9.

Microgranular structure development and biological activity

Analysis of the distribution of the different types of microgranular structure in microscopy at low magnification showed that type C corresponds to portions of channels or cavities walls similar to those described by [START_REF] Eschenbrenner | Contribution des termites à la micro-agrégation des sols tropicaux[END_REF] (Fig. 6). Contours of channels and cavities were particularly visible in L1, L5 and L9 (Fig. 6). Thus, type C originates from soil material kneading by termites resulting in coalescent microaggregates forming mechanically stable and water proof walls. On the other hand, the types A and B originate from an infilling activity with microaggregates resulting from fragmentation and burrowing action of termites and possibly of ants, too. The findings of this study are comparable to those of [START_REF] Miklós | Biodynamique d'une couverture pédologique dans la région de Botucatu, Brésil[END_REF] and [START_REF] Cooper | Origen of microaggregates in soils with ferralic horizons[END_REF]. These authors also worked on Brazilian Latosols and recorded the microgranular structure as originating partly from biological activity. [START_REF] Holt | Termites and soil properties[END_REF] also discussed soil bioturbation by termites during mound-building and cavities and galleries constructions. These constructions, in a long period of time, modified the soil properties and pedogenesis tends as consequence soil turnover over and soil translocation [START_REF] Holt | Termites and soil properties[END_REF]. [START_REF] Schaefer | Brazilian latosols and their B horizon microstructure as long-term biotic constructs[END_REF] proposed that the Bw horizon microgranular structure of Brazilian Latosols would result from long-term biological activity of termites. Based on this model, the structure of the Latosols studied would originate from long-term action of biotic activity, mainly of termite activity as shown in Figure 7. Such a long-term bioturbation might explain the poor vertical differentiation of the Latosols studied and the lack of relationship between Latosol composition and microgranular structure development. The DIRV observed in the soil type studied that also show a microgranular structure of type C would correspond to remnants of mound, cavities or galleries walls made up of coalesced microaggregates (Table 2 and Fig. 6).

Conclusions

The Latosols studied showed a compound strong microgranular structure and weak to moderate medium sub-angular blocky structure independently of the landscape position on the South American Surface or Velhas Surface as well as the parental material, except for the Latosol developed on a quartzite. They were classified as gibbsitic-sexquioxic, kaolinitic nonsesquioxid, or kaolinitic-sesquioxid Latosols.. The study revealed several type of microgranular structure with different proportions according to the Latosol studied but without any relationship with their texture and mineralogy. The types of microgranular structure appeared to be related to soil bioturbation conducted mainly by termites, and secondarily by ants. Thus, the weak macrostructure and strong microgranular structure of the Latosols studied, their little vertical differentiation of the horizons, and the lack of clear relationship between their landscape position and parent material characteristics would result from long-term biotic action. Tables Table 1 General characteristics of the Latosols (L) studied (modified after [START_REF] Reatto | Hydraulic properties of the diagnostic horizon of Latosols of a regional topossequence across the Brazilian Central Plateau[END_REF][START_REF] Reatto | Variation of the kaolinite and gibbsite content at regional and local scale in the Latosols of the Brazilian Central Plateau[END_REF] Table 2 Field morphological descriptions of the Latosols (L) studied Table 3 Physical and chemical characteristics of the diagnostic horizons (Bw) of the Latosols (L) studied Table 4 Micromorphological analyses of the diagnostic horizons (Bw) of the Latosols (L) studied Table 5 Mineralogy of the diagnostic horizons (Bw) of the Latosols (L) studied.

Table 6 Volume of voids resulting from the assemblage of the microaggregates (inter-aggregates voids) for the Bw horizons studied. (g kg -1 ) (g kg -1 ) (g kg -1 ) (g cm -3 ) (g cm -3 ) (g kg -1 ) (cmolc kg -1 ) (cmolc kg -1 ) (cmolc kg -1 ) (cmolc kg -1 ) (cmolc kg -1 ) (cmolc kg -1 ) (g kg -1 ) (g kg -1 ) (g kg -1 ) (g kg Db = bulk density, Dp = particle density, OC = Organic carbon, CEC = Cation exchange capacity.

Table 4 Micromorphological analyses of the diagnostic horizons (Bw) of the Latosols (L) studied (description according to [START_REF] Bullock | Handbook for soil thin section description[END_REF][START_REF] Stoops | Guidelines for analysis and description of soil and regolith thin sections[END_REF].

  Figures 1

Fig. 1 -

 1 Fig. 1 -Distribution of the Latosols on the Brazilian Central Plateau.

Fig. 2 -

 2 Fig. 2 -Location of the ten Latosols studied (L) on each geomorphic surface along a regional

Fig. 3 -

 3 Fig. 3 -Microstructure of the Bw horizon observed in optical microscopy (natural

Fig. 4 -

 4 Fig. 4 -Microstructure of the Bw horizon observed in scanning electron microscopy using the

Fig. 5 -

 5 Fig. 5 -Different types of microstructure (a), particle size distribution (b) and clay

Fig. 6 -

 6 Fig. 6 -Distribution of the different types of microgranular structure (types A, B and C) in the

Fig. 7 -

 7 Fig. 7 -Model based on long-term termite activity for the Latosols studied showing the

  The pH w and pH KCl of the Bw horizons was greater in the Latosols on the SAS (5.20 ≤ pH w ≤ 5.34 and 5.73 ≤ pH KCl ≤ 6.20) than in those of the Latosols located on the VS (4.80 ≤ pH w ≤ 5.20 and 4.00 ≤ pH KCl ≤ 5.50). The organic carbon content ranged from 0.01 to 0.62 g kg -1 and the cation exchange capacity (CEC) from 1.73 to 10.60 cmol c kg -1 . On the SAS, the Bw horizons showed high levels of Al 2 O 3 content, (average

	Al	3 2 O

Table 3

 3 Physical and chemical characteristics of the diagnostic horizons (Bw) of the Latosols (L) studied

	(L) Hor. <2µm 2-50µm 50-2000µm Db	Dp pH w pH KCl OC Ca 2+ + Mg 2+	K +	Na +	Al 3+	H + + Al 3+	CEC	SiO 2 Al 2 O 3 Fe 2 O 3 TiO 2

  highly to moderate separated and well developed pedality and microstructure (60%) ► moderately separated and well developed pedality and granular microstructure (20%) ► weakly separated and weakly developed pedality and massive microstructure (20%) highly separated and well developed pedality and granular microstructure (60%) ► moderately separated and well developed pedality and granular microstructure (30%) ► weakly separated and weakly developed pedality and massive microstructure (10%)

	Latosol	Microstructure	Voids		Groundmass		Micromass
					Fabric elements
				Distribution	Orientation	c/f -Related	Color	b-fabric
				patterns	patterns	Distribution
	L1	► highly separated and well developed pedality and granular microstructure (50%)	► compound packing voids	random	unreferred double-spaced-	red to dark	undifferentiated
		► moderately separated and well developed pedality and granular microstructure (35%)	► complex packing voids			porphyric	red
		► weakly separated and weakly developed pedality and massive microstructure (15%)	► regular and star-shaped vughs and planes				2.5YR 4/8 to
			voids				10R 3/5
	L2	► highly separated and well developed granular microstructure (70%)	► compound packing voids	random	unreferred open-porphyric	red	undifferentiated
		► weakly separated and weakly developed pedality and massive microstructure (30%)	► regular and star-shaped vughs and planes				2.5YR 4/8
			voids
	L3	► highly separated and well developed pedality and granular microstructure (95%)	► compound packing voids and complex	random	unreferred open-porphyric	brownish	undifferentiated
		► weakly separated and weakly developed pedality and massive microstructure (5%)	packing voids				yellow
							10YR 6/6
	L4	► highly to moderate separated and well developed pedality and granular microstructure and	►simple to compound packing voids	random	unreferred	chitonic to	yellow	undifferentiated
		pellicular grain microstructure (100%)			(coarse-monic)	10YR 7/8
	L5	► ► compound packing voids	random	unreferred	open to	dark red	undifferentiated
			► complex packing voids		double-spaced-	2.5YR 3/6
			► channels and planes voids			porphyric
	L6	► highly separated and well developed pedality and granular microstructure (95%)	► simple to compound packing voids	random	unreferred open-porphyric	dark red	undifferentiated
		► weakly separated and weakly developed pedality and massive microstructure (5%)					10R 3/6
	L7	► highly separated and well developed pedality and granular microstructure (65%)	► compound packing voids	random	unreferred open-porphyric	red	undifferentiated
		► moderately separated and well developed pedality and granular microstructure (20%)	► complex packing voids				2.5YR4/8
		► weakly separated and weakly developed pedality and massive microstructure (15%)	► regular and star-shaped vughs and planes
			voids
	L8	► ► compound packing voids	random	unreferred open-porphyric dark red to	undifferentiated
			► complex packing voids				red 10R3/6
			► regular and star-shaped vughs and planes				to 4/8
			voids
	L9	► highly separated and well developed pedality and granular microstructure (55%)	► compound packing voids	random	unreferred open-porphyric	red	undifferentiated
		► moderately separated and well developed pedality and granular microstructure (30%)	► complex packing voids				2.5YR4/8
		► weakly separated and weakly developed pedality and massive microstructure (15%)	► regular and star-shaped vughs and planes
			voids
	L10	► highly separated and well developed pedality and granular microstructure (75%)	► compound packing voids	random	unreferred open-porphyric	red	undifferentiated
		► moderately separated and well developed pedality and granular microstructure (20%)	► complex packing voids				2.5YR4/8
		► weakly separated and weakly developed pedality and massive microstructure (5%)
	c/f: course/fine material in the groundmass

Table 6

 6 Volume of voids resulting from the assemblage of the microaggregates (inter-aggregates voids) for the Bw horizons of the Latosols (L) studied Volume of inter-microaggregates voids (

	inter v V )

Acknowledgements

We thank the Empresa Brasileira de Pesquisa Agropecuária (EMBRAPA) for supporting studies of A. Reatto in Orléans, France. This research is part of the project Embrapa Cerrados -IRD, n°0203205 (Mapping of the Biome Cerrado Landscape and Functioning of Representative Soils). We also like to acknowledge the support of C. Le Lay (INRA) for assistance in the laboratory.

Table 1 General characteristics of the Latosols (L) studied (modified after [START_REF] Reatto | Hydraulic properties of the diagnostic horizon of Latosols of a regional topossequence across the Brazilian Central Plateau[END_REF][START_REF] Reatto | Variation of the kaolinite and gibbsite content at regional and local scale in the Latosols of the Brazilian Central Plateau[END_REF] Red Latosol (1) Rhodic Acrustox (2) Rhodic Ferralsol (3) 1,050 Median 3 < 1

Latosols

Metamorphic Rocks-Complex Goiano -(Anápolis -Itauçu) and metaclastic rocks -Group Araxá -Superior Precambrian Granulite Mafic granulite L2 S 15°37'127'' W 47°45'576'' South America

Red Latosol (1) Typic Acrustox (2) Orthic Ferralsol (3) 1


B Iooss

V Vergès

V Larget

BEPU robustness analysis via perturbed law-based sensitivity indices

Keywords: Computer experiments, Density perturbation, Importance sampling, Nuclear safety, Sensitivity Analysis, Uncertainty

The "best-estimate plus uncertainty" (BEPU) methodology is the term used in the nuclear engineering community when dealing with uncertainty quantification issues in realistic numerical simulation models. One of the most critical hypothesis in these studies is the choice of the probability distributions of uncertain input variables which are propagated through the model. Bringing stringent justifications to the BEPU approach, especially in a safety study, requires quantifying the impact of potential uncertainty on the input variable distribution. To solve this problem, this paper deepens the robustness analysis based on the "Perturbed Law-based sensitivity Indices" (PLI). The PLI quantifies the impact of a perturbation of an input distribution on the quantity of interest (as a quantile of a model output or a safety margin) in the BEPU study. The mathematical formalism of the PLI is applied to two particular quantities of interest: the quantile and the superquantile. For both quantities, the PLI can be easily computed using a unique Monte-Carlo sample containing model inputs and output. Numerical tests are developed in order to define validity criteria of a PLI-based robustness analysis. The practical use of the method is illustrated on thermal-hydraulic computer experiments, simulating a cold leg Intermediate Break Loss Of Coolant Accident (IBLOCA) in a pressurized water nuclear reactor.

Introduction

Uncertainty quantification from numerical simulation models [START_REF]Uncertainty in industrial practice[END_REF][START_REF] Smith | Uncertainty quantification[END_REF][START_REF] Morio | Estimation of rare event probabilities in complex aerospace and other systems[END_REF] is gaining more and more attention in the engineering fields (e.g., nuclear and aerospace). It not only allows to assess the reliability and safety of complex systems (e.g., by estimating statistical risk measures as variance, quantile, probability of threshold exceedance, . . . ), but also to provide robust methodological tools which can complete the deterministic (and sometimes over-conservative) analyzes already performed, in order to better identify the safety margins. As an example, safety analyzes of nuclear power plants involving thermal-hydraulic issues, such as a cold leg Intermediate-Break Loss-Of-Coolant Accident (IBLOCA) in a pressurized water reactor, can be treated using probabilistic and statistical approaches [START_REF] Prosek | The state-of-the-art theory and applications of best-estimate plus uncertainty methods[END_REF][START_REF] Scarrott | Extreme-value-model-based risk assessment for nuclear reactors[END_REF]. In the nuclear engineering community, when the computer code, that is used to simulate the involved physical phenomena, lies on realistic assumptions, the uncertainty quantification framework is known as the "best-estimate plus uncertainty" (BEPU) methodology [START_REF] Wilson | Historical insights in the development of Best estimate Plus Uncertainty safety analysis[END_REF].

It is widely recognized that global sensitivity analysis plays a key role in the BEPU methodology [START_REF] De Crécy | Uncertainty and sensitivity analysis of the LOFT L2-5 test: Results of the BEMUSE programme[END_REF][START_REF] Baccou | SAPIUM: A generic framework for a practical and transparent quantification of thermal-hydraulic code model input uncertainty[END_REF] and, more widely, in systems modeling and policy support [START_REF] Saltelli | Global sensitivity analysis -The primer[END_REF][START_REF] Borgonovo | Sensitivity Analysis -An Introduction for the Management Scientist[END_REF][START_REF] Razavi | The future of sensitivity analysis: An essential discipline for systems modelling and policy making[END_REF]. When the quantity of interest is related to a rare event simulated by the numerical model, reliability-oriented sensitivity analysis methods can be put in action in order to quantify the sources of uncertainty (the model inputs) that drive the rare event [START_REF] Zio | Monte Carlo simulation-based sensitivity analysis of the model of a thermal-hydraulic passive system[END_REF][START_REF] Wei | Global reliability sensitivity analysis of motion mechanisms[END_REF][START_REF] Wei | Structural reliability and reliability sensitivity analysis of extremely rare failure events by combining sampling and surrogate model methods[END_REF][START_REF] Marrel | Statistical developments for target and conditional sensitivity analysis: application on safety studies for nuclear reactor[END_REF][START_REF] Zhang | A novel reliability sensitivity method based on directional sampling and Monte Carlo simulation[END_REF]. However, one of the most critical hypothesis in classical BEPU studies is the choice of the probability distributions of the uncertain inputs which are propagated through the numerical model [START_REF] Baccou | Development of good practice guidance for quantification of thermal-hydraulic code model input uncertainty[END_REF][START_REF] Chabridon | Global Reliability-oriented Sensitivity Analysis under Distribution Parameter Uncertainty[END_REF].

When experimental data are available, the methodology called inverse uncertainty quantification allows to obtain the probabilistic distributions of some physical model inputs [START_REF] Baccou | Development of good practice guidance for quantification of thermal-hydraulic code model input uncertainty[END_REF]. The inverse uncertainty quantification relies on a statistical inference process, involving a frequentist (via maximum likelihood inference) or a Bayesian approach [START_REF] Tarantola | Inverse problem theory and methods for model parameter estimation[END_REF][START_REF] Mendizábal | Bayesian perspective in bepu licensing analysis[END_REF][START_REF] Damblin | Bayesian inference and non-linear extensions of the CIRCE method for quantifying the uncertainty of closure relationships integrated into thermal-hydraulic system codes[END_REF]. Alternative approaches exist, as combining the use of design of experiments and forward uncertainty propagation (see [START_REF] Kovtonyuk | Development and assessment of a method for evaluating uncertainty of input parameters[END_REF][START_REF] Freixa | Testing methodologies for quantifying physical models uncertainties. A comparative exercise using CIRCE and IPREM (FFTBM)[END_REF] for examples in the nuclear field). For inputs without associated experimental data allowing to model their uncertainties, expert judgment is most often used. Whatever the method used for the choice of the probability density functions (pdf) for the model inputs, these pdf are established with various levels of confidence. In this paper, we focus on this particular issue associated to the residual uncertainty about each input pdf (whatever it was obtained). Indeed, facing to regulatory authorities, bringing stringent justifications to the BEPU approach requires quantifying the impact of this second-level uncertainty (i.e. the uncertainty on the input pdf) on the specific quantity of interest (QoI), most often a 2 high-order quantile, related to the BEPU study.

To address such issues, several strategies exist in the literature as a probabilistic modeling of uncertainty in the parameters of the inputs' distributions (see [START_REF] Pasanisi | Estimation of a quantity of interest in uncertainty analysis: Some help from Bayesian decision theory[END_REF] for the uncertainty propagation issues and [START_REF] Chabridon | Reliability-based sensitivity estimators of rare event probability in the presence of distribution parameter uncertainty[END_REF][START_REF] Meynaoui | New statistical methodology for second level global sensitivity analysis[END_REF] for global sensitivity analysis tools in this context) or getting out of the probabilistic framework by dealing with imprecise probabilities theories (see [START_REF] Beer | Imprecise probabilities in engineering analyses[END_REF] for a review and [START_REF] Song | Global sensitivity analysis for model with random inputs characterized by probability-box[END_REF] for global sensitivity analysis in the context of the p-box theory). Recently, the new branch of robustness analysis in uncertainty quantification has emerged in the field of sensitivity analysis of model outputs [START_REF] Iooss | Sensitivity analysis of model outputs: methods and issues for BEPU methodology[END_REF][START_REF] Perrin | Efficient estimation of reliability-oriented sensitivity indices[END_REF][START_REF] Hart | Robustness of the Sobol' indices to marginal distribution uncertainty[END_REF]. It consists of evaluating the impact of the choice of the inputs' distributions and, more precisely, by analyzing the QoI variations with respect to this choice. Two practical interests of this approach are to stay within the probabilistic framework and to avoid the difficult step of probabilistically modeling the input distributions parameters. These two aspects are clearly advantageous for the industrial issues related to our IBLOCA application of interest [START_REF] Larget | How to bring conservatism to a BEPU analysis[END_REF].

One particularly interesting solution for the robustness analysis has been proposed in the context of reliability-oriented sensitivity analysis by [START_REF] Lemaître | Density modification based reliability sensitivity analysis[END_REF] with the Perturbed Law-based sensitivity Indices (PLI). This measure only needs a Monte Carlo sample of model inputs and output, and presents the important advantage to be not limited by the number of model inputs that can be considered. Another advantage of the PLI is that it explicitly focuses on a particular QoI, that has to be specified by the user. The QoI of most BEPU studies is a high-order (e.g. 95%) quantile of the model output variable [START_REF] Prosek | The state-of-the-art theory and applications of best-estimate plus uncertainty methods[END_REF][START_REF] Mousseau | Uncertainty quantification in a regulatory environment[END_REF][START_REF] Larget | How to bring conservatism to a BEPU analysis[END_REF]. However, in simulation-based risk analysis, several works [START_REF] Labopin-Richard | Bregman superquantiles. Estimation methods and applications[END_REF][START_REF] Larget | Increasing conservatism in BEPU IB LOCA safety studies using complementary and industrially cost effective statistical tools[END_REF] have shown the interest of also considering the superquantile (i.e. the mean of the pdf tail exceeding the quantile) as a QoI. This risk measure gives interesting information about what is happening in the distribution-tail above the quantile. Its estimation also presents more regularities than those of the quantile estimation.

In this paper, the PLI-quantile concept [START_REF] Sueur | Perturbed-law based sensitivity indices for sensitivity analysis in structural reliability[END_REF][START_REF] Sueur | Sensitivity analysis using perturbedlaw based indices for quantiles and application to an industrial case[END_REF][START_REF] Delage | Robustness analysis of epistemic uncertainties propagation studies in LOCA assessment thermal-hydraulic model[END_REF] is extended to the PLIsuperquantile one. Via intensive simulation studies, our work also addresses a fundamental question that arises in practice: which level of perturbations can be applied on an input pdf so that the method remains valid? The following section develops the mathematical formalism of the PLI applied to the quantile and superquantile. It shows how, for both the QoI, the PLI can be easily computed using a unique Monte-Carlo sample (ensemble of code runs) containing model inputs and output. Section 3 describes the numerical tests that allow to define validity criteria of a PLI-based robustness analysis. Section 4 illustrates the PLI approach to assess the robustness of the IBLOCA thermal-hydraulic simulation scenario that uses the calculation code CATHARE2. A conclusion section ends this paper.

Principles of PLI-quantiles and PLI-superquantiles 2.1 Introduction to PLI

Given the random vector

X = (X 1 , . . . , X d ) ∈ X ⊆ R d of d independent uncertain input variables of pdf f (x) = f 1 (x 1 )×. . .×f d (x d ), G(•) a numerical model and Y = G(X) ∈ R the model
output, the QoI is a statistical quantity derived from Y . For instance, it can be the mean, the variance, a quantile or a threshold exceedance probability of the model output. In this paper, we limit the PLI analysis in the case of independent inputs and let to further studies the development of PLI in the dependent inputs case.

A PLI aims to measure the impact of the modification of an input density on a QoI [START_REF] Lemaître | Density modification based reliability sensitivity analysis[END_REF]. In order to compute the i-th PLI, we change the density f i of X i into a density f iδ , where δ ∈ R represents the level of the perturbation. The PLI is then simply defined as the relative change in the QoI generated by the perturbation:

S iδ = QoI iδ -QoI QoI . (1)
This definition slightly differs from the one proposed in previous studies [START_REF] Lemaître | Density modification based reliability sensitivity analysis[END_REF][START_REF] Sueur | Sensitivity analysis using perturbedlaw based indices for quantiles and application to an industrial case[END_REF]. Indeed, after several applications of the PLI, it has been found more convenient to compute directly the relative variation of the quantile when submitted to a density perturbation [START_REF] Larget | How to bring conservatism to a BEPU analysis[END_REF]. The interpretation of PLI is then straightforward.

Quantiles and superquantiles

The quantile of order α of Y is

q α = q α (Y ) = inf{t ∈ R, F Y (t) ≥ α} , (2)
where F Y is the cumulative distribution function of the random variable Y . In BEPU studies, a high-order (e.g. 95%) quantile value is most often preferred to a failure probability computation, as it can easily be related to safety margin concepts [START_REF] Prosek | The state-of-the-art theory and applications of best-estimate plus uncertainty methods[END_REF].

However, the quantile does not give any information about what is happening in the distribution-tail above the quantile. Another quantity, called superquantile, giving more information about the distribution-tail has therefore been introduced [START_REF] Rockafellar | Coherent approaches to risk in optimization under uncertainty[END_REF]. A superquantile of order α is defined by

Q α = Q α (Y ) = E G(X)1 G(X)≥q α 1 -α , (3)
where E(•) denotes the expectation operator and 1 A is the indicator function of the set A. If F Y is continuous, we have the equivalence with the following definition [START_REF] Rockafellar | Conditional value-at-risk for general loss distributions[END_REF]:

Q α = E(G(X)|G(X) ≥ q α ) . (4)
Also known under the name of expected-shortfall and conditional-value-atrisk in finance, the superquantile has been proven to be a coherent risk measure (which is not the case for the quantile) [START_REF] Rockafellar | Coherent approaches to risk in optimization under uncertainty[END_REF][START_REF] Labopin-Richard | Bregman superquantiles. Estimation methods and applications[END_REF]. Moreover, its empirical estimator often presents more regularities than those of the quantile empirical estimator.

Reverse importance sampling-based estimation formulas

In a lot of applications, for instance in nuclear safety studies (see, e.g., [START_REF] Iooss | Advanced methodology for uncertainty propagation in computer experiments with large number of inputs[END_REF]), the computer models are costly in terms of CPU time and memory. Only a limited number N of code runs is then available for the estimation of all the PLIs. The sample

X N = {x (n) = (x (n) 1 , ..., x (n)
d )} 1≤n≤N of N independent realizations of X (which follow the pdf f ) is propagated through the model G to produce the sample Y N = {y n } 1≤n≤N of N model outputs. The most standard estimation of a quantile is based on the empirical quantile estimator

denoted q α N = inf{t ∈ R, F N Y (t) ≥ α} where F N Y (t) = 1 N N n=1
1 yn≤t is the empirical estimator of the cumulative density function of Y . The Wilks formula or the bootstrap technique can then be used to add a second level of conservatism (due to the limited size of the sample [START_REF] Zio | Bootstrap and order statistics for quantifying thermal-hydraulic code uncertainties in the estimation of safety margins[END_REF]). In this work, the bootstrap technique will be used. From (2), the perturbed quantile writes:

q α iδ = inf{t ∈ R, F Y,iδ (t) ≥ α} , (5)
with F Y,iδ the cumulative distribution function corresponding to the input variable X i sampled from f iδ . In order to give an estimation of the perturbed quantile (denoted q α N,iδ ) from the same sample X N , we use the so-called reverse importance sampling mechanism [START_REF] Hesterberg | Estimates and confidence intervals for importance sampling sensitivity analysis[END_REF] to compute F N Y,iδ [START_REF] Delage | Robustness analysis of epistemic uncertainties propagation studies in LOCA assessment thermal-hydraulic model[END_REF]:

F N Y,iδ (t) = N n=1 L (n) i 1 (yn≤t) N n=1 L (n) i , (6)
with

L (n) i the likelihood ratio f iδ (x (n) i ) f i (x (n) i )
. We obtain

q α N,iδ = inf{t ∈ R, F N Y,iδ (t) ≥ α} . (7)
The estimator S N,iδ of the PLI-quantile consists in replacing QoI by q α N in Eq. ( 1). Its theoretical properties (asymptotic consistence and central limit theorem) have been given in [START_REF] Gauchy | An information geometry approach for robustness analysis in uncertainty quantification of computer codes[END_REF].

From (3), the perturbed superquantile writes:

Q α iδ = E iδ G(X)1 G(X)≥q α iδ 1 -α , (8)
where E iδ means that the expectation is taken under the perturbed distribution (f 1 , . . . , f iδ , . . . , f d ). We define an alternative version (used in the engineering study of [START_REF] Larget | Increasing conservatism in BEPU IB LOCA safety studies using complementary and industrially cost effective statistical tools[END_REF]) of the perturbed superquantile:

Q α iδ = E G(X)1 G(X)≥q α iδ 1 -α . (9)
This equation gives an alternative estimator because the expectation is taken under the initial distribution instead of the perturbed distribution. Its interest is that its estimation is easier and more stable (because the available sample is more representative of the initial distribution than the perturbed distribution). In the following, the alternative perturbed superquantile refers to Eq. ( 9), while the original perturbed superquantile refers to Eq. ( 8).

In order to estimate these quantities from the sample (X N , Y N ), we use the same mechanism (i.e. the reverse importance sampling procedure) as before by using the discretized-integral formula

Q α N,iδ = 1 N (1 -α) N n=1 G(x (n) )1 G(x (n) )≥ q α N,iδ L (n) i . ( 10
)
The estimator S N,iδ of the original PLI-superquantile consists in replacing QoI (resp. QoI iδ ) by Q α N (resp. Q α N,iδ ) in Eq. ( 1). For the alternative PLI-superquantile estimation, we simply write:

Q α N,iδ = 1 N (1 -α) N n=1 G(x (n) )1 G(x (n) )≥ q α N,iδ . (11)
In the following sections, the PLI-quantile and original PLI-superquantile are considered. Moreover, the estimator S N,iδ will be denoted S iδ .

Input density perturbation method

The density perturbation approach proposed in [START_REF] Lemaître | Density modification based reliability sensitivity analysis[END_REF] (and used in [START_REF] Sueur | Sensitivity analysis using perturbedlaw based indices for quantiles and application to an industrial case[END_REF][START_REF] Iooss | Uncertainty and sensitivity analysis of functional risk curves based on Gaussian processes[END_REF]) consists of replacing the density f i of one input X i by a perturbed one f iδ , where δ ∈ R represents a shift of a moment (e.g. the mean or the variance). Amongst all densities with shifted mean or variance of a δ value, f iδ is defined as the one minimizing the Kullback-Leibler divergence from f i . Figure 1 illustrates some perturbed pdf after a perturbation on the mean of initial pdf. Increasing the mean of a Gaussian pdf just consists in translating it by a constant shift (the perturbed pdf remains Gaussian). Increasing the mean of a uniform pdf turns to a non-uniform pdf (because we force its support to be fixed) which gives more weights to the large values.

In the following, we focus on one type of perturbations used in most of the studies, that is the perturbation of the mean of the random variable. However, even in this simplest perturbation case, minimizing the Kullback-Leibler divergence can be difficult or even impossible for some particular distributions (such as the log-normal one) [START_REF] Lemaître | Density modification based reliability sensitivity analysis[END_REF]. A simpler approach proposed in [START_REF] Perrin | Efficient estimation of reliability-oriented sensitivity indices[END_REF] starts by applying an iso-probabilistic operator (such as the Rosenblatt transform [START_REF] Rosenblatt | Remarks on a multivariate transformation[END_REF]) to the initial input pdf, in order to transform all the input random variables into centered normalized Gaussian ones (denoted N (0, 1)). In case of independent inputs, it just corresponds to the inverse probability distribution transform. Then, the perturbations are applied on Gaussian pdf (as in Figure 1, left) and the solution is trivial as the perturbed pdf is N (δ, 1) where δ is the perturbation level. This method allows to make perturbations comparable when applied in this standard space [START_REF] Delage | Robustness analysis of epistemic uncertainties propagation studies in LOCA assessment thermal-hydraulic model[END_REF][START_REF] Perrin | Efficient estimation of reliability-oriented sensitivity indices[END_REF][START_REF] Larget | How to bring conservatism to a BEPU analysis[END_REF].

In practice, one has to define the variation range of δ as it can be of interest to decrease and increase it. Looking at the PLI evolution as a function of δ will be also informative. To illustrate this, we compute the PLI-quantile at order α = 0.95 on the linear model output Y X i ∼ N (0, 1) for i = 1, 2, 3. The estimations are made by using a Monte Carlo sample (X N , Y N ) of size N = 5000, by applying a mean perturbation on the range δ ∈ [-1, 1] (with a regular discretization). Figure 2 shows the PLI curve for each input. In this linear model case, the PLI curves are linear because the output quantile increases linearly with the increase of each input mean. Moreover, the larger weight of X 1 in the model equation explains its larger impact than that of X 2 and X 3 .

= 2X 1 + X 2 + X 3 /2 with
Thanks to the use of a Monte Carlo sample, applying a bootstrap technique (here with 200 replicas) allows to associate a confidence interval (CI), typically at 95%, to each PLI estimate. Evaluating this estimation uncertainty is crucial in practice as the limited sample size can lead to erroneous result from a certain level of perturbation (e.g. when there are not enough sample values close to the perturbed quantile). For δ = 0, the PLI logically becomes null without any estimation error.

Numerical validation tests

In order to validate the PLI-based robustness analysis, we develop in this section several numerical tests on toy functions with a small number of inputs (three inputs X i ∼ N (0, 1) for i = 1, 2, 3). Note that the method can be applied to large-dimensional models and non-Gaussian inputs as it will be shown in Section 4. The two toy models are the linear model shown in Section 2.4:

Y = 2X 1 + X 2 + X 3 /2 , (12)
and the well-known (in global sensitivity analysis) Ishigami function:

Y = sin(X 1 ) + 7 sin 2 (X 2 ) + 0.1X 4 3 sin(X 1 ) . (13)
The PLI-quantile and (original) PLI-superquantile estimation (which are based on the same principles) are studied by applying a mean perturbation on the range δ ∈ [-3, 3]. In this section, the bootstrap-based 95% PLI-CI are computed with R = 200 replicas. The Figure 3 shows the PLI-quantiles (first and second lines) and PLIsuperquantiles (third and fourth lines) with α = 0.95, for the input X 2 of the linear model 12, computed by using a Monte Carlo sample (X N , Y N ) of size N = 1000 (first and third lines) and N = 10000 (second and fourth lines).

Results for X 1 and X 3 are not shown for sake of conciseness and as they are quite similar than those for X 2 (except than the PLI amplitudes are larger for X 1 and smaller for X 3 as preliminary shown in Fig. 2) due to their similar structural effects in the linear model [START_REF] Zio | Monte Carlo simulation-based sensitivity analysis of the model of a thermal-hydraulic passive system[END_REF].

The left column of Fig. 3 consists in comparing the PLI estimates in black (formulas of Section 2.3) with the reference results based on resimulations in red (i.e. results that are obtained without the reverse importance sampling procedure but by resampling the input following the perturbed pdf). Of course, in a real-case study, only the reverse sampling approach is applicable. Note that the CI of the PLI estimates are obtained by bootstrap (as in Fig. 2) while the CI of the reference results are obtained by 100 Monte Carlo repetitions of the resimulation process (so these CI are no subject to bootstrap-based approximation). We can observe that the PLI (and the PLI-CI) are valid in the central zone of the plots (centered on the value δ = 0) and turn to be invalid at large δ values (when red CI are not inside black CI). Moreover, increasing the sample size extends the PLI-CI validity domain with respect to the δ variation range. Indeed, as the sample size is finite, at a certain level of perturbation, there is not enough sample values to correctly compute the perturbed quantile and its CI.

From these findings, two heuristics are proposed to provide warnings to users when δ is too large:

1. The number of values in X N , smaller or larger of the δ-perturbed mean of an input variable, has to be sufficient. A value of N x = 30 has been chosen (from several numerical tests) as the smallest size. This criterion is easily computed from X N and each δ-value. It is represented by the green vertical lines on each plot. Of course, this value of 30 is quiet subjective and will be subject to further sensitivity studies.

2. The number of values in Y N , smaller or larger of the δ-perturbed quantile of the model output, has to be sufficient. A value of N y = 5 (resp. N y = 10) has been chosen (from several numerical tests) as the smallest size for the PLI-quantile (resp. PLI-superquantile). It is represented by the blue vertical line on each plot. Of course, these values of 5 and 10 are quiet subjective and will be subject to further sensitivity studies.

In Figure 3 (left column), these two criteria are visualized in each plot. It shows their relevance to warn the user that the PLI values (and their CI) no longer make sense beyond these limits (e.g. because PLI-CI do not contain the PLI reference values). Of course, the green and blue vertical lines are not visible in the plots when they are outside of the δ variation range. The right column of Fig. 3 visualizes the uncertainty on the PLI-CI by showing the 95%-CI of the PLI-CI lower bound (in orange) and the 95%-CI of the PLI-CI upper bound (in red). These 95%-CI are obtained by 50 Monte Carlo repetitions of the bootstrap-based estimation of the PLI-CI. It aims to visualize the stability of the PLI-CI, knowing that looking at the value of the PLI should not be done without looking at its associated CI. We verify that, beyond the limits given by the green and blue vertical lines, the CIs of the lower and upper bounds of the PLI-CI explode.

For the Ishigami function ( 13), the Figure 4 provides the results of the validity tests, computed by using a Monte Carlo sample (X N , Y N ) of size N = 1000 and N = 10000. The PLI-quantiles (first and second lines) and PLI-superquantiles (third and fourth lines), with α = 0.95, are given for X 1 (left column) and X 2 (right column). The interest of this toy case is to show that the PLI can also capture some strong non-linear behavior. Inside the bounds given by the green vertical lines (criteria on the sufficient number of points smaller or larger of the δ-perturbed mean of the input variable), we observe the validity of the PLI-CI (the PLI-CI encompasses the red curve corresponding to resimulation which gives reference PLI values, i.e. estimates without bias). The blue vertical line does not appear on these plots, meaning that the number of output variable points smaller or larger of the δ-perturbed quantile or superquantile is sufficient. We conclude that these results show similar behavior than the linear model ones and confirm the relevance of the two heuristic validity criteria.

Application on IBLOCA industrial use-case

As a nuclear facility operator, EDF (the French company of electricity) has to perform safety analysis on his pressurized water nuclear reactors. Those studies include the "Loss of Coolant Accident" resulting in a cold leg intermediatebreak on the primary loop of the reactor (noted IBLOCA). This hypothetical scenario is simulated using system thermal-hydraulic computer codes, which include tens of physical parameters such as condensation or heat transfer coefficients [START_REF] Mazgaj | Comparison of CATHARE results with the experimental results of cold leg intermediate break LOCA obtained during ROSA-2/LSTF test 7[END_REF][START_REF] Sanchez-Saez | Uncertainty analysis of large break loss of coolant accident in a pressurized water reactor using non-parametric methods[END_REF]. Yet, the values of the model parameters are known with a limited precision as they are evaluated through tests performed on smallscale experimental facilities. Some other variables are only observed during periodic inspections, such as the characteristics of hydraulics systems. The use-case of the present work is a 4-loop reactor from the French fleet1 . The transient studied is an IBLOCA, the break considered is a 13-inches break located on the cold leg. The reactor is supposed to be at nominal power at the opening of the break. The main common assumptions for IBLOCA safety analysis are considered. The availability of the systems is assumed to be restricted to the minimum safety systems (ECCS, EFWS, . . . ), other systems are considered to be unavailable. The most limiting additional single failure is postulated.

The IBLOCA transients are computed using the CATHARE2 code (developed by CEA, EDF, Framatome and IRSN [START_REF] Geffraye | CATHARE2 V2.5_2: A single version for various applications[END_REF]), which simulates the time evolution of physical quantities during a thermal hydraulic transient. In our simulation scenario, d ≈ 100 random inputs are considered and the output variable of interest is the second peak cladding temperature. For each input, depending on its nature and the information available on it, a pdf type (e.g. uniform, truncated Gaussian, truncated lognormal, triangular) is chosen, in concordance with its parameters. A Monte Carlo sample of N = 2000 input and output values have then be generated.

The QoI for this study is 95%-quantile and 75%-superquantile of the peak cladding temperature, which is estimated from the Monte Carlo sample at q 0.95 = 737 • C and Q 0.75 = 673 • C. Our PLI approach consists in transforming each input pdf in a N (0, 1) one, then perturbing the mean of each input with δ ∈ [-1.64, 1.64]. This range corresponds to the 5%-quantile and 95%quantile of a N (0, 1) variable. Figure 5 illustrates several perturbed empirical pdf from the Monte Carlo sample for four different variables.

On this high dimensional industrial application case (100 inputs), all the PLI can be easily computed. For a sake of conciseness, only 8 inputs with strong influence (described in Table 1) are shown in Fig. 6 through their PLI-quantile and (original) PLI-superquantile curves. The bootstrap-based 95% CI are computed with R = 50 replicas. Validity criteria do not appear on the plots because the size of the Monte Carlo sample is sufficient with respect to the perturbation values that have been considered.

By comparing the results of PLI-quantile and PLI-superquantile, we observe that both PLI give the same trends. However, differences can be noted on the most influential parameters, depending on the level of perturbation δ reached. Furthermore, PLI clearly determine the penalizing variation di- rection of each input: the QoI (0.95-quantile and 0.75-superquantile of the PCT) increase with an increase of the friction coefficients (X78, X110 and X111) and a decrease of the stable film temperature (X52), the exchange coefficients (X64 and X76) and the bubble-rise velocity X113. Decreasing X111 also leads to an increase of the QoI which allows to emphasize the non-monotonic behavior of the peak cladding temperature with respect to the interfacial friction in the core during blowdown phase. Such information is precious for discussing BEPU results.

In Figure 6, one important observation is the behavior of the confidence interval when δ is high. For some variables, for example X111 on the PLIquantile calculation, a saturation of the CI occur for high values of δ. This is due to the limited number of value of the output exceeding the quantile, which is particularly large for X111. That behavior is not observed in the PLI-superquantile calculation because of the inherent regularity of the su- perquantile, and the fact that the level of superquantile considered is low enough. This phenomenon seems to be fostered by the combination of PLIquantile at high quantile level and strong perturbation. Thus, for further industrial applications, in order to have the most accurate estimation of the PLI, the use of PLI-superquantile associated to limited value of δ (typically δ ≤ 1) should be considered, in addition to the application of the heuristic validity criteria. This recommendation has already been considered in the industrial methodology presented in [START_REF] Larget | Increasing conservatism in BEPU IB LOCA safety studies using complementary and industrially cost effective statistical tools[END_REF].

Conclusion

In the context of uncertainty quantification studies and in particular when applying BEPU methodologies, hypotheses have to be made in order to define the input distributions which are not always easy to justify. This work has focused on the PLI approach which allows to challenge the hypotheses made on the uncertain input pdf. The main goal of PLI is to quantify the robustness of a BEPU QoI with respect to uncertainty on the pdf of the model inputs. However, PLI is also a sensitivity measure quantifying the influence of each input on the QoI. In this paper, PLI has allowed to quantify the impact of each input mean on quantiles and superquantiles. Finally, an important contribution of this work relies in the introduction of new criteria for studying the validity of the PLI analysis. This approach has been applied to an IBLOCA scenario (nuclear reactor thermal-hydraulic phenomena) simulated with the CATHARE2 code. In addition to the rich information that is provided on the inputs' influence and on the QoI robustness, it has enlightened the two important advantages of the PLI approach: it only needs a simple Monte Carlo sample and it is not limited in terms of number of model inputs. For further industrial applications, in order to have the most accurate estimation of the PLI, the use of PLI-superquantile associated to limited value of δ (typically less or equal to 1) should be considered.

To develop more deeply the PLI, further works have to focus on more general perturbation type (see, e.g., a first proposition in [START_REF] Gauchy | An information geometry approach for robustness analysis in uncertainty quantification of computer codes[END_REF]), on developing multivariate PLI (perturbing several inputs at the same time), on considering the dependent inputs case (as in [START_REF] Benoumechiara | Detecting and modeling critical dependence structures between random inputs of computer models[END_REF] in another methodological framework) and on using more efficient sampling technique than Monte Carlo scheme (as the importance sampling as preliminary studied in [START_REF] Sueur | Perturbed-law based sensitivity indices for sensitivity analysis in structural reliability[END_REF] for PLI estimation).

Figure 1 :

 1 Figure 1: Examples of perturbed pdf for Gaussian (top) and uniform (bottom) initial pdf and a perturbation on the mean.

Figure 2 :

 2 Figure 2: PLI-quantile (α = 0.95) for the linear model. The points are the PLI at different discretized values of δ. The lines are the corresponding bootstrap 95%-confidence intervals.

Figure 3 :

 3 Figure 3: PLI validation results applied on X 2 of the linear model.

Figure 4 :

 4 Figure 4: PLI validation results applied on X 1 and X 2 of the Ishigami function.

Figure 5 :

 5 Figure 5: Examples of perturbed pdf for four inputs with different pdf of the IBLOCA case. Each case (a, b, c, d) shows the initial pdf (left) and the perturbed pdf by taking δ = 1.28 (right).

Figure 6 :

 6 Figure 6: PLI-0.95-quantile (upper panels) and PLI-0.75-superquantile (lower panels) of the IBLOCA case.

Table 1 :

 1 Definition of the 8 input variables considered for the PLI calculations. All the log-normal distributions are truncated.

	Input Variable description	Distribution
	X52	Minimum Stable Film Temperature	Uniform
	X64	HTC wall-steam exchange in C7 zone	Log-Normal
	X76	HTC primary-secondary -Natural liquid
		convection in SG	Log-Uniform
	X78	Wall-steam friction coefficient (CV) in HA Log-Uniform
	X110	Interfacial friction in the downcomer	Log-Normal
	X111	Interfacial friction in the core during
		blowdown phase	Log-Normal
	X113	Bubble-rise velocity in the upper plenum	Log-Normal
	X116	Bubble-rise velocity in the upper head	Log-Normal

The input data considered for this case, and thus the obtained peak cladding temperature, do not correspond to actual industrial values

Acknowledgments

We are grateful to Jérôme Stenger for his help on theoretical aspects of PLIsuperquantiles. We also thank Thibault Delage and Roman Sueur for useful remarks, as well as two anonymous reviewers. Numerical computations of the PLI have been realized using the sensitivity package [52] of the R software.