Louis Rivière

Christian Artigues

Azeddine Cheref

Nicolas Jozefowiez
email: nicolas.jozefowiez@univ-lorraine.fr

Marie-José Huguet

Sandra U Ngueveu
email: ngueveu@laas.fr

Vincent Charvillat
email: vincent.charvillat@enseeiht.fr

Multi-product, multi-supplier order assignment and routing for an e-commerce application in the retail sector

Keywords: E-commerce, retail market, order assignment, vehicle routing, genetic algorithms

With the rise of virtualization, the share of e-commerce in the retail market continues to grow in an omnichannel context. We consider an existing software tool, developed by the Devatics company, for pooling inventories in stores to meet online orders. The problem which arises therefore consists in seeking the optimal allocation of a set of customers to stores. In this paper we consider a variant of the offline problem corresponding to an evolution of the existing software, consisting of assigning a set of predefined orders when the transportation cost depends on a delivery tour to the customer locations. We show that the problem corresponds to a vehicle routing problem with additional but standard attributes. A mixed-integer linear programming formulation is given and several heuristics are proposed : a giant tour-based genetic algorithm, a simple cluster-first, route second heuristic and an assignment-based genetic algorithm. Preliminary computational results on a set of realistic problem instances suggest that the assignment-based genetic algorithm better scales as the problem size increases.

INTRODUCTION

With the rise of virtualization, the share of ecommerce in the retail market continues to grow in an omnichannel context. One of the services offered by the Devatics company is Onestock 1 , a tool for pooling inventories to meet online orders. The problem which arises therefore consists in seeking the optimal allocation of a set of customers to stores. Two modes of assignment are possible. Indeed, we can consider the "Online" assignment mode, consisting of the allocation of each order as they are declared, and the "Offline" assignment which consists of assigning all the orders in a single large block. In this paper we consider the offline problem. The Onestock software solves a variant where the transport cost are fixed. In this paper we consider an evolution of the problem towards a variant where the transport costs depend on delivery tours to the customer locations. We propose an mixed-integer linear programming formulation (MILP) of the problem. On realistic data 1 https://www.onestock-retail.com/ instances, we compare several heuristics and metaheuristics. We show that a cluster-first, route-second based genetic algorithm obtains the best results. The problem formulation is given in Section 2. The realistic data extraction approach that we use to generate the data instances is then presented in Section 3. Section 4 first gives a quick state of the art review of efficient metaheuristics for vehicle routing problems and propose adaptations for the considered E-commerce problem. Section 5 provides a computational comparison of the proposed exact and heuristic approaches. Concluding remarks are drawn in Section 6.

MODELING THE PROBLEM

The problem considers a fixed set of online orders D for a set of products P. Each order d ∈ D asks for an amount q d p of product p ∈ P. We have a set M of stores, and each store m ∈ M has a stock s mp of product p ∈ P. The problem is then to meet the demand in products of each order by using the store stocks with the possibility to split the demand across different stores. Then, each store must deliver the ordered packages to the customers in a single tour.

Although in the OneStock software the orders are served online upon receipt, the offline problem has concrete applications. Indeed, most real orders take place in the evening, and therefore cannot be processed before the stores open on next morning. In this case, we have a given number of orders to process at a time. Another application of interest is the estimation of a "regret", i.e. the difference in performance between the immediate allocation of successive orders and the optimal allocation of these orders. In this paper, we wish to minimize the distance traveled on tours. The trucks have no capacity, so each store only has to deliver the products in at most one tour. However, a customer can be served by several stores and therefore by several routes. This model therefore explicitly incorporates the "assignment + routing" aspect.

Input Data

D : set of orders. M : set of stores. P : set of products.

N = M ∪ D : set of nodes. V = {(i, j) | i, j ∈ N 2 } set of arcs. R : set of r i j ∈ R : cost of arc from i to j Q : set of q d p ∈ N : amount of product p ∈ P in order d ∈ D. S : set of s mp ∈ N : stock of product p ∈ P in store m ∈ M.

Variables

X : set of x dmp ∈ N : number of products p sent from store m for order d. Y : set of y dm ∈ {0, 1} : indicates whether a package is sent from store m for order d. Z : set of z m i j ∈ {0, 1} : indicates whether the tour of store m takes arc (i, j). U : set of u m i ∈ N : number assigned to node i in the tour of store m (to avoid subtours).

Objective

min ∑ m∈M ∑ (i, j)∈V z m i j r i j (1) 2.3.1 Constraints ∑ m∈M x dmp ≥ q d p , ∀d, p ∈ D, P (2)
∑ d∈D x dmp ≤ s mp , ∀m, p ∈ M, P

x dmp ≤ y dm min(q d p , s mp ), ∀d, m, p ∈ D, M, P (4)

∑ i∈N z m id ≥ y dm , ∀d, m ∈ D, M (5)
∑ i∈N z m mi ≥ y md , ∀d, m ∈ D, M (6)
∑ i∈N z m i j = ∑ k∈N z m jk , ∀m, j ∈ M, N (7)
u m i -u m j + |D|z m i j ≤ |D| -1, (8) ∀i, j ∈ V \ {m}, i = j, ∀m ∈ M
Objective (1) is this time to minimize the total distance traveled by the store deliverers.

Constraints (2) guarantees that each order is sufficiently supplied. Constraints (3) guarantee that stocks are sufficient for shipments and constraints (4) ensure that the number of products p sent from store m to customer d is set to 0 if no package is sent from m to d. Constraints ( 5) and ( 6) guarantee that the tour of store m passes through customer d and by store m if this allocation has been decided, respectively. Constraints (7) ensure flow conservation while (8) guarantees that only cycles including store m are allowed.

REALISTIC DATA INSTANCE GENERATION

To best stick to reality, we use statistics extracted from real instances of a Devatics customer. An offline problem generator uses these statistics to create datasets with an arbitrary number of orders.

Figures 1 and2 show that for the firm considered the distribution of orders is not egalitarian.Wednesday is the busiest day with around 200 orders on average, with peaks of orders on sales days that do not exceed 1000 orders for a day.

For each product, its popularity is given by the sum of its purchases. In addition, a matrix of cooccurrence indicates for each pair of products the number of orders that include both. We notice that the popularity of the products is distributed in a way approaching a Pareto distribution.

We can therefore generate a semi-realistic n-sized order from the data, first choosing a product randomly based on its popularity, then adding products according to their co-occurrences with the products already present until we have n products. We also have information regarding the size distribution of n orders with an almost binomial distribution of the quantity of product ordered. We notice that the average order gathers 2.8 products, or roughly 1% of the available products (203 in total).

Once the orders are composed, we need to provide stores to ensure the feasibility of the problem. To do this, we are inspired by the overall distribution of the stock of products as well as quantities of stocks of each product in each store. We noticed that more than 30% of the stock is in the warehouse, while all the other stores roughly share the rest.

For each unit of product ordered we select a store to receive a unit of stock of this product. If it is the first unit of this product in stock, a stock margin of two units is added. This margin value, given by Devatics, corresponds to the risk aversion of stores that report a quantity of stock below reality to avoid short-age. We generate a set of problems constituting a test bench on which we vary the values of the number of orders and of the stock margin.

PROPOSED METHODS

The modeling of the routing problem, in particular the constraints that eliminate subtours, do not allow a MILP solver to solve the problem effectively, which justifies the use of heuristics. This problem is very close to several well studied problems such as the traveling salesman problem (TSP) the location and routing problem (LRP) and in particular the vehicle routing problem (VRP). Classically, the latter consists of minimizing the cost of visiting all customers without exceeding a given capacity of the truck. One can consider that our problem corresponds to the VRP with additional classical attributes (no capacity on trucks, multi-products, capacity on stocks, several warehouses, multiple deliveries, one truck per warehouse). The VRP as well as many variations involving some of these attributes are very studied problems for their many practical applications in the field of logistics. In [START_REF] Abdulkader | Vehicle routing problem in omni-channel retailing distribution systems[END_REF][START_REF] Martins | A savings-based heuristic for solving the omnichannel vehicle routing problem with pick-up and delivery[END_REF], heuristics are proposed for a related omnichannel retail problem with the difference that vehicles are not attached to stores but located in a central depot and perform pick up and delivery tours. We first review the efficient metaheuristics designed for the VRP before describing our solution approaches.

The properties of efficient metaheuristics for vehicle routing problems

Vidal et al. [START_REF] Vidal | A unified solution framework for multi-attribute vehicle routing problems[END_REF] proposed a unified heuristic framework to solve ta family of VRP with many attributes. The approach consists in categorizing the attributes (additional constraints) of the VRP to be able to offer solutions adapted to each of these categories. In addition, others works [START_REF] Vidal | Heuristics for multi-attribute vehicle routing problems: A survey and synthesis[END_REF] try to extract from the profusion of the methods the main characteristics which make the success of an approach. Although empirical, this analysis can guide us in creating our method. Two tracks inspired by the methods of LRP solution approach emerge. On the one hand, the giant tour method, which focuses on the "routing" aspect of the problem, and on the other, an approach highlighting the "assignment" aspect of the problem.

Giant tours: Route-first, cluster-second approach

The giant tour method is proposed for the first time for the VRP by Beasley [START_REF] Beasley | Route first -cluster second methods for vehicle routing[END_REF]). Prins [START_REF] Prins | A simple and effective evolutionary algorithm for the vehicle routing problem[END_REF] integrates this technique with a memetic algorithm (genetic algorithm with a local search phase). Vidal et al. [START_REF] Vidal | A unified solution framework for multi-attribute vehicle routing problems[END_REF] among others proposed an efficient variant able to solve a large set of routing problems with various constraints.

In this section we first recall the principle of the giant tour and the SPLIT operator for the VRP. Then, we describe the adaptation of the SPLIT operator to our problem and we finally describe the proposed giant tour based hybrid genetic algorithm that uses a MILP-based repair operator.

The principle of the giant tour

The representation by giant tour consists for the classic VRP to consider, instead of explicit solutions, "giant tours", a kind of concatenation of real tours, representing all the ways to cut this tour in order to respect the capacity constraint of the trucks. The idea is that if we know how to quickly find the best solution for this subset, it is faster to consider "macro solutions". An illustration of the giant tour is given in Fig. 3. In this figure, a giant tour of size 3: [3,1,2] must be cut for a truck capacity of 2. The different possible cuts are therefore [3 -1 -2], [3,1 -2] and [3 -1,2]where "-" represents a return to the depot. Finding an optimal division of the giant tour into subtours is polynomial, as a shortest path algorithm. Figure 4 shows how this cutting operator (commonly called SPLIT) determines the optimal solution. It consists in finding the shortest path in a graph whose nodes are the ordered points visited by the giant tour and whose arcs each represent a grouping. Each arc cannot group more points than the capacity of the trucks and the cost associated with an arc is the cost of the grouping it represents.

Adapting the split operator to our problem

To adapt the method to our problem and its attributes, we first observe the adaptations considered for these attributes taken individually. In the liter- ature, to adapt the method to problems with splitdelivery, [START_REF] Boudia | An effective memetic algorithm with population management for the split delivery vehicle routing problem[END_REF] brings two changes:

0 3 1 2 c [3] c [3,1] c [1] c [1,2] c [2]
• A node can appear several times in a giant tour.

• Each occurrence of a node is associated with the quantity of delivered products. The sum of the delivered quantity must match the request.

The SPLIT operator is also slightly modified since when a sub-tour is considered, visiting several times the same node is meaningless. The sum of the delivered quantities is therefore carried over to a single occurrence of the node. Experience shows that the choice of this node, if it can be difficult to determine optimally, can be done deterministically (the local search carried out thereafter rectifying a possible bad choice). In addition to the loss of optimality, the first of these changes complicates the cutting, in fact, in a giant tour, a node can appear up to M times if it is served by all stores. We switch from a fixed size |D| to a variable size up to |D|.|M|. Especially since in a multi-product context and with capacities on stocks, the information associated with the quantity delivered for each product can be voluminous and, as we will see later, does not by itself guarantee the existence of a solution by the SPLIT operator. Indeed, even if the sum of the delivered products corresponds to the demand for each product, it is still necessary that stores have stocks to deliver these products. This feasibility of a sub tour is itself difficult to determine since it depends on the other sub-routes selected.

Vidal et al. (Vidal et al., 2014a) show how the giant tour can be adapted to multiple depots. The change is restricted to the SPLIT method, in which the costs associated with a sub-tour becomes the cost of allocating the tour to the best depot. In our case with capacities, it is not possible to determine the "best" assignment independently of other assignments. We must therefore create for each assignment of a subtour to a store m, an arc with a corresponding cost. The number of route assignments to be considered becomes O(|M| n ) in the worst case where n is the size of the giant tour since we do not limit the capacity of the trucks.

A SPLIT operator adapted to the case of deliveries with capacities is presented by Duhamel et al. [START_REF] Duhamel | A GRASPxELS approach for the capacitated locationrouting problem[END_REF]. The authors use a Bellman algorithm with several labels per node to keep the stocks available when searching for the shortest path. A dominance rule between labels makes it possible to discard some of them, but the method remains too slow for large instances. Several improvements are proposed to speed it up, in particular, limiting the number of labels considered during the evaluation (a parameter to be set). Despite these complications, one of the properties of our problem, which can simplify the SPLIT procedure, is that each store performs only one tour. In this case, the resource whose labels must report the usage is no longer the stock, but the prior use or not of a store. However, the giant tour cut this way represents a less important set of solutions. So we have a trade-off between speed of cutting and search efficiency.

From this analysis, a giant tour for our problem is thus a sequence S of at most |D|.|M| nodes, each node k corresponding to an order δ(k) associated to a delivered amount λ p (k) such that the sum of delivered amount for the nodes corresponding to the same order equals the total required quantity for each product, which correspond to the following invariant:

|S| ∑ k=1 δ(k)=d λ p (k) = q d p ∀d ∈ D, p ∈ P (9)
The SPLIT operator we propose is synthesized in Algorithmm 1. A solution is represented by a chain of labels. a label L labels is composed of a node L.nodein the giant tour S, a score L.score, a store L.store and a set of available stores L.available m, and a parent label L.parent. A Label represent the subsequence from the successor of L.parent to L.node assigned to store m of cost L.score and such that the set of stores available for delivering the successors of L.node in S is stored in L.available m . The notation µ i j denotes the sub-tour comprising nodes i + 1 up to j in S. To save the propagation of unnecessary labels, the procedure ADDLABELTONODE uses a standard dominance rule. A label L 1 dominates a label L 2 (noted L 1 L 2 ) if it uses fewer resources for a better score. Formally: One must notice that the SPLIT operator may fail in obtaining a feasible set of tours. Indeed, it is possible that all active label at some points fail to satisfy the condition at Line 8. This condition express the fact that the stock in store m must be sufficient to deliver all amounts λ p (k) of each order k in the µ i j subsequence for all products p. But if the λ p (k) atttached to the giant tour satisfy invariant (9), which ensure that the customer demand is satisfied, they are not necessarily compatible with the stock. In this case, the best set of tours computed by the SPLIT operator is incomplete.

L 1 .score < L 2 .score and L 2 .available m ⊆ L 1 .available m or L 1 .score ≤ L 2 .score and L 2 .available m ⊂ L 1 .available m

The giant tour-based hybrid genetic algorithm

Given the SPLIT operator, Algorithm 2 describes the genetic algorithm. Individuals are giant tours. In Line 1, a population of giant tours is initialized randomly via function INITPOPGT. The a standard tournament selection is performed to select two parent individuals I1 and I2 (Line 3, function CHOOSEP-ARENTSGT). An adaptation of the classic one point crossover operator is used on the selected giant tours I1 and I2 at Line 4 with function CROSSOVERGT, which yields offspring I3. Recall that a giant tour is defined by a sequence of nodes k and associated λ p (k) delivered values for each product p ∈ P. For the nodes up to the crossover point, λ p (k) values in I3 are the same as in I1. After the crossover point, the λ p (k) values for the nodes k duplicated from I2 are set so that invariant (9) is satisfied. If there is not enough nodes for a given order to satisfy its demand, nodes are duplicated at the end until the invariant is satisfied. Then, at Line 5, the SPLIT operator is applied to obtain a set T of tours (one for each store).

R T mk d + R dT m,k+1 -R T mk T m,k+1
is the insertion cost of d (in the case where it is the only order inserted in the tour on machine m) and c - md represents symmetrically the profit of removing d from T m . Then function MILPREPAIR, use the standard best insertion algorithm to build the route of each store given the orders assigned to the stores. A two-opt local search algorithm (function LOCALSEARCH at Line 7) further improves tours T . Finally, the giant tour offpring I3 is rebuilt by deriving the giant tour from the sequence of tours on T (function GETGT). The SELECTSURVIVORSGT keeps the best giant tours in the population up to a maximal number α. The DIVERSIFYGT diversification function simply reinitializes part of the population randomly.

Cluster-first, route second approaches

Another common approach, in opposition to the giant tour method, first selects the store assignments to customers, then chooses the best routes for visits.

The idea is that once the allocation fixed, the problem is reduced to an instance of the TSP for each store having to visit customers to whom it was assigned.

In the worst case, each customer should be visited, but it can be expected that in general a store will only visit a small number of customers; and, due to the combinatorial nature of the problem, solving a set of small problems is much faster than solving a single bigger problem.

Cluster-first, route-second heuristic

In an article by Fisher [START_REF] Fisher | A generalized assignment heuristic for vehicle routing[END_REF], a heuristic of this type is proposed for the first time for a variant of the VRP and we adapt this method to our problem. It solves a general assignment problem to determine the assignments and then solves the resulting TSP. Intuitively, the problem of allocation resembles the problem of clustering; indeed, if the sets of customers to be visited are located around the stores, we can hope that the associated routing is of good quality. Among all the existing clustering criteria, we choose to minimize the sum of the distances between customers and their associated stores (potentially several) such that the allocation is valid. This criterion has the merit of being simple and of adapting well to the strong constraints on stocks. Other criteria such as minimizing the greatest distance between a customer and a linked store are not adapted since if an assignment happens to be necessary (for example in the case where only one store has a given product) yielding a high cost A, the assignment of other customers to stores becomes indifferent to the cost as long as it is less than A.

Once the assignment of orders to stores is made, we use the the MILP (1-8) with fixed variables y dm and x pdm , which amounts to solve |M| small TSPs.

Hybrid genetic assignment algorithm

In addition to the previously mentioned heuristic, another genetic algorithm based on solving the underlying assignment problem is proposed (Algorithm 3). The solutions (or individuals) are represented by a set of |M| chromosomes, each corresponding to a store tour. The population is initialized randomly (INITPOP), then, at each iteration, two individuals are selected to be crossed (function CHOOSEPARENTS). The classic one point crossover operator is used on each of the store tour pairs since it has the advantage of keeping part of the parents' route ( CROSSOVER). Each individual offspring is not not always viable because of the product stock in the stores and the MILPREPAIR function is then used to repair the offspring.

The offspring is then inserted into the population. If the population size exceeds a given threshold, a selection procedure (SELECTSURVIVORS) determines which individuals are deleted. To choose the survivors, we compare of course the solution scores, but we also take into account their diversity, using a measure based on the average Hamming distance between an individual and his closest neighbors. When the search stagnates during a given number of iterations, the population is diversified (DIVERSIFY). The operation consists in destroying the current solutions by removing certain visits before applying a random repair of individuals. To do this, we use a random variant of MILPREPAIR.

PRELIMINARY COMPUTATIONAL RESULTS

We compare the four solution methods (MILP using the IBM Cplex solver, Cluster-first heuristic, Assignement-based Genetic algorithm, Giant tourbased Genetic algorithm) on randomly generated problem instances for given dimensions. Figure 5 presents the performances of the methods compared to the best solution among the four. A single 15 min run was performed for each algorithm. Each point correspond to the average score on 3 instances of the given size. So a score of 1 corresponds to the best solution found and a score of 0.5 corresponds to a so- As you would expect, The MILP formulation, even if it finds the optimal solution for small instances, is quickly outperformed when it comes to solving larger instances. The assignment-based genetic (named AG affectation in Fig. 5) algorithm is more efficient than the one based on giant tours, especially for the larger instances with a realistic size of 100 orders, 10 stores and 10 products. In third position, the assignment heuristic obtains surprisingly honorable scores (about 1.25 times higher than the best score). As parameter tuning was handcrafted and only single runs were performed, these preliminary computational experiments should be extended in the future but they suggest good scaling properties of the assignment-based genetic algorithm.

CONCLUSION

We have considered an industrial problem problem consisting in finding the optimal allocation and routing of a set of customer orders to stores, a variant of the location routing problem / vehicle routing problem with complicating constraints. We proposed mixed-integer linear programming formulations and several heuristics for the problem. Our hybrid genetic algorithm based on assigning customers to store first and routing the order second with MILP-based repair operators showed good scaling properties on preliminary computational experiments. Regarding the routing problem, further improvements can be made in the configuration of genetic algorithms. In particular, the Split operator could be improved as in (Duhamel et al., 2011). An idea to experiment would be not to fix the quantities of products, but that each label

Figure 1 :

 1 Figure 1: Number of orders per day over a three-month period

Figure 3 :

 3 Figure 3: Giant tour [3,1,2] and associated solution.

Figure 4 :

 4 Figure 4: Subgraph for cutting by SPLIT the giant tour [3,1,2] with capacity 2.

Algorithm 1 :

 1 The SPLIT procedure 1 //Insert start label; 2 L ← (node = 0; cost = 0; available store = M; parent = nil); 3 ADDLABELTONODE(0, L); 4 for i ∈ {0, ..,

Algorithm 2 :

 2 Giant tour-based genetic algorithm 1 Pop = InitPopGT (n); If tour T is incomplete (see reasons above), a MILP-based repairing operator is used (Line 6). Based on tour T , the repair operator searches to compute a new feasible tour for each store with optimized cost. Let ȳ the current assignment of orders to stores as stated by tour T such that ȳdm = 1 if store m delivers at least a part of order d and ȳdm = 0 if no part of order d is delivered by store m. A new feasible assignment of orders to stores is first computed by solving a variant of multi order, multi product facility location problem issued from MILP (1-8) while replacing the routing constraints by estimated removal profits and insertion costs. Let T m denote the tour on machine m and T mk for k = 1, . . . |T m | the kth order visited by tour T m . The repair MILP has variables y dm and x dmp ∈ {0, 1} for all d ∈ D, m ∈ M, d ∈ D and the following objective and constraints. k=1,...,|T m |-1

Algorithm 3 :

 3 Assignment based genetic algorithm 1 Pop = INITPOP(n);

Figure 5 :

 5 Figure 5: Performance of the different methods compared to the best solution

Acknowledgements

This research was carried out in the context of the "One Stock Performance" project, funded by the EASYNOV program (FEDER / Midi-Pyrénées region), in partnership with DEVATiCS.

We also thank the anonymous referees for their enlightening suggestions.


Pauline Rivoire
email: pauline.rivoire@giub.unibe.ch

Olivia Martius

Philippe Naveau

A comparison of moderate and extreme ERA-5 daily precipitation with two observational data sets

The timing and the intensity of daily precipitation are assessed over Europe

Introduction

Natural hazards related to extreme precipitation (river floods, flash floods, landslides, debris flows and avalanches) cause casualties, damages to infrastructures and buildings and have direct and indirect economic impacts [START_REF] Munichre | A stormy year: Natural Catastrophes 2017[END_REF]. For infrastructure planning and prevention measures, information about rare events, e.g., events that occur on average only once in a hundred years, is important. Such information can be obtained from precipitation data with statistical tools. Assessing the accuracy in high quantiles depends on spatial domain sizes and temporal availability. Different types of global precipitation data sets are available (Sun et al., 2018): global precipitation data sets are based on ground observations, satellite observations, combinations of ground observations and satellite observations and on short term weather model forecasts in reanalyses data sets. Reanalyses combine past observations with weather forecast models to reconstruct past weather. The main advantage of this type of precipitation data set is its regular spatial and temporal coverage. Reanalyses ensure consistency of the precipitation data with the atmospheric conditions, which is important for weather and climate process studies. Here, we focus on ERA-5 precipitation (C3S, 2017). ERA-5 is the latest reanalysis product from the European Centre for Medium-Range Weather Forecasts (ECMWF). ERA-5 precipitation is computed in short-term forecast started from reanalysis initial conditions [START_REF] Hennermann | Era5: data description[END_REF]. The ERA-5 precipitation production process does not include precipitation observation inputs. Hence comparison with observational data makes sense, keeping in mind that observation data have (partly substantial) uncertainties as well (Sun et al., 2018;[START_REF] Kulie | Uncertainties in microwave properties of frozen precipitation: Implications for remote sensing and data assimilation[END_REF][START_REF] Prein | Impacts of uncertainties in European gridded precipitation observations on regional climate analysis[END_REF].

ERA-5 precipitation has already been widely used since its release in 2018, but very few assessments of this data set have been conducted over large regions. Only precipitation over restricted areas and precipitation associated with specific type of events have been assessed (Wang et al., 2018;[START_REF] Hénin | Assessing the use of satellite-based estimates and high-resolution precipitation datasets for the study of extreme precipitation events over the Iberian Peninsula[END_REF][START_REF] Mahto | Does ERA-5 Outperform Other Reanalysis Products for Hydrologic Applications in India[END_REF]Tarek et al., 2020;[START_REF] Amjad | Performance evaluation of satellite-and model-based precipitation products over varying climate and complex topography[END_REF]. The goal of this article is to assess daily precipitation in ERA-5 against observational data sets over large regions: Europe, comparing with the station-based data set EOBS [START_REF] Haylock | A European daily high-resolution gridded data set of surface temperature and precipitation for 1950-2006[END_REF], and the entire globe, comparing with the satellite-based data set CMORPH (Joyce et al., 2004). The verification of daily precipitation will focus on the intensity distribution and on the temporal consistence, i.e., the co-occurrence of events. We use the extended generalized Pareto distribution (Tencaliec et al., 2019) to evaluate the intensity distribution and we calculate co-occurrence hit rates to assess the joint occurrence of precipitation events. This paper is structured as follows. Section 2 describes the data used for this study.

We introduce methods used for the comparison of co-occurrence and intensity in section 3. The results of our analysis are presented in section 4. Finally, the results are summarized and discussed in section 5.

Data

ERA-5 Precipitation

Reanalysis precipitation in this study are extracted from ERA-5 reanalysis data set. ERA-5 is the latest global reanalysis data set provided by the European Center for Medium-Range Weather Forecasts (C3S, 2017;[START_REF] Hersbach | The ERA5 global reanalysis[END_REF]. In this data set, precipitation stem from short-term forecasts and are available at an hourly resolution that we aggregate to daily precipitation. The precipitation data calculation does not rely on observed precipitation. The data is interpolated to a regular grid with 0.25 • resolution.

Observation-based Data Sets

The two gridded observation-based precipitation data sets used in this study are EOBS [START_REF] Haylock | A European daily high-resolution gridded data set of surface temperature and precipitation for 1950-2006[END_REF] that is based on European station observations and CMORPH that is based on satellite observations (Joyce et al., 2004).

The EOBS data set is provided by the European Climate Assessment & Dataset and is a daily gridded data set based on spatially interpolated station data. The version used is 19.0e, with a 0.25 • by 0.25 • grid. The interpolation to a 0.25 • by 0.25 • grid is a combination of monthly precipitation totals and daily anomalies products. Figure 1a in [START_REF] Cornes | An Ensemble Version of the E-OBS Temperature and Precipitation Data Sets[END_REF] displays the station coverage for version 16.0. This coverage is heterogeneous, with a very dense network in Ireland, the Netherlands, Germany, Switzerland, and northern Italy, for example, and very few stations in northern Africa, in the Middle East, in Iceland, in Norway, and in Sweden. EOBS covers land precipitation only, and the comparison with ERA-5 is conducted for the time period between January 1979 and December 2018.

The second observational data, CMORPH, is provided by the National Center for Atmospheric Research (NCAR) (Climate Prediction Center, National Centers for Environmental Prediction, National Weather Service, NOAA, U.S. Department of Commerce, 2011). This gridded precipitation product combines passive microwave satellite scans and geostationary satellite infrared data and provides 3 hour accumulations that we aggregate in daily accumulation. CMORPH stands for climate prediction center morphing method, the name of this combination technique. The precipitation estimation algorithm of this data set is not able to capture snow (Joyce et al., 2004). The spatial resolution of the gridbox is also 0.25 • . The comparison with ERA-5 is conducted for the period 2003-2016, for latitudes between 60 • S and 60 • N.

The two observation-based data sets have the same grid resolution as ERA-5 but a shift of 0.125 • in latitude and longitude is present for the coordinates of the grid points compared to ERA-5.

Data Processing

The study evaluates seasonal precipitation for September, October, November (SON); December, January, February (DJF); March, April, May (MAM); and June, July, August (JJA). Separation between seasons ensures stationarity of the time series. The in-tensity distribution analyses are based on wet days, defined as days with precipitation accumulations exceeding 1 mm. The 1 mm threshold corresponds to standard recommendations for station data [START_REF] Hofstra | Testing E-OBS Eu-teorology[END_REF] and eliminates potential drizzle effect in reanalysis data [START_REF] Maraun | Bias correction, quantile mapping, and downscaling: Revisiting the inflation issue[END_REF]. The co-occurrence analysis is conducted on the entire seasonal time series, including days with precipitation lower than 1 mm.

The precipitation time series are not de-trended in the present study as the response of precipitation to increasing atmospheric CO 2 varies with the precipitation intensity [START_REF] Pendergrass | Changes in the distribution of rain frequency and intensity in response to global warming[END_REF] and trends depend on the length of time series [START_REF] Scherrer | Emerging trends in heavy precipitation and hot temperature extremes in Switzerland[END_REF].

Moreover, [START_REF] Donat | Consistency of temperature and precipitation extremes across various global gridded in situ and reanalysis datasets[END_REF] identified mostly small or insignificant trends for the past thirty years.

Methods

For the sake of simplicity, in the method section OBSER denotes the observationbased data sets. OBSER can be either EOBS or CMORPH, as the comparison procedures between ERA-5 and EOBS and between ERA-5 and CMORPH are identical.

Co-occurrence of Precipitation Events

Binary events are defined here as occurrences of daily precipitation above the P th seasonal percentile with P ∈ {75, 90, 95, 99}. Percentile values can be different between ERA-5 and OBSER. In Figure 1, the 95 th precipitation percentiles in SON is displayed for: (a) ERA-5 over the EOBS domain , (b) EOBS over its entire domain (1979-2018), (c) ERA-5 over the CMORPH domain (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016) and (d) CMORPH over the entire domain (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016).

We define a co-occurrence between two data sets when two exceedances occur either at the same grid point on the same day, at the same grid point with one day lag, or at one of the eight surrounding grid points on the same day. During the spatial and temporal shift, a single event is never used more than once when looking for co-occurrences.

We allow for one day shift to bypass uncertainties that arise having a fixed 24h time window [START_REF] Haylock | A European daily high-resolution gridded data set of surface temperature and precipitation for 1950-2006[END_REF]. The extension to the eight grid points around the centre point addresses potential issues arising from the precipitation interpolation to the different grids.

The hit rate is the ratio between the number of joint events and the total number of events [START_REF] Rhodes | Can reanalyses represent extreme precipitation over England and Wales[END_REF]. The total number of events is the same if computed from ERA-5 or from OBSER.

Intensity Assessment

Extreme value theory is often used in hydrology and climate sciences (e.g., [START_REF] Lamb | Confidence intervals for a spatially generalized, continuous simulation flood frequency model for Great Britain[END_REF][START_REF] Cooley | Bayesian spatial modeling of extreme precipitation return levels[END_REF]Tramblay et al., 2013;[START_REF] Kang | Parameter and quantile estimation for the generalized Pareto distribution in peaks over threshold framework[END_REF]. This approach states that peaks over high thresholds, i.e., amounts of rain exceeding a given threshold u, may be approximated by a generalized Pareto distribution, provided the threshold and the number of observations are large enough and some additional mild conditions are satisfied (see section 3.2 for generalized Pareto distribution definition). However, the generalized Pareto distribution fitting has drawbacks. First, it only captures the upper tail behavior. A distribution combining gamma behavior for small and moderate precipitation amounts with generalized Pareto distribution behavior for high amounts can be a solution. Second, a threshold has to be determined for every station or grid point to separate the upper tail from the rest of the distribution [START_REF] Dupuis | Exceedances over high thresholds: A guide to threshold selection[END_REF]. To overcome these challenges, here we use the extended generalized Pareto distribution (EGPD) [START_REF] Naveau | Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection[END_REF]Tencaliec et al., 2019).

To study wet day precipitation intensity distributions, this section presents a comparison of quantiles and a homogeneity test based on the Kullback-Leibler divergence.

Both parts rely on our EGPD fit.

For the intensity comparison, we discard grid points where the number of wet days is smaller than 500 days for the comparison with EOBS and smaller than 200 days for the comparison with CMORPH. Moreover, auto-correlation can be present in daily time series, for example when two consecutive wet days are fostered by the same weather system (e.g., Lenggenhager et al., 2019a). To address the auto-correlation in time, we consider that two precipitation events separated by two days are independent [START_REF] Barton | Clustering of Regional-Scale Extreme Precipitation Events in Southern Switzerland[END_REF][START_REF] Lenggenhager | Atmospheric blocks modulate the odds of heavy precipitation events in Europe[END_REF][START_REF] Fukutome | Automatic threshold and run parameter selection: a climatology for extreme hourly precipitation in Switzerland[END_REF]. To ensure independence of the time series, the intensity assessment is conducted on one-third of the data that is randomly drawn. This approach is a trade-off between keeping enough data in the sub-samples to ensure robust fitting and best removing the auto-correlation in the data.

Extended General Pareto Distribution

In extreme value theory, one way to model the extremal tail behavior is the so-called peak-over-threshold approach (see, e.g. [START_REF] Coles | An introduction to statistical modeling of extreme values[END_REF][START_REF] Katz | Statistics of Extremes in Climatology and Hydrology[END_REF]. Under this framework, rainfall exceedances above a large threshold u are assumed to follow a Generalized Pareto distribution defined as

H ξ (z) = 1 -(1 + ξz) -1/ξ + if ξ = 0, 1 -e -z otherwise, (1)
where the positive scalar σ represents a scale parameter and the real ξ drives the upper tail behavior. A negative, null and positive ξ corresponds respectively to the "bounded", "light" and "heavy" tail case, i.e. an upper tail that is bounded for ξ < 0, with exponential decay for ξ = 0 or polynomial decay for ξ > 0. The selection of the threshold u is not trivial for large data sets, as each grid point may need a different optimal threshold (e.g., [START_REF] Deidda | A multiple threshold method for fitting the generalized Pareto distribution to rainfall time series[END_REF]. A large threshold implies a small sample size of extremes and consequently, large uncertainties in the estimation of σ and ξ. Conversely, a moderate threshold leads to a possible incorrect approximation by a Generalized Pareto distribution, i.e. a large model error. To bypass this complex threshold selection step, [START_REF] Naveau | Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection[END_REF] proposed a simple scheme to smoothly transition between the main body of the distribution and its upper tail, while keeping the constraint of modeling extremes with a Generalized Pareto distribution. The proposed model can be written as

F (x) = G{H ξ (x/σ)}, for all x > 0, (2)
where G, the transition function, is a continuous cumulative distribution function (cdf) in the unit interval. By imposing the two constraints, lim u s is finite and positive for some s > 0, the new cdf F (.) is bound to be in compliance with extreme value theory for its lower and upper tails. This class of cdf is called extended generalized Pareto distribution (EGPD) family. In this study, the cdf G(.) is estimated using a specific Bernstein polynomial approximation, more information can be found in Tencaliec et al. (2019). The R code is available upon request.

Quantile Confidence Intervals

Confidence intervals for the quantiles of ERA-5 and OBSER precipitation are computed using a semi-parametric bootstrap on EGPD fitting. For each grid point, the following bootstrap procedure is conducted. Two subsamples containing one-third of the time series each are randomly drawn from the initial wet day time series. Each of these two subsamples is bootstrapped 100 times each, and each bootstrapped sample is fitted by a EGPD. Having our disposal 200 bootstrapped estimates of G(.), σ and ξ, quantiles for any given non-exceedance probability can be computed from Eq. ( 2). In particular, 95% confidence intervals of the quantiles are obtained by calculating the empirical 2.5% and 97.5% quantiles of the 200 bootstrapped quantiles values. An important feature to assess the proximity of our different data sets is to check if the confidence intervals from ERA-5 overlap (or not) with the observational data sets.

Kullback-Leibler Divergence Test

The well-known Kullback-Leibler divergence used in various fields "measures" the distance between two probability density functions, say f 1 and f 2 . It is given by equation

E f1 log f 1 (X) f 2 (X) + E f2 log f 2 (Y) f 1 (Y) (3)
with X and Y being random variables following respectively the probability density functions f 1 and f 2 .

Let X ERA-5 =(X i ) i=1,...,n and Y OBSER =(Y j ) j=1,...,m be the time series (after removing the auto-correlation) of wet day precipitation in ERA-5 and OBSER. With f1 and f2 estimated by the EGPD fitting, we obtain the empirical value of the Kullback-Leibler divergence with equation

1 n n i=1 log f1 (X i ) f2 (X i ) + 1 m m j=1 log f2 (Y j ) f1 (Y j ) . (4)
The null hypothesis of our test is "X ERA-5 and Y OBSER have the same distribution", i.e. f1 = f2 . The alternative hypothesis is f1 = f2 .

The distribution of the Kullback-Leibler divergence under the null hypothesis is estimated using 300 values of the divergence between two vectors randomly drawn from a concatenation of X ERA-5 and Y OBSER . The probability of "The Kullback-Leibler divergence between X ERA-5 and Y OBSER is greater than the Kullback-Leibler divergence under the null hypothesis" is the p-value of the test. This p-value is empirically determined from the 300 values of the Kullback-Leibler divergence under the null hypothesis. The null hypothesis is rejected with a confidence level of 5% if the p-value is greater than 0.05.

Difference in Number of Wet Days

The intensity comparison is based on the EGPD, which is fitted to wet days only.

Discrepancies in the number of wet days between ERA-5 and the observational data sets could have an impact when comparing the EGPD fitted to ERA-5 and the observational data sets. To quantify these discrepancies at a fixed grid point, we use two simple measures. The first measure is the ratio of the seasonal number of wet days, defined by :

N wd ERA-5 N wd OBSER (5)
where, for a fixed season, N wd ERA-5 and N wd OBSER are the number of wet days in ERA-5 and OBSER, respectively.

The second measure is the absolute value of the difference in the number of wet days between ERA-5 and OBSER, normalized by the time series length of OBSER, given by:

100 × N wd ERA-5 -N wd OBSER N wd OBSER . (6)
Note that the absolute difference quantifies the distance between the ratio of the number of wet days and 1.

Results

Number of Wet Days

Table 1 presents the mean absolute value of the difference in the number of wet days.

The differences are computed only over grid points retained for the intensity comparison, i.e. with time series longer than 200 days for CMORPH and 500 days for EOBS.

Over Europe, the mean absolute difference is between 11% (SON) and 21% (MAM)

of EOBS number of wet days. In SON and DJF, the number of wet days is lower in ERA-5 than in EOBS in northern Europe and higher in southern Europe. In MAM and JJA the number of wet days is almost always higher in ERA-5 than in EOBS (see Figure A1a in the appendix for a map of the ratio of the number of wet days in SON). The difference in number of wet days between ERA-5 and EOBS can be considered as low and will not have impact on the EGPD fitting.

The global comparison with CMORPH reveals larger discrepancies in the seasonal number of wet days than with EOBS. The mean difference in the number of wet days corresponds to between 66% (DJF) and 76% (SON) of the CMORPH number of wet days.

This difference is mainly due to the ERA-5 wet days being more numerous than in CMORPH.

The number of wet days in ERA-5 is twice as large as in CMORPH for 19% (DJF) to 25% (JJA) of grid points. Figure A1b in appendix displays the map of the ratio of number of wet days in SON. This ratio is greater than 2 over bands at fixed latitudes e.g. over bands close 60 • S, 20 • S, 20 • N and 60 • N.

Co-occurrence of Precipitation Events

In the comparison between ERA-5 and both observation-based data sets, the hit rate decreases with increasing intensity of the events and is similar across the seasons (Table 2). Grid points with a given percentile of less than 1 mm are not considered.

Over Europe the average hit rate between ERA-5 and EOBS for the 75 th percentile is between 73 % (in JJA) and 77 % (in SON), i.e. about three quarters of the events exceeding the 75 th percentiles coincide. For the 95 th percentile, the mean hit rate is between 53% (JJA) and 60% (SON). For the 99 th percentile, the hit rate varies between 39 % and 45 % depending on the season. The global mean hit rate is of the same order of magnitude as for Europe. The mean hit rate between ERA-5 and CMORPH for the 95 th percentiles is above 50%. The mean hit rate associated with the 99 th percentile is between 35% and 37% depending on the season.

Maps of the hit rate for the 95 th percentile can be found in the appendix (Figure B1 and Figure B2). The spatial pattern does not strongly depend on the season or the percentile. For Europe, the hit rate has a large variability near arid regions (Maghreb and Turkey). The rest of Europe is quite homogeneous. A lower hit rate is observed in Iceland and southern Italy. For the global comparison, the best hit rate is reached over the oceans in the mid-latitudes. The hit rate is substantially lower in Eastern China, along the equator, in South America and in tropical Africa.

Intensity Verification

Confidence Intervals on Quantiles

The confidence interval overlap between ERA-5 and EOBS is independent of the probability of non-exceedance, i.e., the intensity of the events. Figure 2a shows the relative position of the 95% confidence intervals for quantiles associated with probability of non-exceedance 0.9 in SON, between ERA-5 and EOBS. A grid point is displayed in yellow if the confidence intervals are overlapping, in orange if the upper boundary of ERA- Quantiles with probability of 0.3 have a larger number of grid points with disagreement in JJA. ERA-5 quantiles for probability 0.3 during JJA are underestimated compared to EOBS quantiles for a major part of Europe (not shown).

5
The global comparison of ERA-5 with CMORPH shows less overlap of the confidence intervals with increasing precipitation intensity, as we will see in section 5. Figure 2b displays the relative position of the 95% confidence intervals for quantiles associated with probability of non-exceedance 0.9 in SON, between ERA-5 and CMORPH and Figure 3b shows the number of seasons with an overlap of the confidence intervals for quantiles with non-exceedance probability 0.9 between ERA-5 and CMORPH. The confidence intervals for non-exceedance probabilities of 0.3 and 0.5 overlap for more than 85% of the grid points. For probabilities 0.75 to 0.95, in a band along the equator the confidence intervals do not overlap for all the seasons for many grid points, see e.g. The Kullback-Leibler test between ERA-5 and CMORPH has a clear signal of agreement in the mid-latitudes and disagreement in the tropics, for all seasons. The summary over all seasons mainly informs about intensity agreement over the oceans, because of the time series length constrain removes most land grid points. Figures 4b and4d present the p-value for JJA and SON, and show over land the same general pattern of disagreement in the tropics and agreement in the mid-latitudes. One exception to this pattern is the disagreement over mountainous regions of the mid-latitudes (western North America, Himalayas, South Chile), in agreement with the results in section 4.3.1.

Summary and Discussion

The analysis of precipitation event co-occurrence between ERA-5 and EOBS and ERA-5 and CMORPH reveals a decreasing agreement with increasing intensity of events, independently of the season.

Key results of the intensity comparison of ERA-5 with EOBS over Europe depend on the season (Table 3). Quantiles in MAM and SON show a good agreement for all nonexceedance probabilities p. Indeed, between 81% (for p = 0.3 in MAM and for p = 0.5 in DJF) and 90% (for p = 0.5 in MAM and SON, and for p = 0.75 in SON) of the grid points have overlapping confidence intervals. The percentages of grid points where the distributions agree (Kullback-Leibler test) are highest for MAM (34%) and SON (39%).

In DJF, the agreement between quantiles is between 69% (p = 0.9) and 82% (p = 0.3).

In JJA agreement is much better for high quantiles (up to 94% for p = 0.95) than low ones, the confidence intervals for quantiles with probability p = 0.3 overlap for only 39% of grid points. This discrepancy for low precipitation intensity has an impact on the Kullback-Leibler test: the null-hypothesis could not be rejected for only 10% of grid points in JJA.

The other seasons show a higher fraction of grid points for which the null-hypothesis was not rejected (between 29% in DJF and 39% in SON).

The study of the wet day precipitation distribution (Kullback-Leibler test) between ERA-5 and EOBS over Europe reveals a robust agreement over regions where the station coverage of EOBS is dense. Areas with the largest differences in the distribution are areas with thin station coverage, e.g. in southern Europe and Russia. [START_REF] Cornes | An Ensemble Version of the E-OBS Temperature and Precipitation Data Sets[END_REF] highlighted that "station coverage is the most important factor in determining the success of the gridded data". In areas with sparse station data, the precipitation is interpolated from distant stations [START_REF] Hofstra | Testing E-OBS Eu-teorology[END_REF]. Additionally, extreme precipitation is smoothed by the spatial interpolation (Hofstra et al., 2010), which justifies our quantiles larger in ERA-5 than in EOBS for extreme precipitation.

The intensity comparison between ERA5 and CMORPH indicates a decreasing agreement between the two data sets with increasing precipitation intensity (Table 3). The percentage of grid points with confidence intervals overlapping is between 92% (JJA) and 94% for p = 0.3, and between 70% (SON) and 75% (DJF) for p = 0.9. One exception is the slightly better agreement of quantiles for extreme precipitation (non-exceedance probability p = 0.95) than for moderately extreme precipitation (p = 0.9), with an overlap rate 1% higher, for all seasons. This can be explained by confidence intervals becoming larger for larger quantiles, and there is thus a higher chance of overlap between ERA-5 and the observational data set. This remark holds also for the EOBS results. The Kullback-Leibler test presents little variation with season, like for the quantile study. Between 52% (SON) and 57% (DJF) of the grid points studied did not reject the null-hypothesis of the test.

The analysis of the entire precipitation distribution reveals proportionally more grid points agreeing between ERA-5 and CMORPH than between ERA-5 and EOBS (see last row of Table 3). This can be due to the longer time series in EOBS leading to a stricter test. Another explanation can be that there are proportionally more challenging regions for a model over Europe, with the Alps for example, whereas globally the largest regions compared are oceans, where the agreement is good in general.

The global comparison of the wet day precipitation distributions between ERA-5 and CMORPH over the period 2003-2016 shows a rather good agreement in the midlatitudes, and a strong disagreement over the tropics. This result is robust over the seasons and does not depend on the method used. The band along the equator where precipitation intensities are lower in ERA-5 compared to CMORPH corresponds to a region with a ratio of the number of wet days rather close to 1 (see Figure A1). The number of wet days does not differ substantially in this region and therefore does not play a role in the robust disagreement between ERA-5 and CMORPH over the tropics. In this region, ERA-5 quantiles are lower than CMORPH ones, especially for non-exceedance probabilities larger than 0.75. This feature has already been observed by [START_REF] Pfahl | Quantifying the Relevance of Cyclones for Precipitation Extremes[END_REF] with ERA-interim, another reanalysis product from ECMWF. They computed the empirical 99 th percentiles of 6-hourly precipitation in ERA-interim and CMORPH for the period 2003-2016 (Figure 2 in their article). The authors showed a strong underestimation of ERA-interim precipitation over the tropics compared to CMORPH. They concluded that the deep convection, a central process in tropical extreme precipitation, was not properly captured by the reanalysis data set. Even though Nogueira (2020) found an improvement of the precipitation simulation over the tropics in ERA-5 compared to ERA-interim the previous reanalysis data set from ECMWF, we assume ERA-5 to still contain an underestimation of the tropical extreme rainfall.

One limit of CMORPH that it is important to emphasize is the limitations of this data set to capture snow (Joyce et al., 2004) and low-intensity precipitation during winter in the mid-latitudes (Sun et al., 2018). This property leads to smaller number of wet days in winter (DJF or JJA depending on the hemisphere) as seen in section 4.1, hence the large areas where our intensity analysis can not be conducted. Some regions at high latitudes (close to 60 • S and 60 • N) and some mountainous regions have a number of wet days large enough for the intensity comparison to be conducted even if snow can be expected. However a substantial difference in the number of wet days is observed. The For the period 2014-2018, [START_REF] Amjad | Performance evaluation of satellite-and model-based precipitation products over varying climate and complex topography[END_REF] showed that ERA-5 overestimates the precipitation observed over Turkey compared to ground based stations, independently of the wetness and slope classes. Our comparison with EOBS indicates the same signal for high quantiles, but also a dry bias for lower quantiles. This can be due to the fact that our study period is longer or that EOBS has a poorer station coverage in this region.

In their study, Tarek et al. ( 2020) compared the mean seasonal precipitation in ERA-5 with station observations over North America between 1979 and 2018. In JJA, they found an underestimation of precipitation in ERA-5 over Florida, and an overestimation along the west coast of Canada, which is in agreement with our comparison of quantiles between ERA-5 and CMORPH for this season and for all probabilities of non-exceedance.

In DJF, they showed an underestimation precipitation in ERA-5 over the west coast of USA and Florida, and an overestimation over the west coast of Canada. Our analysis highlighted that ERA-5 presents larger quantiles than CMORPH over the west coast of North America for all non-exceedance probabilities. Our results are thus in agreement for Canada but not for USA. This can be due to the fact that the time periods studied are different and that CMORPH underestimates precipitation during the cold months (Sun et al., 2018).

Mahto and Mishra (2019) assessed ERA-5 precipitation in India against observation comparing precipitation sums during the monsoon season (June-September) between 1980 and 2018. They found a wet bias over Indo-Gangetic Plain and foothills of Himalaya and a dry bias in semi arid regions of western India. These results are in agreement with our quantile analysis in ERA-5 and CMORPH in JJA for the period 2003-2016.

Conclusion

We compare daily precipitation from the ERA-5 reanalysis data set with daily precipitation from two observation-based data sets, EOBS and CMORPH. The comparison addresses three aspects i) the temporal co-occurrence of moderate to high extreme events in two data sets, ii) the agreement of return values for moderate to extreme nonexceedance probabilities derived from the extended generalized Pareto distribution (EGPD), and iii) a comparison of the full precipitation distribution captured by the EGPD using the Kullback-Leibler divergence. We quantify the co-occurrence of precipitation events with the hit rate. We compare the EGPD distributions between ERA-5 and the observational data sets with confidence intervals for several non-exceedance probabilities and with a test based on the Kullback-Leibler divergence.

Between ERA-5 and EOBS over Europe the hit rate is above 65% for moderate precipitation and approximately 50% for extreme precipitation. Between ERA-5 and CMORPH globally the hit rate is above 60% for moderate precipitation and around 40% for extreme precipitation. Over Europe areas with the least agreement are the southern Mediterranean region and Iceland and for the global comparison areas with the least agreement are land areas between 15 • S and 15 • N, North-West America and Central Asia.

For a majority of grid points confidence intervals for non-exceedance probabilities of 0.3 to 0.95 overlap between ERA-5 and EOBS. We find a disagreement between ERA-5 and EOBS in areas where EOBS uses fewer input stations. We therefore hypothesize that the reanalysis data set might better capture moderate to extreme precipitation in regions where the station coverage is sparse. The analysis also showed that ERA-5 underestimates extreme precipitation compared to CMORPH in the tropics. In general, the magnitudes of the non-exceedance probabilities agree between ERA-5 and the observationbased data sets in the mid-latitudes.

The Kullback-Leibler test on the entire precipitation distributions over Europe shows an agreement of the EGPD distributions in ERA-5 and EOBS over Germany, Ireland, Sweden and Finland. The precipitation distributions differ significantly between in ERA-5 and EOBS in all four seasons in Iceland, Norway, Karelia, Hungary and the Balkan.

The Kullback-Leibler test between ERA-5 and CMORPH shows that precipitation distributions are generally in agreement over the mid-latitudes and differ significantly over the tropics for all seasons, confirming the results of the quantile comparison. ERA-5 should only be used with great care to study extreme precipitation over the tropics. Note: the mean absolute value of the difference in the number of wet days is defined in equation 6, and is computed here for grid points with more than 500 wet days for ERA-5 vs EOBS and for grid points with more than 200 wet days for ERA-5 vs CMORPH. Note: for a given percentile, the mean is computed over all grid points where the precipitation percentile is larger than 1 mm. See section 3.1 for the definition of the hit rate. with non-exceedance probability 0.9 between (a) ERA-5 andEOBS (1979-2018) and between (b) ERA-5 andCMORPH (2003-2016). See section 3.2.2 for computational details. Grid points with an insufficient number of wet days (see section 3.2) are discarded and displayed in white.

-19-manuscript submitted to Earth and Space Science

  confidence interval is lower than the lower boundary of the EOBS confidence interval, and in blue if the lower boundary of ERA-5 confidence interval is larger than the upper boundary of the EOBS confidence interval. Figure3ashows the number of seasons with a confidence intervals overlap for quantiles with non-exceedance probability 0.9 between ERA-5 and EOBS. The confidence intervals overlap during all the seasons for a major part of Europe. The exceptions are Iceland, Norway and Western Russia, Romania, the Adriatic sea coast and some grid points in the Alps. Non-overlapping confidence intervals correspond primarily to an underestimation of the quantiles by ERA-5 for low precipitation intensities (not shown), and overestimation for large intensities (Figure2a.).

  . For all seasons and between 15 • S and 15 • N, ERA-5 quantiles are smaller compared to CMORPH. Another disagreement that deserves to be highlighted is the higher ERA-5 quantiles compared to CMORPH over the mountainous regions of the west coast of North America, the south of Chile, in Papua New Guinea, the Himalayas and the Alps, regardless of the non-exceedance probabilities and seasons.4.3.2 Comparison of the Full Distributions Using the Kullback-Leibler TestThe Kullback-Leibler divergence of the full precipitation distribution points to regions of agreement and disagreement independent of the seasons, for both the comparison between ERA-5 and EOBS, and ERA-5 and CMORPH, see Figure 4. Figures 5a and 5b display the number of seasons for which the p-value of the Kullback-Leibler test is greater than 0.05, i.e. where the EGPD distributions fitted to ERA5 and to EOBS and CMORPH do not differ significantly. ERA-5 and EOBS wet day precipitation intensities agree best over Germany, Ireland, Sweden and Finland. Wet day precipitation intensity follows the same distribution in ERA-5 and EOBS for most grid points in these countries. Regions with the least agreement, i.e. where the null hypothesis is rejected for all the seasons, are Iceland, Norway, Hungary and the Balkans. The area with at least one season where the null hypothesis is rejected is rather large. The Kullback-Leibler test gives weight to differences in the entire distribution, thus the low intensity precipitation disagreement in JJA mentioned previously in section 4.3.1 has an impact on the Kullback-Leibler divergence (see Figure 4a and Figure 4c for the plots of the p-value in JJA and SON). Note the pattern following the border between Norway and Sweden, and Finland and Russia: a very good agreement is observed in Sweden and Finland, whereas the null hypothesis is rejected for almost all seasons in Norway and Karelia.

  higher ERA-5 precipitation compared to CMORPH over mountainous regions might be related to snow.Timmermans et al. (2019) revealed a disagreement between CMORPH and gauge-based product for extreme precipitation in mountainous regions of western USA in DJF and interpreted it as a consequence of the post-processing performed in CMORPH leading to lots of missing data in winter(Xie et al., 2017).[START_REF] Hénin | Assessing the use of satellite-based estimates and high-resolution precipitation datasets for the study of extreme precipitation events over the Iberian Peninsula[END_REF] assessed ERA-5 daily accumulated precipitation during extreme precipitation events over the Iberian Peninsula for the period 2000-2008 against precipitation from a ground based gridded data set. They found an overestimation of daily sums for moderate extreme events and and underestimation for the most extreme events. Our study of quantiles in ERA-5 and EOBS for non-exceedance probabilities greater than 0.75 reveals a moderate signal of overestimation of ERA-5 precipitation in the same region. One exception is the southern Basque Country where an underestimation of ERA-5 quantiles for these probabilities in DJF. Our comparison of moderate and large extremes over the period 1979-2018 is therefore only partially in agreement with their study over the period2000-2008.

Figure 2 .

 2 Figure 2. Relative position of the confidence intervals (CIs) for SON quantiles associatedwith non-exceedance probability 0.9 between (a)ERA-5 and EOBS (1979-2018) and between (b)ERA-5 and CMORPH (2003-2016). See section 3.2.2 for computational details. Grid points with an insufficient number of wet days (see section 3.2) are discarded and displayed in white.

Figure 3 .

 3 Figure 3. Number of seasons with overlapping confidence intervals for quantiles associated

Figure 4 .

 4 Figure 4. p-value of the Kullback-Leibler divergence test (as defined in section 3.2.3) between ERA-5 and EOBS (1979-2018) in JJA (a) and SON (b) and between ERA-5 and CMORPH (2003-2016) in JJA (c) and SON (d). p-values ¡0.05 indicate that the distributions differ significantly. Grid points with an insufficient number of wet days (see section 3.2) are discarded and displayed in white.

Figure 5 .

 5 Figure 5. Number of seasons where the distributions are similar, i.e., without rejection of the null hypothesis of the Kullback-Leibler test (as defined in section 3.2.3) (a) between ERA-5 and EOBS (1979-2018) and (b) between ERA-5 and CMORPH (2003-2016). Grid points with an insufficient number of wet days (see section 3.2) are discarded and displayed in white.

Figure A1 .

 A1 Figure A1. Ratio of the number of wet days as defined with Eq. (5) in SON between (a)ERA-5 and EOBS (1979-2018) and between (b)ERA-5 and CMORPH (2003-2016). Grid points with an insufficient number of wet days (see section 3.2) are discarded and displayed in white.

Figure B1 .

 B1 Figure B1. Hit rate for events greater than the 95 th percentile between ERA-5 and EOBS in (a) SON, (b) DJF, (c) MAM and (d) JJA.

Figure B2 .

 B2 Figure B2. Hit rate for events greater than the 95 th percentile between ERA-5 and CMORPH in (a) SON, (b) DJF, (c) MAM and (d) JJA. Grid points with an insufficient number of wet days (see section 3.2) are discarded and displayed in white.

Table 1 .

 1 Sun, Q., Miao, C., Duan, Q., Ashouri, H., Sorooshian, S., & Hsu, K.-L. Interim near surface air temperature and precipitation over Arctic sea ice: Effects on sea ice thermodynamics and evolution. The Cryosphere Discussions, 1-28. doi: 10.5194/tc-2018-245 Xie, P., Joyce, R., Wu, S., Yoo, S. H., Yarosh, Y., Sun, F., & Lin, R. (2017). Reprocessed, bias-corrected CMORPH global high-resolution precipitation es-Mean Absolute value of the Difference in the Number of Wet Days

	ropean high-resolution gridded data set of daily precipitation and surface
	temperature.	Journal of Geophysical Research Atmospheres, 114 (21).	doi:
	10.1029/2009JD011799
	Hofstra, N., New, M., & McSweeney, C.	(2010).	The influence of interpola-
	tion and station network density on the distributions and trends of climate
	variables in gridded daily data.	Climate Dynamics, 35 (5), 841-858.	doi:
	10.1007/s00382-009-0698-1

Joyce, R. J., Janowiak, J. E.,

Arkin, P. A., & Xie, P. (2004)

. CMORPH: A Method that Produces Global Precipitation Estimates from Passive Microwave and Infrared Data at High Spatial and Temporal Resolution. Journal of Hydrome-

Table 2 .

 2 Mean Hit Rate ERA-5 vs EOBS and ERA-5 vs CMORPH

	Percentile		EOBS				CMORPH
		DJF	MAM	JJA	SON	DJF	MAM	JJA	SON
	75 th	76%	75%	73%	77%	74%	73%	72%	73%
	90 th	66%	65%	61%	67%	62%	61%	60%	60%
	95 th	59%	58%	53%	60%	53%	52%	51%	52%
	99 th	44%	45%	39%	45%	37%	35%	35%	36%

manuscript submitted to Earth and Space Science

The strengths of ERA-5 daily precipitation data are the regular spatial and temporal resolution and the consistency with the large-scale circulation and there is generally a good agreement with observation-based data sets in the extra-tropics. The reanalysis data set provides valuable complementary information to observational data in regions where observational data sets are sparse, e.g. in areas where the EOBS station coverage is poor or for CMORPH in regions and seasons where snow is prevalent. In the tropics, an observational data set should be preferred over ERA-5. The authors declare that they have no conflict of interest.

The authors acknowledge the E-OBS data set from the EU-FP6 project UERRA (http://www.uerra.eu) and the Copernicus Climate Change Service, and the data providers in the ECA&D project (https://www.ecad.eu). EOBS daily precipitation data at resolution 0.25°over Europe lands are available on this link: https://surfobs.climate .copernicus.eu/dataaccess/access eobs.php Variable rr, ensemble mean. One must be a registered EOBS user and be signed in to download this dataset.

For ERA-5 we use global total daily precipitation data at resolution 0.25°. Hourly values were downloaded from the ECMWF MARS server (a valid ECMWF account required): https://apps.ecmwf.int/data-catalogues/era5/?type=fc&class=ea&stream= oper&expver=1 Forecast steps 6 to 17, variable: Total precipitation (228.128). The MARS / EMOSLIB interpolation library has been used.

CMORPH global daily precipitation rate data at resolution 0.25°have been downloaded from https://rda.ucar.edu/datasets/ds502.0/ (registration on https://rda .ucar.edu/ required) Variable: CMORPH precipitation estimate.

The data sets to reproduce the figures and tables can be found in Zenodo (https:// doi.org/10.5281/zenodo.4443804).

The codes for the intensity assessment are available from GitHub (https://github .com/PauRiv/characterization ERA-5 daily precipitation). The codes for the cooccurrence verification are available upon request.

manuscript submitted to Earth and Space Science Note: for a given non-exceedance probability p, the percentage denotes the proportion of grid points for which the confidence intervals are overlapping. For the whole distribution, the percentage denotes the proportion of grid points where the null-hypothesis can not be rejected, i.e. where the distributions are similar.

Appendix B Hit Rate