Philippe Roudier
email: philippe.roudier@locean.upmc.fr

Jafet C M Andersson

Chantal Donnelly

Luc Feyen

Wouter Greuell

Fulco Ludwig

Projections of future floods and hydrological droughts in Europe under a +2°C global warming

We present an assessment of the impacts of a +2°C global warming on extreme floods and hydrological droughts (1 in 10 and 1 in 100 year events) in Europe using eleven bias-corrected climate model simulations from CORDEX Europe and three hydrological models. The results show quite contrasted results between northern and southern Europe. Flood magnitudes are expected to increase significantly south of 60 o N, except for some regions (Bulgaria, Poland, south of Spain) where the results are not significant. The sign of these changes are particularly robust in large parts of Romania, Ukraine, Germany, France and North of Spain. North of this line, floods are projected to decrease in most of Finland, NW Russia and North of Sweden, with the exception of southern Sweden and some coastal areas in Norway where floods may increase. The results concerning extreme droughts are less robust, especially for drought duration where the spread of the results among the members is quite high in some areas. Anyway, drought magnitude and duration may increase in Spain, France, Italy, Greece, the Balkans, south of the UK and Ireland. Despite some remarkable differences among the hydrological models' structure and calibration, the results are quite similar from one hydrological model to another.

Finally, an analysis of floods and droughts together shows that the impact of a +2°C global warming will be most extreme for France, Spain, Portugal, Ireland, Greece and Albania. These results are particularly robust in southern France and northern Spain.

Introduction

In Europe, freshwater resources are of crucial importance for many sectors including agriculture, hydropower generation, cooling water for power plants and domestic and industrial water supply. At the same time, water can have a direct impact on safety and livelihoods through floods that can lead to disastrous human and economic losses (e.g. the 2013 floods in Central Europa resulted in a loss of more than €12bn1 ). During the last decades, several studies have underlined that water resources and especially river flows have had strong variations across Europe [START_REF] Kovats | Europe[END_REF] due to climate, water extractions and land use change [START_REF] Sterling | The impact of global land-cover change on the terrestrial water cycle[END_REF]. Even if the relative share of these three driving factors is difficult to assess over the past, it is clear that the strong climate changes expected in Europe for the 21 st century will have a significant impact on river flows [START_REF] Cisneros | Freshwater resources[END_REF].

Thus, a first step in order to make relevant mitigation and adaptation policies is to develop a clear picture of the potential future stream flows extremes.

Changes in river flow extremes at a +2°C global warming are currently of central interest as this is the global target defined by policymakers to lower international greenhouse gases emissions (European Commission, 2007). Therefore, in this study, we do not select a specific time period; rather we focus on defining the hydrological impacts in a world with a +2°C global warming relative to pre-industrial levels (Vautard et al., 2014). Describing the impacts of a +2°C global warming on topics such as water resources, agriculture or infrastructures is the main aim of the FP7 project IMPACT2C under which this study was conducted.

Future changes in hydrological extremes are still highly uncertain. There is a general consensus that in most part of the world climate change will result in more rainfall extremes (IPCC, 2012). However, to what extent this will affect hydrological extremes is still highly uncertain and differs from region to region. To address this uncertainty many previous studies have used multiple climate models to force a single hydrological model (e.g. [START_REF] Dankers | Flood hazard in Europe in an ensemble of regional climate scenarios[END_REF]; Feyen et al. (2009)). However recent work, mainly at global scale, has shown that not only the choice of the climate models affect future change in hydrological extremes: different hydrological models give sometimes very different results too [START_REF] Haddeland | Multimodel Estimate of the Global Terrestrial Water Balance: Setup and First Results[END_REF]; [START_REF] Schewe | Multimodel assessment of water scarcity under climate change[END_REF]; [START_REF] Prudhomme | Hydrological droughts in the 21st century, hotspots and uncertainties from a global multimodel ensemble experiment[END_REF]; [START_REF] Dankers | First look at changes in flood hazard in the Inter-Sectoral Impact Model Intercomparison Project ensemble[END_REF]). Many multiple climate and hydrological models impact assessments have used global climate models. These global climate models tend to have significant biases in the representation of precipitation extremes.

The aim of this study is to improve the assessment of hydrological extreme impacts at the European continent scale including a better description of uncertainties through the use of multiple hydrological models and state-of-the art climate projections from the high-resolution CORDEX project. We first focus on the skill of each hydrological model to simulate extreme floods and hydrological drought (section 3.1) and we next assess the impact that a +2°C warming would have on meteorological variables which are relevant for flood and drought generation (section 3.2) and on extreme flows (return periods are 10 and 100 years, section 3.3 and 3.4). We finally summarize the results about flood magnitude, droughts magnitude and drought duration for twelve European cities and with a final assessment on which European region will be the most affected by extreme flows at +2°C warming.

Material and data

Data

Forcing data

This study uses a sub-ensemble from the latest (as of 2014) ensemble of high-resolution, dynamically downscaled daily climate simulations from CMIP5, 5 th Coupled Model Intercomparison Project and CORDEX, A Coordinated Regional Downscaling Experiment [START_REF] Jacob | EURO-CORDEX: new high-resolution climate change projections for European impact research[END_REF]. The 11 ensemble members were chosen to be representative of the larger ensemble and consist of 5 GCM/RCM combinations and 3 Representative Concentration Pathways (RCPs) (Moss et al. ( 2010)) (Table 1). One of the main differences with previous impact studies is that rather than focusing on a future time slice, change is quantified at the 30-year period when each driving GCM reaches +2°C in global mean temperature relative to pre-industrial levels (1881( -1910( , Vautard et al. (2014))). Thus, the ensemble is a representation of a world where greenhouse gas emissions have caused a +2°C global warming. The climate variables were regridded to a 0.25° resolution grid. As climate models are known to be affected by some biases that can clearly affect the results [START_REF] Chen | Finding appropriate bias correction methods in downscaling precipitation for hydrologic impact studies over North America[END_REF], daily climate variables (precipitation, maximum, minimum and average temperature, dew point temperature, shortwave and longwave downward radiations) were bias-corrected using quantile mapping [START_REF] Themeßl | Empirical-statistical downscaling and error correction of daily precipitation from regional climate models[END_REF][START_REF] Themeßl | Empirical-statistical downscaling and error correction of regional climate models and its impact on the climate change signal[END_REF] and [START_REF] Wilcke | Multi-variable error correction of regional climate models[END_REF]). Bias-correction was made using the E-OBS gridded observational dataset [START_REF] Haylock | A European daily highresolution gridded data set of surface temperature and precipitation for 1950-2006[END_REF] as a reference. These bias corrected data were then used to run the three hydrological models . Climatic changes were subsequently assessed by comparing the +2°C period with the baseline period 1971-2000 (Table 1).

Observed discharge data

In order to assess the skills of the hydrological models in representing specific floods and droughts, we used 428 discharge stations over Europe selected from the Joint Research Center database which gathers several sources like the Global Runoff Data Centre (2013) data or other publicly available datasets (e.g., HYDROBanque, CEDEX, Norwegian Water Resources and Energy Directorate, National river flow archive, Waterinfo.be , eHyd). The original database was filtered to include only stations with data for the time period 1971-2000 and with no day with missing data in order to not to miss any extreme discharges. Model performance was validated using the median of the values from each hydrological model forced by the ensemble of climate models over the 30 years control period.

Summary of models

Three pan-European hydrological models were used to simulate daily discharge: Lisflood [START_REF] Burek | LISFLOOD Distributed Water Balance and Flood Simulation Model[END_REF], E-Hype (Donnelly et al., 2015) and VIC [START_REF] Liang | A simple hydrologically based model of land surface water and energy fluxes for general circulation models[END_REF]. The models differ in both complexity of process description, input data and setup. For example, differences include spatial resolution (one model is subbasin rather than grid based), description of evapotranspiration and snow processes, the number of soil layers and depth assumptions and the calibration procedure. A summary of these models can be found in the Appendix and more details are available in [START_REF] Greuell | Evaluation of five hydrological models across Europe and their suitability for making projections underof climate change[END_REF].

Because each of these models has a different output resolution (and in the case of E-HYPE, a highresolution subbasin rather than grid output), the output from each model was regridded to a comparable 0.5°*0.5° grid.

Methodology for extreme value analysis

High flows

We focus here on the discharge magnitude of the 1 in 10 and 1 in 100 year floods (QRP10 and QRP100). These two return periods were chosen as the first one represents an extreme event occurring often enough that it is remembered by individuals and communities and the other one is a standard value used in some countries to design flood protection. In order to compute QRP10 and QRP100 for each pixel of the grid, we followed a typical extreme value analysis fitting methodology (see e.g. [START_REF] Roudier | Calculation of design rainfall and runoff on the Bani basin (Mali) : a study of the vulnerability of hydraulic structures and of the population since the drought[END_REF]). First, the daily maximum discharge was selected for each year of the 30 years long period, then fitted a Generalized Extreme Value (GEV) distribution using the L-moments and finally, we calculated QRP10 and QRP100 using the fitted function. The goodness-of-fit was also checked using the Anderson-Darling test (at 5%), as recommended by [START_REF] Meylan | Hydrologie Frequentielle. Science et ingénierie de l'environnement. PPUR Moss RH et al[END_REF]. For each of the 3 hydrological models, and each of the 11 climate runs, the relative QRP10 and QRP100 change between the +2°C period and the baseline was computed, resulting in a set of 33 relative changes. To describe this set, in this paper the median of all the ensemble members is used (the combined climatological and hydrological model ensemble) rather than the mean, in order to avoid giving excessive weight to potential outliers. The significance of changes is also assessed using a Wilcox test (5% threshold) between future period and baseline. Therefore, in all relative change plots, we set as missing values pixel that do not pass the test. Moreover, those that do not pass the goodness-of-fit test were also set as missing value.

Low flows

We focus for low flows on the same return periods as for the high flow analysis, following the methodology used by Feyen and Dankers (2009). First the daily discharge time-series is smoothed applying a seven day moving average in order to remove the day-to-day variations and then for the magnitudes of low flows (QlowRP10 and QlowRP100), the lowest smoothed discharge event is selected, every year. For the duration of low flows a threshold approach is used with the 20 th flow percentile of the smoothed flow duration curve as threshold. After computing the 20 th percentile (for each pixel, projection and hydrological model) we select all the days that have a smoothed discharge value below this 20 th percentile and then we compute the duration of each event below that threshold and select the maximum drought duration for each year. Finally, for both drought magnitude and duration, the rest of the methodology is equal to section 2.3.1 (fitting a GEV distribution on each set of 30 values in order to compute QlowRP10 and QlowRP100 magnitude and duration; we also set as missing values pixel that do not pass the Wilcox or Anderson-Darling test)

Results and discussion

Hydrological model validation and selection for ensemble projections

A broad and detailed validation of the hydrological models focusing on average conditions is presented in [START_REF] Greuell | Evaluation of five hydrological models across Europe and their suitability for making projections underof climate change[END_REF]. The aim of the validation presented here is to assess the models' skill to simulate specific indicators used in this study. For each hydrological model we therefore assess if the median QRP100 (QlowRP100) computed based on the 11 bias corrected climate runs is close to the QRP100 (QlowRP100) computed with observed discharge data. Results are shown in Figure S1 A likely explanation for these differences in performance could be the way these models are tuned to observation data. The Lisflood model is calibrated in individual catchments using a highresolution (5 km) interpolated observation data set (EFAS-meteo, see [START_REF] Ntegeka | EFAS-Meteo: A European daily high-resolution gridded meteorological data set for 1990 -2011[END_REF]) and default parameters in ungauged catchments. The VIC model uses a general set of parameters applicable anywhere in the model domain, but linked to soil-type and landuse. This parameter tuning was done using a specific forcing dataset described in [START_REF] Nijssen | Global Retrospective Estimation of Soil Moisture Using the Variable Infiltration Capacity Land Surface Model, 1980-93[END_REF] .The E-HYPE model also uses a general set of parameters, linked to soil-type and landuse. For E-HYPE, these parameters were originally calibrated to a small set of representative gauged basins for each soil and land use using a corrected ERA-INTERIM forcing data set (see Donnelly et al. (2015)), but for this study an evapotranspiration parameter was slightly adjusted to better balance the model performance with the E-OBS data set. This was deemed important by the E-HYPE modeling team as the climate scenarios were bias-corrected to the E-OBS data set and recent studies have shown that bias-correction and climate impact results can vary between reference observation data sets (e.g. [START_REF] Gutmann | An intercomparison of statistical downscaling methods used for water resource assessments in the United States[END_REF]). It is therefore not surprising that overall E-HYPE performs well when forced with data bias-corrected to E-OBS, at least regarding mean discharge. Lisflood is likely to perform better in regions where the E-OBS forcing is not significantly different to the EFAS-meteo data as the other models compromise performance in individual catchments for relatively good performance in multiple catchments using the same parameters. For hydrological extremes, again Lisflood can be expected to perform well in calibration catchments. E-HYPE may outperform VIC for extremes because the representative gauged basins used to calibrate E-HYPE are generally smaller catchments (1000 to 5000 km 2 ) where runoff generating processes dominate over routing and lake processes.

Given the variability in model performance for extremes, and perhaps inability of some models to successfully represent these extremes, we suggest to remove for the initial selection the hydrological models that have a NSE and a NSE(log) below 0.5 for the index studied, i.e. QRP100 and QlowRP100.

Based on the NSE threshold, an ensemble consisting of results from E-HYPE, Lisflood, and VIC will be used for floods, and Lisflood and E-HYPE for low flows.

Summary of projected changes in meteorological variables pertaining to floods and droughts

We first aim at studying here the future evolution of some of the drivers of low and high flow changes. We therefore plot in Figure 1 the projected change in maximum annual snowpack (outputs from the 3 hydrological models) and intense daily rainfall (return period is 10 years, computed using a GEV distribution using the aforementioned methodology for discharge) at +2°C warming. Results show that there is generally a clear decrease in maximum snowpack, reflecting the future warming over Europe. The regions that are less affected like southern Spain or southern Italy are those with currently marginal snowfall. The snowpack decrease is particularly important in Fenno-Scandinavia and the Alps.

Intense rainfall events are projected to increase significantly over the whole continent, with no particular spatial pattern, which is consistent with previous studies [START_REF] Madsen | Review of trend analysis and climate change projections of extreme precipitation and floods in Europe[END_REF]; [START_REF] Rajczak | Projections of extreme precipitation events in regional climate simulations for Europe and the Alpine Region[END_REF]).

Floods (QRP10 and QRP100)

For changes in flood magnitude there is a clear North to South gradient, with a strong increase with a strong increase in flood magnitudes south the 60°N line, except for some regions in Bulgaria, Czech Republic, Poland, the western Balkans, the Baltic countries, and southern Spain where no significant changes can be detected (Figure 2). Almost everywhere the increase in 100 year floods (QRP100) is stronger than the 10 year floods (QPR10). Floods are even increasing in areas such as southern Mediterranean where the average discharge is projected to decrease [START_REF] Greuell | Evaluation of five hydrological models across Europe and their suitability for making projections underof climate change[END_REF].

However flood changes are consistent with the extreme rainfall changes south of 60°N (Figure 1). Above the 60 o N line, the situation is more heterogeneous with a relatively strong decrease in flood magnitude in parts of Finland, NW Russia and North of Sweden with the exception of southern Sweden and some coastal areas in Norway where increases in floods are projected. Projections of decreasing flood magnitudes are mainly due to the decreases in snowpack in areas where most of the floods are caused by spring snowmelt in combination with rainfall. Increases in flood magnitude in Scandinavia are mainly seen in coastal areas where the rain-fed floods will increase (similar results were found by [START_REF] Vormoor | Climate change impacts on the seasonality and generation processes of floods -projections and uncertainties for catchments with mixed snowmelt/rainfall regimes[END_REF]). Note also that the date of occurrence of annual maximum discharge is expected to be earlier in spring for these areas (Figure S8) while for the rest of Europe, changes are quite limited.

Following this reduction in snow-melt floods in Fenno-Scandinavia and according to the decreasing snow pattern shown in Figure 1, it is expected a similar decrease in floods would be expected for the Alps. However, our results shown an increase in floods which is also reported in several other studies [START_REF] Rojas | Assessment of future flood hazard in Europe using a large ensemble of bias-corrected regional climate simulations[END_REF], [START_REF] Gobiet | 21st century climate change in the European Alps-A review[END_REF] or [START_REF] Köplin | Seasonality and magnitude of floods in Switzerland under future climate change[END_REF]). For Switzerland, [START_REF] Köplin | Seasonality and magnitude of floods in Switzerland under future climate change[END_REF] showed that there are large contrasts in the flood regime in the Alps, and snowmelt-floods only occur in a limited part of country. Thus, at a 0.5°*0.5° resolution, we are not able to capture some local decreases of QRP10 and QRP100 due to earlier snowmelt. This hypothesis is confirmed by the same analysis performed only for Lisflood at 5*5 km grid that shows for some areas a decrease of QRP100.

Secondly, we think that at the 0.5°*0.5° resolution over the Alps, there is for most of the pixels a change of flood regime (from snowmelt to rain-fed or mixed snowmelt-rain-fed) rather than only an earlier occurrence of snowmelt floods, leading to a different situation compared to Fenno-Scandinavia. This hypothesis, is supported by Figure S8 which shows (i) a clear reduced number of high flows in June (southern Alps) and July (rest of the Alps) but (ii) no strong increase in earlier spring. A deeper analysis, at a finer scale, is needed to fully understand all the changes in the Alps but beyond the scope of the present paper.

Finally, we also provide model-by-model details of these results in the Appendix. Generally, all three hydrological models show the same pattern. However, in Western Europe and the Mediterranean area, Lisflood projects a stronger flood increase than the other models.

Even if it is difficult to compare with studies that do not focus on exactly the same set of parameters (climate models, resolution, time period), these results are mainly in line with the recent literature as reviewed by [START_REF] Madsen | Review of trend analysis and climate change projections of extreme precipitation and floods in Europe[END_REF] (table 3 in their paper). Out of 22 studies dealing with future flood projections in Europe, 4 of them are not directly comparable (they differentiate winter and spring floods), 14 have a global agreement with our findings and 4 give different results. Those four studies include (i) the UK, where [START_REF] Reynard | Regionalised Impacts of Climate Change on Flood Flows[END_REF] find few catchments with changes in flood frequency above +20% and [START_REF] Kay | RCM rainfall for UK flood frequency estimation. I. Method and validation[END_REF], with a limited ensemble of driving climate runs, depict a decrease in flood magnitude in south and east England; (ii) France (Seine and Somme rivers) where [START_REF] Ducharne | Climate change impacts on water resources and hydrological extremes in northern France[END_REF] cannot conclude on any robust change while we find a significant increase of flood events and (iii) eastern Germany, Poland and southern Sweden where [START_REF] Rojas | Assessment of future flood hazard in Europe using a large ensemble of bias-corrected regional climate simulations[END_REF] find decreases in flood magnitude (this latter study which covers the whole continent agrees however with our results for the rest of Europe). Finally, [START_REF] Dankers | First look at changes in flood hazard in the Inter-Sectoral Impact Model Intercomparison Project ensemble[END_REF] find more contrasted changes for 1 in 30 year flood change in Europe using several GCMs and hydrological models: according to these results, floods would decrease in large parts of Europe including Greece, Italy and eastern Europe.

Low flows: magnitude and duration

Low floods (QlowRP10) are expected to decrease for many countries mainly located in the southern part of Europe: Spain, France, Italy, Greece, the Balkans and also south of the UK and Ireland. This is mainly due to less rainfall [START_REF] Rajczak | Projections of extreme precipitation events in regional climate simulations for Europe and the Alpine Region[END_REF] and also higher potential evapotranspiration in some regions like Italy (see Figure S4 and Van Vliet et al.,2015). The duration of these droughts is also increasing (Figure 3, top), especially in Spain. For the rest of Europe, the projections show generally a decrease of drought magnitude and duration. Moreover, changes are not significant in some areas like western Germany or southern Sweden: areas with in-significant changes are larger for low flows than for floods. This reduction of low flow duration and magnitude is mainly caused (i) by less snowfall and more precipitation for areas with low flows in winter and (ii) by a general increase of rainfall for areas with low flows in summer (Vautard et al., 2014). Both hydrological models depict generally the same pattern (see Figure S5 for magnitude). However, for southern Sweden, E-HYPE predicts a small decrease in low-flow magnitudes and Lisflood an increase. Lisflood also tends to predict a significant increase of QlowRP10 across Eastern Europe while for E-HYPE this change is also positive, but not significant (white pixels). The same spatial pattern is found for the 100 year return period (see Figure S6) but the results are generally less significant, see for example southern France.

Our results on drought magnitude change are similar to those presented by [START_REF] Forzieri | Ensemble projections of future streamflow droughts in Europe[END_REF].

Focusing on the QlowRP20 they also find for the 2080s a decrease for the Mediterranean and the UK but their findings for Sweden and Norway show an almost uniform increase of low flows although we observed a North/South difference that may reflect the different set of climate models used. [START_REF] Prudhomme | Hydrological droughts in the 21st century, hotspots and uncertainties from a global multimodel ensemble experiment[END_REF], using several climate and hydrological models find a general increase of hydrological droughts over Europe, but they focus on less extreme droughts, and they use RCP 8.5, at the end of the century. At a more local scale, other studies agree with our results. For example, in France, [START_REF] Chauveau | Quels impacts des changements climatiques sur les eaux de surface en France a l'horizon 2070 ?[END_REF] also predict a decrease in low flow magnitude; in Germany Huang et al.

(2014) find a decrease of QlowRP50 in 2080 for some areas like the Rhine basin and uncertainty elsewhere.

Uncertainties and summary of the results

We analyzed the spread in future change of hydrological extremes using three different parameters (QRP10, QlowRP10 and QlowRP10 duration) and focusing on the 1 st and 3 rd quartiles of relative changes distribution (Figure 4 and Figure S7). The spread among results is the largest for QlowRP10 duration, depicted in Figure 4 by the inter-quartile range. Despite this spread, the sign of QlowRP10 duration change (positive or negative) is the same for the first and third quartiles for areas with a very large spread (e.g. eastern Europe), i.e. the ensemble of projections somewhat agree on the direction of change, but uncertainty is high. Second, for QlowRP10 magnitude the spread between both quartiles is generally smaller especially in areas like France, the UK, Italy, Portugal and Greece where it is generally below 20%. However, both quartiles do not have the same sign in all these areas (e.g. the UK and Italy) thus depicting projections less robust than for Portugal, south of France or south of Spain.

Third, for floods (QRP10) the quartiles agreement is slightly better than for QlowRP10 magnitude except in southern Fenno-scandinavia This general larger uncertainty for low flows is mainly due to (i) differences between the two models in soil moisture and evapotranspiration calculation which are directly related to low flows [START_REF] Greuell | Evaluation of five hydrological models across Europe and their suitability for making projections underof climate change[END_REF] and (ii) the number of selected models is different for floods (VIC, E-HYPE, Lisflood) and droughts (E-HYPE, Lisflood). Figure S1 and Figure S2 show clearly that the three hydrological models perform similar for floods while for droughts, Lisflood tends to overestimate QlowRP10 and E-HYPE to underestimate it, thus resulting in a larger spread. Moreover, with this kind of assessment using quartiles, the areas with robust results are generally the same than the ones using the Wilcox test, except for some areas like southern Spain where this latter test gives significant changes.

In order to detect hotspot regions that will be subject to negative changes in several extreme flow indicators, we combined the floods, drought magnitude and drought duration changes in a single analysis (Figure 5). In most parts of France, Spain, Portugal, Ireland, Greece and Albania, the projected changes under +2°C are generally more extreme. We mean by "more extreme" regions where there is a consistent worsening (of at least 5%) in all the extreme indicators considered for a 10-years return period: more intense floods (QRP10 >+5%), more intense hydrological droughts (QlowRP10<-5%) and longer droughts (QlowRP10 duration > +5%). In parts of Norway, Sweden, Finland and western Russia future warming will see a reduction in both streamflow floods and droughts.

Conclusion

Our aim is to make a robust assessment of the impact that a +2°C global warming would have on hydrological extremes (floods, droughts magnitude and duration) in Europe by using an ensemble of eleven high-resolution RCM outputs and three pan-European hydrological models. Results show that such a warming could increase flood magnitudes (10 years and 100 years return period) significantly in most parts of Europe (e.g. about +20% close to London; and Warsaw for QRP10), even for areas where the annual rainfall is expected to decrease in the future, e.g. Spain (about +10% QRP10 close to Sevilla).

However, in Fenno-Scandinavia the situation is more contrasted with (i) a large area that is expected to have less intense snowmelt floods, occurring earlier in spring except and (ii) the southern and coastal areas of Fenno-Scandinavia where we predict an increase of rain-fed flood magnitude. In the Alps, even though snowpack is also projected to reduce, floods are expected to increase generally due to a change of flood regime from snowmelt to rain-fed and possibly because the spatial resolution of this study (0.5°*0.5°) potentially hides local decreases of intense snowmelt floods. Moreover, despite some significant differences among the hydrological models' structure and calibration, the results are quite similar from one model to another and consistent with other studies.

Future changes in hydrological drought magnitude and duration show contrasting patterns across Europe. Our projections show that for large areas of Italy, France, Spain, Greece, the Balkans, Ireland and the UK, droughts will become more intense and longer mainly due to less rainfall and higher evapotranspiration, in some areas. The sign of these changes is particularly robust in southern France, parts of Spain, Portugal and Greece. For the rest of Europe changes in droughts are not significant or there is a reduction of droughts length and magnitude, especially in northern Fenno-Scandinavia and Western Russia where the sign of the changes is also very robust.

Our results show that for a significant part of Europe there will be a clear intensification of the hydrological cycle resulting in both increases in droughts and floods. Extreme flows will be particularly harmful in Spain, Greece, France, Ireland and Albania: it is thus urgent to integrate these future changes for policy making in water resources management and flood protection design. (bypass of soil layer), exchange of soil moisture between the two soil layers and drainage to the groundwater, sub-surface and groundwater flow, and flow through river channels. Runoff produced for every grid cell is routed through the river network using a kinematic wave approach.

• VIC

The Variable Infiltration Capacity model (VIC) performed calculations on a regular 0.5x0.5 degrees lat-lon grid. In VIC each grid cell is subdivided into an arbitrary number of tiles, each covered by a particular vegetation type. Evapotranspiration is computed for each tile, whereas the soil is uniform within each grid cell. The soil consists of 3 layers with a total depth of ~3 m for most cells. The energy balance approach of Penman-Monteith [START_REF] Shuttleworth | A land-surface hydrology parameterization with subgrid variability for general circulation models[END_REF] is used to compute evapotranspiration, requiring a forcing consisting of atmospheric temperature and humidity, wind speed and incoming radiation. Penman-Monteith is also employed to simulate snow melt. VIC owes its name to the treatment of surface runoff, which is generated from net precipitation under the assumption that the infiltration capacity varies within each grid cell according to a one-parameter distribution function (Wood et al., 1992). Routing of the runoff is taken into account with the model described in [START_REF] Lohmann | A large-scale horizontal routing model to be coupled to land surface parametrization scheme[END_REF]. This model first transports the runoff generated in a grid cell to the main river leaving the grid cell with a unit hydrograph and then routs the water through the river network connecting the cells.

A2. Validation

  (see supplementary material) for floods and Figure S2 for low flows. The skills of the three hydrological models are generally better for high than low flows. Lisflood performs slightly better than the other models, according to the Nash-Sutcliffe Efficiency coefficient (NSE = 0.82) and the Root Mean Squared Error (RMSE = 961 m3/s). However, when focusing on the NSE of the logarithmic modeled and observed values (NSE(log), used for extreme values, see Krause et al., (2005)), E-Hype is slightly better than Lisflood. For low flows, Lisflood somewhat overestimates QlowRP100 while the other two models underestimate it. E-Hype has the best performance for low flows (NSE = 0.72, NSE(log)=0.68, RMSE = 75 m3/s) compared to Lisflood (NSE=0.54, NSE(log)=0.51, RMSE=95 m3/s) and VIC (NSE=0.38, NSE(log)=0.38, RMSE=112 m3/s).

Figures

  Figures

Figure 1 :

 1 Figure 1: Change in maximum annual snowpack (all three hydrological models, top) and change in intense rain (RP10, all three hydrological models, bottom)

Figure 3 :

 3 Figure 3: characteristics of low flows (RP10): duration (top) and magnitude (bottom). The median is computed over 22 ensemble members. Only significant changes (i.e. passing the Wilcox test at 5%) are shown here. When QlowRP10 is zero for the baseline period, we set the relative change as missing value.

FigureFigure 5 :

 5 Figure 4: quartiles of relative change. Left column: absolute difference between the 1 st and 3 rd quartiles of relative change (in %) for QRP10 (top), QlowRP10 magnitude (middle) and QlowRP10 duration (bottom). Right column: agreement on the sign of the 1 st and 3 rd quartiles (e.g. the green area means that both quartiles are positive or negative)

Figure S1 :

 S1 Figure S1: Observed QRP100 428 stations, years 1971/2000) vs. modelled one (median QRP100 over 5 climate members), for each of the three hydrological models. The dashed line is y=x, NSE is Nash-Sutcliff Efficiency coefficient, NSE(log) is Nash-Sutcliff Efficiency coefficient of log(observed QRP100) vs. log(modeled QRP100) and RMSE is the root mean squared error.

Figure S5 :

 S5 Figure S5: QlowRP10 magnitude relative change, for both hydrological models used for low flows: Lisflood (top) and E-HYPE (bottom). The median is computed over 11 members. Only significant changes are shown here. When QlowRP10 is zero for the baseline period, we set the relative change as missing value

Figure

  Figure S7: 25 th percentile (left column) and 75 th percentile (right one) relative change for QRP10 (top row), QlowRP10 (middle row) and QlowRP10 duration (bottom row). For each pixel, the percentiles are computed over 33 members (QRP10) or 22 members (QlowRP10 and QlowRP10 duration).

http://www.munichre.com/en/media-relations/publications/press-releases/2013/2013-07-09-press-release/index.html

Acknowledgments

The authors would like to thank the FP7 project IMPACT2C and all the contributing members for funding this study and providing climate data. Moreover, we thank Goncalo Gomes from JRC for his help with observed discharges and Alessandra Bianchi for GIS support. We finally think three anonymous reviewers for their helpful comments.

Supplementary material

A1. Summary of models

• E-HYPE

The E-HYPE model (Donnelly et al., 2015) is a pan-European application of the Hydrological Predictions for the Environment (HYPE, [START_REF] Lindström | Development and testing of the HYPE (Hydrological Predictions for the Environment) water quality model for different spatial scales[END_REF]) which simulates hydrological variables in more than 35000 sub-basins (median size 200 km 2 ) across Europe. A modified Hargreaves-Semani equation using daily minimum and maximum temperature and a time constant radiation that is assumed to vary with latitude is used to compute evapotranspiration. Snow melt is calculated using a degree day method. In HYPE evapotranspiration, snow and runoff generation calculations are made for hydrological response units (HRUs) representing unique land use and soil type combinations and consisting of up to three soil layers down to 2.5 m. Each sub-basin may consist of any number of HRUs.

The model is forced by daily precipitation and temperature and then calculates flow paths in the soil based on snow melt, evapotranspiration, surface runoff, infiltration, percolation, macropore flow, tile drainage, and lateral outflow to the stream from soil layers with water content above field capacity.

HRUs are connected directly to the stream and act in parallel. The groundwater level in each HRU is fluctuating, may saturate the soil layers and water may percolate between sub-basins. Routing of the runoff is made along local and main rivers within each sub-basin and includes delay and dampening processes. Lakes, where they exist, cause delay of flow using the weir equation and may be local (off the main river) or main (on the main river).


Jordan Ninin

Global Optimization based on Contractor Programming: an Overview of the IBEX library

IBEX is an open-source C++ library for constraint processing over real numbers. It provides reliable algorithms for handling non-linear constraints. In particular, roundoff errors are also taken into account. It is based on interval arithmetic and affine arithmetic. The main feature of IBEX is its ability to build strategies declaratively through the contractor programming paradigm. It can also be used as a black-box solver or with an AMPL interface. Two emblematic problems that can be addressed are: (i) System solving: A guaranteed enclosure for each solution of a system of (nonlinear) equations is calculated; (ii) Global optimization: A global minimizer of some function under non-linear constraints is calculated with guaranteed and reliable bounds on the objective minimum.

Kernel of IBEX

. These methods and algorithms do not consider single numerical values, or floating-point numbers, but manipulate sets. The interval arithmetic offers a solid theoretical basis to represent and to calculate with subsets of R n .

Interval Arithmetic

An interval is a closed connected subset of R. A non-empty interval [x] can be represented by its endpoints: [x] = [x, x] = {x : x ≤ x ≤ x} where x ∈ R ∪ {-∞}, x ∈ R ∪ {+∞} and x ≤ x. The set of intervals is denoted by IR.

In IBEX, three external implementations of the Interval Arithmetic can be linked: Filib++ [START_REF] Nehmeier | FILIB++, Expression Templates and the Coming Interval Standard[END_REF], Profil-Bias [START_REF] Knuppel | PROFIL/BIAS-A fast interval library[END_REF], Gaol [START_REF] Goualard | Gaol: NOT Just Another Interval Library[END_REF]. To improve the portability and the compatibility of IBEX, a homemade interval arithmetic is available without using low-level functionality which can be dependent on the architecture of the CPU. All the arithmetic has been patched to comply with the new IEEE 1788-2015 Standard for Interval Arithmetic [START_REF]IEEE Standard for Interval Arithmetic[END_REF].

Affine Arithmetic

The affine arithmetic is a technique to compute lower and upper bounds of functions over an interval. It is based on the same principle as interval arithmetic excepted that the quantities are represented by an affine form, see [START_REF] Stolfi | Self-validated numerical methods and applications[END_REF].

As in interval arithmetic, the usual operations and functions are extended to deal with affine forms. For the non-affine operations and some transcendental functions, such as the square root, the logarithm, the inverse and the exponential, several algorithms exist depending on the use: if you focus on the computation of the bounds, or on performance, or on linear approximation, or on the reliability, etc. Indeed, seven different versions are available to satisfy all the needs of the user [START_REF] Stolfi | Self-validated numerical methods and applications[END_REF].

Contractor Programming

The concept of contractor is directly inspired by the ubiquitous concept of filtering algorithms in constraint programming [START_REF] Chabert | Contractor programming[END_REF]. The strength of IBEX lies mainly in this concept. Every algorithm in IBEX is included as a Contractor.

Definition 1 Let X ⊆ R n be a feasible region. The operator C X : IR n → IR n is a contractor for X if: ∀[x] ∈ IR n , C X ([x]) ⊆ [x], (contraction) C X ([x]) ∩ X ⊇ [x] ∩ X. (completeness)
This definition means that: (i) Filtering gives a subdomain of the input domain [x]; (ii) the resulting subdomain C X ([x]) contains all feasible points. No solution is lost. A contractor is defined by a feasible region X, and its purpose is to eliminate a part of a domain which is not in X.

All set operators can be extended to contractors. For example, the intersection of two contractors creates a contractor for the intersection of these two sets. In the same way, the hull of two contractors creates a contractor for the disjunction of these constraints.

Definition 2 Let X and Y ⊆ R n be two feasible regions.

CtcCompo:

(

C X ∩ C Y )([x]) = C X ([x]) ∩ C Y ([x]) CtcUnion: (C X ∪ C Y )([x]) = C X ([x]) ∪ C Y ([x]) CtcFixPoint: C ∞ = C • C • C • . . .
Using these properties, interacting, combining and merging heterogeneous techniques becomes simple.

Lists of Contractors

In this section, a small part of all contractors available in IBEX is described:

CtcFwdBwd and CtcHC4: the atomic contractors Forward-backward is a classical algorithm in constraint programming for contracting quickly with respect to one equality or inequality constraint. See, e.g., [START_REF] Benhamou | Continuous and interval constraints[END_REF], [START_REF] Collavizza | A Note on Partial Consistencies over Continuous Domains[END_REF]. However, the more occurrences of variables in the expression of the (in)equality, the less accurate the contraction. Hence, this contractor is often used as an atomic contractor embedded in a higher-level operator like Propagation or Shaving.

HC4 is another classical algorithm of constraint programming. It allows to contract according to a system of constraints. The basic idea is to calculate the fix point of a set of

n contractors C 1 , . . . , C n , i.e. (C • • • • • C) ∞ ,
without calling a contractor when it is unnecessary.

Ctc3BCid and CtcAcid: the shaving contractors Ctc3BCID is a shaving operator. It is an implementation of the 3BCID algorithm defined in [START_REF] Trombettoni | Constructive interval disjunction[END_REF]. The shaving operator applies a contractor C on sub-parts (slices) of the input box. If a slice is entirely eliminated by C, the input box can be contracted by removing the slice from the box. This operator can be viewed as a generalization of the SAC algorithm in discrete domains [START_REF] Bessiere | Theoretical analysis of singleton arc consistency[END_REF]. This concept with continuous constraint was first introduced in [START_REF] Lhomme | Consistency techniques for numeric CSPs[END_REF] with the "3B" algorithm. In [START_REF] Lhomme | Consistency techniques for numeric CSPs[END_REF], the sub-contractor C was CtcHC4. In IBEX, the idea was extended and Ctc3BCID can be combined with every contractor.

CtcAcid is an adaptive version of the 3BCID contractor. The handled number of variables for which a shaving will be performed is adaptively tuned. The ACID algorithm alternates: (i) small tuning phases (during e.g. 50 nodes) where the shaving is called on a number of variables double of the last tuned value (all variables in the first tuning phase); statistics are computed in order to determine an average number of interesting calls, The number of variables to be handled in the next running phase is computed at the end of the tuning phase; (ii) and large running phases (during e.g. 950 nodes) where 3BCID is called with the number of variables determined during the last tuning phase.

CtcPolytopeHull: Contractors based on Linear relaxation Considering a system of linear inequalities, CtcPolytopeHull gives the possibility to contract a box to the hull of the polytope (the set of feasible points). This contractor calls a linear solver linked with IBEX (CPLEX, Soplex or CLP) to calculate for each variable x i , the following bounds: min

Ax≤b∧x∈[x]
x i and max

Ax≤b∧x∈[x]
x i , where [x] is the box to be contracted.

If some constraints are nonlinear, Linearization procedures can automatically linearize the non-linear constraints. There exists some built-in linearization techniques in IBEX : (i) LinearRelaxXTaylor: a corner-based Taylor relaxation [START_REF] Araya | A contractor based on convex interval Taylor[END_REF]; (ii) LinearRelaxAffine2: a relaxation based on affine arithmetic [START_REF] Ninin | A reliable affine relaxation method for global optimization[END_REF]; (iii) LinearRelaxCombo: a combination of the two previous techniques (the polytope is basically the intersection of the polytopes calculated by each technique).

CtcQInter: the q-relaxed intersection If a set of constraints is based on physical data, it is not uncommon that some of this data is wrong. In this situation, the q-relaxed intersection of contractors can be applied to this problem.

The q-relaxed intersection of m subsets X 1 , . . . , X m of R n is the set of all x ∈ R n which belong to at least (m -q) X i . We denote it by X {q} = {q} X i . Since the q-relaxed intersection is a set operator, we have extended this notion to contractors:

{q} C Xi ([x]) = {q} (C Xi ([x]
)). This contractor allows modeling the possibility of invalid constraints: it can also be used for robust optimization.

In [START_REF] Carbonnel | Qintersection Algorithms for Constraint-Based Robust Parameter Estimation[END_REF], Carbonnel et al. found an algorithm with a complexity θ(nm 2 ) to compute a box which contains the q-relaxed intersection of m boxes of R n .

CtcExist and CtcForAll: the contractors with quantifiers Another possibility is to project a subset of R n over one or more dimensions. For example, if a constraint needs to be satisfied for all values of a parameter in a given set, such as {x ∈ R n : ∀t ∈ X ⊆ R m , g(x, t) ≤ 0}, few solvers are available to deal with it. Another example is when a constraint needs to be satisfied for at least one value of the parameter, such as {x ∈ R n : ∃t ∈ X ⊆ R m , g(x, t) ≤ 0}.

Two operators are defined as contractors. The first one is CtcForAll and the second one is CtcExist. CtcForAll contracts each part of

Optimization Strategies

To find the global optimum of a problem in a reliable way, IBEX included several global optimization strategies. The principle of these algorithms is based on a branch-and-bound technique [START_REF] Ninin | Optimisation Globale basée sur l'Analyse d'Intervalles : Relaxation Affine et Limitation de la Mémoire[END_REF]. At each iteration, the domain under study is bisected to improve the computation of bounds. Boxes are eliminated if and only if it is certified that no point in the box can produce a better solution than the current best one, or that at least one constraint cannot be satisfied by any point in such a box. To accelerate convergence, contractors are used at each iteration to prune the width of boxes.

The default optimization strategy is based on a mathematical model of the optimization problem that needs to be solved. This model can be constructed directly using the symbolic kernel of IBEX or using the AMPL interface. The default contractor inside is the following: CtcAcid∩(CtcPolytopeHull ∩ CtcHC4)

∞ . The performance of the default optimizer is comparable to the global optimizer BARON. However, our approach is completely reliable, and can deal with more general problems (with trigonometric function).

Moreover, a general pattern is also available [START_REF] Ninin | Global optimization based on contractor programming[END_REF]. This pattern only requires a contractor defined on the feasible set X of the problem and another contractor on the unfeasible set X. Indeed, IBEX can address more complex real-life problems with disjunction constraint, quantifiers, outliers, non-linearity, trigonometric functions, etc.

  [x] which is contracted by C([x] × {y}) for any y ∈ Y. Indeed, each part [a] of [x], such as ∃y ∈ Y, ([a], y) / ∈ Z, can be removed. Thus, each part [b] of [x], such as ∀y ∈ Y, ([b], y) ∈ Z, is kept. A similar algorithm is used in CtcExist.

* ENSTA-Bretagne, Lab-STICC, IHSEV team, 2 rue Francois Verny, 29806 Brest, France. jordan.ninin@ensta-bretagne.fr


Dominique Monnet
email: dominique.monnet@ensta-bretagne.fr

Jordan Ninin
email: jordan.ninin@ensta-bretagne.fr

Benoit Clement
email: benoit.clement@ensta-bretagne.fr

Global Optimization of H ∞ problems: Application to robust control synthesis under structural constraints

In this paper, a new technique to compute a synthesis structured Robust Control Law is developed. This technique is based on global optimization methods using a Branch-and-Bound algorithm. The original problem is reformulated as a min/max problem with non-convex constraint. Our approach uses interval arithmetic to compute bounds and accelerate the convergence.

Context

Controlling an autonomous vehicle or a robot requires the synthesis of control laws for steering and guiding. To generate efficient control laws, a lot of specifications, constraints and requirements have been translated into norm constraints and then into a constraint feasibility problem. This problem has been solved, sometimes with relaxations, using numerical methods based on LMI (Linear Matrix Inequalities) or SDP (Semi Definite Program) [START_REF] Arzelier | Multi-objective H 2 /H ∞ /impulseto-peak control of a space launch vehicle[END_REF][START_REF] Boyd | Linear Matrix Inequalities in System and Control Theory[END_REF]. The main limitation of these approaches is the complexity of the controller for implementation in an embedded system. However, if a physical structure is imposed on the control law in order to make the implementation easier, the synthesis of this robust control law is much more complex. And this complexity has been identified as a key issue for several years. An efficient first approach based on local non-smooth optimization was given by Apkarian and Noll [START_REF] Apkarian | Nonsmooth H ∞ synthesis. Automatic Control[END_REF].

In this talk, we will present a new approach based on global optimization in order to generate robust control laws.

H ∞ control synthesis under structural constraints

We illustrate our approach with an example of the control of a periodic second order system G with a PID controller K subjected to two frequency constraints

W 1 W 2 K G - + constraint1 constraint2

Controler

Dynamic System ? ?

w z 1 z 2 y Figure 1: 2-blocks H ∞ problem
on the error e and on the command u of the closed-loop system, see Figure 1.

The objective is to find k = (k p , k i , k d ) to stabilize the closed-loop system while minimizing the H ∞ norm of the controlled system to ensure robustness. The H ∞ norm of a dynamic system P is defined as follows:

||P || ∞ = sup ω (σ max (P (jω))),
with σ max the greatest singular value of the transfer function P and j the imaginary unit.

In our particular case, the closed-loop system can be interpreted as two SISO systems (Single In Single Out). The H ∞ norm of a SISO system is the maximum of the absolute value of the transfer function. Indeed, to minimize the H ∞ norm of our example, we need to solve the following min/max problem:

       min k max sup ω W 1 (jω) 1 + G(jω)K(jω) , sup ω W 2 (jω)K(jω) 1 + G(jω)K(jω) , s.t.
The closed-loop system must be stable.

(

) 1
The stability constraint of a closed-loop is well-known: the roots of denominator part of the transfer function 1 1+G(s)K(s) must have a non-positive real part [START_REF] Petersen | Robust control of uncertain systems: Classical results and recent developments[END_REF]. Using Routh-Hurwitz stability criterion [START_REF] Hurwitz | Ueber die Bedingungen, unter welchen eine Gleichung nur Wurzeln mit negativen reellen Theilen besitzt[END_REF], this constraint can be reformulated as a set of non-convex constraints.

Proposition 1 Let us consider a polynomial Q(s) = a n s n + a n-1 s n-1 + • • • + a 1 s + a 0 .
The real parts of its roots are negative if the entries in the first column of the following table are positive:

v1,1 = an v1,2 = an-2 v1,3 = an-4 v1,4 = an-6 . . . v2,1 = an-1 v2,2 = an-3 v2,3 = an-5 v2,4 = an-7 . . . v3,1 = -1 v 2,1 v1,1 v1,2 v2,1 v2,2 v3,2 = -1 v 2,1 v1,1 v1,3 v2,1 v2,3 v3,3 = -1 v 2,1 v1,1 v1,4 v2,1 v2,4 . . . . . . v4,1 = -1 v 3,1 v2,1 v2,2 v3,1 v3,2 v4,2 = -1 v 3,1 v2,1 v2,3 v3,1 v3,3 . . . . . . . . . v5,1 = -1 v 4,1 v3
(s) = 1 + G(s)K(s)
, the H ∞ control synthesis under structural constraint is reformulated as a min/max problem with non-convex constraints.

Global optimization of min/max problems

In order to solve Problem (1), our approach is based on a Branch-and-Bound technique [START_REF] Ninin | Optimisation Globale basée sur l'Analyse d'Intervalles : Relaxation Affine et Limitation de la Mémoire[END_REF]. At each iteration, the domain under study is bisected to improve the computation of bounds. Boxes are eliminated if and only if it is certified that no point in the box can produce a better solution than the current best one, or that at least one constraint cannot be satisfied by any point in such a box.

The non-convex constraint can be handled with constraint programming techniques. In our approach, we use the ACID algorithm [START_REF] Trombettoni | Constructive interval disjunction[END_REF] which reduces the width of the boxes and so accelerates the convergence of the branch-and-bound.

The key point of our approach concerns the computation of the bounds of the objective function. In our example, the objective function can be reformulated as the following expression, with x = (k p , k i , k d ):

f (x) = sup ω∈[ωmin,ωmax] g(x, ω). ( 2
)
At each iteration, Algorithm 1 is used to compute a lower bound of this function over a box [x]. This algorithm is also a branch-and-bound algorithm based on Interval Arithmetic. But, for not wasting time, we limit the maximum number of iterations for computing faster lower bounds. Each element

([ω], ub ω ) stored in L is composed of: (i) [ω] a sub-interval of [ω min , ω max ] and (ii) ub ω an upper bound of g over [x] × [ω].

Algorithm 1 Computation of bounds of f over a box [x]

Require: g: the function under study (see Equation 2); x: an initial box; L: the list of boxes; nbIter: the maximal number of iterations. 1: Initialization: (lb out , ub out ) = (-∞, ∞). 2: for nb := 1 to nbIter do 3:

Extract an element (ω, ub ω ) from L.

4:

Bisect ω into two sub-boxes ω 1 and ω 2 .

5:

for i:=1 to 2 do 6:

Compute lb ωi and ub ωi a lower and an upper bound of g(x, ω) over [x] × [ω i ] using Interval Arithmetic techniques [START_REF] Jaulin | Applied Interval Analysis, with Examples in Parameter and State Estimation, Robust Control and Robotics[END_REF].

Application

In our example, we consider a second-order system and weighting functions W 1 and W 2 penalizing the error signal and control signal respectively:

G(s) = 1 s 2 + 1.4s + 1 , K(s) = k p + k i s + k d s 1 + s . W 1 (s) = s + 100 100s + 1 , W 2 (s) = 10s + 1 s + 10 .
We want to find k p , k i and k d the coefficients of the structured controller K such that the closed-loop system respects the constraints:

max || W 1 (jω) 1 + G(jω)K(jω) || ∞ , || W 2 (jω)K(jω) 1 + G(jω)K(jω) || ∞ ≤ 1
The control is bounded by [-2, 2], and we limit the interval of ω to [10 -2 , 10 2 ].

Our algorithm gives the following result:

max sup ω W 1 (jω) 1 + G(jω)K(jω) , sup ω W 2 (jω)K(jω) 1 + G(jω)K(jω) = 2.1414 with k p = -0.0425, k i = 0.4619, k d = 0.2566
Unfortunately, the value of the solution of the min/max problem is greater than 1. So, the constraints are not respected as shown in Figure 2 (solid lines are above dotted lines of the same color at some frequencies). In this example, the main advantage of our global optimization approach is that unlike classical method based on non-smooth optimization, we can certify that no robust solution of our problem exists.

Figure 2 :

 2 Figure 2: Weighting functions and singular values of the solution

7 :

 7 if lb ωi > lb out then 8:lb out := lb ωi , {Update the best lower bound} Remove from L all the elements j such as ub ωj < lb out , return (lb out , ub out ): a lower and an upper bound of f over x.Thanks to Interval Analysis, at the end of Algorithm 1, we can ensure that the value of the maximum of f over [x] is included in [lb out , ub out ].

	9:
	10:	end if
	11:	if ub ωi > lb out then
	12:	Add (ω i , ub ωi ) in L,
	13:	end if
	14:	end for
	15: end for
	16: ub out :=	max (ωi,ubω i )∈L	ub ωi
	17: