Matteo Ponzano
email: matteo.ponzano@umr-cnrm.fr

Bruno Joly

Laurent Descamps

Philippe Arbogast

Systematic error analysis of heavy-precipitation-event prediction using a 30-year hindcast dataset

The western Mediterranean region is prone to devastating flash floods induced by heavy-precipitation events (HPEs), which are responsible for considerable human and material losses. Quantitative precipitation forecasts have improved dramatically in recent years to produce realistic accumulated rainfall estimations. Nevertheless, there are still challenging issues which must be resolved to reduce uncertainties in the initial condition assimilation and the modelling of physical processes. In this study, we analyse the HPE forecasting ability of the multi-physics-based ensemble model Prévision d'Ensemble ARPEGE (PEARP) operational at Météo-France. The analysis is based on 30-year (1981-2010) ensemble hindcasts which implement the same 10 physical parameterizations, one per member, run every 4 d. Over the same period a 24 h precipitation dataset is used as the reference for the verification procedure. Furthermore, regional classification is performed in order to investigate the local variation in spatial properties and intensities of rainfall fields, with a particular focus on HPEs. As grid-point verification tends to be perturbed by the double penalty issue, we focus on rainfall spatial pattern verification thanks to the feature-based quality measure of structure, amplitude, and location (SAL) that is performed on the model forecast and reference rainfall fields. The length of the dataset allows us to subsample scores for very intense rainfall at a regional scale and still obtain a significant analysis, demonstrating that such a procedure is consistent to study model behaviour in HPE forecasting. In the case of PEARP, we show that the amplitude and structure of the rainfall patterns are basically driven by the deep-convection parametrization. Between the two main deep-convection schemes used in PEARP, we qualify that the Prognostic Condensates Microphysics and Transport (PCMT) parametrization scheme performs better than the B85 scheme. A further analysis of spatial features of the rainfall objects to which the SAL metric pertains shows the predominance of large objects in the verification measure. It is for the most extreme events that the model has the best representation of the distribution of object-integrated rain.

Introduction

Episodes of intense rainfall in the Mediterranean affect the climate of western Europe and can have an important societal impact. During these events, daily rainfall amounts associated with a single event can reach annual equivalent values. These rainfall events coupled with a steep orography are responsible for associated torrential floods, which may cause considerable human and material losses. In particular, southern France is prone to devastating flash-flood events such as those at Aude [START_REF] Ducrocq | The extreme flash flood case of November 1999 over Southern France[END_REF], Gard [START_REF] Delrieu | The Catastrophic Flash-Flood Event of 8-9 September 2002 in the Gard Region, France: A First Case Study for the Cévennes-Vivarais Mediterranean Hydrometeorological Observatory[END_REF], and Vaison-La-Romaine [START_REF] Sénési | The Vaison-La-Romaine Flash Flood: Mesoscale Analysis and Predictability Issues[END_REF], which occurred on 12-13 November 1999, 22 September 1992, and 8-9 September 2002, respectively. For instance, in the Gard case more than 600 mm was observed locally during a 2 d, event and 24 people were killed during the associated flash flooding. Extreme rainfall events generally occur in a synoptic environment favourable for such events [START_REF] Nuissier | A Statistical Downscaling to Identify the Large-Scale Circula-tion Patterns Associated with Heavy Precipitation Events over Southern France[END_REF].

A detailed list of the main atmospheric factors which contribute to the onset of HPEs are reported by [START_REF] Lin | Some Common Ingredients for Heavy Orographic Rainfall[END_REF]:

(1) a conditionally or potentially unstable airstream impinging on the mountains, (2) a very moist low-level jet, (3) a steep mountain, and (4) a quasi-stationary convective system that persists over the threat area. However, not all these factors necessarily need to be present at the same time to

Published by Copernicus Publications on behalf of the European Geosciences Union.

produce HPEs. In southeastern France, the Mediterranean Sea acts as a source of energy and moisture which is fed to the atmospheric lower levels over a wide pronounced orography above the Massif Central, Pyrenees, and southern Alps [START_REF] Delrieu | The Catastrophic Flash-Flood Event of 8-9 September 2002 in the Gard Region, France: A First Case Study for the Cévennes-Vivarais Mediterranean Hydrometeorological Observatory[END_REF]. Extreme rainfall amounts are enhanced, especially along the southern and eastern foothills of mountainous chains [START_REF] Frei | A Precipitation Climatology of the Alps from High-Resolution Rain-Gauge Observations[END_REF][START_REF] Nuissier | A Numerical Study of Three Catastrophic Precipitating Events over Southern France. I: Numerical Framework and Synoptic Ingredients[END_REF], in particular the southeastern part of the Massif Central (Cévennes). [START_REF] Ehmele | The influence of Sardinia on Corsican rainfall in the western Mediterranean Sea: A numerical sensitivity study[END_REF] emphasized the important role played by complex orography, the mutual interaction between two close mountainous islands in this case, in heavy rainfall under strong synoptic forcing conditions. Nevertheless, other regions are also affected by rainfall events with a great variety of intensity and spatial extension. [START_REF] Ricard | A Climatology of the Mesoscale Environment Associated with Heavily Precipitating Events over a Northwestern Mediterranean Area[END_REF] studied this regional spatial distribution based on a composite analysis and showed the existence of mesoscale environments associated with heavy-precipitation events. Considering four subdomains, they found that the synoptic and mesoscale patterns can greatly differ as a function of the location of the precipitation.

Extreme rainfall events are generally associated with coherent structures slowed down and enhanced by the relief, whose extension is often larger than a single thunderstorm cell. At some point, this mesoscale organization can turn into a self-organization process, leading to a mesoscale convective system (MCS) when interacting with its environment, which in turn leads to high-intensity rainfall [START_REF] Nuissier | A Numerical Study of Three Catastrophic Precipitating Events over Southern France. I: Numerical Framework and Synoptic Ingredients[END_REF].

Among the list of factors contributing to HPE creation, some are clearly only within the scope of high-resolution convection-permitting models. Indeed, vertical motion and moisture processes need to be explicitly solved to get realistic representation of convection. On the other hand, as we have just highlighted, some other factors linked with synoptic circulations or orography representations can be well estimated in global models, in particular when horizontal resolution gets close to 15-20 km. Consequently, the corresponding predictability of such factors can reach advantageous lead times for early warnings, i.e. longer than the standard 48 h that the limited area model may be expected to achieve. Indeed, if long-term territorial adaptations are necessary to mitigate the impact of HPEs, a more reliable and earlier alert would be beneficial in the short term. Weather forecasting coupled with hydrological impact forecasting is the main source of information for triggering of weather warnings. Severe weather warnings are issued for the 24 h forecast only. However, in some cases, the forecast process could be issued some days prior to the severe weather warnings. A better understanding of the sources of model uncertainty for such a time range may provide a major source of improvement for early diagnosis.

Forecast uncertainties can be related to initialization data (analysis) or lateral boundary conditions, and they have been investigated with both deterministic models [START_REF] Argence | Impact of Initial Condition Uncertainties on the Predictability of Heavy Rainfall in the Mediterranean: A Case Study[END_REF] and ensemble models [START_REF] Vié | Cloud-Resolving Ensemble Simulations of Mediterranean Heavy Precipitating Events: Uncertainty on Initial Conditions and Lateral Boundary Conditions[END_REF]. Several pre-vious studies showed that predictability associated with intense rainfall and flash floods decreases rapidly with the event scale (Walser et al., 2004;Walser and Schär, 2004;[START_REF] Collier | Flash Flood Forecasting: What Are the Limits of Predictability?[END_REF]. Several studies based on ensemble prediction systems have shown the general ability of such models to sample the sources of uncertainty in HPE probabilistic forecasting [START_REF] Du | Short-Range Ensemble Forecasting of Quantitative Precipitation[END_REF][START_REF] Petroliagis | Potential Use of the ECMWF Ensemble Prediction System in Cases of Extreme Weather Events[END_REF][START_REF] Stensrud | Using Ensembles for Short-Range Forecasting[END_REF][START_REF] Schumacher | Ensemble-Based Forecast Uncertainty Analysis of Diverse Heavy Rainfall Events[END_REF]World Meteorological Organization, 2012). In ensemble forecasting, the uncertainty associated with the forecast is usually assessed by taking into account initial and model error propagation. As for the initial uncertainty, major meteorological centres implement different methods, the most common of which are singular vectors [START_REF] Buizza | The Singular-Vector Structure of the Atmospheric Global Circulation[END_REF][START_REF] Molteni | The ECMWF Ensemble Prediction System: Methodology and Validation[END_REF], bred vectors (Toth andKalnay, 1993, 1997), and perturbed observation in the analysis process [START_REF] Houtekamer | A System Simulation Approach to Ensemble Prediction[END_REF][START_REF] Houtekamer | Data Assimilation Using an Ensemble Kalman Filter Technique[END_REF]. The model error is related to grid-scale unsolved processes in the parametrization scheme and is assessed in the models with two main techniques. Some models use stochastic perturbations of the inner-model physics scheme (Palmer et al., 2009); others use different parametrization schemes in each forecast member [START_REF] Charron | Toward Random Sampling of Model Error in the Canadian Ensemble Prediction System[END_REF][START_REF] Descamps | Representing model uncertainty using the multiparametrization method[END_REF].

The global ensemble model Prévision d'Ensemble ARPEGE (PEARP; [START_REF] Boisserie | Calibrated Forecasts of Extreme Windstorms Using the Extreme Forecast Index (EFI) and Shift of Tails (SOT)[END_REF] implemented at Météo-France is based on the second technique, also known as a multi-physics approach. Compared to the stochastic perturbation, the error model distribution cannot be explicitly formulated in the multi-physics approach. It is then difficult to know a priori the influence of the physics scheme modifications on the forecast ability of the model. This is even more the case when highly non-linear physics with processes of a high order of magnitude are considered. In order to improve the understanding and interpretation of ensemble forecasts in tense decision-making situations as well as for model development and improvement purposes, it would be of great interest to have a full and objective analysis of the model behaviour in terms of HPE forecasting. This is one of the main aims of this study.

In order to achieve such a systematic analysis, standard rainfall verification methods can be used. They are usually based on grid-point-based approaches. These techniques, especially when applied to intense events, are subject to time or position errors, leading to low scores [START_REF] Mass | Does Increasing Horizontal Resolution Produce More Skillful Forecasts?[END_REF], also known as the double-penalty problem [START_REF] Rossa | Overview of Methods for the Verification of Quantitative Precipitation Forecasts[END_REF]. To counteract this problem, spatial verification techniques have been developed with the goal of evaluating forecast quality from a forecaster standpoint. Some of these techniques are based on object-oriented verification methods [START_REF] Ebert | Verification of Precipitation in Weather Systems: Determination of Systematic Errors[END_REF]Davis et al., 2006a;[START_REF] Wernli | SAL-A Novel Quality Measure for the Verification of Quantitative Precipitation Forecasts[END_REF][START_REF] Davis | The Method for Object-Based Diagnostic Evaluation (MODE) Applied to Numerical Forecasts from the 2005 NSSL/SPC Spring Program[END_REF]AghaKouchak et al., 2011;[START_REF] Mittermaier | Feature-Based Diagnostic Evaluation of Global NWP Forecasts[END_REF]. The feature-based quality measure SAL [START_REF] Wernli | SAL-A Novel Quality Measure for the Verification of Quantitative Precipitation Forecasts[END_REF][START_REF] Wernli | Spatial Forecast Verification Methods Intercomparison Project: Application of the SAL Technique[END_REF] is used in this study. Another element required to achieve such an analysis is the availability of fore-cast datasets long enough to get a proper sampling of the events to verify.

In our study, we profit from a reforecast dataset based on a simplified version of the PEARP model available over a 30year period. Such reforecast datasets have been previously shown to be relevant for calibrating operational models in various ways. In [START_REF] Hamill | Probabilistic Quantitative Precipitation Forecasts Based on Reforecast Analogs: Theory and Application[END_REF], [START_REF] Hamill | Probabilistic Forecast Calibration Using ECMWF and GFS Ensemble Reforecasts[END_REF], [START_REF] Hamill | Verification of TIGGE Multimodel and ECMWF Reforecast-Calibrated Probabilistic Precipitation Forecasts over the Contiguous United States[END_REF][START_REF] Boisserie | Calibrated Forecasts of Extreme Windstorms Using the Extreme Forecast Index (EFI) and Shift of Tails (SOT)[END_REF], the reforecast is used as a learning dataset to fit statistical models to calibrate forecast error corrections that are then applied to operational forecasting outputs. [START_REF] Boisserie | Calibrated Forecasts of Extreme Windstorms Using the Extreme Forecast Index (EFI) and Shift of Tails (SOT)[END_REF] and [START_REF] Lalaurette | Early detection of abnormal weather conditions using a probabilistic extreme forecast index[END_REF] have shown the possibility of using a reforecast dataset as a statistical reference of the model to which the extremeness of a given forecast is compared. In this paper, we analyse the ensemble model PEARP forecast predictability at lead times between day 2 and day 4 of daily rainfall amounts. This analysis is performed on the long reforecast 30-year dataset. One aim is to determine whether a multi-physics approach could be considered as a model error sampling technique appropriate for a good representation of HPEs in the forecast at such lead times. In particular, the behaviour of the different physics schemes implemented in PEARP has to be estimated individually. One main side aspect of this work focuses on developing a methodology suitable for evaluating the performances of an ensemble reforecast in the context of intense precipitation events using an object-oriented approach. In particular, we focus on intense precipitation over the French Mediterranean region. In addition to the analysis of diagnostics from the SAL metric, a statistical analysis of 24 h rainfall objects identified in the forecasts and the observations is performed in order to explore the spatial properties of the rainfall fields.

The data and the methodology are presented in Sect. 2. Section 2.1 describes the reforecast ensemble dataset, and Sect. 2.2 details the creation of the daily rainfall reference, the HPE statistical definition, and the regional clustering analysis. Results arising from the spatial verification of the overall reforecast dataset are presented in Sect. 3.1. Section 3.2 presents SAL diagnostics divided into all different physical parametrization schemes of the ensemble reforecast and for the spatial properties of individual objects. Conclusions are given in Sect. 4.

Data and methodology

PEARP hindcast

The PEARP reforecast dataset consists of a 10-member ensemble computed daily from 18:00 UTC initial conditions, covering 4 months (from September to December), every year of a 30-year period . This period has been chosen since HPE occurrence in the considered region is largest during the autumn season (see Fig. 3 from [START_REF] Ricard | A Climatology of the Mesoscale Environment Associated with Heavily Precipitating Events over a Northwestern Mediterranean Area[END_REF]. It uses ARPEGE (Action de Recherche Petite Echelle Grande Echelle; [START_REF] Courtier | The ARPEGE project at Météo-France[END_REF], the global operational model of Météo-France with a spectral truncation T798, 90 levels in the vertical, and a variable horizontal resolution (mapping factor of 2.4 with a highest resolution of 10 km over France). One ensemble forecast is performed every 4 d of the 4-month period up to 108 h lead time. Our initialization strategy follows the hybrid approach described in Boisserie et al. (2016), in which first the atmospheric initial conditions are extracted from the ERA-Interim reanalysis [START_REF] Dee | The ERA-Interim Reanalysis: Configuration and Performance of the Data Assimilation System[END_REF] available at the European Centre for Medium-Range Weather Forecasts. Second, the land-surface initialization parameters are interpolated from an offline simulation of the land-surface SURFEX model (Masson et al., 2013) driven by the 3-hourly near-surface atmospheric fields from ERA-Interim. The 24 h accumulated precipitation forecasts are extracted on a 0.1 • × 0.1 • grid that defines domain D (see Fig. 1c), which encompasses southeastern France (Fig. 1a). The reforecast dataset does not have any representation of initial uncertainty, but it implements the same representation of model uncertainties (multi-physics approach) as in the PEARP operational version of 2016.

Nine different physical parameterizations (see Table 1) are added to the one that corresponds to the ARPEGE deterministic physical package. This set of parameterizations is the same as the one implemented in PEARP. Two turbulent diffusion schemes are considered: the turbulent kinetic energy scheme (TKE; [START_REF] Cuxart | A turbulence scheme allowing for mesoscale and large-eddy simulations[END_REF][START_REF] Bazile | The Turbulent Kinetic Energy (TKE) scheme in the NWP models at Meteo France[END_REF] and the Louis scheme (L79; Louis, 1979). TKE mod is a slightly modified version of TKE, in which horizontal advection is ignored. For shallow convection, different schemes are used: a mass flux scheme introduced by [START_REF] Kain | Convective Parameterization for Mesoscale Models: The Kain-Fritsch Scheme[END_REF] and modified by [START_REF] Bechtold | A Mass-Flux Convection Scheme for Regional and Global Models[END_REF], hereafter called the KFB approach; the Prognostic Condensates Microphysics and Transport scheme (PCMT; [START_REF] Piriou | An Approach for Convective Parameterization with Memory: Separating Microphysics and Transport in Grid-Scale Equations[END_REF]; the eddy diffusivity and Kain-Fritsch scheme (EDKF); and the PMMC (Pergaud, Masson, Malardel, Couvreux) scheme (Pergaud et al., 2009). The deep-convection component is parameterized by either the PCMT scheme or the [START_REF] Bougeault | A Simple Parameterization of the Large-Scale Effects of Cumulus Convection[END_REF] scheme (hereafter B85). Closing the equation system used in these two schemes means relating the bulk mass flux to the in-cloud vertical velocity through a quantity γ qualifying the convection area coverage. Two closures are considered: the first one (C1) is based on the convergence of humidity, and the second one (C2) is based on the CAPE (convective available potential energy). The B85 scheme originally uses the C1 closure, while PCMT alternatively uses the closure (C1 or C2) which maximizes the γ parameter. Physics package 2 uses a modified version of the B85 scheme in which deep convection is triggered only if cloud top exceeds 3000 m (B85 mod in Table 1). The same trigger is used in physics package 3 in which deep convection is parameterized using the B85 scheme along with a CAPE closure (CAPE in Table 1). Finally the oceanic flux is solved by means of the ECUME (Exchange Coefficients from Unified Multicampaigns Estimates) scheme [START_REF] Belamari | Report on uncertainty estimates of an optimal bulk formulation for surface turbulent fluxes[END_REF]. In ECUME mod evaporation fluxes above sea surfaces are enhanced. Control member and member 9 are characterized by the same parametrization set-up, but member 9 differs for the modelling of orographic waves.

Daily rainfall reference

The 24 h accumulated precipitation is derived from the in situ Météo-France rain-gauge network, covering the same period as the reforecast dataset. The 24 h rainfall amounts collected from 14 French departments within the reforecast domain D are used (Fig. 1b). In order to maximize the rain-gauge network density within the region, all daily available validated data covering the period have been used. Rain-gauge observations are used to build gridded precipitation references by a statistical spatial interpolation of the observations. The aim of this procedure is to ensure a spatial and temporal homogeneity of the reference, as well as the same spatial resolution as the reforecast dataset. [START_REF] Ly | Different methods for spatial interpolation of rainfall data for operational hydrology and hydrological modeling at watershed scale: a review[END_REF] provided a review of the different methods for spatial interpolation of rainfall data. They showed that kriging methods outperform deterministic methods for the computation of daily precipitation. However, both types of methods were found to be comparable in terms of hydrological modelling results. For the interpolation, we use a mixed geo-statistical and deterministic algorithm, which implements ordinary kriging (OK; [START_REF] Goovaerts | Geostatistics for natural resources evaluation[END_REF] and inversedistance-weighting methods (IDW; [START_REF] Shepard | A Two-Dimensional Interpolation Function for Irregularly-Spaced Data[END_REF]. For the kriging method, three semi-variogram models (exponential, Gaussian, and spherical) are fitted to a daily sample semivariogram drawn from raw and square-root-transformed data [START_REF] Gregoire | Regression Estimation Following the Square-Root Transformation of the Response[END_REF][START_REF] Erdin | Data Transformation and Uncertainty in Geostatistical Combination of Radar and Rain Gauges[END_REF]. This configuration involves the use of six different geo-statistical interpolation models. In addition, four different IDW versions are used, by varying the geometric form parameter d used for the estimation of the weights (see Eq. ( 2) in Ly et al., 2011) and the maximum number n of neighbour stations involved in the IDW computation. Three versions are defined by fixing parameter d = 2 and alternatively assigning n values equal to 5, 10, and N (with N being the total number of stations available for that specific day). In the fourth version we set n = N and d = 3. For each day, a different interpolation method is used, and its selection is based on the application of a crossvalidation approach. We select 55 rain gauges as a training dataset (see the red diamonds in Fig. 1c) in order to have sufficient coverage over the domain, especially on the mountainous area. Root-mean-square error (RMSE) is used as a criterion of evaluation. For each day, the method which minimizes the RMSE computed within the rain gauges of the training dataset is selected and the spatial interpolation is then performed on a regular high-resolution grid of 0.05 • . The highest-resolution estimated points are then upscaled to the 0.1 • grid resolution of domain D, by means of a spatial average. This upscaling procedure aims at reproducing the filtering effect produced by the parameterizations of the model on the physical processes that occur below the grid resolution.

HPE database

We implement a methodology in order to select the HPEs from the daily rainfall reference. Anagnostopoulou and Tolika (2012) have examined parametric and non-parametric approaches for the selection of rare events sampled from a dataset. Here we adopt a non-parametric peak-over-threshold approach, on the basis of World Meteorological Organization (WMO) guidelines (World Meteorological Organization, 2016). The aim is to generate a set of events representative of the tail of the rainfall distribution for a given region and season. Following the recommendation of Schär et al. (2016), an all-day percentile (P 0≤n≤1 ) formulation is applied. A potential weakness of the research methodology based on the gridded observation reference is that a few extreme pre- cipitation events affecting a smaller area than the grid resolution may not be identified. However, this approach has been preferred to a classification using rain gauges because spatial and temporal homogeneity are ensured.

We proceed as follows: first the domain is split into two subregions based on the occurrence of climatological intense precipitation during the 30-year period. The subregion A includes all the points whose climatological percentile 99.5 is lower than or equal to a threshold T , and subregion B includes all the other points. Threshold T , after several tests, has been set to 85 mm. This choice was made in order to separate the domain into two regions characterized by different frequency and intensity of HPEs. Subregion A designates a geographical area where a large number of cases of intense precipitation are observed. Subregion B primarily covers the plain area, where HPE frequency is lower. For this reason, two different level threshold values are selected to define an event, depending on the subregion. More specifically, a day is classified as an HPE if one point of subregion A accumulated rainfall is greater than 100 mm or if one point of subregion B rainfall is greater than its percentile 99.5. The selection led to a classification of 192 HPEs, corresponding to a climatological frequency of 5 % over the 30-year period. The 24 h rainfall amount maxima within the HPE dataset range from 100 to 504 mm. It is worth mentioning that since we consider daily rainfall, rainfall events that would have high 48 h or 72 h accumulated rainfall may be disregarded. Figure 2 shows for each point of the domain the number of HPEs as well as the composite analysis of HPEs. The composite analysis involves computing the grid-point average from a collection of cases. The signal is enhanced along the Cévennes chain and in the Alpine region. It should be noted that some points are never taken into account for the HPE selection (white points of Fig. 2) because the required conditions have not been met. The analysis of the rainfall fields across the , percentage of HPEs, and fraction of HPEs. HPEs (%) refers to the ratio between the number of HPEs within the cluster and the total number of HPEs. Fraction of HPEs (%) refers to the ratio between the number of HPEs within the cluster and the total number of dates included in the corresponding cluster.

Cluster

Total (%) HPEs (%) Fraction of HPEs (%) HPE database exhibits the presence of patterns of different shape and size, revealing potential differences in terms of the associated synoptic and mesoscale phenomena (not shown).

1

Clustering analysis

Clustering analysis methods can be applied to daily rainfall amounts in order to identify emergent regional rainfall patterns. This classification is largely used for assessing the between-day spatial classification of heavy rainfall [START_REF] Romero | Daily Rainfall Patterns in the Spanish Mediterranean Area: An Objective Classification[END_REF][START_REF] Peñarrocha | Classification of Daily Rainfall Patterns in a Mediterranean Area with Extreme Intensity Levels: The Valencia Region[END_REF]Little et al., 2008;[START_REF] Kai | A Spatial Cluster Analysis of Heavy Rains in China[END_REF]. We applied a cluster analysis, as an exploratory data analysis tool, in order to assess geographical properties of the precipitation reference dataset. The size of the dataset is first reduced, and the signal is filtered out by means of a principal component analysis [START_REF] Morin | Use of Principal Component Analysis to Identify Homogeneous Precipitation Stations for Optimal Interpolation[END_REF][START_REF] Mills | Principal Component Analysis of Precipitation and Rainfall Regionalization in Spain[END_REF][START_REF] Teo | Principal Component Analysis of Observed and Modeled Diurnal Rainfall in the Maritime Continent[END_REF]. The first 13 principal components (PCs), whose projection explains 90 % of the variance, are retained. Then the K-means clustering method is applied. It is a non-hierarchical method based on the minimization of the intraclass variance and the maximization of the variance between each cluster. A characteristic of k-means method is that the number of clusters (K) into which the data will be grouped has to be a priori prescribed. Consequently, we first have to implement a methodology to find the number of clusters which leads to the most classifiable subsets. The analysis is applied to the full reference dataset, including rainy and dry days. We run 2000 tests for a range of a priori cluster numbers K that lies between 3 and 13, by varying a random initial guess each time. Then, for a given K, an evaluation of the stability of the assignment into each cluster is performed. The number of clusters is considered stable if each cluster size is almost constant from one test to another. K = 5 is retained as the most stable number of clusters and because it suggests a coherent regional stratification of the daily rainfall data. The final classification within the 2000 tests is selected by minimizing the sum of the dis-tance between the cluster centroids from each test and the geometric medians of cluster centroids computed from all the tests. The test which minimizes this quantity has been selected as the reference classification. The results from the cluster classification are summarized in Table 2. The clusterization shows large differences in terms of cluster size: more than three-fourths of the dataset is grouped in cluster 4, which mostly collects the days characterized by weak precipitation amounts or dry days. The percentage of HPEs within the clusters shows that the most intense events are represented in clusters 2, 3, and 5, among which cluster 5 shows the largest proportion of HPEs (65 % of HPEs within this cluster). Clusters 2, 3, and 5 together account for 86 % of the HPEs.

The same composite analysis as the one previously applied to HPE class is now computed for each cluster class (Fig. 3). It shows significant differences between clusters. The relative intensity of events and the location are different for each of the clusters. Rainfall range is weak for cluster 1 and close to zero for cluster 4. Cluster 2 includes some moderate 24 h rainfall amounts related to generalized precipitation events and a few HPEs. For cluster 1, composite values are slightly higher in the northwestern area of the domain, while for cluster 2 rainfall amount values are more significant on the eastern side of domain D. Clusters 3 and 5 together account for 63 % of the HPEs of the whole period, but rainfall events seem to affect different areas. Cluster 3 includes most of the events impacting the Cévennes mountains and the eastern departments on the southern side of the Alps. Cluster 5 average rainfall is enhanced along the southern side of the Cévennes, especially the Languedoc-Roussillon region.

The bottom right panel of Fig. 3 shows the density distributions computed from the maximum daily rainfall for each cluster. It is worth noting that cluster rainfall distributions cover different intervals of maximum daily rainfall amounts. Cluster 4 includes all the dry days. As this paper focuses on the most severe precipitation events, results will only be shown for clusters 2, 3, and 5 for the remainder of the paper.

2.3

The SAL verification score

The SAL score definition

The SAL score is an object-based quality measure introduced by [START_REF] Wernli | SAL-A Novel Quality Measure for the Verification of Quantitative Precipitation Forecasts[END_REF] for the spatial verification of numerical weather prediction (NWP). It consists in computing three different components: structure S is a measure of volume and shape of the precipitation patterns, amplitude A is the normalized difference of the domain-averaged precipitation fields, and location L is the spatial displacement of patterns on the forecast/observation domains.

Different criteria for the identification of the precipitation objects could be implemented: a threshold level [START_REF] Wernli | SAL-A Novel Quality Measure for the Verification of Quantitative Precipitation Forecasts[END_REF][START_REF] Wernli | Spatial Forecast Verification Methods Intercomparison Project: Application of the SAL Technique[END_REF], a convolution threshold (Davis et al., 2006a, b), or a threshold level conditioned to a cohesive min- imum number of contiguous connected points [START_REF] Nachamkin | Application of the Composite Method to the Spatial Forecast Verification Methods Intercomparison Dataset[END_REF][START_REF] Lack | An Object-Oriented Multiscale Verification Scheme[END_REF]. The threshold level approach needs only one estimation parameter, so it has been preferred to the other methods for its simplicity and interpretability. Since we focus on the patterns associated with the HPEs, we decided to adapt the threshold definition given by T f = x max ×f , where x max is the maximum precipitation value of the points belonging to the domain and f is a constant factor (= 1/15, in the paper of [START_REF] Wernli | SAL-A Novel Quality Measure for the Verification of Quantitative Precipitation Forecasts[END_REF]. Here the coefficient f has been raised to one-fourth because a smaller value results in excessively large objects spreading out over most of domain D. Choosing a higher f factor enables us to obtain more realistic features within the domain considered. Threshold levels T f are computed daily for the reforecast and the reference dataset. Although objects are smaller than the domain for most of the situations, a few objects extending outside the domain are consequently limited by the boundaries of the region concerned.

If we consider domain D, the amplitude A is computed as follows:

A = R for D -R obs D 0.5 ( R for D + R obs D ) ∈ [-2, 2] , (1)
where D denotes the average over domain D. R for and R obs are the 24 h rainfall amounts over D associated with the forecast and the observation, respectively. A perfect score is achieved for A = 0. The domain-averaged rainfall field is overestimated by a factor of 3 if A = 1. Similarly, it is underestimated by a factor of 3 if A = -1. The amplitude is maximal

(A = 2) if R for D R obs D → +∞ and minimal (A = -2) if R for D R obs D
→ 0. The two other components require the definition of precipitation objects (thereafter {Obj}), also called features, which represent contiguous grid points belonging to domain D, characterized by rainfall values exceeding a given threshold. The location L is a combined score defined by the sum of two contributions, L1 and L2. L1 measures the magnitude of the shift between the centre of mass of the whole precipitation field for the forecast ( xfor ) and observation ( xobs ):

L1 = | xfor -xobs | d ∈ [0, 1] , ( 2
)
where d is the largest distance between two boundary points of the considered domain D. The second metric L2 takes into account the spatial distribution of the features inside the dohttps://doi.org/10.5194/nhess-20-1369-2020 Nat. Hazards Earth Syst. Sci., 20, 1369-1389, 2020 main, that is, the scattering of the objects:

r = N n=1 M n | x -x n | N n=1 M n , (3)
where M n is the integrated mass of the object n, x n is the centre of mass of the object n, N is the number of objects, and x is the centre of mass of the whole field.

L2 = 2 |r for -r obs | d ∈ [0, 1] (4)
L = L1 + L2 ∈ [0, 2] (5)
L2 aims at depicting object differences between observed and forecasted scattering of the precipitation objects. We can notice that the scattering variable (Eq. 3) is computed as the weighted distance between the centre of total mass and the centre of mass of each object. Therefore L is a combination of the information provided by the global spatial distribution of the fields (L1) and the difference in scattering of the features over the domain (L2). The location score is perfect if L1 = L2 = 0, so if L = 0 all the centres of mass match each other.

The S component is based on the computation of the integrated mass M k of one object k, scaled by the maximum rainfall amount of the object k:

V k = M k maxR(x; x ∈ Obj k ) . (6)
Then, the weighted average V of all features is computed, in order to obtain a scaled, weighted total mass:

V = N n=1 M n V n N n=1 M n , (7)
S = V for -V obs 0.5(V for + V obs ) ∈ [-2, 2] . (8)
Then, S represents the difference of both forecasted and observed volumes, scaled by their half-sum. It is important to scale the volume so that the structure is less sensitive to the mass, meaning that it relates more to the shape and extension of the features rather than their intensities. In particular S < 0 means that the forecast objects are large and/or flat compared to the observations. Inversely, peaked and/or smaller objects in the forecast give positive values of S. We refer to [START_REF] Wernli | SAL-A Novel Quality Measure for the Verification of Quantitative Precipitation Forecasts[END_REF] for the exploration of the behaviour of SAL for some idealized examples.

On the basis of the definition of the score, it can be noticed that A and L1 components are not affected by the object identification and depend only on the total rainfall fields.

A selected example of the application of SAL

An example of the SAL score applied to an HPE, which occurred on 28 October 2004, is shown in Fig. 4 (60 h lead time forecast run using the physical package n.8). For the rainfall reference, a 24 h rainfall maximum value (121.3 mm) was registered in the southeastern coastal region. Therefore the threshold level T f is set to 30.3 mm. For the forecast, the maximum value is 123.1 mm (T f = 30.8 mm), and, in contrast with the reference, it is located in Cévennes. The number of objects, three, is equivalent in both fields. The value of A is 0.08, which means that the domain-averaged precipitation field of the forecast is nearly similar to the reference one.

The structure S component is positive (0.28), which could be explained by the larger forecast object over the Cévennes area, while the object along the southeastern coast is smaller and less intense. The contribution of the third object is negligible for the computation of S. The L component is equal to 0.23, with L1 = 0.13 and L2 = 0.10. The location error L1 means that the distance between the centres of total mass (see diamonds in Fig. 4) is 13/100 of the largest distance between two boundary points of the considered domain. This error is mostly due to the fact that the most intense rainfall patterns are far apart from each other in the observations and the forecast.

Analysis of the reforecast HPE representation

An SAL verification score has been applied to the reforecast dataset to perform statistical analysis of QPF (quantitative precipitation forecast) errors. The reforecast dataset is considered to be a test-bed model in order to study sources of systematic errors in the forecast. The overall reforecast performance is first examined for HPEs and non-HPEs, then according to the clusters. In a second step, the behaviour of the different physics schemes is analysed by separately considering the SAL results of each reforecast member. Similarly, the analysis is again allocated to HPEs and non-HPEs and subsequently to each cluster. For both the reforecast and the reference, we set all the days with at least one grid point beyond 0.1 mm as a rainy day. In order to facilitate the comparison between the parameterizations, SAL verification is only performed when all the members and the reference are classified as rainy days. Table 3 shows the contingency table of the rainy and dry days. Therefore 84 false alarms, 226 missed cases, and 62 correctly forecast dry days are not involved in the SAL analysis. No HPEs belong to the misses, and no simulated HPEs belong to the false alarms. The SAL measure is then applied to the 3258 rainy days.

SAL evaluation of the HPE forecast

HPEs and non-HPEs

First the relationship between the A component of SAL and the maximum grid-point error is investigated (Fig. 5). The 36 and 60 h lead times (LT12 hereafter) and 84 and 108 h lead times (LT34 hereafter) are grouped together. Maximum daily absolute errors range between -250 and 250 mm. Rare higher values are observed, which are likely related to strong double-penalty effects that often occur in grid-pointto-grid-point verification. Points are mostly scattered along the amplitude axis, showing that the error dependence on the A component is weak. Concerning HPEs, the scatter plot shows A component values under 1, which means that the scaled average precipitation in the forecast never exceeds 3 times the observation. In contrast, A-component negative values are predominant, in particular at LT34, in relation to strong underestimations of the domain-averaged rainfall field. Some cases of significant maximum grid-point errors in conjunction with a moderate negative A component must be related to strong location errors. In these cases, the domainaveraged field may be similar to the observed one while the maximum rainfall is spatially deviated. For the non-HPEs, we can see that, especially for LT34, the model could significantly overestimate both the A component and the maximum grid-point error.

The relationship between the different SAL components might help to understand sources of model error. In Fig. 6 the S and A components are drawn for the HPEs only. Perfect scores are reached for the points located on the origin O of the diagram. Very few points are located in the top left-hand quadrant. This indicates that an overestimation of precipitation amplitude associated with too small rainfall objects is rarely observed. The points, especially for LT34, are globally oriented from the bottom left-hand corner to the top righthand corner. This suggests a linear growth of the A com-ponent as a function of the S component, which means that the average rainfall amount is roughly related to the structure of the spatial extension. For the two diagrams, it can also be noticed that many of the points are situated in the lower right quadrant, suggesting the presence of too large and/or flat rainfall objects compared to the reference while the corresponding A component is negative. This is supported by the values of the medians of the distribution of the two components (dashed lines) and the quartile values (respective limits of the boxes). The positive bias in the S component is even stronger for the most extreme HPEs (red triangles). The distortion of S-component error compared to the A component shows that the model has more difficulties reproducing the complex spatial structure than simulating the average volume of a heavy rainfall. This deficiency may be related to the convection part not represented in the parametrization scheme. It may also be related to the representation of orography at a coarse resolution. As shown by [START_REF] Ehmele | The influence of Sardinia on Corsican rainfall in the western Mediterranean Sea: A numerical sensitivity study[END_REF], an adequate representation of topographic features and local dynamic effects is required to correctly describe the interaction between orography and atmospheric processes. Furthermore, initial conditions have been shown to have a significant influence on rainfall forecasting [START_REF] Kunz | The severe hailstorm in southwest Germany on 28 July 2013: characteristics, impacts and meteorological conditions[END_REF][START_REF] Khodayar | Multi-scale observations of atmospheric moisture variability in relation to heavy precipitating systems in the northwestern Mediterranean during HyMeX IOP12[END_REF]Caldas-Álvarez et al., 2017).

For each point of the diagram in Fig. 6, we compute its distance from the origin (perfect score -A = 0, S = 0). The dotted circles respectively contain the 25 %, 50 %, and 75 % points with the smallest distance. The radii of the circles are much larger for LT34, confirming a degradation of the scores for longer lead times.

Clusters

We use our clustering procedure (as defined in Sect. 2.2.2) to analyse the characteristics of the forecast QPF errors along with the regional properties. SAL components are stated for each day of each cluster associated with HPEs, i.e. C2, C3, https://doi.org/10.5194/nhess-20-1369-2020

Nat. Hazards Earth Syst. Sci., 20, 1369-1389, 2020 4). A positive correlation is found for all three clusters, which corroborates the results from Fig. 6 where HPEs correspond to the highest S-component values. The maximum correlation is found for cluster 3. Although correlations are statistically significant, it is worth noting that values are quite weak (in particular for cluster 5).

Sensitivity to physical parameterizations

The SAL measure is analysed separately for the 10 different physical packages to study corresponding systematic errors. More specifically, we raise the following questions. Do the errors based on an object-quality measure and computed for the different physics implemented in an ensemble system show different rainfall structure properties? Which physical packages are more sensitive to the intense rainfall forecast errors? As in Sect. 3.1, we first distinguish the results for the HPE group before the cluster ones.

HPEs

Probability density distributions for each SAL component are separately computed for each physics reforecast (Fig. 8), considering only the HPEs. Colours correspond to four categories, depending on the parametrization of the deep convection. The figure highlights that members from each of the two main parametrization schemes (B85 and PCMT) have similar behaviours. Considering the A component, PCMT members are more centred around zero than B85 at LT12. This effect is higher at LT34, for which B85 and PCMT density distributions are more shifted. At LT34, more events with a positive A component are associated with PCMT, whereas negative values are more recurrent in B85. The A component never exceeds +1, but significant underestimations are observed. This range of values stems from the fact that the forecast verification is applied to a subsample of the observation limited to the most extreme events. For these specific events, a model underestimation is more frequent than an overestimation. At short lead times, the separation between the two deep-convection schemes is also well established for the S component (Fig. 8), but it becomes mixed up at LT34 (Fig. 8). One reason for this behaviour could be that predictability decreases at LT34, so that discrepancies in spatial rainfall structure assigned to the physics families become less identifiable. The S component is positively skewed in all cases (in particular for the B85 physics at LT12 lead time). This supports the previous analysis of the S component (Figs. The results for the S component also highlight better skills for PCMT schemes for HPEs, especially at short lead times.

Focusing on high values of S, B85 exhibits a stronger distribution tail at LT12, while both schemes seem comparable for LT34.

For the L component, the maxima of the density distributions are higher for PCMT at lead time LT12, implying a more significant number of good estimations of pattern location. Regarding the tail of the L-component PDF, it is globally more pronounced at LT34 than LT12. This means that the location of HPEs is poorly forecasted at long lead times. Concerning the behaviour of the forecasts that use the CAPE or B85 mod schemes, their A-component PDFs are close to the B85 PDFs. This is not observed for the other components. For the S component, the CAPE distribution follows the PCMT one at LT12. For the L component, the B85 mod PDF is close to the B85 ones, while CAPE shows different behaviour from all the other physics. The use of a closure based on CAPE, rather than on the convergence of humidity, seems to modulate the location of precipitation produced by this deep-convection parametrization scheme. Moreover, at LT34 CAPE is characterized by a lower number of strong location errors, compared to the other physics.

Clusters

According to the results of the previous section, which show that the predictability of intense rainfall events is sensitive to the parametrization of the deep convection, we have continued to analyse the model behaviour for the four different deep-convection schemes: B85, B85 mod , CAPE, and PCMT. The link between the behaviour of the physical schemes and belonging to a particular cluster is statistically assessed through the SAL component differences between the schemes.

Any parametric goodness-of-fit tests, which assume normality, have been discarded, because SAL values are not normally distributed. We choose the k-sample Anderson-Darling (AD) test [START_REF] Scholz | K-Sample Anderson-Darling Tests[END_REF][START_REF] Mittermaier | Feature-Based Diagnostic Evaluation of Global NWP Forecasts[END_REF], in order to evaluate whether differences between two given distributions are statistically significant. It is an extension of the two-sample test (Darling, 1957), originally developed starting from the classic Anderson-Darling test (Anderson and Darling, 1952). The k-sample AD test is a nonparametric test designed to compare continuous or discrete subsamples of the same distribution. In this case the test is implemented for the evaluation of the pairs of distributions.

The tests are performed for the comparison of each pair of PDFs combined from the four deep-convection families and from the three clusters' classification. For the A component, PCMT physics distributions depart significantly from B85 schemes at all lead times, while B85 mod and CAPE perform as B85, meaning that the modified versions of B85 weakly affect physics behaviour (not shown).

With respect to the S-component distributions, k-sample AD tests show significant differences between B85 and PCMT physics for LT12, but not for the longest lead times (not shown). At LT34 we observe a convergence of the physics scheme towards a homogeneous distribution, meaning that the differences between physics are negligible.

The test applied to the location component does not reveal significant differences between the PDFs. We suppose that the limited dimensions of the domain employed in this study, as well as its irregular shape, may lead to a less coherent estimation of the location, resulting in a degradation of the score significance. Since the L-component result is not informative about HPEs, it is ignored hereafter.

Once the statistical differences between the PDFs of the physics have been examined, it is interesting to compare the relative error on the amplitude and structure components. S and A component errors are estimated by comparing the shapes of their distributions. Empirical cumulative density functions (ECDFs) of S and A components are computed separately for each cluster and lead time (LT12 and LT34). We show an example of an ECDF for cluster 2 at LT34 (Fig. 9). Forecasts are perfect when the ECDF tends towards a Heaviside step function, which means that the distribution tends towards the Dirac delta function centred on zero. These functions are estimated over a bounded interval, corresponding to the finite range of S and A components. The deviation from the perfect score was quantified by estimating the area under the ECDF curve on the left side and the area above the ECDF curve on the right side:

err -= 0 -2 F (x)dx - 0 -2 H (x)dx = 0 -2 F (x)dx -0 = 0 -2 F (x)dx, (9) err
+ = 2 0 H (x)dx - 2 0 F (x)dx = 2 - 2 0 F (x)dx, (10) err
= err -+ err + = 2 - 2 0 F (x)dx + 0 -2 F (x)dx, (11)
where F (x) is the ECDF computed for A or S, H (x) is the Heaviside step function, and err is the forecast error for a given component. The lower and upper boundaries of the integrals are equal to -2 and +2, because A and S components range between these two values by construction. Since the previous k-sample AD test highlighted significant differences within the two main classes B85 and PCMT, the evaluation of the errors is limited to these two specific classes.

The results of the error diagnostic err for the A component are shown in Fig. 10a. Errors increase with lead time. We note that the negative errors are always at least twice as large as the positive ones. Forecasted averaged rainfall amounts are almost always underestimated. PCMT produces overall better A-component statistics, except for cluster 3 at LT34. It is interesting to observe that the weakest errors are associated with cluster 3, which is the most extreme one. Since cluster 3 collects a large number of precipitation events impacting the Cévennes chain, we may suppose that the domain-averaged rainfall amounts are more predictable in situations of precipitation driven by the orography. Concerning the S-component evaluation (see Fig. 10b), structures of rainfall patterns are better forecasted for heavy-rainfall events (clusters 3 and 5) than for the remaining classes of events. In contrast to the A component, the S component exhibits the highest err + for the B85 scheme for most of the cases (majority of + signs in Fig. 10b), whereas this trend is not systematic for PCMT physics. PCMT globally performs better than B85, except for cluster 2. As with the amplitude A, the S component gets worse for longer lead times, resulting in a shift to larger err -for both B85 and PCMT physics (more -sign for LT34 in Fig. 10a,b). The lowest errors of the S component are achieved for cluster 5. Cluster 5 HPEs are known to have specific regional properties whose influence on S-component results should be studied with further diagnostics.

Rainfall object analysis

We now analyse the physical properties of the objects, i.e. the number of objects from a rainfall field and the objectintegrated volumes, according to the different clusters. All the statistics are applied separately to the B85, PCMT physics, and observations. For each day of the dataset period, the thresholds defined in Sect. 2.3.1 lead to the identification of a certain number of precipitating objects. The frequency of the number of objects per day is plotted by means of normalized histograms for the three clusters (Fig. 11). Clusters 2 and 3 show maximum frequency for one-and three-object ranges, whereas cluster 5 is dominated by one object per day. This specific property of cluster 5 can explain the best result obtained for the S component (Sect. 3.2.2). Indeed, we may assume that S-component estimation is more accurate for a one-to-one object comparison. The other clusters frequently display rainfall accumulation bands split over the domain, typically over the Cévennes and Alpine regions. Object identification for PCMT forecast shows that there is an overestimation of single object days compared to the observation and to the B85 physics scheme, a behaviour emphasized in clusters 3 and 5.

More details about the magnitude of the objects can be produced by computing the integrated mass per object, M k (see Sect. 2.3.1). First, for each day, objects are sorted from the largest to the smallest integrated mass. Integrated mass distribution of the two heaviest objects (noted O 1 and O 2 ) is then dispatched as a function of the number of objects for each cluster in Fig. 12. First, the range value of M is highly variable from one cluster to another. Maximum values are observed for cluster 3, while the magnitude for clusters 2 and 5 is comparable. The decrease in the mass for O 1 is clearer for cluster 3, meaning that a high number of objects over the domain leads to a natural decrease in the M value of the heaviest ones. We think that a part of the total integrated mass is then redistributed to the other objects. This is confirmed by O 2 curves since its mass increases with the number of objects. Conversely, for cluster 5, O 1 mass increases with the number of the objects, while O 2 is almost stable. The gap between O 1 and O 2 masses is at a maximum in the most extreme clusters (3 and 5). This suggests that when computing the volume V (see Eq. 7) and L2 (see Eq. 4), the weighted average is dominated by the object O 1 . This implies that the verification could be considered a single-to-single object metric.

We now examine the ratio between the daily maximum rainfall of objects O 1 and O 2 . This ratio ranges between 1.5 and 3, which means that O 1 represents the essential contribution of the daily rainfall peak. Since O 1 base area tends to be significantly larger that O 2 , the information related to the inner object maximum rainfall is diluted in the large base area, resulting in a flat weak mean intensity of the object. This last result appears to support the fact that the SAL metric gives more weight to the object that contains the most intense rainfall.

The comparison between the model reforecast physics and the observations is addressed using the whole distribution of daily mass M from the objects O i identified across the full reforecast dataset, where i ranges between 1 and the total number N of objects. We proceed separately for each physical package. For a given scheme and cluster, the quantile values corresponding to the selected dataset are sorted in ascending order and then plotted versus the quantiles calculated from observations (Fig. 13). Half of the quantile distributions are not visible as they correspond to very weak pattern masses. For cluster 2 and PCMT physics most of the distribution of object mass is close to the observations. However all other physics distributions are skewed to the right compared to the observations for values below 10 000 mm. This behaviour is also observed for cluster 5 and it involves PCMT physics as well, for values between percentile 0.5 and percentile 0.7. Overall, in the quantile-quantile plot for cluster 5, the PCMT outperforms B85. In cluster 3, discrepancies between PCMT, B85, and the observations are of opposite sign, with PCMT being slightly above the observations, while B85 shows a weak underestimation. The CAPE physics distribution is left skewed compared to the observations and to the other physics. These results highlight some interesting properties of the models in predicting the rainfall objects. Except for some deviation concerning a few extreme cases of cluster 2 and a small portion of distributions of cluster 5, object mass distribution of physics is similar to the distribution drawn from the observation, especially for cluster 3. This means that the forecast is able to reproduce

err + ≥ 2; -if 1.1 ≤ err - err + < 2; ≈ if 0.9 < err - err + < 1.1; + if 0.5 < err - err + ≤ 0.9; + (bold) if err - err + ≤ 0.5.
the same proportion of rainfall amounts inside a feature as the observations, even concerning the extreme right tail of the distributions, which corresponds to the major events of the series.

Summary and conclusions

In this study we have characterized the systematic errors of 24 h rainfall amounts from a reforecast ensemble dataset, covering a 30-year fall period. A 24 h rainfall observation reference has been produced on a regular grid with a resolution identical to the model in order to run point-to-point verification. We applied an object-based quality measure in order to evaluate the performance of the forecasts of any kind of HPE. Then, we took advantage of a rainfall clustering to analyse the dependence of systematic errors on clusters.

The selection of the HPEs within the reference dataset was based on a peak-over-threshold approach. The spatial regional discrepancies between HPEs are highlighted on the basis of the k-means clustering of the 24 h rainfall. Finally, we analysed the rainfall object properties in the model and in the observation to underline the rainfall field object properties for which the model acts distinctly.

The peak-over-threshold criterion leads to the selection of 192 HPEs, confirming that the most impacted regions are the Cévennes area and part of the Alps. The composite analysis for the five clusters shows that each cluster is associated with a specific class and location of 24 h precipitation events. It was found that 86 % of the number of HPEs are included in clusters 2, 3, and 5. Cluster 2 and 3 HPEs predominantly impact the Cévennes and Alps areas, while cluster 5 HPEs are located over the Languedoc-Roussillon region. Moreover, clusters 3 and 5 include the most extreme HPEs. Only diagnostics for clusters 2, 3, and 5 are considered.

The SAL object-quality measure has been applied distinctly to the 10 physics schemes (one per member) of the reforecast dataset and compared to the rainfall reference. It shows that the model's overall behaviour for HPE forecasting is characterized by negative A components and positive S components. As in grid-point rainfall verification, all the SAL components get worse as a function of lead time. The model HPE rainfall objects tend to be more extended and less peaked. Even though their corresponding domainaverage amplitude is weaker, it does not mean that the event maximum intensity is always weaker. This result is important, showing modellers that, even for intense rainfall events when orography interaction and quasi-stationarity mesoscale systems play a great role, the model tends to reproduce rainfall patterns with greater extension, rather than both smaller extension and weaker-intensity patterns.

In order to show regional disparities in the model behaviour, the SAL diagnostics have been divided according to the clusters and it shows interesting results. First, the A component negative contribution for the whole sample is higher, showing that on average more underestimation than overestimation is observed for the amplitude SAL component. It is notably the case for the most extreme clusters (over Cévennes and over Languedoc-Roussillon). However, when considering both positive and negative contributions to the integrated A component, the most extreme cluster (cluster 3) leads to better scores. This could mean that the variability of the A component is positively reduced for the most intense events. This is quite surprising and could reinforce the role of orography in this error decrease. As for the S-component distribution, we showed it is slightly positively skewed for clusters 2 and 3, while for cluster 5 the distribution of the S component is more centred. Likewise for the A component the integrated balance of positive and negative S-component contributions leads to better results for clusters 3 and 5. It is even more remarkable for cluster 5, for which the S compo- nent reaches the best score. Though it is difficult at this point to determine whether this characterizes an actual contrast in the model behaviour or if it is due to the physical properties of the cluster 5 events. One hypothesis could be related to the large number of single objects characterizing this cluster. forecasting until lead time ranges longer than 3 d, after which no significant differences appear. This difference is clearly in favour of the PCMT scheme, which performs better than B85 for both SAL A and S components and in the majority of the subsampled scores considering the HPEs or the regional clusters. However, this PCMT asset is not huge, and both physics schemes can contribute to good or bad forecasts. The main significant difference is for the S component for the most in-tense rainfall, which shows that PCMT better approximates the structure of the rainfall patterns in these cases.

In light of the ability of our method to produce significant results even after several subsampling steps, we decided to study further statistical characterization of the SAL rainfall objects. It has been shown that in most cases, one large object stands out among other smaller objects, which often gathers the biggest part of the rain signal. For cluster 5, characterized by the Languedoc-Roussillon HPEs, the rainfall distribution could even be considered a single-object rainfall field. Then we focused on the ranked distributions (quantile-quantile analysis) of the object masses to compare the overall rainfall climatology of the model with the reference. First, this analysis showed that in particular the weakest precipitation is overestimated by all physics schemes. However, looking at the object mass distributions for the whole period, we find they are relatively close between all the physics schemes and the observation for most extreme rainfall events, especially for the PCMT deep-convection scheme. This statistical result implies that a global model should be able to reproduce a reliable distribution of rainfall objects along a long time period, e.g. the climate of the model and of the observations are close to each other. Therefore, in the case of PEARP, most of the forecast errors are mainly related to a low consistency between observed and forecasted fields, rather than to an inability of the prediction system to produce intense precipitation amounts.

This last result, objectively quantified for high-rainfall event thresholds (around 100 to 500 mm) on a long enough period, is important for two reasons. The first one concerns atmospheric modellers, showing that the physics schemes are able to reproduce climatological distributions of the most challenging rainfall events. On this basis, future research could investigate other sources of uncertainties like from the analysis set-up and implement ensuing model improvements. The model physics perturbation technique should then play a greater role in the control of the ensemble dispersion. From this perspective, the novel reanalysis ERA5 would be interesting to use, in particular its perturbed members, to improve the uncertainty from initial conditions in the reforecast. The second lesson to be learned from this study is that it is worth focusing on the study of model behaviour on intense event forecasting as it provides important learning to ensemble model end-users, in particular in the context of decisionmaking based on weather forecast. Quantifying systematic errors could also be used to favourably improve their inclusion in nested forecast tool processes.

In terms of methodology, this study also highlights that the combination of SAL verification and clustering is a relevant approach to show systematic errors associated with regional features for intense precipitation forecasting. This achievement is only enabled by the availability of a long reforecast dataset. This methodology could be further extended to a different model and another geographic region, on the condition of sampling a large number of HPEs.

The inter-comparison between some model physics deepconvection schemes and their role in HPEs predictability shows it is of course very sensitive for designing multiphysics ensemble forecasting systems. While the sensitivity to the initial perturbations was not studied in this work, the forecast of intense rainfall seems to be mainly driven by the classes of deep-convection parameterizations. Since the physical parametrization set-up is built by replicated schemes, the model error representation might lack an exhaustive sampling of the forecasted trajectories. Using more than two deep-convection parametrization schemes may improve the representation of model errors, at least for heavyprecipitation events.

Figure 1 .

 1 Figure 1. Panel (a) shows a situation map of the investigated area (rectangle with red edges) with respect to western Europe and the Mediterranean Sea. Panel (b) shows the rain-gauge network used for the study. Red diamonds represent the rain gauges selected for crossvalidation testing, and blue dots represent the rain gauges selected for cross-validation training. Panel (c) shows the 0.1 • × 0.1 • model grid (in blue), along with the location of three key areas. Domain D is located within the borders of the model grid (c).

Figure 2 .

 2 Figure 2. Annual average of HPE occurrence per grid point (in green). The composite of daily rainfall amounts (mm d -1 ) of the HPE dataset is represented by the blue isohyets.

Figure 3 .

 3 Figure 3. Rainfall composites (mm d -1 ) for the five clusters selected by the K-means algorithm. The bottom right panel shows the probability density distribution of the maximum daily rainfall (mm) for each cluster class.

Figure 4 .

 4 Figure 4. SAL pattern analysis for the case of 28 October 2004, applied to the observation data (a) and one 60 h lead time forecast (b). Base contours of the identified objects are shown as red lines. Grey points stand for the rain barycentre of each pattern, and the grey diamond depicts the rain barycentre for the whole field. The size of the barycentre points is proportional to the integrated mass of the associated object.

Figure 5 .

 5 Figure 5. Relationship between the daily rainfall grid-point maximum algebraic error and the A component of the SAL score. HPE days are plotted in red, while other days are in black. Panel (a) is for LT12 lead time, and (b) shows LT34 lead time. Linear regression analysis is added to the plot.

Figure 6 .

 6 Figure 6. Relationship between the A component and the S component of the SAL score (SAL diagrams) for HPEs only, for lead times LT12 (a) and LT34 (b). Blue triangles represent HPEs with grid-point maximum rainfall under 200 mm d -1 , and red triangles show rainfall amounts above 200 mm. Triangles are proportional to the rainfall value. Some main characteristics of the component distribution are plotted. The median value is shown as dashed lines, and percentiles 25 and 75 delimitate the boxes. Circles represent the limits 25th, 50th, and 75th percentiles to the best score (A = 0, S = 0).

Table 4 .Figure 7 .

 47 Figure 7. A-component (a) and S-component (b) normalized histograms and probability density functions for clusters 2, 3, and 5. Results for lead time LT12 are plotted as black lines, and results for lead times LT34 are in grey.

  6 and 7), showing that for intense rainfall the https://doi.org/10.5194/nhess-20-1369-2020 Nat. Hazards Earth Syst. Sci., 20, 1369-1389, 2020 model mostly produces a larger and flatter rainfall signal.

Figure 8 .

 8 Figure 8. Probability density functions of the three SAL components for the HPEs and for each physics of the reforecast system (coloured lines). Physics scheme are gathered in four categories depending on the parametrization of the deep convection: PCMT (blue), B85 (orange), B85 mod (green), and CAPE (purple). The left column corresponds to lead time LT12, and the right column relates to lead time LT34.

Figure 9 .

 9 Figure 9. Empirical cumulative distribution function of the A component computed from cluster 2 at lead time LT34 for the four classes of physics schemes.

Figure 10 .

 10 Figure 10. Dumbbell plot of integrated error diagnostics computed using Eq. (11). Colours refer to B85 (orange) and PCMT (blue) deepconvection parametrization schemes. Results are stratified on the basis of the clusters and lead times. Symbols denote whether positive or negative errors dominate. These signs are defined using the following definition: -(bold) if err -

Figure 11 .

 11 Figure 11. Normalized histograms of the daily number of SAL patterns, for the B85 physics scheme (red), PCMT (blue), and observation (green). Panels correspond to the three-cluster classification.

Figure 12 .

 12 Figure 12. Distribution of the SAL first-pattern O 1 rain amount according to the number of patterns per day. Curves stand for the median of the distribution, and shaded areas range between the 25th and 75th percentiles. The dashed lines correspond to the second ranked SAL pattern O 2 rain amount.

Figure 13 .

 13 Figure 13. Quantile-quantile plot between SAL pattern rain amounts from the model (Y axis) and from the observation (X axis). Physics schemes are gathered into four classes (B85, PCMT, B85 mod , CAPE). Observation deciles correspond to the vertical dashed lines.

Table 1 .

 1 Physical parameterizations used in the ensemble reforecast.

		Turbulence Shallow convection Deep convection Oceanic flux
	Ref TKE	KFB	B85	ECUME
	1	TKE	KFB	B85	ECUME mod
	2	L79	KFB	B85 mod	ECUME
	3	L79	KFB	CAPE	ECUME
	4	TKE mod	KFB	B85	ECUME
	5	TKE	EDKF	B85	ECUME
	6	TKE	PMMC	PCMT	ECUME
	7	TKE	KFB	PCMT	ECUME
	8	TKE	PCMT	PCMT	ECUME
	9	TKE	KFB	B85	ECUME

Table 2 .

 2 Classification of days computed from 24 h rainfall amounts in southern France

Table 3 .

 3 Contingency table computed for rainy and dry days. Contingency table Obs. rainy day Obs. dry day

	Model rainy day	3258	84
	Model dry day	226	62

Nat. Hazards Earth Syst. Sci., 20, 1369-1389, 2020 https://doi.org/10.5194/nhess-20-1369-2020

https://doi.org/10.5194/nhess-20-1369-2020Nat. Hazards Earth Syst. Sci., 20, 1369-1389, 2020

Acknowledgements. We wish to thank the Occitanie region and Météo France for co-funding the PhD programme during which this work has been done. The authors are grateful to the two reviewers for their remarks and careful reading.

Review statement. This paper was edited by Joaquim G. Pinto and reviewed by two anonymous referees.

Data availability. Research data can be accessed by contacting Matteo Ponzano (matteo.ponzano@meteo.fr) and the other authors.

Author contributions. MP, BJ, and LD conceived and designed the study. MP carried out the formal analysis, wrote the whole paper, carried out the literature review, and produced the observation reference dataset. BJ built the hindcast dataset. BJ, LD, and PA reviewed and edited the original draft.

Competing interests. The authors declare that they have no conflict of interest.


Leila Amgoud
email: amgoud@irit.fr

Jonathan Ben-Naim
email: bennaim@irit.fr

Weighted Bipolar Argumentation Graphs: Axioms and Semantics

The paper studies how arguments can be evaluated in weighted bipolar argumentation graphs (i.e., graphs whose arguments have basic weights and may be supported and attacked). It introduces principles that an evaluation method (or semantics) would satisfy, analyzes existing semantics with respect to them, and finally proposes a new semantics for the class of non-maximal acyclic graphs. 1

Introduction

Argumentation is a form of common-sense reasoning, which consists of justifying claims by arguments. The latter have generally basic strengths, and may be both attacked and supported by other arguments, leading to the so-called weighted bipolar argumentation graphs.

Two families of semantics were proposed in the literature for the evaluation of arguments in such graphs: extension semantics [Cayrol and Lagasquie-Schiex, 2005b;[START_REF] Oren | Semantics for evidence-based argumentation[END_REF][START_REF] Brewka | [END_REF][START_REF] Boella | [END_REF][START_REF] Nouioua | Bipolar argumentation frameworks with specialized supports[END_REF], and gradual semantics [Cayrol and Lagasquie-Schiex, 2005a;[START_REF] Evripidou | Valentinos Evripidou and Francesca Toni. Quaestio-it.com: a social intelligent debating platform[END_REF][START_REF] Baroni | [END_REF][START_REF] Rago | [END_REF][START_REF] Pazienza | [END_REF]. The former evaluate mainly sets of arguments (extending Dung's 1995 semantics) while the latter focus on individual arguments.

This paper extends our previous works on axiomatic foundations of semantics for unipolar graphs (support graphs [Amgoud and Ben-Naim, 2016b] and attack graphs [Amgoud and Ben-Naim, 2016a]). It defines principles that a semantics would satisfy in a bipolar setting. Such principles are useful for defining reasonable semantics, for a better understanding of the design choices or foundations of each semantics, and for comparing pairs of semantics. As a second contribution, the paper analyzes existing semantics against the principles. The main conclusion is that extension semantics do not harness the potential of support relations. For instance, when there is no attack in a graph, the existing semantics declare all (supported, non-supported) arguments of the graph as equally strong. Gradual semantics take into account supporters in this particular case, however they violate some key principles. The third contribution of the paper is the definition of a novel gradual semantics for the sub-class of non-maximal acyclic bipolar graphs. We show that it satisfies all the proposed principles. Furthermore, it avoids a big jump problem that may impede the relevance of existing gradual semantics in some practical applications.

Basic Concepts

Let us now introduce weighted bipolar argumentation graphs.

Definition 1 A weighted bipolar argumentation graph is a quadruple

A = A, w, R, S , where A is a finite set of argu- ments, w a function from A to [0, 1], R ⊆ A×A, S ⊆ A×A.
Let WAG denote the set of all such graphs.

Given two arguments a and b, aRb (resp. aSb) means a attacks (resp. supports) b, and w(a) is the intrinsic strength of a. The latter may be the certainty degree of the argument's premises, trustworthiness of the argument's source, . . ..

A semantics is a function assigning an overall strength from the unit interval [0, 1] to each argument of a weighted bipolar graph. Arguments that get value 1 are extremely strong whilst those that get value 0 are worthless.

Definition 2 A semantics is a function S transforming any

A = A, w, R, S ∈ WAG into a function Deg S A from A to [0, 1]. For a ∈ A, Deg S A (a)
is the overall strength of a. Let us recall the notion of isomorphism between graphs.

Definition 3 Let

A = A, w, R, S , A ′ = A ′ , w ′ , R ′ , S ′ ∈ WAG. An isomorphism from A to A ′ is a bijective function f from A to A ′ such that: i) ∀ a ∈ A, w(a) = w ′ (f (a)), ii) ∀ a, b ∈ A, aRb iff f (a)R ′ f (b), iii) ∀ a, b ∈ A, aSb iff f (a)S ′ f (b).

Notations: Let

A = A, w, R, S ∈ WAG and a ∈ A. Att A (a) = {b ∈ A | bRa} and sAtt A (a) = {b ∈ Att A (a) | Deg S A (b) = 0}. Similarly, Supp A (a) = {b ∈ A | bSa} and sSupp A (a) = {b ∈ Supp A (a) | Deg S A (b) = 0}. Let A ′ = A ′ , w ′ , R ′ , S ′ ∈ WAG such that A ∩ A ′ = ∅. Let A ⊕ A ′ be A ′′ , w ′′ , R ′′ , S ′′ ∈ WAG such that A ′′ = A ∪ A ′ , R ′′ = R ∪ R ′ , S ′′ = S ∪ S ′ , ∀x ∈ A ′′ , the following holds: w ′′ (x) = w(x), if x ∈ A; w ′′ (x) = w ′ (x), if x ∈ A ′ .

Principles for Semantics

In what follows, we propose principles that shed light on foundational properties behind semantics. In other words, properties that help us to better understand the underpinnings of semantics, and that facilitate their comparisons. The first nine principles are simple combinations of properties proposed for graphs with only one type of interactions (support in [Amgoud and Ben-Naim, 2016b] and attack in [Amgoud and Ben-Naim, 2016a]). The three last principles are new.

The first basic principle, Anonymity, states that the degree of an argument is independent of its identity.

Principle 1 (Anonymity) Semantics S satisfies anonymity iff, for all A = A, w, R, S , A ′ = A ′ , w ′ , R ′ , S ′ ∈ WAG, for any isomorphism f from A to A ′ , the following property holds: ∀ a ∈ A, Deg S A (a) = Deg S A ′ (f (a)). Bi-variate independence principle states the following: the overall strength of an argument a should be independent of any argument b that is not connected to it (i.e., there is no path from b to a, ignoring the direction of the edges).

Principle 2 (Bi-variate Independence) Semantics S satisfies bi-variate independence iff, for all A = A, w, R, S ,

A ′ = A ′ , w ′ , R ′ , S ′ ∈ WAG such that A ∩ A ′ = ∅, the fol- lowing property holds: ∀ a ∈ A, Deg S A (a) = Deg S A⊕A ′ (a)
. Bi-variate directionality principle states that the overall strength of an argument should depend only on its incoming arrows, and not on the arguments it itself attacks or supports.

Principle 3 (Bi-variate Directionality) Semantics S satisfies bi-variate directionality iff, for all

A = A, w, R, S , A ′ = A ′ , w ′ , R ′ , S ′ ∈ WAG such that A = A ′ , w = w ′ , R ⊆ R ′ ,
and S ⊆ S ′ , the following holds: for all a, b, x ∈ A, if R ′ ∪ S ′ = R ∪ S ∪ {(a, b)} and there is no path from b to x, then Deg S A (x) = Deg S A ′ (x). A path can mix attack and support relations, but the edges must be directed from b to x.

Bi-variate equivalence principle ensures that the overall strength of an argument depends only on its basic strength and the overall strengths of its attackers and supporters.

Principle 4 (Bi-variate Equivalence) Semantics S satisfies bi-variate equivalence iff, for any

A = A, w, R, S ∈ WAG, for all a, b ∈ A, if: • w(a) = w(b),
• there exists a bijective function

f from Att A (a) to Att A (b) s.t ∀x ∈ Att A (a), Deg S A (x) = Deg S A (f (x)), • there exists a bijective function f ′ from Supp A (a) to Supp A (b) s.t ∀x ∈ Supp A (a), Deg S A (x) = Deg S A (f (x)), then Deg S A (a) = Deg S A (b).
Stability principle states the following: if an argument is neither attacked nor supported, its overall strength should be equal to its intrinsic strength.

Principle 5 (Stability) Semantics S satisfies stability iff, for any

A = A, w, R, S ∈ WAG, for any a ∈ A, if Att A (a) = Supp A (a) = ∅, then Deg S A (a) = w(a)
. Neutrality principle states that worthless attackers and worthless supporters have no effect on their targets.

Principle 6 (Neutrality) Semantics S satisfies neutrality iff for any

A = A, w, R, S ∈ WAG, for all a, b, x ∈ A, if: • w(a) = w(b), • Att A (a) ⊆ Att A (b), • Supp A (a) ⊆ Supp A (b), • Att A (b) ∪ Supp A (b) = Att A (a) ∪ Supp A (a) ∪ {x} such that Deg S A (x) = 0, then Deg S A (a) = Deg S A (b).
If an argument a is equally or less attacked than an argument b, and equally or more supported than b, then Bi-variate Monotony ensures that a should be at least as strong as b.

Principle 7 (Bi-variate Monotony) Semantics S satisfies bi-variate monotony iff, for any

A = A, w, R, S ∈ WAG, for all a, b ∈ A such that: • w(a) = w(b), • Att A (a) ⊆ Att A (b), • Supp A (b) ⊆ Supp A (a),
the following holds:

• Deg S A (a) ≥ Deg S A (b); (Monotony) • if (Deg S A (a) > 0 and sAtt A (a) ⊂ sAtt A (b)) or (Deg S A (b) < 1 and sSupp A (b) ⊂ sSupp A (a)), then Deg S A (a) > Deg S A (b). (Strict Monotony)
Bi-variate Reinforcement principle concerns the quality of attackers and supporters. It states that any argument becomes stronger if the quality of its attackers is reduced and the quality of its supporters is increased.

Principle 8 (Bi-variate Reinforcement) Semantics S satisfies bi-variate reinforcement iff, for any A = A, w, R, S ∈ WAG, for all C, C ′ ⊆ A, for all a, b ∈ A, for all x, x ′ , y, y ′ ∈ A \ (C ∪ C ′ ) such that:

• w(a) = w(b), • Att A (a) = C ∪ {x} Att A (b) = C ∪ {y}, • Supp A (a) = C ′ ∪ {x ′ } Supp A (b) = C ′ ∪ {y ′ }, • Deg S A (x) ≤ Deg S A (y) Deg S A (x ′ ) ≥ Deg S A (y ′ ), the following holds: • Deg S A (a) ≥ Deg S A (b); (Reinforcement) • if (Deg S A (a) > 0 and Deg S A (x) < Deg S A (y)) or (Deg S A (b) < 1 and Deg S A (x ′ ) > Deg S A (y ′ )), then Deg S A (a) > Deg S A (b). (Strict Reinforcement)
Our next principle is not mandatory. Its suitability depends on the nature of arguments (opinion-based, analogical, etc.). It combines Imperfection property from [Amgoud and Ben-Naim, 2016b] with Resilience one from [Amgoud and Ben-Naim, 2016a]. Imperfection states that an argument whose basic strength is less than 1 cannot be fully rehabilitated by supports. This property prevents irrational behaviors, like fully accepting fallacious arguments that are supported. Resilience states that an argument whose basic strength is positive cannot be completely destroyed by attacks.

Principle 9 (Resilience) Semantics S satisfies resilience iff, for any A = A, w, R, S ∈ WAG, for all a ∈ A, if 0 < w(a) < 1, then 0 < Deg S A (a) < 1. The next principles are new and answer the question: how the overall strengths of attackers should be aggregated with those of supporters? To answer this question, it is important to specify which of the two types of interactions is more important. We argue that a support is weaker than an attack since an attack is sufficient for weakening a whole argument while a support may only strengthen a part of it. Hence, in what follows we consider the case where an attack is more important than a support (Franklin principle) and the case where both relations are equally important (Strict Franklin principle). The latter ensures that an attacker and a supporter of equal strength counter-balance each other. This principle is suitable in multiple criteria decision making contexts.

Principle 10 (Franklin) A semantics S satisfies franklin iff, for any

A = A, w, R, S ∈ WAG, for all a, b, x, y ∈ A, if • w(a) = w(b), • Att A (a) = Att A (b) ∪ {x}, • Supp A (a) = Supp A (b) ∪ {y},
• Deg S A (x) = Deg S A (y), then the following hold: Weakening states that if the attackers of an argument overcome its supporters, the argument looses weight. The idea is that supporters are not sufficient for counter-balancing attackers, however, they may mitigate the global loss due to attacks.

• Deg S A (a) ≤ Deg S A (b), (Franklin)
Principle 11 (Weakening) Semantics S satisfies weakening iff, for any A = A, w, R, S ∈ WAG, for any a ∈ A, if w(a) > 0 and there exists an injective function f from Supp A (a) to Att A (a) such that:

• ∀x ∈ Supp A (a), Deg S A (x) ≤ Deg S A (f (x)); and • sAtt A (a) \ {f (x) | x ∈ Supp A (a)} = ∅ or ∃x ∈ Supp A (a) such that Deg S A (x) < Deg S A (f (x)), then Deg S
A (a) < w(a). Strengthening states that if the supporters of an argument overcome its attackers, the argument gains weight. Indeed, attacks are not sufficient for counter-balancing supports, however, they may mitigate the global gain due to supports.

Principle 12 (Strengthening) Semantics

S satisfies strengthening iff, for any A = A, w, R, S ∈ WAG, for any a ∈ A, if w(a) < 1 and there exists an injective function f from Att A (a) to Supp A (a) such that:

• ∀x ∈ Att A (a), Deg S A (x) ≤ Deg S A (f (x)); and • sSupp A (a) \ {f (x) | x ∈ Att A (a)} = ∅ or ∃x ∈ Att A (a) such that Deg S A (x) < Deg S A (f (x)),

Formal Analysis of Existing Semantics

There are two families of semantics in case of bipolar argumentation graphs: extension semantics [Cayrol and Lagasquie-Schiex, 2005b;[START_REF] Oren | Semantics for evidence-based argumentation[END_REF][START_REF] Brewka | [END_REF][START_REF] Boella | [END_REF][START_REF] Nouioua | Bipolar argumentation frameworks with specialized supports[END_REF] and gradual semantics [Cayrol and Lagasquie-Schiex, 2005a;[START_REF] Baroni | [END_REF][START_REF] Rago | [END_REF].

Extension semantics extend Dung's 1995 ones for accounting for supports between arguments. They take as input an argumentation graph A, w, R, S whose arguments have all the same intrinsic strength, and return sets of arguments, called extensions. From the extensions, a three-valued qualitative degree is assigned to every argument. Indeed, an argument is accepted if it belongs to all extensions, undecided (or credulously accepted) if it belongs to some but not all extensions, and rejected if it does not belong to any extension. When the support relation is empty, the semantics proposed in [Cayrol and Lagasquie-Schiex, 2005b;[START_REF] Oren | Semantics for evidence-based argumentation[END_REF][START_REF] Brewka | [END_REF][START_REF] Boella | [END_REF][START_REF] Nouioua | Bipolar argumentation frameworks with specialized supports[END_REF] coincide with Dung's ones. Thus, they violate the principles that are violated by Dung's semantics (see [Amgoud and Ben-Naim, 2016a] for a detailed analysis of Dung's semantics). For instance, stable semantics violates Bi-variate Independence, Bi-variate Equivalence, Stability, Resilience, and Strict Monotony. When the attack relation is empty, the approaches from [Cayrol and Lagasquie-Schiex, 2005b;[START_REF] Oren | Semantics for evidence-based argumentation[END_REF][START_REF] Boella | [END_REF]] return a single extension, which contains all the arguments of the graph at hand. Thus, all arguments are equally accepted. This shows that the support relation does not play any role, and a supported argument is as strong as a non-supported one. Formally, these approaches violate strengthening principles, which capture the role of supports. The approaches developed in [START_REF] Brewka | [END_REF][START_REF] Nouioua | Bipolar argumentation frameworks with specialized supports[END_REF]] also return a single extension when the attack relation is empty. In case of acyclic graphs, the extension is the set of all arguments. Thus, these approaches violate strengthening and the support relation may not be fully exploited in the evaluation of arguments.

Gradual semantics were introduced for the first time in [Cayrol and Lagasquie-Schiex, 2005a]. The authors presented some properties that such semantics should satisfy (like a particular case of strengthening). However, they did not define concrete semantics. To the best of our knowledge, the first gradual semantics was proposed in [START_REF] Evripidou | Valentinos Evripidou and Francesca Toni. Quaestio-it.com: a social intelligent debating platform[END_REF]. It was extended to QuAD in [START_REF] Baroni | [END_REF] for evaluating arguments in acyclic graphs.

Definition 4 Let A = A, w, R, S ∈ WAG. A is acyclic iff the following holds: for any non-empty finite sequence a = a 1 , a 2 , . . . , a n of elements of A, if ∀i ∈ {1, 2, . . . , n -1}, a i , a i+1 ∈ R ∪ S, then a n , a 1 ∈ R ∪ S.

We define next restricted semantics. All principles for semantics are straightforwardly adapted to restricted semantics.

Definition 5 A restricted semantics is a function S transforming any acyclic

A = A, w, R, S ∈ WAG into a function Deg S A from A to [0, 1].
QuAD is a restricted semantics assigning a numerical value to every argument on the basis of its intrinsic strength, and the overall strengths of its attackers and supporters. It evaluates differently supporters and attackers before aggregating them. QuAD was recently extended to the restricted semantics DF-QuAD in [START_REF] Rago | [END_REF]. Unlike QuAD, it uses the same function for aggregating supporters and attackers separately. It satisfies Strict Franklin axiom, thus it treats equally attacks and supports. It violates Strengthening and Weakening in presence of attackers/supporters of degree 1. However, the semantics avoids the irrational behavior of QuAD.

Proposition 5 DF-QuAD violates Strict Monotony, Strict Reinforcement, Resilience, Weakening, and Strengthening. It satisfies all the remaining principles.

Both (QuAD, DF-QuAD) suffer from a big jump problem. Consider the graph depicted in Figure 1. For x being QuAD or DF-QuAD, Deg x A (i) = 0.991. Note that the value of i makes a big jump from 0.10 to 0.991 and i becomes even stronger than its supporter j. There are two issues with such jump: First, the gain is enormous. Assume that i is a fallacious argument. The supporter may increase slightly its strength but cannot fully rehabilitate it. Second, such jump impedes the discrimination between different cases where w(i) > 0.1 since whatever the value of w(i), the overall strength is almost 1.

Euler-based Graded Semantics

We have seen that no existing semantics satisfies all our principles. We define a novel restricted semantics satisfying them all for the subclass of acyclic non-maximal graphs. Without loss of generality, the basic strengths of arguments are considered less than 1. Note that few arguments are intrinsically perfect. The probability of false information, exceptions, etc., is rarely 0. In contrast, the loss of cyclic graphs is important. But, we consider that the class of all acyclic non-maximal graphs is expressive enough to deserve attention. Definition 6 A restricted semantics is a function S transforming any acyclic non-maximal A = A, w, R, S ∈ WAG into a function Deg S A from A to [0, 1]. Before presenting our semantics, we need to introduce a relation between arguments based on the longest paths to reach them (mixing support and attack arrows). Definition 7 (Well-founded relation) Let A = A, w, R, S be an acyclic BAG and a ∈ A. A path to a in A is a non-empty finite sequence a = a 1 , a 2 , . . . , a n of elements of A such that a n = a and ∀i ∈ {1, 2, . . . , n -1}, a i , a i+1 ∈ R ∪ S. We denote by Rel(A) the well-founded binary relation ≺ on A such that ∀x, y ∈ A, x ≺ y iff max{n | there exists a path to x of length n} < max{n | there exists a path to y of length n}. Since A is acyclic, those maximum lengths are well-defined, so is Rel(A).

We are ready to define the Euler-based restricted semantics. The general idea is to take into account supporters and attackers in an exponent E of e (the Euler's number). More precisely, the stronger or more-numerous the supporters, the greater and more-likely-positive that exponent. Obviously, the inverse is true with the attackers. Then, the overall strength of an argument a is naturally defined as w(a)e E . Finally, we need certain tweaking (including a double polarity reversal) to make our function a restricted semantics in the first place, and to have it satisfy all our axioms. Definition 8 We denote by Ebs the restricted semantics such that for any acyclic non-maximal A = A, w, R, S ∈ WAG, Ebs(A) is the function f from A to [0, 1] recursively defined with Rel(A) as follows: ∀a ∈ A,

f (a) = 1- 1 -w(a) 2 1 + w(a)e E , E = x∈Supp(a) f (x)- x∈Att(a) f (x).
Theorem 1 Let A = A, w, R, S ∈ WAG be acyclic nonmaximal and a ∈ A. The following holds: Note also that being supported by a very strong argument does not cause a weak argument to become extremely strong, which shows that Ebs does not suffer from the big jump problem. Indeed, Deg Ebs A (i) = 0.22 and thus the jump is not big. Theorem 2 Ebs satisfies all the 12 principles.

Deg Ebs A (a) = 1 - 1 -w(a) 2 1 + w(a)

Conclusion

The paper presented for the first time principles that serve as guidelines for defining semantics for weighted bipolar argumentation graphs. It also analyzed existing semantics with regard to the principles. The results revealed that extension semantics violate key principles, and do not fully exploit support relations. Existing gradual semantics satisfy more but not all the principles. We proposed a novel semantics which satisfies all the principles. This work can be extended in several ways. An urgent future work consists of defining a semantics that satisfies the principles while dealing with any typology of graphs including thus cyclic ones. Another perspective consists of defining semantics that take strict precedence to attacks over supports. Finally, new principles can be defined for capturing other features of semantics like those proposed very recently in [Baroni et al., 2018;[START_REF] Mossakowski | Bipolar weighted argumentation graphs[END_REF].

Proposition 1

 1 Let S be a semantics that satisfies Bi-variate Independence, Bi-variate Directionality, Stability, and Strict Franklin. For any A = A, w, R, S ∈ WAG, for any a ∈ A, if there exists a bijective function f from Att A (a) to Supp A (a) such that for any x ∈ Att A (a), Deg S A (x) = Deg S A (f (x)), then Deg S A (a) = w(a).

Proposition 4

 4 QuAD violates Strict Monotony, Strict Reinforcement, Resilience, (Strict) Franklin, Weakening, and Strengthening. It satisfies all the remaining principles. As a consequence of violating Weakening and Strengthening, QuAD may behave irrationally. Consider a weighted bipolar graph where A = {a, b 1 , b 2 , b 3 }, w(b 1 ) = w(b 2 ) = 0.8, w(b 3 ) = 0.9, R = {(b 2 , a), (b 3 , a)}, and S = {(b 1 , a)}. Thus, a has an attacker and a supporter of equal strengths (Deg S A (b 1 ) = Deg S A (b 2 ) = 0.8), and an additional attacker b 3 . If w(a) = 0.2, then Deg S A (a) = 0.422 meaning that the single supporter is privileged to the two attackers. However, if w(a) = 0.7, Deg S A (a) = 0.477 meaning that attacks are privileged to support. More generally, we can show that if w(a) ≥ 0.5, then Deg S A (a) < w(a), else Deg S A (a) > w(a). Hence, choosing which of support and attack should take precedence depends on the intrinsic strength of an argument.

Example 1

 1 Consider the graph depicted in Figure 1. Every circle contains [argument name]:[intrinsic strength] and below [overall strength].The neutrality principle can be checked with g and e, stability with d, g, h, and j, bivariate monotony with a and b, bivariate reinforcement with b and c, Imperfection with i, Strict Franklin with a, weakening with e.g. b, and strengthening with i.

  Figure 1: Bipolar graph

This is an abridged version of a paper titled "Evaluation of arguments on weighted bipolar graphs" which won a best-paper award at ECSQARU-17 conference.

Acknowledgments

This work was supported by ANR-13-BS02-0004 and ANR-11-LABX-0040-CIMI.