J M Baudoin
email: baudouin@univ-metz.fr

F Guérold

V Felten

E Chauvet

P Wagner

P Rousselle

Elevated Aluminium Concentration in Acidified Headwater Streams Lowers Aquatic Hyphomycete Diversity and Impairs Leaf-Litter Breakdown

Aquatic hyphomycetes play an essential role in the decomposition of allochthonous organic matter which is a fundamental process driving the functioning of forested headwater streams. We studied the effect of anthropogenic acidification on aquatic hyphomycetes associated with decaying leaves of Fagus sylvatica in six forested headwater streams (pH range, 4.3-7.1). Non-metric multidimensional scaling revealed marked differences in aquatic hyphomycete assemblages between acidified and reference streams. We found strong relationships between aquatic hyphomycete richness and mean Al concentration (r=-0.998, p<0.0001) and mean pH (r=0.962, p<0.002), meaning that fungal diversity was severely depleted in acidified streams. By contrast, mean fungal biomass was not related to acidity. Leaf breakdown rate was drastically reduced under acidic conditions raising the issue of whether the functioning of headwater ecosystems could be impaired by a loss of aquatic hyphomycete species.

Introduction

Running-water ecosystems that are subjected to anthropogenic acidification are typically first-to third-order streams, i.e. headwater streams often draining forested watersheds. Due to several ecological factors, such as light limitation by the riparian vegetation, the primary production in these streams is very low. Conversely, forested headwater streams receive by far most of their energy from the riparian vegetation. Therefore, the metabolism of headwater streams is generally heterotrophic [START_REF] Cummins | Shredders and riparian vegetation[END_REF][START_REF] Vannote | The river continuum concept[END_REF]. Leaf litter, which represents an important fraction of the allochthonous input entering the ecosystems, is used by decomposers and detritivores. Micro-organisms and leaf-shredding macroinvertebrates play an essential role in allochthonous organic matter processing and its incorporation into the trophic webs [START_REF] Gessner | A perspective on leaf litter breakdown in streams[END_REF][START_REF] Webster | Vascular plant breakdown in freshwater ecosystems[END_REF]. Aquatic fungi, particularly aquatic hyphomycetes which form a phylogenetically diverse group of anamorphous species (mainly ascomycetes and few basidiomycetes), contribute the most to the microbial biomass associated with leaves in streams [START_REF] Findlay | A cross-system comparison of bacterial and fungal biomass in detritus pools of headwater streams[END_REF][START_REF] Gulis | Effect of inorganic nutrients on relative importance of fungi and bacteria in carbon flow from submerged decomposing leaf litter[END_REF][START_REF] Weyers | Fungal and bacterial production during the breakdown of yellow poplar leaves in two streams[END_REF]. These filamentous fungi commonly occur in freshwaters and produce tetraradiate, variously branched or sigmoid asexual spores (i.e. conidia) adapted to dispersal in flowing water and adhesion to decaying plant tissues [START_REF] Bärlocher | Research on aquatic hyphomycetes: historical background and overview[END_REF]. A rapid conidial germination and mycelial invasion of the leaf matrix, typically within a few weeks, allow aquatic hyphomycetes to quickly colonize their ephemeral resource [START_REF] Gessner | Fungal decomposers of plant litter in aquatic ecosystems[END_REF]. In addition, the mycelial development coincides with or is closely followed by a very high release of conidia which constitutes a major part of the fungal production [START_REF] Gessner | Fungal decomposers of plant litter in aquatic ecosystems[END_REF]. Aquatic hyphomycetes are responsible for the initial breakdown of leaves and for enhancing the palatability of leaf detritus, thus providing a more suitable food source for shredders which further contribute to the conversion of leaf litter into fine particles and other decomposition products [START_REF] Arsuffi | Selective feeding by caddisfly (Trichoptera) detritivores on leaves with fungal colonized patches[END_REF][START_REF] Arsuffi | Selective feeding by shredders on leaf-colonizing stream fungi: comparison of macroinvertebrate taxa[END_REF][START_REF] Bärlocher | Freshwater fungal communities[END_REF][START_REF] Chamier | Effect of pH on microbial degradation of leaf litter in seven streams of the English Lake District[END_REF][START_REF] Covich | The role of benthic invertebrate species in freshwater ecosystems[END_REF][START_REF] Gessner | A perspective on leaf litter breakdown in streams[END_REF][START_REF] Golladay | Factors affecting food utilization by a leaf shredding aquatic insect: leaf species and conditioning time[END_REF][START_REF] Suberkropp | Microorganisms and organic matter processing[END_REF]. Consequently, due to their function, aquatic hyphomycetes occupy an essential place in the functional role that leaf breakdown plays in forested headwater streams.

Several studies have clearly shown that chemical alteration of aquatic ecosystems can have deleterious effects on fungal communities. A decline in aquatic hyphomycete diversity has been reported in streams contaminated with either heavy metals [START_REF] Bermingham | Effects of a coal mine effluent on aquatic hyphomycetes. I. Field Study[END_REF][START_REF] Duarte | Effects of zinc on leaf decomposition by fungi in streams: studies in microcosms[END_REF] or organic compounds [START_REF] Au | Fungi and cellulolytic activity associated with decomposition of Bauhinia purpurea leaf litter in a polluted and unpolluted Hong Kong waterway[END_REF][START_REF] Raviraja | Breakdown of ficus and eucalyptus leaves in an organically polluted river in India: fungal diversity and ecological functions[END_REF]. Similarly, leaf breakdown rates are often affected in rivers polluted with either mineral or organic pollutants (for review, see [START_REF] Gessner | A case for using litter breakdown to assess functional stream integrity[END_REF]). In particular, several studies have clearly demonstrated that leaf-litter breakdown is severely reduced in acidified ecosystems [START_REF] Dangles | Impact of stream acidification on litter breakdown: implication for assessing ecosystem functioning[END_REF][START_REF] Mulholland | Relathionships between stream acidity and bacteria, macroinvertebrates, and fish: a comparison of north temperate mountain streams[END_REF].

During the last three decades, a large number of studies have focused on the effects of acidification on aquatic biota (fish, macroinvertebrates, zooplankton, macrophytes, algae). Most of them have revealed deleterious effects leading to the loss of biodiversity. In particular, increased concentrations of aqueous aluminium arising from freshwater acidification are considered a major environmental problem due to its high toxicity to aquatic organisms. Aqueous aluminium has been recognized as a main toxicant for aquatic animals [START_REF] Gensemer | The bioavailability and toxicity of aluminum in aquatic environments[END_REF], and several studies have also suggested a direct effect on microbial metabolism [START_REF] Chamier | Effects of aluminium in acid streams on growth and sporulation of aquatic hyphomycetes[END_REF][START_REF] Myrold | Effect of acid rain on soil microbial processes[END_REF][START_REF] Pina | Microbial interactions with aluminium[END_REF]. Surprisingly and despite the key role they play in headwater stream functioning, little is known about the impact of anthropogenically acidified running waters on the diversity of microorganisms in general and on aquatic hyphomycetes in particular. Results from the few studies dealing with acidic water and the diversity of aquatic hyphomycetes are contradictory. In some studies, aquatic hyphomycetes have been shown to benefit from lower pH [START_REF] Bärlocher | Aquatic hyphomycete spora in 10 streams of New Brunswick and Nova Scotia[END_REF][START_REF] Wood-Eggenschwiller | Aquatic hyphomycetes in sixteen streams in France, Germany and Switzerland[END_REF], whereas in other studies, the richness was lower [START_REF] Iqbal | Aquatic hyphomycete spora of some Dartmoor streams[END_REF][START_REF] Shearer | Aquatic hyphomycete communities in the River Teign. I. Longitudinal distribution patterns[END_REF] or not related to pH [START_REF] Chamier | Effect of pH on microbial degradation of leaf litter in seven streams of the English Lake District[END_REF]. Obviously, additional data are needed to more fully understand the effects of low-pH waters on fungi colonizing leaves in streams [START_REF] Suberkropp | Fungal growth, production, and sporulation during leaf decomposition in two streams[END_REF].

To answer the question of whether acidified soft waters with high aluminium concentrations affect the diversity of aquatic hyphomycete assemblages associated with decaying leaves, a litter-breakdown experiment was conducted in six streams characterized by different acidification levels. After immersion in the streams, the leaf mass remaining was determined, and the sporulation of the leaf detritus was induced in the laboratory, allowing the identification of fungal species through released conidia. The relationship between diversity and breakdown rate was investigated to link the potential effect of species loss on leaf-litter processing in headwater streams.

Materials and Methods

Study Sites

The study was conducted in the Vosges Mountains (Northeastern France) where anthropogenic acidification has adversely affected surface waters (Fig. 1) [START_REF] Dambrine | Evidence of current soil acidification in spruce stands in the vosges mountains, north-eastern France[END_REF][START_REF] Probst | Evidence of springwater acidification in the Vosges Mountains (North-East of France): Influence of bedrock buffering capacity[END_REF]. Based on previous studies [START_REF] Dangles | Impact of stream acidification on litter breakdown: implication for assessing ecosystem functioning[END_REF][START_REF] Dangles | Influence of shredders in mediating breakdown rates of beech leaves in circumneutral and acidic forest streams[END_REF][START_REF] Felten | Short-term physiological responses to a severe acid stress in three macroinvertebrate species: a comparative study[END_REF][START_REF] Tixier | Plecoptera response to acidification in several headwater streams in the Vosges Mountains (northeastern France)[END_REF]

Leaf-Litter Breakdown

We used beech (F. sylvatica) leaves because this species is by far the most common deciduous riparian tree species and provides almost all leaf-litter input. Beech leaves were collected from trees just before abscission in fall 2002. Leaf bags were created by placing 5 g (±0.05) of air-dried leaves in 10 mm plastic mesh bags (15×20 cm). Each exact mass of leaf material was recorded and each bag identified with a plastic label. Weighed leaves were moistened with distilled water to avoid breakage during placement into mesh bags. Leaf bags were closed with nylon line so that they took the shape of a tetrahedron thus ensuring a good circulation flow and a natural arrangement of leaves. About 24 bags were submerged at each study site in zones expected to naturally accumulate leaves. They were individually secured to the bank with plastic-coated wire, which was anchored to the stream bottom with large boulders. Four replicate bags were randomly retrieved from the six streams after [START_REF] Farmer | The effects of lake acidification on aquatic macrophytes-a review[END_REF][START_REF] Griffith | Fungal biomass and leaf-litter processing in streams of different water chemistry[END_REF][START_REF] Legendre | Numerical Ecology[END_REF][START_REF] Thioulouse | ADE-4: a multivariate analysis and graphical display software[END_REF]83 and 124 days of exposure, immediately placed in individual sealed plastic bags with stream water and transported to the laboratory in a cool box. For each sampling date, 7 l of stream water was collected for laboratory processing of the leaf bags.

The leaves were rinsed individually with water from the respective stream to remove attached fine particles and invertebrates. Two leaf disks were cut avoiding the central vein with a 12 mm diameter cork borer from five randomly selected leaves in each bag. Five of them were stored at -80°C in a zip-lock plastic bag until ergosterol determination, and the five others were immediately processed to allow sporulation. The remaining leaf material was oven-dried to constant mass (105°C, 48 h) and weighed to the nearest 0.1 mg, and subsamples (500 mg) were ignited in a muffle furnace (550°C, 4 h) to relate dry mass to ash-free dry mass (AFDM). Ten additional leaf bags were kept in the laboratory before the beginning of the experiment. Four replicates were used to determine leaching loss by submerging them in a bowl under a slow flow of dechlorinated water for 24 h and six others to estimate the initial oven-dried mass and AFDM of all leaf bags.

Fungal Biomass

Ergosterol was extracted from leaf disks and determined as previously described [START_REF] Gessner | Use of solid-phase extraction to determine ergosterol concentrations in plant tissue colonized by fungi[END_REF]. Briefly, the leaf disks were lyophilised, weighed to the nearest 0.1 mg and heated in 5 ml hot alkaline ethanol (KOH, 8 g l -1 ) for 30 min at 80°C. The extract was purified by solid-phase extraction on cartridges (Waters Oasis HLB, 60 mg, 3 cc). Ergosterol was separated by reversed phase high performance liquid chromatography on C 18 and quantified by measuring absorbance at 282 nm. Ergosterol was converted to fungal biomass using a conversion factor of 5.5 mg ergosterol g -1 mycelial dry mass [START_REF] Gessner | Ergosterol-to-biomass conversion factors for aquatic hyphomycetes[END_REF].

Fungal Diversity

The sporulation rate and species composition of the leafassociated fungal assemblages were determined in the laboratory [START_REF] Gessner | Qualitative and quantitative analyses of aquatic hyphomycetes in streams[END_REF]. Five leaf discs from each replicate bag were placed in 100-ml Erlenmeyer flasks with 25 ml filtered water (Whatman glass fiber filter GF/F) from the stream in which the leaves were retrieved. Then, the flasks were placed on an orbital shaker (100 rpm) for 48 h at 10°C. After incubation, the suspension was poured into 50-ml centrifuge tube without the discs, which were rinsed in the flask with distilled water (3×2 ml) to dislodge remaining attached conidia. Rinse water was combined in the tube, and the volume was adjusted to 35 ml with distilled water and 2 ml of 37% formalin. Conidial suspensions were stored in the dark until analysis. The five leaf discs from each leaf pack were oven-dried (105°C, 48 h) and weighed to the nearest 0.1 mg, and AFDM was determined (550°C, 4 h).

The conidial suspension was transferred to a beaker with 1 ml Triton X-100 solution (0.5%) and shaken on a magnetic stirrer for 10 min, and an aliquot (4-8 ml) was filtered (membrane filter, 5-mm pore size, 25 mm diameter), stained with 0.1% Trypan blue in 60% lactic acid and examined microscopically. Conidia were identified by relying on specific literature [e.g. 39] and counted by scanning the entire filter. For each species, the sporulation rate (conidia mg -1 AFDM day -1 ) was determined.

Data Analysis

Principal components analysis (PCA) was carried out to ordinate streams with respect to physical and chemical variables. PCA was performed using the following variables: pH, ANC, conductivity, total Al, [NO 3 +SO 4 ] and [Ca+Mg].

Non-metric multidimensional scaling (NMDS) analysis of sporulation data was used to assess differences among sites in aquatic hyphomycete assemblages. This ordination method is a robust procedure for analysing ecological data [START_REF] Minchin | An evaluation of the relative robustness of techniques for ecological ordination[END_REF]. We used the Bray-Curtis coefficient to quantify the dissimilarity among sites based on joint occurrence and abundance of taxa [START_REF] Clarke | Change in marine communities: an approach to statistical analysis and interpretation[END_REF]. NMDS attempts to maximise the fit between measured dissimilarities and distance between resulting data points within a predefined number of spatial dimensions [START_REF] Legendre | Numerical Ecology[END_REF][START_REF] Legendre | Distance-based redundancy analysis: testing multispecies responses in multifactorial ecological experiments[END_REF]. The goodness of fit was estimated with a stress function, which ranges from 0 to 1, with values close to zero indicating a good fit. Axes from the NMDS analysis were correlated (Spearman rank correlation) with physical and chemical data to identify variables most strongly corresponding to among-site differences in aquatic hyphomycete assemblages [START_REF] Hawkins | Channel morphology, water temperature, and assemblage structure of stream insects[END_REF].

Leaf breakdown rates (k) were estimated by fitting the AFDM data with the linear model m t = m 0kt where m t is the AFDM remaining at time t, and m 0 is the initial AFDM at the beginning of the experiment. The linear model was used because it fitted the data better than the exponential model proposed by Petersen and Cummins [START_REF] Petersen | Leaf processing in a woodland stream[END_REF]. Daily decay rates were compared using analysis of covariance (ANCOVA) followed by a multiple-comparison (Tukey's test, [START_REF] Zar | Biostatistical analysis[END_REF]).

A one-way analysis of variance (ANOVA) was used to test for differences in fungal biomass, and Pearson correlations were performed to investigate possible relationships between diversity and physical and chemical parameters.

Statistical analyses were performed using ADE4 (PCA, [START_REF] Thioulouse | ADE-4: a multivariate analysis and graphical display software[END_REF]), SPSS for windows (NMDS) and Statistica (ANOVA, ANCOVA, Spearman rank and Pearson correlations).

Results

Water Chemistry

The six headwater streams ranged from near neutral [mean pH = 7.1 and 6.6 in La Maix (LM) and Bihet (BH), respectively] to strongly acidified (mean pH=4.6 and 4.3 in Gentil Sapin (GS) and Basse Escaliers (BE), respectively).

The first factorial plane of the PCA explained 91.4% of the total variance (Fig. 2). The first axis of the PCA was defined primarily by the ANC and Ca + Mg (Fig. 2a). It explained 68.6% of the total variance and strongly separated the circumneutral streams (LM and BH) from the acidified ones [GS, BE, Wassongoutte (WA), Longfoigneux (LG); Fig. 2b). An additional 22.8% of the total variance was explained by the second axis, which discriminated among streams draining catchments underlain by sandstone vs granite bedrock. Streams draining granite catchments (BH, WA, LG) were mainly characterized by lower concentrations of strong acid anions. In addition, acidified streams draining sandstone (BE, GS) exhibited higher total Al concentrations than acidified streams draining granite (WA, LG).

Fungal Diversity and Mycelial Biomass

A total of 37 species of aquatic hyphomycetes was observed during the leaf-litter breakdown experiment (Table 2). The higher fungal diversity was found in the two circumneutral streams BH and LM (27 and 25 species, respectively) and the lowest in BE (seven taxa) which was also the most acidic one. Approximately 17 species were never recorded in acidic streams. Among these species, eight were shared by the two circumneutral streams BH and LM: Alatospora pulchella, Culicidospora gravida, Lemonniera aquatica, Tricladium chaetocladium, Tricladium splendens, Tumularia aquatica, Tumularia tuberculata and Ypsilina graminea. The nine other species were observed only in one of the circumneutral streams, LM or BH. On the other hand, only one species, Casaresia sphagnorum, was recorded in the four acidic streams and never in the circumneutral ones.

During beech-leaf processing, Clavariopsis aquatica was the dominant species in the two circumneutral streams (i.e. 83.4% of the conidia produced in LM and 49.5% in BH), and Flagellospora curvula was dominant in the four acidic streams (i.e. from 83.1% in WA to 98.9% in BE).

The NMDS ordination described well the overall differences in aquatic hyphomycete assemblages between the six streams (Fig. 3). Stress was very low (0.003), indicating a high degree of representation. Axis 1 of the NMDS was strongly correlated (Spearman rank correlation) with pH (r= 0.89; p<0.05), ANC (r=0.89; p<0.05) and total aluminium concentration (r=-0.89; p<0.05). Position along this axis reflected the acidification status of the streams and strongly separated the acidic streams with low fungal richness (Table 2) from the circumneutral streams with markedly higher richness. Axis 2 correlated best with [SO 2À 4 þ NO À 3 ] concentrations (r=0.83; p<0.05) and separated streams according to the bedrock of their catchments.

Pearson correlation analyses between physical and chemical parameters and fungal species richness revealed significant relationships (Fig. 4) only with mean total aluminium concentration (log 10 transformed; r=-0.998, p<0.0001) and mean pH (r=0.962; p<0.002).

When considering fungal biomass, no clear tendency was observed, and mean values did not differ significantly between streams (one-way ANOVA; F=1.352, p=0.271), demonstrating that fungal biomass associated with decaying leaves was not affected by acidification (Fig. 5).

Breakdown Rates

Breakdown rates of beech leaves were significantly different among streams (ANCOVA; F=41.9, p<0.001) and ranged from 0.0045 day -1 in the most acidic site (BE) to 0.0285 day -1 in the circumneutral stream BH (Fig. 6; Table 3). The breakdown rate in the circumneutral stream draining granite (BH) was significantly faster than that observed in the circumneutral stream draining sandstone (LM). However, both streams exhibited faster breakdown rates than acidic streams, while no significant differences occurred among the four acidic streams. After 124 days of exposure in these acidic streams, the mass loss did not exceed 16%, indicating that leaf-litter breakdown in these headwater streams had virtually stopped.

Discussion

Acidification of ecosystems resulting from atmospheric pollution is one of the most revealing demonstrations that human activities can strongly impair terrestrial and aquatic ecosystems located in remote areas. Among the most striking consequences of freshwater acidification, the loss of biodiversity has been well documented for different taxonomic groups such as fish [START_REF] Heard | Episodic acidification and changes in fish diversity in Pennsylvania headwater streams[END_REF][START_REF] Leivestad | Fish kill at low pH in a Norwegian river[END_REF], invertebrates [START_REF] Driscoll | Acidic deposition in the Northeastern United States: Sources and inputs, ecosystem effects, and management strategies[END_REF][START_REF] Guérold | Macroinvertebrate community loss as a result of headwater stream acidification in the Vosges Mountains (N-E France)[END_REF], algae [START_REF] Meegan | Periphyton communities in headwater streams of different water chemistry in the Central Appalachian Mountains[END_REF][START_REF] Muller | Effects of artificial acidification on the growth of periphyton[END_REF] and macrophytes [START_REF] Farmer | The effects of lake acidification on aquatic macrophytes-a review[END_REF]. By comparison, the effect of anthropogenic acidification of aquatic ecosystems on the diversity of micro-organisms in general and aquatic hyphomycetes in particular has received much less attention despite their importance at multiple ecological scales. The few studies that have attempted to link aquatic hyphomycete richness and acidity have provided very contradictory results. For instance, in a study of aquatic hyphomycetes in ten streams of New Brunswick and Nova Scotia, Bärlocher [START_REF] Bärlocher | Aquatic hyphomycete spora in 10 streams of New Brunswick and Nova Scotia[END_REF] concluded that "their tolerance of low pH values makes them one of the rare groups of stream organisms that may actually benefit from the effects of acid rain". This conclusion contrasts markedly with the results of the previous studies of Iqbal and Webster [START_REF] Iqbal | Aquatic hyphomycete spora of some Dartmoor streams[END_REF] and Shearer and Webster [START_REF] Shearer | Aquatic hyphomycete communities in the River Teign. I. Longitudinal distribution patterns[END_REF] in Devon who reported impoverished fungal diversity in upland acidic water compared to lowland circumneutral sites, whereas Chamier [START_REF] Chamier | Effect of pH on microbial degradation of leaf litter in seven streams of the English Lake District[END_REF] found low species richness in upland streams regardless of pH (4.9-6.8) in seven streams of the English Lake District.

In our studies, we observed a severe depletion of species in acidified streams (from 44 to 75%), and more than 45% of the species associated with decaying leaves in circumneutral streams were never recorded in acidified streams. The species richness was strongly correlated with aluminium concentration and pH. However, because mineral acidity is the main factor regulating Al concentrations in acid surface waters, Al is more likely to control the diversity of aquatic hyphomycetes than H + . In experiments to test for the effect of Al on aquatic hyphomycetes, Chamier and Tipping [START_REF] Chamier | Effects of aluminium in acid streams on growth and sporulation of aquatic hyphomycetes[END_REF] found that deleterious effect of Al depend on the species considered, the effects being more marked on species common in circumneutral streams such as T. splendens. Similarly, differences in sensitivity to other metals among aquatic hyphomycete species have also been reported from experiments performed to assess the toxicity of Cd [START_REF] Abel | Effects of cadmium on aquatic hyphomycetes[END_REF].

Several hypotheses can be proposed to explain these discrepancies. Chemical parameters such as Al were not analysed in earlier studies, and the type of acidity (natural, i.e. organic vs anthropogenic) was not considered. Therefore, the potential toxicity of water can markedly differ even in streams with similar pH. Furthermore, the results provided by these previous studies of species richness in relation to pH are difficult to compare as several of them have been performed in upland and lowland streams and/or in forested streams and streams without riparian trees.

In our study, the acidity of headwater streams did not originate from organic acidity (DOC < 2.3 mg l -1 ), all the streams were forested (mainly A. alba, P. abies and F. sylvatica) and exhibited very similar environmental characteristics, allowing us to compare fungal diversity and species composition and to investigate the effects of anthropogenic acidification. Our results show that in the Vosges mountains, aquatic hyphomycete communities are severely impoverished, as are macroinvertebrate communities, as established in previous studies [START_REF] Guérold | Macroinvertebrate community loss as a result of headwater stream acidification in the Vosges Mountains (N-E France)[END_REF][START_REF] Tixier | Plecoptera response to acidification in several headwater streams in the Vosges Mountains (northeastern France)[END_REF]. In these streams, a combined toxic effect of high protons and aluminium concentrations has been demonstrated for different species of invertebrates [START_REF] Felten | Haemolymph [Na + ] and [Cl -] loss in Gammarus fossarum exposed in situ to a wide range of acidic streams[END_REF][START_REF] Felten | Short-term physiological responses to a severe acid stress in three macroinvertebrate species: a comparative study[END_REF] as commonly reported for other aquatic organisms such as fish [START_REF] Masson | Use of blood parameters in fish to assess acidic stress and chloride pollution in French running waters[END_REF], crayfish [START_REF] Appelberg | Changes in haemolymph ion concentrations of Astacus astacus L. and Pacifastacus leniusculus (Dana) after exposure to low pH and aluminum[END_REF], macrophytes [START_REF] Farmer | The effects of lake acidification on aquatic macrophytes-a review[END_REF] and algae [START_REF] Claesson | The toxicity of aluminium to two acido-tolerant green algae[END_REF][START_REF] Havens | The role of aluminium contamination in determining phytoplankton and zooplankton responses to acidification[END_REF].

It is interesting to note that fungal biomass associated with decaying leaves was not affected by acidification, thus corroborating the results from Dangles and Chauvet [START_REF] Dangles | Effects of stream acidification on fungal biomass in decaying beech leaves and leaf palatability[END_REF] in five streams from the Vosges Mountains but contrasting markedly with the findings of Griffith and Perry [START_REF] Griffith | Fungal biomass and leaf-litter processing in streams of different water chemistry[END_REF], who found significantly reduced fungal biomass associated with decomposing oak leaves in a strongly acidified West Virginia stream. Such a discrepancy illustrates the fact that leaf-mass loss is no straightforward indicator of the fungal activities. A possible explanation for the significant differences in leaf breakdown rates despite similar fungal biomasses in our study lies in the sensitivity of pectin lyases to water chemistry. The enzymatic degradation of pectin represents a key process in leaf-mass loss [START_REF] Jenkins | The influence of water chemistry on the enzymatic degradation of leaves in streams[END_REF], and different studies have shown that the activities of pectin lyases are higher in neutral or alkaline waters than in acid or soft waters [START_REF] Chamier | Water chemistry[END_REF][START_REF] Griffith | Leaf litter processing and exoenzyme production on leaves in streams of different pH[END_REF][START_REF] Jenkins | The influence of water chemistry on the enzymatic degradation of leaves in streams[END_REF]. Moreover, Chamier and Tipping [START_REF] Chamier | Effects of aluminium in acid streams on growth and sporulation of aquatic hyphomycetes[END_REF] have demonstrated that treatment with monomeric Al decreased pectinase production by four aquatic hyphomycete species.

Not surprisingly, our results showed that leaf-litter breakdown was severely depressed under acidic conditions, as previously reported in different studies [START_REF] Chamier | Effect of pH on microbial degradation of leaf litter in seven streams of the English Lake District[END_REF][START_REF] Dangles | Impact of stream acidification on litter breakdown: implication for assessing ecosystem functioning[END_REF][START_REF] Dangles | Influence of shredders in mediating breakdown rates of beech leaves in circumneutral and acidic forest streams[END_REF][START_REF] Griffith | Leaf litter processing and exoenzyme production on leaves in streams of different pH[END_REF][START_REF] Mulholland | Relathionships between stream acidity and bacteria, macroinvertebrates, and fish: a comparison of north temperate mountain streams[END_REF]. In a recent study of 25 woodland first-and secondorder streams along an acidification gradient in the Vosges Mountains (France), Dangles et al. [START_REF] Dangles | Impact of stream acidification on litter breakdown: implication for assessing ecosystem functioning[END_REF] found that breakdown rates of beech (F. sylvatica) leaves varied more than 20-fold between the most acidified and circumneutral sites (k=0.0002-0.0055 day -1 ) consistently with the associated fungal biomass and microbial respiration. Total abundance, biomass and species richness of leaf-shredding invertebrates associated with decomposing leaves were, however, not related to stream acidity (the diversity of aquatic hyphomycetes was unfortunately not investigated in this study) which suggests that fungi played a prominent role in acidified streams. That the slow breakdown rates we have measured in acidified headwater streams is linked to the decreased aquatic hyphomycete diversity is questionable. Several experimental studies have provided evidence that the loss of species can affect the efficiency with which resources are processed within an ecosystem [START_REF] Hector | Plant diversity and productivity experiments in European grasslands[END_REF][START_REF] Loreau | Biodiversity and ecosystems functioning: synthesis and perspectives[END_REF][START_REF] Tilman | Diversity and productivity in a long-term grassland experiment[END_REF] although the mechanisms underlying the response are not always well understood [START_REF] Zak | Plant diversity, soil microbial communities and ecosystem function: are there any links?[END_REF]. Whether diversity of aquatic hyphomycetes interferes with their functional role in leaf breakdown in the context of the present study should deserve special attention, e.g. through complementary in vitro and field experiments [e.g. 64]. Potential implications of such relationships are of major importance for higher trophic levels within the stream ecosystem.

Acknowledgment The present research was financially supported by the Office National des Forêts, the Conseil Général des Vosges and by the French ministry of ecology and sustainable development. Additional support was provided by the EU Commission (RivFunction: contact EVK1-CT-2001-00088).

Figure 2

 2 Figure 2 PCA performed on the physico-chemical variables in the six headwater streams during the leaf breakdown experiment. a Correlation circle on the F1xF2 factorial plane. b Ordination of the six sites on the F1×F2 factorial plane. Small squares represent the sample

Figure 3 Figure 4

 34 Figure 3 NMDS plot based on hyphomycete assemblages associated with decaying leaves in the six headwater streams. Circumneutral streams are shown in white and acidic streams in black. Squares represent streams draining granite bedrock and circles streams draining sandstone bedrock

1 )Figure 5 Figure 6

 156 Figure 5 Mean (+SD) fungi biomass in each stream. Streams are arranged following decreasing pH from the left to the right

Table 1

 1 Mean values (min-max) of the physico-chemical variables for the six streams over the study period (n=9)

	Al tot (μg l -1 )	55 (10-	148)	49 (10-	150)	221 (124-	376)	212 (92-	407)	379 (102-	540)	652 (316-	960)
	K +		1.60 (1.50-	1.67)	0.43 (0.29-	0.68)	0.11 (0.08-	0.15)	0.19 (0.10-	0.50)	1.08 (0.90-	1.14)	0.76 (0.67-	0.91)
	Na +		1.04 (0.98-	1.17)	3.73 (1.79-	6.64)	1.16 (0.98-	1.41)	1.32 (1.07-	1.84)	0.77 (0.65-	0.85)	0.61 (0.52-	0.73)
	Mg 2+		3.73 (3.13-	4.41)	1.37 (0.81-	2.25)	0.22 (0.19-	0.29)	0.29 (0.19-	0.94)	0.61 (0.43-	0.68)	0.26 (0.22-	0.27)
	Ca 2+	(mg l -1 )	6.98 (6.15-	8.24)	2.74 (1.73-	3.52)	0.56 (0.38-	1.00)	0.52 (0.41-	0.70)	1.35 (0.90-	1.60)	0.68 (0.61-	0.75)
	SO 2À 4		6.53 (5.55-	7.33)	2.78 (2.52-	3.15)	2.76 (2.58-	3.05)	2.69 (2.60-	2.82)	4.81 (4.23-	5.06)	5.15 (4.33-	5.51)
	NO À 3		4.04 (3.51-	4.35)	2.50 (1.11-	3.44)	0.89 (0.26-	1.46)	0.99 (0.35-	1.42)	5.28 (4.79-	5.83)	3.15 (2.30-	3.63)
	Cl -		1.36 (1.27-	1.48)	5.92 (1.37-	12.5)	0.78 (0.73-	0.84)	0.78 (0.66-	0.90)	1.09 (0.99-	1.19)	0.77 (0.69-	0.92)
	Conductivity ANC	(μS cm -1 ) (μeq l -1 )	81 (73-100) 379 (260-	590)	51 (29-75) 89	(32-194)	15 (13-18) -6	(-14-3)	15 (13-16) -7	(-13-0)	35 (31-39) -18	(-48-8)	33 (28-39) -33	(-58-21)
	Longitude pH	E	7°03′17″ 7.1 (6.4-	7.7)	6°52′32″ 6.6 (5.8-	7.6)	6°51′57″ 5.1 (4.2-	5.8)	6°53′13″ 4.9 (4.0-	5.4)	7°03′45″ 4.6 (3.9-	4.9)	7°05′45″ 4.3 (3.8-	4.5)
	Stream Acronym Latitude	N	La Maix LM 48°29′	24″	Bihet BH 47°57′	56″	Longfoigneux LG 47°57′	56″	Wassongoutte WA 47°58′	23″	Gentil Sapin GS 48°26′	58″	Basse BE 48°27′	Escaliers 59″

Table 2

 2 List of aquatic hyphomycete species associated with decaying beech leaves in each stream

	Species

Table 3

 3 Breakdown rates (day -1 ) determined by fitting a linear model

	Stream	-k	r 2
	LM	0.0143 a	0.979
	BH	0.0285 c	0.837
	LG	0.0045 b	0.873
	WA	0.0054 b	0.914
	GS	0.0054 b	0.884
	BE	0.0051 b	0.870
	All regressions were significant (p<0.005). Rates with different
	superscript letters are significantly different (Tukey's test multiple
	comparisons of slopes, p<0.005)


Paul Lemaître
email: lemaitre.paul.andre@gmail.com

Ekaterina Sergienko

Aurélie Arnaud

Nicolas Bousquet

Fabrice Gamboa

Bertrand Iooss

Arnaud

Density modification based reliability sensitivity analysis

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

Introduction

In the context of structural reliability, computer models are used in order to assess the safety of industrial systems relying on complex physical phenomena. For instance, an electric operator would like to predict the level of a potential river flood in order to determine the height of a dyke preventing any disaster. In this example, the computer model (simulating the hydraulic model) has some uncertain input variables (flow rate, river length, water height, etc.), that are modelled by random variables. In this paper, the computer code is a "black-box" deterministic numerical model and the study focuses on one of its output. Due to the randomness of the model inputs, this output is a random variable more or less sensitive to the uncertainty of the input variables.

Sensitivity analysis (SA) is a tool used to explore, understand and (partially) validate numerical models. It aims at explaining the outputs regarding the input uncertainties ( [START_REF] Saltelli | Sensitivity analysis for importance assessment[END_REF]). We use the "global SA" definition given by Saltelli et al. [START_REF] Saltelli | Sensitivity analysis in practice: A guide to assessing scientific models[END_REF] wherein the whole variation range of the inputs is considered. The application of such an approach can be model simplification (by removing irrelevant modelling elements), input variables ranking or research prioritization. There is a wide range of SA techniques, regarding what type of problem the experimenter faces with ( [START_REF] Iooss | Revue sur l'analyse de sensibilité globale de modèles numériques[END_REF]). For instance, screening methods are to be applied when there is a large number of inputs, and few models assumptions. From a quantitative point of view, the most popular techniques are variance-based methods and the so-called Sobol' indices ( [START_REF] Saltelli | Sensitivity analysis in practice: A guide to assessing scientific models[END_REF][START_REF] Sobol | Sensitivity analysis for non-linear mathematical models[END_REF]). These are based upon Hoeffding decomposition of L 2 function and functional variance decomposition [START_REF] Antoniadis | Analysis of variance on function spaces[END_REF].

It should be noticed that most SA methods focus on real-valued continuous numerical output variables. When the output is a binary value (e.g. when the numerical model returns "faulty system" or "safe system"),

• the practitioner wants to determine which are the inputs that impact the most the failure event the inputs distributions being set and supposed to be perfectly known. This amounts to an absolute ranking objective.

• P will be impacted by the choice of the input distributions; the engineer wants to assess the influence of this choice on the FP. Therefore the objective here is to quantify the sensitivity of the model output to the family or shape of the inputs.

• In practice, input distributions are estimated from data, thus leading to uncertainty on the values of the distribution parameters. The practitioner wants to assess the influence of the distribution parameters on P . Therefore the objective here is to test the sensitivity of the model to the parameters of the inputs.

In most studies, sensitivity indices for failure probabilities are defined in strong correspondence with a given method of estimation (e.g. [START_REF] Lemaire | Structural reliability[END_REF][START_REF] Zuniga | Adaptive directional stratification for controlled estimation of the probability of a rare event[END_REF]). Their interpretation is consequently limited. We propose in this article to define new generic sensitivity indices. Our sensitivity index is based upon density modification, and is adapted to failure probabilities. A methodology to estimate such indices is derived. For simplicity reasons, a classical Monte Carlo framework is considered in the following. Additionally, the sensitivity index can be computed using the sole set of simulations that has already been used to estimate the failure probability P , thus limiting the number of calls to the numerical model.

The outline of the article is the following: first we define a generic strategy of input perturbation in Section 2, based upon maximum entropy rules. We then present our index and its theoretical properties in Section 3, altogether with the estimation methodology. The behaviour of the indices is examined in Section 4 through numerical simulations in various complexity settings, involving toy examples and a realistic case-study. Comparisons with two reference sensitivity analysis methods (FORM indices and Sobol' indices) highlight the relevance of the new indices in most situations. The main advantages and remaining issues are finally discussed in the last section of the article, that introduces avenues for future research.

Methodologies of input perturbation

Our sensitivity analysis method requires to define a perturbation for each input. In general, and especially in preliminary reliability studies, there is no prior rule allowing to elicit a specialized perturbation for each input variable. We thus would like to propose a simple perturbation methodology, allowing the practitioner to answer the questions itemized in the Introduction. Furthermore, we make the implicit hypothesis that the extreme values of the inputs lead to the (rare) failure event.

Given a unidimensional input variable X i with pdf f i , let us call X iδ ∼ f iδ the corresponding perturbed random input. This perturbed input takes the place of the real random input X i , in a sense of modelling error : what if the correct input were X iδ instead of X i ?

More precisely, we suggest to define a perturbed input density f iδ as the closest distribution to the original f i in the entropic sense and under some constraints of perturbation. Information-theoretical arguments ( [START_REF] Cover | Elements of information theory 2nd edition (Wiley series in telecommunications and signal processing[END_REF]) led us to choose the Kullback-Leibler (KL) divergence between f iδ and f i as a measure of the discrepancy to minimize under those constraints. Given the hypotheses and the needs, we focus on linear constraints that can be interpreted in term of moments perturbation. This will lead to a quantification of the impact on P of each variable. This will also provide results on the sensitivity of P to the choice of input distributions. We will later present the perturbations corresponding to a mean shift and a variance shift.

Recall that between two pdf p and q we have

KL(p, q) = +∞ -∞ p(y) log p(y) q(y) dy if log p(y) q(y) ∈ L 1 (p(y)dy). (1)
Let i = 1, • • • , d, the constraints are expressed as follows in function of the modified density f mod :

g k (x i )f mod (x i )dx i = δ k,i (k = 1 • • • K) . (2)
Here, for k = 1, • • • , K, g k are given functions and δ k,i are given real. These quantities will lead to a perturbation of the original density. The modified density f iδ considered in our work is:

f iδ = argmin f mod |(2) holds KL(f mod , f i ) (3)
and the result takes an explicit form ( [START_REF] Csiszár | I-divergence geometry of probability distributions and minimization problems[END_REF]) given in the following proposition.

Proposition 2.1 Let us define, for λ = (λ 1 , • • • , λ K ) T ∈ R K , ψ i (λ) = log f i (x) exp K k=1 λ k g k (x) dx , (4)
where the last integral can be finite or infinite (in this last case ψ i (λ) = +∞). Further, set Dom

ψ i = {λ ∈ R K |ψ i (λ) < +∞}.
Assume that there exists at least one pdf f m satisfying [START_REF] Auder | Global sensitivity analysis based on entropy[END_REF] and that Dom ψ i is an open set. Then, there exists a unique λ * such that the solution of the minimisation problem ( 3) is

f iδ (x i ) = f i (x i ) exp K k=1 λ * k g k (x i ) -ψ i (λ * ) . (5)
The theoretical technique to compute λ is provided in appendix A. Here are presented two kinds of perturbations used further on.

Mean shifting

The first moment is often used to parameterize a distribution. Thus the first perturbation presented here is a mean shift, that is expressed with a single constraint:

x i f mod (x i )dx i = δ i . (6)
In term of SA, this perturbation should be used when the user wants to understand the sensitivity of the inputs to a mean shift -that is to say "what if the mean of input X i were δ i instead of E [X i ]?". Proposition 2.2 Considering the constraint [START_REF] Cover | Elements of information theory 2nd edition (Wiley series in telecommunications and signal processing[END_REF], under the assumptions of Proposition 2.1 the expression of the optimal perturbed density is

f iδi (x i ) = exp(λ * x i -ψ i (λ * ))f i (x i ) (7)
where λ * is such that equation ( 6) holds.

Note that Equation (4) becomes

ψ i (λ) = log f i (x i ) exp(λx i )dx i = log (M Xi (λ)) (8)
where M Xi (u) is the moment generating function (mgf) of the i-th input. With this notation, λ * is such that

x i exp (λ * x i -log (M Xi (λ * ))) f i (x i )dx i = δ i , which leads to x i exp (λ * x i ) f i (x i )dx = δ i M Xi (λ * ) .
This can be simplified to:

M ′ Xi (λ * ) M Xi (λ * ) = δ i . (9)
This equation may be easy to solve when the expression of the mgf of the input X i and of its derivative is known.

Variance shifting In some cases, the expectation of an input may not be the main source of uncertainty.

One might be interested in perturbing its second moment. This case may be treated considering a couple of constraints. The perturbation presented is a variance shift, therefore the set of constraints is:

x i f mod (x i )dx i = E [X i ] , x 2 i f mod (x i )dx i = V per,i + E [X i ] 2 . (10)
The perturbed distribution has the same expectation E [X i ] as the original one and a perturbed variance

V per,i = Var [X i ] ± δ i .
Proposition 2.3 Under the assumptions of Proposition 2.1, for the constraint [START_REF] Iooss | Revue sur l'analyse de sensibilité globale de modèles numériques[END_REF], the expression of the optimal perturbed density is:

f iδi (x i ) = exp(λ * 1 x + λ * 2 x 2 -ψ i (λ * ))f i (x i ) where λ *
1 and λ * 2 are so that equation ( 10) holds. As an example, the two kind of perturbations previously presented are provided for two families of inputs (Gaussian and Uniform) in figure 1. The perturbations are respectively a mean and variance increasing. It is noticeable (and will be proved further on) that the shape is conserved for the Gaussian distribution when shifting the mean or the variance. On the other hand, when increasing its mean, the Uniform distribution is packed down on the right-hand boundary of its support. When increasing its variance, the density is packed down on both boundaries of its support.

Perturbation of Natural Exponential Family

In general, when perturbating the input densities, the shape is not conserved. However in the specific case of Natural Exponential Family (NEF), the following proposition can be derived. Proposition 2.4 Assume that the original random variable X i belongs to the NEF, i.e. its pdf can be written as:

f i,θ (x i ) = b(x i ) exp [x i θ -η(θ)]
where θ is a parameter from a parametric space Θ, b(.) is a function that depends only of x i and

η(θ) = log b(x) exp [x i θ] dx i
is the cumulant distribution function. Considering the assumptions of Proposition 2.1, the optimal pdfs proposed respectively in Proposition 2.2 and Proposition 2.3 are also distributed according to a NEF.

The proof comes from theorem 3.1 in [START_REF] Csiszár | I-divergence geometry of probability distributions and minimization problems[END_REF]. The details of computation are given for a mean shift and a variance shift in Appendix B. 3 Definition, estimation and properties of a sensitivity index Given a unidimensional input variable X i with pdf f i and the corresponding perturbed random input X iδ ∼ f iδ . The perturbed failure probability becomes:

P iδ = 1 {G(x)<0} f iδ (x i ) f i (x i ) f (x)dx (11)
where x i is the i th component of the vector x. Independently of the mechanism chosen for the perturbation (see previous section for proposals), a good sensitivity index S iδ should have intuitive features that make it appealing to reliability engineers and decision-makers. We believe that the following definition can fulfil these requirements.

Definition 3.1 Define the Density Modification Based Reliability Sensitivity Indices (DMBRSI) the quantity S iδ :

S iδ = P iδ P -1 1 {P iδ ≥P } + 1 - P P iδ 1 {P iδ <P } = P iδ -P P • 1 {P iδ ≥P } + P iδ • 1 {P iδ <P } .
Firstly, S iδ = 0 if P iδ = P , as expected if X i is a non-influential variable or if δ expresses a negligible perturbation.

Secondly, the sign of S iδ indicates how the perturbation impacts the failure probability qualitatively. It highlights the situations when P iδ > P i.e. if the remaining (epistemic) uncertainty on the modelling X i ∼ f i can increase the failure risk. In this case, the uncertainty on the concerned variable should be more accurately analysed. Conversely, if if P iδ < P , P can be interpreted as a conservative assessment of the failure probability, with respect to variations of X i . In such a case, deeper modelling studies on X i appear less essential.

Thirdly, given its sign, the absolute value of S iδ has simple interpretation and provides a level of the conservatism or non-conservatism induced by the perturbation. A value of α > 0 for the index means that

P iδ = (1 + α)P . If S iδ = -α < 0 then P iδ = (1/(1 + |α|))P .
The postulated ability of S iδ to enlighten the sensitivity of P to input perturbations must be tested in concrete cases, when an estimator PN of P can be computed using an already available design of N numerical experiments. In this paper, N is assumed to be large enough such that statistical estimation stands within the framework of asymptotic theory. Besides, we assume for simplicity a standard Monte Carlo design of experiments, according to which PN = N n=1 1 {G(x n )<0} /N where the x 1 , • • • , x N are independent realisations of X. The strong Law of Large Numbers (LLN) and the Central Limit Theorem (CLT) ensure that for almost all realisations PN ----→

N →∞ P and N/[P (1 -P )]( PN -P ) L ----→ N →∞ N (0, 1).
The Monte Carlo framework allows P iδ to be consistently estimated without new calls to G, through a "reverse" importance sampling mechanism:

PiδN = 1 N N n=1 1 {G(x n )<0} f iδ (x n i ) f i (x n i )
.

This property holds in the more general case when P is originally estimated by importance sampling rather than simple Monte Carlo, which is more appealing when G is time-consuming [START_REF] Beckman | Monte-Carlo estimation under different distributions using the same simulation[END_REF][START_REF] Hesterberg | Estimates and confidence intervals for importance sampling sensitivity analysis[END_REF]. This generalization is discussed further in the text (Section 5). The following lemma ensures the asymtotic behaviour of such an estimator.

Lemma 3.1 Assume the usual conditions [START_REF] Saltelli | Variance based sensitivity analysis of model output. design and estimator for the total sensitivity index[END_REF]. It can be consistently estimated by

(i) Supp(f iδ ) ⊆ Supp(f i ), (ii) Supp(fi) f 2 iδ (x) f i (x) dx < ∞, then PiδN ----→ N →∞ P iδ and √ N σ -1 iδN PiδN -P iδ L ----→ N →∞ N (0, 1). The exact expression of σ -1 iδN is given in Appendix C, equation
σ2 iδN = 1 N N n=1 1 {G(x n )<0} f iδ (x n i ) f i (x n i ) 2 -P 2 iδN .
The proof of this Lemma is given in Appendix C.1.

The asymptotic properties of any estimator of S iδ will depend on the correlation between PN and PiδN . The next proposition summarizes the features of the joint asymptotic distribution of both estimators.

Proposition 3.1 Under assumptions (i) and (ii) of Lemma 3.1,

√ N PN PiδN - P P iδ L ----→ N →∞ N 2 (0, Σ iδ )
where Σ iδ is given in Appendix C, Equation ( 16) and can be consistently estimated by

Σiδ = PN (1 -PN ) PiδN (1 -PN ) PiδN (1 -PN ) σ2 iδN .
The proof of this Proposition is given in Appendix C.2.

Given ( PN , PiδN ), the plugging estimator for S iδ is

ŜiδN = PiδN PN -1 1 { PiδN ≥ PN } + 1 - PN PiδN 1 { PiδN < PN } . (12)
In corollary of Proposition 3.1, applying the continuous-mapping theorem to the function s(x, y) = y x -1 1 y≥x + 1 -x y 1 y<x , ŜiδN converges almost surely to S iδ . The following CLT results from Theorem 3.1 in [START_REF] Van Der | Asymptotic statistics[END_REF].

Proposition 3.2 Assume that assumptions (i) and (ii) of Lemma 3.1 hold and further that

P = P iδ , we have √ N ŜiδN -S iδ L ----→ N →∞ N 0, d T Σd ( 13
)
with d = ∂s ∂x (P, P iδ ), ∂s ∂y (P, P iδ )

T for x = y, and

∂s ∂x (x, y) = -y1 {y≥x} /x 2 - 1 y 1 {y<x}, ∂s ∂y (x, y) = 1 x 1 {y≥x} + x1 {y<x} /y 2 .
This holds when P = P iδ . Indeed, one has for x * = 0 :

lim y ≥ x x y → x * y * ∇s(x, y) = lim y < x x y → x * y * ∇s(x, y) = - 1 x * , 1 x * T .

Numerical experiments

In this section, the methodology presented throughout the article is tested on two academic cases and a more realistic industrial model. The new indices are compared to the results provided by two reference methods, FORM indices (or importance factors) and Sobol' indices. First, we briefly present these two methods.

4.1 Presentation of the two reference SA methods

FORM importance factors

The First Order Reliability Method (FORM) [START_REF] Lemaire | Structural reliability[END_REF] is an estimation technique for a failure probability based upon approximating the failure domain and solving an optimization problem. In practice, it is considered a standard solution in structural reliability, mainly because of its low cost. This method proceeds in four steps:

• transformation of the input variables space into the standard space (a space where all the variables are standard independent Gaussian);

• search of the closest failure point to the origin of the standard space (also referred as the design point) via an optimisation algorithm;

• approximation of the failure surface by an hyperplane;

• failure probability estimation based on the geometry of the failure domain.

Importance factors are byproducts of this method. These sensitivity measures aims at quantifying the importance of a variable on the failure probability. From the design point u * one writes:

u * = β HL α *
where β HL is the distance between the origin of the standard space and u * and α * is the normalised vector of direction. Then for each variable U i , one can obtain the importance factor α * 2 i .

Sobol' sensitivity indices

In the SA framework, let us have X = (X 1, ..., X d ), a random vector where the variables are mutually independent and Y = G(X), the output of a deterministic code G(). Thus a functional decomposition of the variance is feasible, often referred as functional ANOVA:

Var[Y ] = d i=1 V i (Y ) + d i<j V ij (Y ) + • • • + V 12...d (Y ) where V i (Y ) = Var[E(Y |X i )], V ij (Y ) = Var[E(Y |X i , X j )] -V i (Y ) -V j (Y
) and so on for higher order interactions. The so-called Sobol' indices or sensitivity indices [START_REF] Sobol | Sensitivity analysis for non-linear mathematical models[END_REF] are obtained as follows:

S i = V i (Y ) Var[Y ] , S ij = V ij (Y ) Var[Y ] , • • •
These indices express the share of variance of Y that is due to a given input or a set of inputs. The number of indices growths in an exponential way with the number d of dimension: there are 2 d -1 indices. For computational time and interpretation reasons, the practitioner rarely estimate indices of order higher than 2. Therefore Homma & Saltelli [START_REF] Homma | Importance measures in global sensitivity analysis of nonlinear models[END_REF] introduced the so-called total Sobol index or total effect that writes as follow:

S Ti = S i + j =i S ij + j =i,k =i,j<k S ijk + ...
This total sensitivity index contains all the effects due to X i and the interactions between X i and other inputs.

Testing methodology

Importance Factors and Sobol' indices (SI) are computed using the methodologies given in [START_REF] Lemaire | Structural reliability[END_REF] and [START_REF] Sobol | Sensitivity analysis for non-linear mathematical models[END_REF], respectively. The Sobol' indices are computed using two initial samples of size 10

Hyperplane failure surface

For the first example, X is set to be a 4-dimensional vector, with d = 4 independent marginal distributions normally distributed with parameters 0 and 1. Therefore f Xi ∼ N (0, 1) for i = 1, .., 4. The failure function is defined as:

G(X) = k - 4 i=1 a i X i
where k and a = (a 1 , a 2 , a 3 , a 4 ) are the parameters of the model. For this numerical example, parameters are set with values k = 16 and a = (1, -6, 4, 0). An explicit expression for P can be given:

P = φ   -k/ 4 i=1 a 2 i   = φ -16 √ 53 ≃ 0.014
where φ(.) is the standard normal cumulative distribution function.

It is expected that the influence of X i on P uniquely depends on |a i |. The larger the absolute value of the coefficient, the larger the expected influence. The aim for choosing one non-influential variable X 4 (because a 4 = 0) is to assess if the SA methods can identify this variable as non-influential on the failure probability.

FORM

In this ideal hyperplane failure surface case, FORM provides an approximated value PF ORM = 0.01398, which is as expected ( [START_REF] Lemaire | Structural reliability[END_REF]) close to the exact value. 39 model calls have been required. The importance factors, given in Table 1, provide an accurate variable ranking for the failure function.

Variable X 1 X 2 X 3 X 4
Importance factor 0.018 0.679 0.302 0

Sobol' indices

The first-order and total indices are displayed in Table 4.2.2. The interpretation of the results is that X 2 and X 3 concentrate most of the variance of the indicator function. At first order, 25% of its variance is explained by X 2 without any interaction. It should be noted that the total index for X 4 is null, assessing that this variable does not impact the failure probability. The C.o.v. of the total indices estimators are small, meaning that this method is reproducible and that 6 × 10 5 points are enough to estimate in an efficient way the indices S T i . On the other hand, some C.o.v. values for small first order indices are quite high. This effect is well-known in SI estimation problems and can be corrected by improved formulas [START_REF] Saltelli | Variance based sensitivity analysis of model output. design and estimator for the total sensitivity index[END_REF].

Sobol' Index S 1 S 2 S 3 S 4 S T 1 S T 2 S T 3 S T

Density modification based reliability indices

The method presented throughout this article is applied on the hyperplane function. As explained in section 2, several ways to perturb the input distributions exist. For this case, a mean shifting is first applied, then a variance shifting with fixed mean. A simple calculus gives that the perturbed pdf are Gaussian, respectively with the constraint mean and variance 1 for the mean shifting perturbation (see Table 7), and with mean 0 and the constraint variance for the variance shifting perturbation.

The MC estimation gives P = 0.01446 with 10 5 function calls. For the mean shifting (see ( 6)), the domain variation for δ ranges from -1 to 1 with 40 points, reminding that δ = 0 cannot be considered as a perturbation. For the variance shifting (see [START_REF] Iooss | Revue sur l'analyse de sensibilité globale de modèles numériques[END_REF]), the variation domain for V per ranges from 1/20 to 3 with 28 points, where V per = 1 is not a perturbation. The estimated indices are plotted respectively in Figure 2 for mean shifting and in Figure 3 for variance shifting. The 95% confidence intervals are plotted around the indices, using the presented asymptotic formulas in Section 3.

Mean perturbation indices

The indices S iδ behave in a monotonic way given the importance of the perturbation. The slope at the origin is directly related to the value of a i . For influential variables (X 2 and X 3 ), the increasing or the decreasing is faster than linear, whereas the curve seems linear for the slightly influential variable (X 1 ). Modifying the mean with a positive amplitude slightly rises the failure probability for X 1 , highly decreases it for X 2 and increases it for X 3 (Figure 2). The effects are reversed with similar amplitude for negative δ. It can be seen that X 4 has no impact on the failure probability for any perturbation. Those results are consistent with the expression of the failure function. One can see that the confidence intervals (CI) associated to X 2 and X 3 are fairly well separated, except for the small absolute value of δ. On the other hand, the CI associated to X 1 and X 4 are not separated until absolute value of δ higher than 0.6. It can be concluded from the observation of the CI is that the impact of variable X 1 and X 4 cannot be differentiated, unless a broad change of the mean occurs.

Variance perturbation indices Increasing the variance of input X 2 and X 3 increases the failure probability, whereas it decreases when decreasing the variance (Figure 3). Modifying the variance of X 1 and X 4 have no effect on the failure probability. The increasing of the indices is linear for X 2 and X 3 , and the decreasing of the indices is faster than linear, especially for X 2 . Considering the CI, one can see that they are well separated for variable X 2 and X 3 , assessing the relative importance of these variables. On the other hand, the CI associated to X 1 and X 4 are not separated and contain 0.

Conclusion

The DMBRSI has brought the following conclusions: when shifting the mean (that is to say the central tendency in this case), the most influential variable is X 2 , followed by X 3 . X 1 is slightly influential while X 4 is not influential at all. When shifting the variance, variable X 2 is more influential than variable X 3 . Variables X 1 and X 4 have no impact when shifting the variance. We argue that all these information are much richer than the ones provided by importance factors and by Sobol' indices. Indeed, the information are provided about regions of the input space leading to failure event. This is, in our opinion, more of interest to the practitioner than a "simple" variable ranking.

Thresholded Ishigami function

The Ishigami function is a common test case in SA since it presents high degree of non-linearity with interactions between the variables. A modified version of the Ishigami function will be considered in this paper. One has G(X) = sin (X 1 ) + 7 sin (X 2 ) 2 + 0.1X 4 3 sin (X 1 ) + 7 where X is a 3-dimensional vector of independent marginals uniformly distributed on [-π, π] . In Figure 4, the failure points (where G(x) < 0) are plotted in a 3-d scatterplot.

There are 614 failure points on a MC sample of 10 5 points. Therefore the failure probability here is roughly P = 6.14 × 10 -3 . The complex repartition of the failure points can be noticed. Those points lay in a zone defined by the negative values of X 1 , the extremal and mean values of X 2 (around -π, 0 and π), and the extremal values of X 3 (around -π and π).

FORM

The algorithm FORM converges to an incoherent design point (6.03, 0.1, 0) in 50 function calls, giving an approximate probability of PF ORM = 0.54. The importance factors are displayed in Table 3. The bad performance of FORM is expected given that the failure domain consists in six separate domains and that the function is highly oscillant, leading to optimization difficulties. The design point is aberrant, therefore the FORM results of SA are incorrect.

Variable X 1 X 2 X 3
Importance factor 1e -17 1 0

Table 3: Importance factors for Ishigami function

Sobol' indices

The first-order and total indices, computed with 5 × 10 5 function calls are displayed in Table 6. The small values of first order indices show that no variable has impact on the variance of the indicator of failure on its own. The three total indices have relatively high and similar values. This states that all the variables highly interact with each other to cause system failure. The SI method is thus non-discriminant in this case. The low C.o.v. show that the method is reproducible.

Density modification based reliability indices

The method presented throughout this article is applied on the thresholded Ishigami function. As for the hyperplane test case, a mean shifting and a variance shifting are applied. The modified distribution when a mean shift is applied on a uniform distribution is given in Table 7. The modified pdf when shifting the variance and keeping the same expectation is proportional to a truncated Gaussian when decreasing the variance. When increasing the variance, the perturbed distribution is a symmetrical distribution with 2 modes close to the endpoints of the support (see figure 1). As previously, the same MC sample of size 10 5 (also used to produce Figure 4) is used to estimate the indices with both perturbations. For the mean shifting (see ( 6)), the variation domain for δ ranges from -3 to 3 with 60 points -numerical consideration forbidding to choose a shifted mean closer to the endpoints. For variance shifting, the variation domain for V f ranges from 1 to 5 with 40 points. Let us recall that the original variance is Var[X i ] = π 2 /3 ≃ 3.29. The estimated indices are plotted respectively in Figure 5 for mean shifting and in Figure 6 for variance shifting.

Mean perturbation indices A perturbation of the mean for X 2 and X 3 will increase the failure probability, though the impact for the same mean perturbation is stronger for X 3 ( S 3,-3 and S 3,3 approximately equal respectively 9.5 and 10, Figure 5). On the other hand, the indices concerning X 1 show that a mean shift between -1 and -2 increases the failure probability, whereas an increasing of the mean or a large decreasing strongly diminishes the failure probability ( S 1,3 approximatively equals -7.10 11 ). Therefore, Figure 5 leads to two conclusions. First, the failure probability can be strongly reduced when increasing the mean of the first variable X 1 (this is also provided by Figure 4 wherein all failure points have a negative value of X 1 ). Second, any change in the mean for X 2 or X 3 will lead to an increase of the failure probability. The CI are well separated, except in the -1 to 1 zone. One can notice that the CI associated to X 2 contains 0 between values of δ from -1.5 to 1.5, thus the associated indices might be null in these case. This has to be taken into account when assessing the relative importance of X 2 .

Variance perturbation indices Figure 6 (upper part) shows that a change in the variance has little effect on X 2 and X 1 , though the change is of opposite effect on the failure probability. However, considering that the indices S 2,Vper,i and S 1,Vper,i lie between -0.4 and 0.4, one can conclude that the variance of theses variables are not of great influence on the failure probability. On the other hand, Figure 6 (lower part) shows that any reduction of Var [X 3 ] strongly decreases the failure probability, and that an increase of the variance

Industrial case : flood case

The goal of this test case is to assess the risk of a flood over a dyke for the safety of industrial installations [START_REF] Bernardara | Uncertainty analysis in flood hazard assessment: hydrological and hydraulic calibration[END_REF]. This comes down to model the level of a flood. As a function of hydraulical parameters, many of them being randomized to account for uncertainty. From a simplification of the Saint-Venant equation, a flood risk model is obtained. The quantity of interest is the difference between the level of the dyke and the height of water. If this quantity is negative, the installation is flooded. Hydraulical parameters are the following: Q the flow rate, L the watercourse section length studied, B the watercourse width, K s the watercourse bed friction coefficient (also called Strickler coefficient), Z m and Z v respectively the upstream and downstream bottom watercourse level above sea level and H d the dyke height measured from the bottom of the watercourse bed.

The water level model is expressed as:

H =   Q K s B Zm-Zv L   3 5
.

Therefore the following quantity is considered:

G = H d -(Z v + H).
Among the model inputs, the choice is made that the following variables are known precisely: L = 5000 (m), B = 300 (m), H d = 58 (m), and the following are considered to be random. Q (m 3 .s FORM assesses that K s is of extremely high influence, followed by Q that is of medium influence. Z v has a very weak influence and Z m is negligible. It can be noticed that the estimated failure probability is twice as small as the one estimated with crude MC, but remains in the same order of magnitude.

Sobol indices

The first-order and total indices are displayed in Table ??. It can be seen that the estimates of some indices are negative despite the fact that Sobol's indices are theoretically positive. The estimation can indeed produce negative results for values close to 0. Considering the first order indices, Z v and Z m are of null influence on their own. Q is considered to have a minimal influence (1% of the variance of the indicator function) by itself, and K s explains 24% of the variance on its own. When considering the total indices, it can be noticed that both Z v and Z m have a weak impact on the failure probability. On the other hand, Q has a major influence on the failure probability. K s total index is close to one, therefore K s explains (with or without any interaction with other variables) almost all the variance of the failure function.

Let us compare the informations provided by the Sobol' indices with the information provided by the importance factors. One cannot conclude from the total Sobol' indices that Z m is not influent whereas the importance factors assess that this variable is of negligible influence. Additionally, the total Sobol' index associated to K s and Q state that both these variables are of high influence whereas the importance factors state that K s is of high influence and Q is of medium influence.

Density modification based reliability indices

The method presented throughout this article is applied on the flood case. Only the mean shifting will be applied here. The modified pdf are given in Table 7 (Appendix D)and a numerical trick is used to deal with truncated distributions, as stressed in Appendix E. One can notice that the different inputs follow various distributions (unlike the other examples), thus the question of "equivalent" perturbation arises. It will be discussed further in Section 5. Here the choice has been made to shift the mean relatively to the standard deviation, hence including the spread of the various inputs in their respective perturbation. So for any input, the original distribution is perturbed so that its mean is the original's one plus δ times its standard deviation, δ ranging from -1 to 1 with 40 points.

Figure 7 assesses that an increasing of the mean of the inputs increases the failure probability slightly for Z v , strongly for Q, and diminishes it slightly for Z m and strongly for K s . This goes the opposite way when decreasing the mean. In terms of absolute modification, K s and Q are of same magnitude, even if K s has a slightly stronger impact. On the other hand, the effects of mean perturbation on Z m and Z v are negligible. The CI associated to Q and K s are well separated from the others, except in a δ = -.3 to .3 zone. The CI associated to Z v and Z m overlap. Thus even though the indices seem to have different value, it is not possible to conclude with certainty about the influence of those variables.

Conclusion on the method

The method presented in this paper gives relevant complementary information in addition of traditional SA methods applied to a reliability problem. Traditional SA methods provide variable ranking, whereas the proposed method provides an indication on the variation in the probability of failure given the variation of parameter δ. This is useful when the practitioner is interested on which configurations of the problem lead to an increase of the failure probability. This might also be used to assess the conservatism of a problem, if every variations of the input lead to decrease in the probability of failure. Additionally, it has two advantages:

• The ability for the user to set the most adapted constraints considering his/her problem/objective,

• The MC framework allowing to use previously done function calls, thus limiting the CPU cost of the SA, and allowing the user to test several perturbations.

We argue that with an adapted perturbation, this method can fulfil the three presented engineers objective.

Equivalent perturbation

The question of "equivalent" perturbation arises from cases where all inputs are not identically distributed. Indeed, problems may emerge when some inputs are defined on infinite intervals and when other inputs are defined on finite intervals (such as uniform distributions). For instance, consider a two-dimensional model with one Gaussian distribution and one uniform distribution as inputs. Thus, a mean shift will be a translation for the first input, whereas it will lead to a Dirac distribution in one endpoint for the other input. Hence, a mean shift cannot be considered as an "equivalent" perturbation. A "relative mean shift" seems promising idea. But if we consider a model with two Gaussian inputs of equal variance 1 and of mean respectively 1 and 10000. Then, a relative mean shift of 10% will result in Gaussian distributions with mean respectively 1.1 and 11000, and still variance 1. This counter-example shows that relative mean shift might not be an adequate perturbation in terms of "equivalence".

Support pertubation

In the examples given throughout this paper, the perturbations of the inputs left the support of those variables unperturbed. However, the practitioner might be interested in the sensitivity on the boundaries of the support. The proposed method will be applied with support perturbations in further tests. However, we stress that given the estimation method (reverse importance sampling), it is mandataory that the support of the perturbed density is included in the support of the original density. Thus one cannot perturb the inputs so that the perturbed support is wider than the original one.

Further work

Finally, these first results provide some avenues for future research:

• Adapting the estimator of the indices S iδ in term of variance reduction and of number of function calls. Further work will be made with importance sampling methods, and possibly subset methods. The use of sequential methods [START_REF] Bect | Sequential design of computer experiments for the estimation of a probability of failure[END_REF] may also be tested.

• Second, there is a need to find a way to perturb "equivalently" several distributions of different natures.

To answer this problem, perturbations that are not based upon moment constraints can be proposed.

For instance, quantile constraints might be considered. As far as we noticed, in most cases the values of the input leading to the failure event comes from the tails of the input distributions. What if these tails were badly modeled ? Therefore a perturbation based on the quantiles is proposed. More precisely, suppose that the weigth of the left tail is increased. That is to say that the value q 5% becomes for the modified density, for instance the δ quantile. This writes: 1 ]-∞;q 5% ] (x)f mod (x)dx = δ.

We believe that this kind of perturbation is equivalent with respect to inputs of different natures.

A perturbation based on an entropy constraint might also be proposed. The differential entropy of a distribution can be seen as a quantification of uncertainty [START_REF] Auder | Global sensitivity analysis based on entropy[END_REF]. Thus an example of (non-linear) constraint on the entropy can be:

-f X iδ (x) log f X iδ (x)dx = -δ f Xi (x) log f Xi (x)dx.
Yet further computations have to be made to obtain a tractable solution of the KL minimization problem under the above constraint.

Thus, using the results of [START_REF] Csiszár | I-divergence geometry of probability distributions and minimization problems[END_REF], one has λ * = arg min H(λ).

The expression of the gradient of H with respect to the j th variable is

∇ j H(λ) = g j (x)f i (x) exp( K k=1 λ k g k (x))dx exp ψ i (λ) -δ j .
Similarly, the expression of the second derivative of H with respect to the h th and the j th variables is

D hj H(λ) = g h (x)g j (x)f i (x) exp( K k=1 λ k g k (x))dx exp ψ i (λ) - g j (x)f i (x) exp( K k=1 λ k g k (x))dx exp ψ i (λ) g h (x)f i (x) exp( K k=1 λ k g k (x))dx exp ψ i (λ)
.

This method has been used in this paper for computing the optimal vector λ * when a variance shifting was applied. The integrals were evaluated with Simpson's rule.

B Proofs of the NEF properties

In this Appendix, the details of the calculus for the Proposition 2.4 are provided.

NEF specificities : If the original density f i (x) is a NEF, then under a set of K linear constraints on f (x), one has :

f (x) = b(x) exp [xθ -η(θ)] ,
thus :

f δ (x) = f (x) exp K k=1 λ k g k (x) -ψ(λ)
The regularization constant from (4) can be written as:

ψ(λ) = log b(x) exp xθ + K k=1 λ k g k (x) -η(θ) dx (14)
If the integral on ( 14) is finite, f δ exists and is a density.

Mean shifting With a single constraint formulated as in ( 6), ( 14) becames :

ψ(λ) = log b(x) exp [xθ + λx -η(θ)] dx = log b(x) exp [x (θ + λ) -η(θ) + η(θ + λ) -η(θ + λ)] dx if η(θ + λ) is well defined. ψ(λ) = (η(θ + λ) -η(θ)) + log b(x) exp [x (θ + λ) -η(θ + λ)] dx = η(θ + λ) -φ(θ) since b(x) exp [x (θ + λ) -η(θ + λ)] = f θ+λ (x)
with notation from (2.4), is a density of integral 1. Thus

f δ (x) = b(x) exp [xθ -φ(θ)] exp [λx -η(θ + λ) + η(θ)] = b(x) exp [x [θ + λ] -η(θ + λ)] = f θ+λ (x)
Thus the mean shifting of a NEF of CDF η(.) results in another NEF with mean η ′ (θ + λ) = δ (constraint) and variance η ′′ (θ + λ).

Variance shifting With a single constraint formulated as in [START_REF] Iooss | Revue sur l'analyse de sensibilité globale de modèles numériques[END_REF], using [START_REF] Saltelli | Sensitivity analysis for importance assessment[END_REF], the new distribution has for density:

f δ (x) = b(x) exp xθ + xλ 1 + x 2 λ 2 -ψ(λ) -η(θ)
Since λ is known or computed, and θ is also known, consider the variable change z = √ λ 2 x assuming λ 2 is strictly positive (the variable change is z = √ -λ 2 x if λ 2 is strictly negative). Thus,

f δ (x) = b( z √ λ 2 ) exp z 2 exp z √ λ 2 (θ + λ 1 ) -ψ(λ) -η(θ) = exp η (θ + λ 1 ) √ λ 2 -η(θ) -ψ(λ) c(z) exp z (θ + λ 1 ) √ λ 2 -η (θ + λ 1 ) √ λ 2 with c(z) = b( z √ λ 2 ) exp z 2 .
By (4),

ψ(λ) = log b(x) exp xθ + xλ 1 + x 2 λ 2 -η(θ) dx = log b( z √ λ 2 ) exp z 2 exp (θ + λ 1 ) √ λ 2 z -η(θ) + η (θ + λ 1 ) √ λ 2 -η (θ + λ 1 ) √ λ 2 dx = η (θ + λ 1 ) √ λ 2 -η(θ) + log c(z) exp (θ + λ 1 ) √ λ 2 z -η (θ + λ 1 ) √ λ 2 dx = η (θ + λ 1 ) √ λ 2 -η(θ)
Thus one has :

f δ (x) = c(z) exp z (θ + λ 1 ) √ λ 2 -η (θ + λ 1 ) √ λ 2
thus the variance shifting of a NEF results in another NEF parameterized by (θ+λ1) √ λ2 .

C Proofs of asymptotic properties

C.1 Proof of Lemma 3.1

Under assumption (i), we have

Supp(f iδ ) 1 {G(x)<0} f iδ (x i ) f i (x i ) f (x) dx ≤ Supp(f iδ ) f iδ (x i ) dx i = 1.
So that, the strong LLN may be applied to PiδN . Defining

σ 2 iδ = Var 1 {G(X)<0} f iδ (X i ) f i (X i ) , (15)
one has

σ 2 iδ = Supp(fi) 1 {G(x)<0} f 2 iδ (x i ) f i (x i ) d j =i f j (x j ) dx -P 2 iδ < ∞ under Condition (ii).
Therefore the CLT applies:

√ N σ -1 iδ PiδN -P iδ L -→ N (0, 1) .
Under assumption (ii), the strong LLN applies to σ2 iδN . So that, the final result is straightforward using Slutsky's lemma.

C.2 Proof of Proposition 3.1

First, note that

E P P iδ -P P iδ = E 1 N 2 N n=1 1 {G(x n )<0} N n=1 1 {G(x n )<0} f iδ (x n i ) f i (x n i ) -P P iδ = 1 N 2 E   N n=1 1 {G(x n )<0} 2 f iδ (x n i ) f i (x n i ) + N n=1 N j =i 1 {G(x n )<0} 1 {G(x j )<0} f iδ (x j i ) f i (x j i )   -P P iδ = 1
N 2 [N P iδ + N (N -1) P P iδ ] -P P iδ = 1 N (P iδ -P P iδ ) .

Assuming the conditions under which Lemma 1 is true, the bivariate CLT follows with Σ iδ = P (1 -P ) P iδ (1 -P )

P iδ (1 -P ) σ 2 iδ .
Each term of this matrix can be consistently estimated, using the results in Lemma 1 and Slutsky's lemma.

D Summary Table with modified distributions for mean shift E Numerical trick to work with truncated distribution

In the case where a mean shifting is considered on a left truncated distribution. We present a tip that can help to compute λ * . The studied trucated variable Y T has distribution f Y T . Let us denote Y ∼ f Y the corresponding nontruncated distribution. The truncation occurs for some real value a. This truncation may happen for some physical modelling reason. One has:

f Y T (y) = 1 1 -F (a)
1 [a,+∞[ (y)f Y (y).

The formal definition of M Y T (λ) the mgf of Y T for some λ is:

M Y T (λ) = 1 1 -F Y (a)
When the expression does not take a practical form, one can use numerical integration to estimate the integral terms. Unfortunately, for some heavy tailed distribution (for instance Gumbel distribution), this numerical integration might be complex or not possible. This is due to the multiplication by an exponential of y. The following tip helps to avoid such problems. Denoting M Y (λ) the mgf of the non-truncated distribution, one can remark that: The integral term is much smaller in the left heavy tailed distribution case. Therefore the numerical integration (for instance using Simpson's method) is much more precise or became possible. The same goes for M ′ Y T (λ) which has alternative expression:

M Y (λ) =
M ′ Y T (λ) = 1 1 -F Y (a) M ′ Y (λ) - a -∞
yf Y (y) exp [λy] dy .

Finally, another form of 16 is:

δ = M ′ Y (λ) -
This alternative expression may lead to more precise estimations of λ * when M Y (λ) and M ′ Y (λ) are known (which is the case for most usual distribution) since the integral term are much smaller than in the first expression . A reference to this Appendix is made in the summary table 7.

Figure 1 :

 1 Figure 1: Mean shifting (left) and variance shifting (right) for Gaussian (upper) and Uniform (lower) distributions. The original distribution is plotted in solid line, the perturbed one is plotted in dashed line.

Figure 2 :Figure 3 :

 23 Figure 2: Estimated indices S iδ for hyperplane function with a mean shifting

Figure 4 :

 4 Figure 4: Ishigami failure points from a MC sample

Figure 5 :

 5 Figure 5: Estimated indices S iδ for the thresholded Ishigami function with a mean shifting

  -1 ) follows a positively truncated Gumbel distribution of parameters a = 1013 and b = 558 with a minimum value of 0. K s (m 1/3 s -1 ) follows a truncated Gaussian distribution of parameters µ = 30 and σ = 7.5, with a minimum value of 1. Z v (m) follows a triangular distribution with minimum 49, mode 50 and maximum 51. Z m (m) follows a triangular distribution with minimum 54, mode 55 and maximum 56. A 10 5 MC sample gives an estimation of the failure probability P = 8.6 × 10 -4 .

Figure 6 :

 6 Figure 6: Estimated indices S i,Vper for the thresholded Ishigami function with a variance shifting. The upper figure is a zoom where the S i,Vper axis lies into [-0.5, 0.5]. The lower figure shows almost the whole range variation for S i,Vper . The curves cross for the value of V per that corresponds to the original variance, namely π 2 /2

Figure 7 :

 7 Figure 7: Estimated indices S iδ for the flood case with a mean perturbation

  +∞ a f Y (y) exp [λy] dy.Let us recall that we are looking for λ * such as:δ = M ′ Y T (λ * ) M Y T (λ * ) = +∞ a yf Y (y) exp [λy] dy +∞ a f Y (y) exp [λy] dy .

  Y (y) exp [λy] dy = a -∞ f Y (y) exp [λy] dy + +∞ a f Y (y) exp [λy] dy Thus another expression for M Y T (λ) is:M Y T (λ) = 1 1 -F Y (a) M Y (λ) -a -∞f Y (y) exp [λy] dy .

1 σ √ 2π exp - 1 2 x i -δ σ 2 λ * = δ-µ σ 2 1 2 x i -µ-σ 2 λ * σ 2 δ

 112212 pdf f iδ Link between λ * and δ NEF(θ)N EF (θ + λ * ) f iδ (x i ) = b(x i ) exp [x i [θ + λ * ]η(θ + λ * )] η ′ (θ + λ * ) = δ Special case of NEF: N (µ, σ) N (δ, σ) f iδ (x i ) = Uniform distribution: U [a,b] ∝ truncated exponential f iδ (x i ) = λ * e λ * b -e λ * a 1 [a,b] (x i )e λx i δ = 1 (b-a) e λ * b (λ * b-1)+e λ * a (1-λ * a) λ * (e λ * b -e λ * a ) Left Tr Gaussian N T (µ, σ, a) N T (µ + σ 2 λ * , σ, a) f iδ (x i ) = 1 [a,+∞[ (x i ) 1= µ + σ 2 λ *σ φ a-(µ+σ 2 λ * ) σ 1-Φ a-(µ+σ 2 λ * ) σ Triangle T (a, b, c) f iδ (x i ) = exp(x i λ *ψ(λ * ))f (x i ) δ = (a-1 λ * )e λ * a (b-c)+(b-1 λ * )e λ * b (c-a)+(c-1 λ * )e λ * c (a-b) e λ * a (b-c)+e λ * b (c-a)+e λ * c (a-b) Left Tr Gumbel G T (µ, β, a) f iδ (x i ) = exp(x i λ *ψ(λ * ))f (x i ) δ = M ′ Y (λ * )-a -∞ yf Y (y) exp[λ * y]dy M Y (λ * )-a -∞ f Y (y) exp[λ * y]dy with : M Y (λ * ) = Γ (1β) exp [λ * µ] M ′ Y (λ * ) = Γ (1β) exp [λ * µ] µ -β̥ (0) (1λ * )

a-∞

  yf Y (y) exp [λy] dy M Y (λ) -a -∞ f Y (y) exp [λy] dy .

Table 1 :

 1 Importance factors for hyperplane function

Table 2 :

 2 Sobol' indices for hyperplane function

	4

Table 4 :

 4 Sobol' indices for Ishigami function

	Sobol' Index	S 1	S 2	S 3	S T 1	S T 2	S T 3
	Mean	0.0234	0.0099	0.0667	0.8158	0.6758	0.9299
	C.o.v.	0.0072	0.0051	0.0095	0.0156	0.0216	0.0094

Table 5 :

 5 The algorithm FORM converges to a design point (1.72, -2.70, 0.55, -0.18) in 52 function calls, giving an approximate probability of PF ORM = 5.8 × 10 -4 . The importance factors are displayed in Table5. Importance factors for flood case

	4.4.1 FORM
	Variable	Q	Ks	Zv	Zm
	Importance factor	0.246	0.725	0.026	0.003

Table 6 :

 6 Sobol' indices for Ishigami function

	Sobol' Index	S 1	S 2	S 3	S T 1	S T 2	S T 3
	Mean	0.0234	0.0099	0.0667	0.8158	0.6758	0.9299
	C.o.v.	0.0072	0.0051	0.0095	0.0156	0.0216	0.0094

Table 7 :

 7 Modified distributions for mean shifting.Note that Φ(.)is the cdf of the standard normal distribution, and φ(.) is its pdf.

Acknoweldgements

Part of this work has been backed by French National Research Agency (ANR) through COSINUS program (project COSTA BRAVA noANR-09-COSI-015). We thank Dr. Daniel Busby (IFP EN) for several discussions and Emmanuel Remy (EDF R&D) for proofreading. We also thank an anonymous reviewer and the associated editor for their helpful comments. All the statistical parts of this work have been performed within the R environment, including the mistral and sensitivity packages.

Appendices

A Computation of Lagrange multipliers

Let H be the Lagrange function: