F Palacios-Rodríguez
email: fatima.palacios@ucm.es.

Keywords: extreme-value theory, precipitation, risk analysis, space-time Pareto processes, stochastic simulation, hydrology

To better manage the risks of destructive natural disasters, impact models can be fed with simulations of extreme scenarios to study the sensitivity to temporal and spatial variability. We propose a semi-parametric stochastic framework that enables simulations of realistic spatio-temporal extreme fields using a moderate number of observed extreme space-time episodes to generate an unlimited number of extreme scenarios of any magnitude. Our framework draws sound theoretical justification from extreme value theory, building on generalized Pareto limit processes arising as limits for event magnitudes exceeding a high threshold. Specifically, we exploit asymptotic stability properties by decomposing extreme event episodes into a scalar magnitude variable (that is resampled), and an empirical profile process representing space-time variability. For illustration on hourly gridded precipitation data in Mediterranean France, we calculate various risk measures using extreme event simulations for yet unobserved magnitudes, and we highlight contrasted behavior for different definitions of the magnitude variable.

Introduction

Extreme events of geophysical processes such as precipitation extend over space and time, and they can entail devastating consequences for human societies and ecosystems. Flash floods in Southern France constitute highly destructive natural phenomena causing material damage and threatening human lives [START_REF] Vinet | Flash flood-related mortality in southern France: first results from a new database[END_REF], as for example during two relatively recent catastrophic flash-flood events in the Gard department in September 2002 [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 near Montpellier in October 2014 [START_REF] Brunet | Retour d'expérience sur les crues des 6 et 7 octobre 2014 à montpellier-grabels (hérault, france) : caractéristiques hydro-météorologiques et contexte historique de l'épisode[END_REF]. Since damage and costs of floods have been increasing over the last decades, the understanding of temporal and spatial variability of rainfall patterns generating such floods receives considerable attention from the authorities (European Environment Agency, 2007). To help with this understanding, we develop a method to stochastically simulate realistic spatio-temporal extreme scenarios, which can be fed to impact models. Examples of impact models are urban flood models, such as the shallow water models of [START_REF] Guinot | Flux and source term discretization in two-dimensional shallow water models with porosity on unstructured grids[END_REF] and [START_REF] Guinot | Flux closures and source term models for shallow water models with depth-dependent integral porosity[END_REF], which produce hydrological variables such as water height or water speed, based upon which experts make decisions about flood risk. Spatial extreme-value theory (EVT) has become a very active field of research and proposes databased stochastic modeling of such extreme events for predicting probabilities, risks and uncertainty behavior [START_REF] Coles | An Introduction to Statistical Modeling of Extreme Values[END_REF][START_REF] De Haan | Extreme Value Theory. An Introduction[END_REF]Davison et al., 2012;[START_REF] Ferreira | The generalized Pareto process; with a view towards application and simulation[END_REF][START_REF] Davison | Statistics of extremes[END_REF]. Throughout this paper, the notion of space will refer to planar geographic space (R 2 ) if not stated otherwise, in contrast to space augmented by the time dimension (i.e., R 2 × R). In physical processes, very complex deterministic patterns (e.g., physical laws, or spatio-temporal trends in the statistical sense) and probabilistic patterns (variability at relatively small spatial and temporal scales) may arise. In combination with the high dimension of datasets, realistic spatio-temporal modeling is challenging. In this work, we develop a data-driven non-parametric spatio-temporal resampling approach by transforming observed quantiles in a spatially and temporally coherent way, i.e., by appropriately accounting for spatio-temporal extremal dependence. We illustrate our method on a high-dimensional data set of gridded hourly reanalysis data. Our procedure draws sound justification from asymptotic theory for threshold exceedances with a strong probabilistic interpretation. We will explain how it allows us to flexibly define extreme episodes in space-time data based on different ways of aggregating marginal return periods over space and time. Recall that return levels are quantiles at a given probability level q ∈ (0, 1) and correspond to an event magnitude that is exceeded with probability 1 -q. Then, 1/(1 -q) is the associated return period.

The block-maxima and peaks-over-threshold (POT) methods are two widely known strategies in univariate EVT to identify extreme events in a data set. While the block-maxima method is based on the partition of the observation period into equally-sized blocks (for instance months or years) to extract the maximum observation in each of these subperiods [START_REF] Ferreira | On the block maxima method in extreme value theory: PWM estimators[END_REF], the POT method consists in the study of exceedances above a given high threshold [START_REF] Pickands | Statistical inference using extreme order statistics[END_REF][START_REF] Embrechts | Modelling extremal events for insurance and finance[END_REF][START_REF] Beirlant | Statistics of Extremes: Theory and Applications[END_REF]. Max-stable processes, introduced by de Haan (1984), are the natural infinite-dimensional generalization of the univariate generalized extreme value (GEV) distribution, which constitutes the only limiting distribution of block maxima when block sizes increase to infinity. [START_REF] Ferreira | The generalized Pareto process; with a view towards application and simulation[END_REF] and [START_REF] Dombry | Functional regular variations, Pareto processes and peaks over threshold[END_REF] showed that generalized Pareto processes are the only possible asymptotic limits for threshold exceedances of stochastic processes. Both approaches are closely linked through theoretical tail stability properties. In practice, the application of block-maxima or POT methods may yield very different approaches and results, as illustrated by [START_REF] Wi | Non-stationary frequency analysis of extreme precipitation in South Korea using peaks-over-threshold and annual maxima[END_REF]. Using stochastic process models defined over continuous space allow us to capture fine-scale behavior of extreme events, and it bypasses the implementation of ad-hoc partitioning methods of space prior to applying univariate or multivariate extreme-value models in alternative approaches; see [START_REF] Mornet | Wind storm risk management: sensitivity of return period calculations and spread on the territory[END_REF]; [START_REF] Carreau | Partitioning into hazard subregions for regional peaksover-threshold modeling of heavy precipitation[END_REF] for examples of space-partitioning. [START_REF] Davison | Geostatistics of extremes[END_REF] remark that max-stable models with appropriately chosen dependence specification enable successful spatial modeling of extremes. Indeed, max-stable processes have been widely applied [START_REF] Davison | Geostatistics of extremes[END_REF][START_REF] De Fondeville | High-dimensional peaks-over-threshold inference[END_REF][START_REF] Tyralis | Estimation of intensitydurationfrequency curves using max-stable processes[END_REF]Shin et al., 2019, and references therein), and several approaches to stochastic simulation of spatial max-stable fields have been proposed [START_REF] Dombry | Conditional simulation of max-stable processes[END_REF][START_REF] Dombry | Exact simulation of max-stable processes[END_REF][START_REF] Oesting | Spatial modeling of drought events using max-stable processes[END_REF]Oesting et al., 2018a,b). Since max-stable processes are linked to the block-maxima approach, their realizations often aggregate information of several of the underlying original events, which may limit the physical interpretation of the simulated processes. Consequently, modeling based on max-stable processes fitted to block maxima data, and the resulting simulations, seem appropriate for studying longterm events. On the other hand, generalized Pareto processes represent the original events that satisfy a threshold exceedance condition. For specific choices of marginal distributions, they can be represented constructively by multiplying a random scaling variable with a so-called spectral process, the latter characterizing the spatial variation in the extreme events [START_REF] Ferreira | The generalized Pareto process; with a view towards application and simulation[END_REF][START_REF] Dombry | Functional regular variations, Pareto processes and peaks over threshold[END_REF][START_REF] Thibaud | Efficient inference and simulation for elliptical Pareto processes[END_REF][START_REF] Opitz | Semi-parametric resampling with extremes[END_REF].

In current practice, one usually first fits a parametric model for the spectral processes, and the estimated values of dependence parameters are then plugged in to define the model used for simulation. de Fondeville and Davison (2018) propose high-dimensional parametric POT inference for data with margins normalized in a first step; they develop estimation based on proper scoring rules instead of likelihood expressions, since the normalizing constants in the latter become unwieldy in moderately high dimension (say, more than 50 spatial locations). They illustrate the potential and flexibility of their approach by modelling extreme rainfall with different purely spatial risk functionals (also called cost functionals in the following); a spatio-temporal variant of the risk functional is used by the same authors in de Fondeville and [START_REF] De Fondeville | Functional peaks-over-threshold analysis[END_REF].

In contrast, we here develop an algorithm for extracting observed spectral processes from data, and we then combine them with newly sampled scaling variables to generate new realizations of the extreme events. Since extreme events are frequently spatio-temporal in nature (i.e., they span over several time steps, often with stochastic temporal dependence), their spatial extension and duration have to be accounted for properly. In recent work, [START_REF] Chailan | A semiparametric method to simulate bivariate space-time extremes[END_REF] build on the semi-parametric approach of [START_REF] Caires | On the determination of the temporal and spatial evolution of extreme events[END_REF] and the simulation technique of [START_REF] Ferreira | The generalized Pareto process; with a view towards application and simulation[END_REF] and propose a semiparametric method with low computational cost to simulate extreme bivariate spatio-temporal fields of wave heights in the Gulf of Lions along the French Mediterranean coast. The notion of "simulation" that we adopt here refers to resampling of observed events, i.e., we "recycle" patterns found in the observed data. [START_REF] Chailan | A semiparametric method to simulate bivariate space-time extremes[END_REF] use this idea to assess event-scale coastal hazards by conducting a risk analysis for some simulated storms from a reanalysis data set. While their approach provides an appropriate simulation technique, several aspects call for improvement.

The method developed in this paper generalizes the approach of [START_REF] Chailan | A semiparametric method to simulate bivariate space-time extremes[END_REF] to overcome some of its drawbacks. It uses probabilistic representations of extreme event episodes, similar to de Fondeville and Davison (2018), and related to the framework of functional extremes as proposed by de Fondeville and [START_REF] De Fondeville | Functional peaks-over-threshold analysis[END_REF]. With respect to [START_REF] Chailan | A semiparametric method to simulate bivariate space-time extremes[END_REF], we provide three major novelties. First, our procedure allows for an infinite number of stochastic simulations, instead of lifting to deterministically specified return levels. Second, we embed our semi-parametric resampling scheme in the framework of generalized Pareto processes. Third, a flexible general procedure is presented to identify extreme events and quantify their magnitude by accounting for space-time aggregation through homogeneous cost functionals that encapsulate operations such as averaging or taking maxima. Our space-time aggregation provides a generalization of spatial cost functionals used in the application of de Fondeville and Davison (2018), but here we resort to a nonparametric dependence specification, thus bypassing the issue of model mis-specification. With respect to multivariate extreme value analysis (i.e., when the observation domain consists of a finite and small number locations), our approach is closely related to empirical spectral measures, which have become a standard tool for estimating extremal dependence (e.g., [START_REF] Beirlant | Statistics of Extremes: Theory and Applications[END_REF][START_REF] Opitz | The spectrogram: A threshold-based inferential tool for extremes of stochastic processes[END_REF].

The remainder of this paper is structured as follows. Section 2 recalls the theory for generalized Pareto processes, and we specifically highlight the space-time setting where the temporal dimensionwhose role is fundamentally different from that of geographic space-must be made explicit. Techniques to practically implement and validate the spatio-temporal generalized Pareto framework are proposed in Section 3. Our algorithm to generate extreme space-time scenarios is developed in Section 4. We illustrate our approach on hourly rainfall reanalysis data available on a 1 km 2 grid in Southern France over a 10-year period from 1997 to 2007 in Section 5. In this case study, we perform a comparative analysis based on two commonly used risk measures using simulated extreme scenarios. Conclusions and and outlook to future research are given in Section 6.

Space-time generalized Pareto processes

We write S for a compact subset of R d to denote the area of interest and T for a compact subset of R + to denote the time dimension, and we denote by C(S × T ) the space of continuous functions on S × T , equipped with the supremum norm. The restriction of C(S × T ) to non-negative functions is written C + (S × T ). Similarly, we define the space of non-negative continuous functions in S as C + (S).

In multivariate EVT, a generalized Pareto limit was introduced in [START_REF] Rootzén | Multivariate generalized pareto distributions[END_REF] by conditioning on an exceedance event in at least one component. The aforementioned idea was extended to infinite-dimensional spaces by the definition of generalized Pareto process in [START_REF] Ferreira | The generalized Pareto process; with a view towards application and simulation[END_REF] where the condition is based on exceedances of the supremum over the compact study domain. To gain flexibility in the definition of the conditioning extreme events, [START_REF] Dombry | Functional regular variations, Pareto processes and peaks over threshold[END_REF] provided the notion of -Pareto processes by considering more general exceedances defined in terms of a homogeneous cost functional denoted . Our focus here is on the spatial and temporal dimensions for the extent of extreme events. With the aim of modeling phenomena that exceed a certain extreme threshold, we start by defining and characterizing space-time generalized -Pareto processes. The following constructive definition and theoretical results adapt notations and theory of [START_REF] Dombry | Functional regular variations, Pareto processes and peaks over threshold[END_REF] to the specific context of spatio-temporal modeling, where we highlight the role of the time dimension, which usually plays a different role than space.

Construction

We define a risk functional (also called cost functional ) : C + (S × T ) → [0, +∞) as a continuous nonnegative function that is homogeneous, i.e. (tf ) = t (f ) for t ≥ 0. Examples of such are the functions of maximum, minimum, average, or the value at a specific point (s 0 , t 0 ) ∈ S × T . Definition 2.1 (Standard space-time -Pareto process). Let W * = {W * (s, t)} s∈S,t∈T be a stochastic process in C + (S × T ). We call W * a standard space-time -Pareto process if it can be represented as

W * (s, t) d = RY (s, t) (1)
where 1. Y is a stochastic process in C + (S × T ) satisfying (Y ) = 1;

2. R has Pareto distribution with scale 1 and shape γ R , i.e., P(R > r) = r -γ R , r > 1;

3. Y and R are stochastically independent.

The above definition is equivalent to the definition through the POT stability property: for any u ≥ 1, the distribution of the renormalized threshold-exceeding process {u -1 W * | (W * ) ≥ u} is equal to the distribution of W * ; see Theorem 2 of [START_REF] Dombry | Functional regular variations, Pareto processes and peaks over threshold[END_REF]. By construction, we get

Y d = W * / (W * ) and R d = (W * ).
A generalized version of such Pareto processes is given in Definition 2.2 by allowing for flexibility in the marginal distributions according to the location-scale-shape parametrization commonly used in univariate EVT.

Definition 2.2 (Generalized space-time -Pareto process). Given an -Pareto process W * (s, t) constructed according to Definition 2.1 and continuous real functions σ(s, t) > 0, µ(s, t) and γ(s, t) in C(S × T ), a generalized space-time -Pareto process is any process constructed as

W (s, t) d = µ(s, t) + σ(s, t){W * (s, t) γ(s,t) -1}/γ(s, t), γ(s, t) = 0, µ(s, t) + σ(s, t) log W * (s, t), γ(s, t) = 0.
(2)

Asymptotic results

We shortly recall the two main asymptotic results for characterizing extremes of stochastic processes : max-stable processes and Pareto processes. We refer the reader to the literature for technical details [START_REF] Lin | On convergence toward an extreme value distribution in c[0,1][END_REF][START_REF] De Haan | Extreme Value Theory. An Introduction[END_REF][START_REF] Ferreira | The generalized Pareto process; with a view towards application and simulation[END_REF][START_REF] Thibaud | Efficient inference and simulation for elliptical Pareto processes[END_REF][START_REF] Dombry | Functional regular variations, Pareto processes and peaks over threshold[END_REF]. We use the symbol "⇒" to represent variants of weak convergence of random elements from the univariate, multivariate or functional domain.

Consider independent copies X 1 , . . . , X n of a stochastic space-time process X = {X(s, t)} s∈S,t∈T with continuous trajectories. We say that the process X is in the functional maximum domain of attraction of a max-stable process Z = {Z(s, t)} s∈S,t∈T with continuous trajectories if there exist continuous functions a n > 0 and b n such that max

1≤i≤n X i (s, t) -b n (s, t) a n (s, t) s∈S,t∈T ⇒ {Z(s, t)} s∈S,t∈T . (3)
Further details about space-time max-stable processes can be found in Davis et al. (2013a,b).

The convergence in (3) can be decomposed into the convergence of marginal distributions and the convergence of the process on a normalized marginal scale (i.e., the convergence of the dependence structure), such that the limiting marginal and dependence behavior can be studied separately; see de Haan and Ferreira (2006, Section 9.2). A standardised process X * = {X * (s, t)} can be defined by X * (s, t) = H -1 (F (s,t) (X(s, t))), s ∈ S, t ∈ T , where H -1 denotes the inverse function of the standard Pareto distribution function H, and F (s,t) denotes the distribution of X(s, t). If X has continuous marginal distributions F (s,t) , then X * has marginal standard Pareto distributions. For a n ≡ n, b n ≡ 0, the max-stable limit for X * in (3) is a standard max-stable process Z * = {Z * (s, t)} s∈S,t∈T with unit Fréchet marginal distributions; see de Haan and Ferreira (2006, Definition 9.2.4).

If X * is in the maximum domain of attraction of a max-stable process Z * and the cost functional is continuous at 0, we get the convergence of -exceedances on the standard scale:

u -1 X * (s, t)| (X * (s, t)) > u ⇒ {W * (s, t)} , u → ∞, (4)
where W * (s, t) is a standard space-time -Pareto process as in Definition 2.1 (Dombry and Ribatet, 2015, Theorem 3). Conversely, if the convergence in (4) holds for chosen as the maximum norm, then convergence in (3) of the max-stable process X * to Z * follows.

Practice

In practice, we use the asymptotic theory exposed in Section 2 for conducting statistical analyses on extreme events based on finite-sample data, which poses a number of practical challenges. In this section, we propose solutions for three issues: the standardisation of marginal distributions (Section 3.1), the definition of extreme space-time episodes (Section 3.2), the analysis and verification of asymptotic stability properties (so-called threshold-stability, see Section 3.3).

Marginal transformations

We first discuss suitable marginal transformations of X such that X * satisfies convergence with respect to -exceedances in (4). In theory, values of X * (s, t) close to 0 are pushed to 0 when u → ∞ in ( 4), but in practice the use of a high but finite threshold u leads to non-zero values in u -1 X * (s, t). Therefore, a certain ambiguity persists in practice to define the standardisation for relatively small, non extreme values of X(s, t). In particular, if the minimum value of the data process X arises with positive and non negligible probability, such as the value 0 for the absence of precipitation in our application study, then this minimum value should be mapped to 0 in the standardised process X * . Here, we develop the general idea of such transformations and a more specific transformation for precipitation data is proposed in Section 5. We choose a distribution function G : R → [0, 1] whose survival function Ḡ satisfies: x Ḡ(x) → 1, x → ∞, and Ḡ(0) = 1; we write G ← for the (generalized) inverse function of G. We then define the transformation T = T (s,t) : R → [0, ∞) towards the standardised process X * as follows:

X * (s, t) = T (X(s, t)) = G ← (F (s,t) (X(s, t))) (5)
where F (s,t) : R → [0, 1] denotes the distribution of X(s, t). The (generalized) inverse transformation of T can be defined as

T ← (f ) = F ← (s,t) (G(f )) for f ∈ C + (S × T ), with F ← (s,t) the (generalized) inverse function of F (s,t) .
Regarding marginal modeling, it is natural to use a tail representation motivated by univariate EVT, whose parametrization corresponds directly to the generalized Pareto process in Definition 2.2. For a fixed high threshold function u(s, t), we assume that

P(X(s, t) > x) = 1 -F (s,t) (x) = 1 + γ(s, t) x -µ(s, t) σ(s, t) -1/γ(s,t) + (6)
for x > u(s, t), with parameter functions for position µ(s, t) < u(s, t), for scale σ(s, t) > 0 and for shape γ(s, t), such that the right-hand side of ( 6) is less than 1 [START_REF] Thibaud | Efficient inference and simulation for elliptical Pareto processes[END_REF]. For data values X(s, t) below u(s, t), we may use appropriately chosen empirical distribution functions or any other useful model, where the probability mass below u(s, t) should amount to F (s,t) (u(s, t)) with F (s,t) defined in (6). The selection of a good threshold u(s, t) is a recurrent issue in extreme-value analysis; for instance, recently, Silva Lomba and Fraga Alves (2020) introduce an automated threshold selection procedure based on the method of L-moments.

From a practical point of view, the transformation that we introduce in (5) does not have a direct physical meaning in terms of a volume of water in the case of precipitation, but it can be interpreted as follows. Due to standardisation towards (approximate) standard Pareto tails, we have P(T (X(s, t)) > T (x)) ∼ 1 T (x) for large x. Therefore, given an independent copy X of X, we obtain the conditional probability

P(T (X (s, t)) > T (X(s, t)) | X(s, t) = x(s, t)) = P(T (X (s, t)) > T (x(s, t))) ∼ 1 T (x(s, t))
.

This means that the transformed value T (x(s, t)) can be interpreted as the return period of observations larger or equal to x(s, t). Therefore, at high quantiles (where we have the standard Pareto tails) we can interpret X * as the space-time process of (marginal) return periods.

Definicion of extreme episodes

For the purpose of simulating realistic spatio-temporal extreme scenarios, we have to define what "extreme" means. With environmental data, we often have only a single observation of the space-time process X, and very high values typically tend to cluster temporally within relatively short sub-periods. We consider such sub-periods as extreme space-time events. If it is realistic to assume that temporal dependence of extremes becomes negligible for relatively large time lags, theoretical results based on independent processes as in Section 2 can be used. In the space-time generalized Pareto process framework, the value of (X) quantifies the magnitude of events. In practice, we apply to a large collection of candidate episodes to extract the most extreme ones. Our extraction algorithm is designed to avoid temporal intersection of the selected extreme episodes.

There is no unique definition of an extreme event, i.e. of the cost functional , rather it depends on the nature of the considered phenomenon, on the data set, on the objective of the study, and also on the structure of the model [START_REF] Mcphillips | Defining extreme events: A cross-disciplinary review[END_REF]. For instance, [START_REF] Carreau | Multivariate density model comparison for multi-site flood-risk rainfall in the French Mediterranean area[END_REF] define rainfall susceptible of causing flash floods in small Mediterranean catchments as rainfall with a spatial average larger than 50 mm. Expert knowledge may suggest how to measure the extreme nature of an event, where the question of how to combine criteria related to duration, spatial extent and magnitude is recurrent. For instance, [START_REF] French | Quantifying the risk of heat waves using extreme value theory and spatio-temporal functional data[END_REF] develop new visualizations of extreme heat waves by composing a temporal and spatial cost functional. [START_REF] Chailan | A semiparametric method to simulate bivariate space-time extremes[END_REF] extract extreme wave heights based on spatio-temporal maxima in sliding time windows.

In the following, we use the idea of sliding space-time windows and specify the support of the cost functional introduced in Section 2.1 as a neighborhood N (s, t) at location s ∈ S and at time t ∈ T . In practice, the window size defines the maximal time duration and spatial extent of extreme events. The space index s may be missing if we consider the full study area for extracting extreme events. This neighborhood could be defined through an event duration δ in time, and the spatial support could be the full study area or a sub-region such as a catchment or a certain distance buffer around a specific site s 0 . To indicate the local support of the cost functional defined as a neighborhood around (s, t), we use the notation s,t (X * ) = ({X * (s , t ), (s , t ) ∈ N (s, t)}).

We propose to define N (s, t) as the product of a spatial neighborhood N (s) (for instance, {s ∈ S | ss ≤ h km} with h ∈ R + ) and a temporal neighborhood

N (t) (for instance, {t ∈ T | |t -t | ≤ k hours} with k ∈ N ), N (s, t) = N (s) × N (t).
The above choice of the spatial extent and temporal duration of the neighborhood allows taking into account the spatial and temporal dependence of extreme episodes in the studied dataset; the resulting space-time process s,t (X * ) can be seen as a local smoothing of X * . For instance, in the application in Section 5 we choose h = 15 km and k = 12 hours. Useful cost functionals for space-time episodes are obtained by composing a spatial functional S with a temporal functional T , the latter applied to the values of S observed over a number of consecutive time steps :

s,t (X * ) = T ( S s,t-(δ-1) (X * ), . . . , S s,t (X * )), (7)
with S s,t (X * ) = S ({X * (s , t) | s ∈ N (s)}) and δ the duration of the episode. Moreover, based on S s,t (X * ) we can define cost functionals that combine the values obtained for all spatial neighborhoods N (s) by taking their maximum value (or again, any other spatial aggregation value). In this case, we define:

t (X * ) = T max s∈S S s,t-(δ-1) (X * ), . . . , max s∈S S s,t (X * ) . ( 8
)
If X * satisfies the functional domain of attraction condition (3), then

P( (X * ) > u) ∼ θ /u, u → ∞, (9)
where θ is the -extremal coefficient (for details, see [START_REF] Engelke | Extremal behaviour of aggregated data with an application to downscaling[END_REF]. When s,t corresponds to the maximum function over N (s, t) (i.e., T = max and S s,t = max), the -extremal coefficient θ s,t defines the classical extremal coefficient of the domain N (s, t) (see Example 4 of [START_REF] Engelke | Extremal behaviour of aggregated data with an application to downscaling[END_REF].

Using (9), we can calculate approximate return levels for extreme episodes characterized as -exceedances above a large threshold u. The simplest case arises for θ = 1, i.e., when θ is known beforehand and we do not have to estimate it from data. For instance, if (s 0 , t 0 ) ∈ S × T is a fixed space-time point, we can define the cost functional value (X * ) as X * (s 0 , t 0 ), and θ = 1. Moreover, θ = 1 if is the average, i.e. s,t (x) = 1 |N (s,t)| N (s,t) x(s , t ) d(s , t ); see Ferreira et al. (2012, Proposition 2.2). When θ = 1, an estimator of θ can be plugged into (9), such as an empirical estimate or alternative estimators [START_REF] Engelke | Extremal behaviour of aggregated data with an application to downscaling[END_REF]. Finally, since P( ((X ) * ) > (X * ) | X * = x * ) ∼ θ / (x * ) at high quantiles of (X * ) for an independent copy X of X, we can interpret (x * )/θ as the return period of an extreme event x * . In addition, by taking into account (5), the cost functional (approximately) aggregates marginal return periods X * (s, t) into return periods (X * ) for space-time episodes.

Techniques to analyze asymptotic dependence properties

The functional domain of attraction condition in (3) is the theoretical basis for using generalized Pareto processes. It requires that a relatively strong type of extremal dependence, known as asymptotic dependence, prevails in the data-generating process X, at least for small distances in space and time. With asymptotic dependence between two points (s, t) and (s , t ) = (s + ∆s, t + ∆t), we observe a strictly positive limit of the probability P F (s ,t ) (X(s , t )) > u | F (s,t) (X(s, t)) > u as u → 1. With asymptotic dependence, so-called threshold stability holds when moving towards higher quantiles, such that the typical spatial and temporal extent of clusters of extreme values does not depend on event magnitude. In practice, we should verify that data exhibit such asymptotic dependence, and we need tools to assess the spatial and temporal range of asymptotic dependence, e.g., for choosing appropriate cost functionals . We describe two approaches: the study of empirical bivariate tail correlation coefficients corresponding to the above conditional probabilities (also known as the χ-measure, see [START_REF] Coles | Dependence measures for extreme value analyses[END_REF]), and the assessment of the independence of observed scale variable (X * ) and profile process X * / (X * ) for when the scale variable exceeds a high threshold.

Spatial and temporal tail correlation coefficients

Let (X 1 , X 2 ) be a bivariate random vector. Following [START_REF] Coles | Dependence measures for extreme value analyses[END_REF], the extremal dependence of (X 1 , X 2 ) can be analysed through the behavior of the conditional probability

P (X 1 > F ← 1 (q)|X 2 > F ←
2 (q)) as q tends to 1, where F ← i , i = 1, 2, denotes the generalized inverse function of X i . The distribution of the random vector (X 1 , X 2 ) is said to be asymptotically dependent if a positive limit value χ arises for the conditional probability, i.e.,

χ(q) = P (X 1 > F ← 1 (q), X 2 > F ← 2 (q)) P (X 2 > F ← 2 (q)) → χ > 0 (q → 1 -), ( 10
)
where χ is also known as the tail correlation coefficient. The case χ = 0 represents asymptotic independence.

Therefore, we can use the χ-measure to obtain a summary of extremal dependence with respect to distance in space and lag in time. For a given finite level q < 1, we define χ spa (h; q) to measure extremal dependence between sites separated by spatial distance h at a given time. Similarly, we define χ tim (k; q) as the χ-index measuring extremal dependence for a time lag k at a given site. We calculate χ spa (h; q) empirically based on observation pairs with structure (X(s, t i ), X(s+∆s, t i )) where ∆s = h, and we empirically calculate χ tim (k; q) based on observation pairs of spatial maxima, (max s∈S X(s, t i ), max s∈S X(s, t i + k)). The estimators for these two χ-functions are obtained by plugging in empirical survival functions and empirical quantiles in (10). From now on, we refer to these two estimators as the empirical χ spa (h) and χ tim (k) functions, respectively.

Independence of scale and profile

The POT stability manifests itself through the (approximate) independence between the profile process Y = X * / (X * ) and the random scale R = (X * ) for (X * ) > u. In practice, the threshold u should be high enough for this property to hold approximately, such that the limit process in (4) becomes a useful approximation to data. Due to the very high dimension of the profile process in the space-time setting, it is difficult to check this independence directly based on observed scales and profiles. Instead, we propose to check for the absence of strong trends in summary statistics of Y with respect to the event magnitude R, which would indicate dependence between Y and R.

In our application, we will focus on checking the scale-profile independence in space by considering the set of extreme spatial episodes W * t satisfying S t (W * t ) > u, and we use two summary statistics calculated from the profile processes

Y t = W * t / t (W * t ) in C + (S)
. First, we consider f u (Y t ) defined as the proportion of sites s where Y t (s) ≤ u : useful values of u are relatively small or large quantiles of Y t , to check for trends in the magnitude of Y t with respect to t (W * t ). Second, we consider the empirical standard deviation sd(Y t ) of Y t (s) = Y t (s): if there are trends with respect to event magnitude, we usually find trends of sd(Y t ). The square root transformation ensures finite standard deviation values.

Another related check consists in ensuring that R/u has standard Pareto distribution; otherwise, the structure of the Pareto process is compromised. Simple checks consist in drawing QQ-plots and associated confidence bounds of observed versus theoretical quantiles, e.g., by adopting a uniform scale by transforming to 1-u/R. Theoretical developments around formal statistical tests to check the Paretodistribution of the scale variables have been conducted by [START_REF] Falk | Testing for a multivariate generalized Pareto distribution[END_REF] for multivariate Pareto distributions (i.e., for the finite-dimensional case), and we here adopt their idea for stochastic processes through the analysis of QQ-plots.

Several empirical studies on climatic data show that extremal dependence may weaken when the event magnitude increases [START_REF] Opitz | The spectrogram: A threshold-based inferential tool for extremes of stochastic processes[END_REF][START_REF] Huser | Modeling spatial processes with unknown extremal dependence class[END_REF][START_REF] Le | Dependence properties of spatial rainfall extremes and areal reduction factors[END_REF][START_REF] Tawn | Modelling spatial extreme events with environmental applications[END_REF]. Then, asymptotic independence may ultimately arise, or the dependence strength may stabilize at very high but unobserved magnitudes. We cannot check this stability behavior with absolute certainty in finite samples. If the extremal dependence strength continues to weaken in data above the selected threshold u, we acknowledge that the generalized Pareto process framework leads to rather conservative probability estimates for observing concomitant high values. A tendency towards weakening extremal dependence would imply that, at higher quantile levels, extreme episodes tend to have smaller extent in space, and tend to stretch less far in time, i.e., have smaller temporal extent.

Uplifting observed extreme episodes

We now describe the general procedure for the extraction of extreme space-time episodes (Section 4.1) and the algorithm to resample new space-time scenarios (Section 4.2). A probabilistic interpretation of this resampling scheme is given in Section 4.3. Throughout and without loss of generality, we here use the same notation for the single observation of the space-time process X(s, t) and the stochastic process itself.

Selection of extreme episodes

Algorithm 1 describes the extraction of extreme episodes from standardised data X * . To start, we define the space-time neighborhoods N (s, t) whose intensities are assessed by applying the cost functional . In particular, it is possible to restrict extraction to subsets S of S and T of T ; for instance, sites s whose neighborhood N (s, t) is not fully contained in the study domain should be removed in S . If the neighborhood is the full study region, we may drop the index s and simply write N (t). We choose a threshold u for the cost functional for which the asymptotic stability properties underpinning our approach are (approximately) satisfied. There must be at least one exceedance of the cost functional above the threshold in the data set. The first step of the algorithm is to compute the values of for each neighborhood N (s, t). We select as the first extreme episode the neighborhood N (s 1 , t 1 ) where s,t reaches its maximum value 1 . We aim at extracting a collection of extreme episodes that are at most weakly dependent; therefore, the algorithm needs a mechanism to "decluster" extreme episodes. The second extracted extreme episode corresponds to the maximum value of s,t (X * ) arising in the data set X * (s, t) with t in the set of reduced time steps after removal of time steps that intersect with N (s 1 , t 1 ) or, more generally, with a larger temporal buffer zone N buffer (t 1 ) around t 1 involving a buffer parameter β ≥ 0 to remove more time steps. We then iterate this procedure of episode extraction and data set reduction. The stopping criterion for the extraction of extreme episodes is two-fold: either a fixed target number m of extreme episodes is reached, or the extreme condition s,t (X * ) > u for a fixed high threshold u cannot be fulfilled any longer in the reduced data set.

If the maximum of s,t (X * ) is not unique and is realized at several coordinates (s, t), we must define a rule to extract a single (s, t) that identifies the corresponding extreme space-time episode. In particular, if we find several consecutive time steps t where t (X * ) in Equation ( 8) is equal to the maximum, we fix the anchor time step t of the extreme episode as follows. In this situation, typically δ consecutive values are equal, and we then set t to the closest value below or equal to the median of these time steps. This rule will tend to center the extreme space-time episode on the strongest values in X * . That is, if the maximum arises at time steps t 0 , . . . , t 0 + δ -1, we fix t = t 0 + δ 2 as the anchor time step of the extreme space-time episode.

Algorithm 1: Algorithm for selecting extreme episodes defined over space-time neighborhoods N (s, t). In Step 8, instead of extracting only the extreme neighborhood N (s i , t i ), we may sometimes want to extract the full study domain N (t i ) × S.

Input:

• {X * (s, t), s ∈ S, t ∈ T }, space-time observations on a standardised scale;

• S ⊆ S sites of interest and T ⊆ T time steps of interest.

• m the maximum number of extreme episodes to select;

• u threshold on s,t (X * ) for the selection of extreme episodes;

• δ > 0 the duration of extreme episodes defining temporal neighborhoods N (t) = [t -(δ -1), t];

• β ≥ 0 buffer time step to ensure independent extreme episodes defining extended temporal neighborhoods N buffer (t) = [t -(δ -1) -β, t + (δ -1) + β];

• N (s) spatial neighborhood for s ∈ S , such that N (s, t) = N (s) × N (t).

Output:

• m: the number of selected extreme episodes (m ≤ m );

• X * [1] , X * [2] , . . . , X * [m]
, {s 1 , s 2 , . . . , s m }, {t 1 , t 2 , . . . , t m }, { 1 , 2 , . . . , m }: collection of extreme episodes; observation sites and times; aggregation values related to extreme episodes.

1 begin 2 Set I = T . 3 Calculate s,t (X * ) for all t ∈ T , s ∈ S with N (s, t) ⊂ S × T . 4 i ← 1. 5 while i ≤ m and max s∈S ,t∈I s,t (X * ) > u do 6 (s i , t i ) ← arg max t∈I,s∈S s,t (X * ) 7 i ← si,ti (X * ) 8 X * [i] ← {X * (s , t ), (s , t ) ∈ N (s i , t i )} I ← I \ N buffer (t i ) i = i + 1 9 m ← i -1 10 return m, X * [1] , X * [2] , . . . , X * [m] , {s 1 , s 2 , . . . , s m }, {t 1 , t 2 , . . . , t m }, { 1 , 2 , . . . , m }

Simulation method

The simulation procedure detailed below is an extension of the simulation method of Ferreira and de Haan (2014) to the spatio-temporal setting. To sample new extreme space-time scenarios, we proceed as follows:

1. Standardisation: Estimate marginal tail parameter functions γ(s, t), σ(s, t) and µ(s, t) in ( 6), and denote by X * = {T (X(s, t))} s∈S,t∈T the resulting standardised process (5).

2. Selection of extreme episodes: Fix the maximum number of extreme episodes m . Use Algorithm 1 to extract the collection of m ≤ m extreme episodes X * [i] , i = 1, . . . , m. 3. Lifting: Sample R i , i = 1, . . . , m according to a Pareto distribution with shape 1 and scale α > 0,

i.e. P(R i > x) = α/x, x ∈ [α, ∞), and generate lifted extreme episodes as

V i (s, t) = R i X * [i] (s, t) i = R i Y i (s, t), (s, t) ∈ N (s i , t i ). ( 11
)
4. Back-transformation to original scale: Lifted extreme episodes are transformed back to the original marginal scale by

W i (s, t) = T ← (V i (s, t)), (s, t) ∈ N (s i , t i ).
When fixing the value m of the number of extreme episodes to extract, we aim for a representative sample of spatio-temporal extremal patterns in the data, but have to keep in mind that for a large value of m the POT stability property may not be satisfied.

Interpretation

Our methodology draws justification from the framework of generalized space-time Pareto processes as given in Definition 2.2. For each extreme episode, the above procedure allows us to sample new realizations of generalized space-time Pareto processes, i.e., there is an appropriate probabilistic interpretation of extreme events resampled through our lifting procedure.

We outline the mathematical derivation of this property. Suppose that we have fixed a threshold u > 0 corresponding to a relatively high quantile of the cost functional . Our modeling assumption in (4) can be informallys written as

W * (s, t) d ≈ X * (s, t) u conditional to (X * (s, t)) ≥ u,
where W * represents a standard -Pareto process. Therefore, the extreme episodes on the standardized marginal scale, X * (s, t) conditional to (X * (s, t)) ≥ u, have distribution according to a generalized Pareto process with µ(s, t) = u, σ(s, t) = u and γ(s, t) = 1, using notations from Definition 2.2.

Next, we recall that the structure of the marginal model with parameters for shape, scale and location in ( 6) is equivalent to

X(s, t) = σ marg (s, t) (X * (s, t)) γmarg(s,t) -1 γ marg (s, t) + µ marg (s, t),
given that we consider the (univariate) tail of X(s, t). Combining the two points above leads to X(s, t) ≈ σ marg (s, t)u γmarg(s,t) (W * (s, t)) γmarg(s,t) -1 γ marg (s, t)

+ µ marg (s, t) + σ marg (s, t)(u γmarg(s,t) -1) γ marg (s, t) , conditional to (X * (s, t)) ≥ u, where the σ marg , µ marg , γ marg -functions are those from the univariate distributions. Therefore, the extracted and lifted extreme episodes possess (approximately) the distribution of a generalized Pareto process as given in Definition 2.2, with parameters γ(s, t) = γ marg (s, t), σ s,t = u γmarg(s,t) σ marg (s, t), and µ(s, t) = µ marg (s, t) + σ marg (s, t)(u γmarg(s,t) -1) γ marg (s, t) .

We do not use the marginal location-scale-shape parametrization below the marginal thresholds, which means that the structure of our processes differs from the generalized Pareto process below the threshold in the margins, i.e., the generalized Pareto approximation may be cruder below the threshold. However, asymptotically (for increasingly high thresholds), we can expect that the transformation of values far below the threshold has little influence.

The following property holds for the return periods of resampled extreme episodes. Using the extracted observed episodes, the lifting procedure in Section 4.2 samples new realizations V i of a space-time Pareto process with support N (s i , t i ) for each extreme episode i. Let X be an independent copy of X, and drop the indices s i , t i of in the following equations. Since P( ((X ) * ) > x) ∼ θ /x for large x, we can write

P ( ((X ) * ) > (RX * / (X * ))|R = r, X * = x * ) = P ( ((X ) * ) > (rx * / (x * ))|R = r, X * = x * ) = P ( ((X ) * ) > r) ∼ θ /r.
Since we resample scale variables R i that are larger than a lower bound α, we obtain that the minimum return period is α/θ for resampled extreme episodes V , and larger α leads to longer return periods. By construction, uplifted scenarios have the same spatial patterns of variability as observed scenarios, but they correspond to longer return periods if α is chosen larger than (x * ).

To go further, it is possible to theoretically show that, after marginal standardisation, the resampled and backtransformed episodes W i (s, t) have approximately the same probability distribution as the observed extreme episodes in X(s, t). Our procedure therefore generates a threshold-stable stochastic process at a higher level than the observed one when α > i .

Application

We use our resampling algorithm to produce large numbers of realistic spatio-temporal extreme precipitation scenarios in a region in Mediterranean France where flash floods are frequent. Furthermore, we show how to calculate two risk measures for the most extreme observed space-time episodes before and after lifting them to longer return periods.

Data set

Our semi-parametric approach does not provide a mechanism to spatially interpolate observations. Therefore, precipitation measurements should be available over a sufficiently dense network of sites. We use hourly precipitation reanalysis data over a 1 km 2 grid, constructed by merging radar signals and observed hourly precipitation totals [START_REF] Tabary | A 10-year (1997-2006) reanalysis of quantitative precipitation estimation over france: methodology and first results[END_REF]. The grid has 10, 914 cells covering a 133.2 km × 104.3 km area in Mediterranean France, see Figure 1, with 87, 642 hourly time steps covering the 10-year period from 1997 to 2007. The unit of measurement is mm/h. This data set was provided by Météo-France (http://www.meteofrance.com). The large dimension of the data set allows us to disregard restrictive parametric assumptions in favour of a nonparametric approach for the extremal dependence model.

Empirical return levels of rainfall intensities at the 98% level (i.e., of strictly positive observations) is reported for each grid cell in Figure 1.

Marginal distributions

The first step of our lifting procedure is the definition of a marginal transformation T , appropriate for extreme hourly precipitation data, to obtain the standardised process X * in (5). We first discuss our choice of the target distribution G. Due to the hourly temporal resolution, zero values occur with very high frequency in the data. Therefore, we include a discrete mass p 0 at 0 to represent the absence of precipitation. Following [START_REF] Opitz | Modeling asymptotically independent spatial extremes based on laplace random fields[END_REF], we construct G to have a mass p 0 ≥ 0 at 0, a uniform density on (0, x 0 ), and a standard Pareto distribution for x > x 0 where x 0 > 1. The junction point x 0 is chosen to ensure the continuity of the density of G for x > 0:

G(x) =        0, x < 0, p 0 , x = 0, p 0 + (1-p0) 2 4 x, 0 < x ≤ 2/(1 -p 0 ), 1 -1/x, x > 2/(1 -p 0 ). ( 12
)
An illustration of G for p 0 = 0.7 is provided in Figure 2. By keeping the value 0 in the distribution G for dry cells (i.e., for values 0 in the original data), we impose that the lifting procedure does not modify dry cells since a value 0 will always remain at 0. This may not always be realistic, but our focus is on extreme rainfall, such that very accurate modeling of dry cells during extreme episodes goes beyond the scope of this work. However, our approach ensures that observed 0 values will not be lifted to positive precipitation values; otherwise, one could generate unrealistically large regions with positive precipitation in the lifting simulations. Next, we choose the distribution function F (s,t) of X(s, t) as the empirical distribution function F (s) (i.e., at each grid cell s) when X(s, t) ≤ u(s, t), and according to (6) when X(s, t) > u(s, t). We use spatial models for the marginal tail parameters, whose estimators μ(s), σ(s) and γ(s) in F (s) are obtained by estimating the marginal distributions site by site using a threshold u(s) chosen as a high empirical quantile for fixed s; here, we choose the 0.95-quantile of hourly rainfall intensities. Thanks to the consistency of these estimators and the continuity of T , we can apply the continuous mapping theorem such that the transformation T (with estimators plugged in) provides a consistent estimate of T .

Spatio-temporal cost functionals

We here consider two cost functionals. Our first choice (1) s,t is the spatio-temporal median, i.e., the median of observations X * (s, t) from the spatio-temporal neighborhood N (s, t) given as a space-time cylinder, calculated for all (s, t) whose neighborhood is fully contained in the study region. Regarding spatial neighborhood (i.e., the base of the cylinder), we specify it through a 15 km disc centered at s, where the choice of 15 km is further motivated in the following subsection presenting empirical extremal dependence properties; regarding temporal neighborhood (i.e., the height of the cylinder), we set it to extend backward from t, such that N (t) = {t-(δ-1), . . . , t} with duration δ = 12 hours, again motivated by considerations detailed in the following section. Our choice of a diameter of 30 km for spatial discs can also be explained by hydrological considerations. The dimension of catchments in the study area is heterogeneous, but the 30 km-diameter is close to the diameter of some important catchments in the study area, especially the catchment of the Lez river flowing through Montpellier, the largest city in the study area, which is often affected by urban flooding. With respect to the temporal length of episodes, we follow the twin aims of representing heavy rainfall in a relatively short time period (e.g., convective precipitation episodes) and moderate rainfall over a larger time period (e.g., cyclonic precipitation episodes). To this end, we here consider 12 hours a good compromise for the temporal window in both cases.

The second cost functional

(2) s,t is chosen as the spatio-temporal maximum and does not depend on spatial location as the spatio-temporal median above; its value

(2) t (X * ) corresponds to the maximum over the whole study area at time t. Using notation from Equation (7), we have N (s) = S, S s,t = max and T = max.

To achieve declustering of extreme episodes in the time domain, we set β = 1 to separate extreme Figure 3: Empirical tail correlation functions for quantiles at probability levels 0.98 and 0.99. Left: empirical χ spa (h), based on a subsample of 1500 pairs of grid cells, with a smoothed curve based on local polynomial regression (turquoise and orange lines for 0.98 and 0.99, respectively). Right: empirical χ tim (k), based on pairs of spatial maxima separated by a time lag k (same color code). Block bootstrap confidence intervals at 95% (dashed lines; blue for 0.98, red for 0.99).

episodes by at least 1 hour in Algorithm 1.

While the median in the first cost functional puts focus on the central tendency of the sample of (local) values, its behavior may be very different from the maximum in the second cost functional, which does not take into account the variation of values around the maximum. For the median, the local specification makes sense, since otherwise the value of the cost functional would be too strongly dominated by mostly small or zero values arising far from the high-intensity region of the precipitation episode.

For both choices, we would need an estimate of the extremal coefficient θ (i) s,t , i = 1, 2, to set or obtain return periods of the lifted events. For this purpose, we propose to use a maximum censored likelihood for estimating the scale parameter θ (i) s,t of a Pareto distribution with fixed shape 1, using observed magnitudes (i) s,t (X * ), i = 1, 2, censored below a high threshold u.

Extremal dependence properties

Using techniques proposed in Section 3.3, we first illustrate pairwise empirical tail correlation functions with respect to spatial distance and temporal lags. We also check if threshold stability is a valid assumption for the data set at high quantiles.

Figure 3 shows the estimates of spatial (left panel) and temporal (right panel) tail correlation functions, using empirical quantiles at 0.98 and 0.99 probability level. Pointwise block bootstrap confidence intervals at 95% are constructed using variable size blocks with block length following a geometric distribution with mean 300 hours [START_REF] Politis | The stationary bootstrap[END_REF][START_REF] Davis | Estimating extremal dependence in univariate and multivariate time series via the extremogram[END_REF]. Lower confidence bounds remain clearly separated away from 0 for spatial distances below 100 km and for time lags shorter than 12 hours, hinting at substantial extremal dependence at finite, observed quantile levels. There is a certain sensitivity of the estimated curves with respect to the probability used for fixing the threshold u to the corresponding empirical quantile, with a slight tendency towards decreasing dependence strength at higher levels; see Figure 3.

We complement these findings by checking threshold stability based on the independence of scales and profiles for high event magnitudes, and based on the Pareto distribution of the cost functionals; see Subsection 3.3.2 for methodological background. Certain calculations for extreme episodes were quite sensitive to the high proportion of 0 values (i.e., absence of precipitation) in the data set, which amounts to around 92 %. With p 0 = 0.92 in the marginal standard distribution G, the standard Pareto tail holds only for values above 2/(1 -p 0 ) = 25, i.e., at very levels, which may be the reason. Therefore, we add a preprocessing step where we remove hourly time steps t i from the data set if the precipitation totals in a sliding 24 hour-window centered at t i , cumulated over all grid cells, are smaller than 550 mm, corresponding to a spatially averaged precipitation total per grid cell of 0.05 mm over 24 hours. The resulting data subset contains only around 23 % of 0 values, and we work with this subset for our lifting procedure. Moreover, by applying this prepocessing step, we avoid to consider too heterogeneous weather patterns, i.e., we avoid to pool together very dry periods and relatively wet periods.

For checking the scale-profile independence at high levels, we first calculate (i) s,t , i = 1, 2, for all (s, t) ∈ S × T , and then use our selection algorithm 1 to extract the declustered extreme space-time episodes based on setting the maximum number of extracted episodes (m ) to ∞, and we set the threshold u of the cost functional above which episodes are considered as extreme to the 99.5%-quantile of all Denote by Y i = {Y i (s, t)} the extracted profile processes, i = 1, . . . , m. Figure 4 shows exploratory plots for the local spatio-temporal median (first row) and the spatio-temporal maximum (second row). In the leftmost and second-to-left displays, the proportion of profile process values Y i (s, t) below or equal to a threshold u ≥ 0, denoted by f u (Y i ), is plotted for u = 0, and for u fixed to the empirical 0.95-quantile of all episodes Y i taken together. In particular, the value u = 0 corresponds to the proportion of dry pixels. The empirical standard deviation sd(Y i ) of the square root Y i (s, t) of profile process values Y i (s, t) is depicted in the third displays (from the left) of Figure 4. For easier visual interpretation, the summary statistics are plotted against the observed cost functionals on uniform scale, i.e. 1 -u/ * i . Finally, we recall that a distribution of transformed cost functionals 1 -u/ * i that deviates from uniformity on [0, 1] may cast doubt on the threshold stability property. A QQ-plot of 1 -u/ * i is shown in the fourth display of Figure 4 with pointwise confidence bounds.

Judging from the shape of the local regression curves in Figure 4, we do not discern strong systematic trends in profile summaries with respect to event magnitude, although some slight trends are perceptible. The QQ-plots hint at slight deviation from uniformity. Overall, there are only weak signs of some nonstationary behavior, and we do not detect strong signals for asymptotic independence. We conclude that our lifting framework based on threshold-stability is appropriate for the data, and for the lifting procedure we work with the extreme episodes extracted for the current analysis. If data were truly asymptotically independent, we acknowledge that our resampling procedure would lead to rather conservative estimates of aggregated extreme risks.

Risk analysis

We first highlight differences in the ranking of the most extreme episodes according to the two cost functionals, and we provide some maps to contrast the original and the standarized marginal scale. Ending times t i for the 6 most extreme precipitation episodes are reported in Table 1. There is no overlap between the extreme episodes according to the two cost functionals, which underscores their fundamentally different approach to measuring the overall magnitude -either by considering only the strongest value in the case of the spatio-temporal maximum, or by using a statistical measure of central tendency of all the observations recorded during the episode in the case of the spatio-temporal median. We also see that the most extreme precipitation scenarios tend to arise during the autumn season.

Figures 5 and6 show the original precipitation data X(s, t) and a corresponding uplift W (s, t) for several time steps with extreme episodes extracted using the spatio-temporal median or the spatiotemporal maximum, respectively. The new magnitude variable R in the lifting step is respectively fixed to twice the maximum value of observed magnitudes, such that we obtain a clear increase of precipitation values.

Risk is a complex notion, in particular for space-time resolved processes. The most commonly used risk measure in hydrology corresponds to the univariate return level (RL) at a probability level q ∈ [0, 1], denoted as Q q . However, the return level fails to give any information about the thickness of the tail of the distribution above this level. An alternative risk measure was proposed in actuarial science, the s,t . First row: most extreme episode, here shown for t =1999-09-14, 04:00:00; second row: second most extreme episode, here shown for t =2005-09-06, 15:00:00. Grey and black contour indicate altitude. so-called Conditional Tail Expectation (CTE, [START_REF] Denuit | Actuarial Theory for Dependence Risks: Measures, Orders and Models[END_REF]. Information about the thickness of the tail of the distribution is included in the CTE, defined for a given level q ∈ [0, 1] and for a random variable X by CT E q (X) = E(X|X > Q q (X)).

We perform a risk analysis that aims at exploring differences in lifted extreme episodes that can be imputed to the choice of cost functionals, and of the Pareto scale parameter α used for sampling new scale variables R i in the lifting step; see Equation 11. We here consider the 50 largest declustered space-time episodes extracted for each of the two cost functionals (2) s,t . For each scenario, we generate 500 lifted simulations, i.e., we sample 500 values of the scale variable R i and combine each one of them with one of the extreme episodes selected at random. We consider two choices for the Pareto scale parameter α. Our first choice α 1 is the 20th largest value of the respective cost functional. Given the 10-year observation period of data, this value approximately corresponds to a return level for 6 months, if we assume temporal stationarity. For the second choice, we set α 2 = 10α 1 , i.e., we increase the return level (and therefore the return period) by 10.

We compute RL and CTE with q = 0.98 and q = 0.99 for the values of each lifted episode W i , where we first aggregate values W (s, t) for each spatial grid cell by taking pixel-wise temporal means over the 12 time steps. Therefore, the return level measure reports the precipitation level that is exceeded in 100(1 -q)% of all pixels (i.e., in 2% or 1% of the surface of the study region with q = 0.98 and q = 0.99, respectively), while the CTE measure reports the average precipitation value over the pixels where this return level is exceeded.

Figure 7 presents boxplots of the RL and CTE values calculated for the 500 simulations and different scenarios: different cost functionals; q = 0.98 and q = 0.99; and lifting parameters α 1 and α 2 . By construction, the CTE values are systematically higher than the RL values. All risk measures appear to be relatively heavy-tailed, which we impute to the heavy tails in the precipitation data. Moreover, the ten times higher scale parameter α 2 leads to substantially higher risk than α 1 , and this increase trends to be stronger than the one observed between risk measures when changing from q = 0.98 to q = 0.99. Furthermore, the aggregation through the spatio-temporal median leads systematically to higher risk measures than the spatio-temporal maximum. Indeed, the spatio-temporal maximum may tend to select episodes with highly localized peaks, i.e., there may be a large majority of zeros or small values with a few spatially strongly confined clusters of very large precipitation intensities. On the other hand, risk measures based on the spatio-temporal median better account for the persistence of moderate to high precipitation intensities through space and time. We underline that the mechanism of spatio-temporal cost functionals allows for a flexible choice according to the modeling context.

Conclusion and outlook

In this work, we set up a general framework for space-time generalized Pareto process. It allows developing a semi-parametric method to simulate extreme space-time scenarios of phenomenona such as precipitation. The extremal dependence structure is fully data-driven, and we require parametric assumptions only for the univariate tails, based on asymptotic theory. A crucial component is the cost functional defined over a sliding space-time window. It characterizes extreme episodes as episodes whose return period, appropriately aggregated over space and time, exceeds a high threshold. The application of our method to a gridded precipitation data set in Mediterranean France was used for a relatively simple risk analysis. It illustrates how spatio-temporal cost functionals can be defined, how they affect the selection of extreme episodes, and how the magnitude of the newly sampled scale variables impacts the magnitude of the lifted Figure 7: Boxplots for spatial RL and CTE at q = 0.98 at q = 0.99. First row: spatio-temporal median cost functional. Second row: spatio-temporal maximum cost functional. With α 1 approximately a RL for 6 months (red boxplots), and α 2 = 10α 1 (blue boxplots). extreme episodes on the original marginal scale. The proposed methodology requires densely gauged networks or gridded data as spatial interpolation is currently not enabled. While model fitting and choice with parametric models for gridded data is hard due to the large sample size and further bears risks of model mis-specification, semi-parametric techniques such as ours are appealing alternatives for gridded data and scale well with the size of the dataset.

In practice, it is difficult to find extreme value data with long observation periods to empirically study extreme value properties for long return periods without strong modeling assumptions. For practitioners, we provide a methodology that allows them to create extreme scenarios where they can control return levels or periods for aggregated data without any need to explicitly model dependence at extreme quantiles. Since our method identifies extreme episodes based on local conditions (e.g., for spatial discs of diameter 30 km), it would be possible to "trade space for time" in order to obtain a satisfactory number of approximately independent replicates of extreme episodes, given the following conditions: the process should be close to stochastic independence at larger distances in space, the study region should be relatively large, and the assumption of spatial stationarity should be realistic. Indeed, while it may be difficult to extract a large number of extreme episodes occurring locally around a fixed location, we may "borrow" extreme episodes from other, relatively distant locations and just shift them in space towards the location where we wish to obtain more observed extreme episodes. The assumption of spatial stationarity may not be realistic in practice for processes such as precipitation, but we could work with milder assumptions, e.g., by borrowing extremal episodes between locations or zones for which hydrological expertise suggests similar behavior of extremes. Assume that we want to extract empirical extreme profile processes with spatial support A, where A is compact in R 2 . We can extract such observed extreme episodes at s + A, centered at some pixel s, and then simply shift them by -s to obtain the desired support A. Through this approach, we may attain a moderately large number of observed extreme episodes (say, around 50 to 100) even if the length of the observation period is short. The lifting procedure then further allows us to substantially increase the number of available extreme episodes for impact studies, and we can draw reliable statistical inferences on summary statistics related to the extreme episodes, such as statistical summaries (mean, variance etc.), risk measures in the actuarial context, or output variables of impact models (e.g., models to study flooding in urbanized areas).

Besides precipitation reanalyses, other types of interesting applications include resampling from regional or global climate model reanalyses or projections. The work of [START_REF] Opitz | Semi-parametric resampling with extremes[END_REF] focused on spatial resampling of heatwaves in France; they combined a lifting step with nonparametric resampling techniques such as Multiple-Point Statistics [START_REF] Mariethoz | Multiple-point geostatistics: stochastic modeling with training images[END_REF], and Direct Sampling in particular [START_REF] Mariethoz | The direct sampling method to perform multiplepoint geostatistical simulations[END_REF], such that new profile processes inheriting spatial dependence patterns from the observed extreme episodes were sampled. The application of such resampling techniques allows for a strong increase in the variety of newly sampled extreme episodes and could be extended to our spatio-temporal setting, but it also implies additional (although rather mild) assumptions on the spatial and temporal structure of extremes.

In future work, space-time distance metrics other than the Euclidean distance could be used to define the space-time neighborhoods N (s, t). To account for orographic structures, the crossing distance could be used, which includes a vertical component related to the crossing of crests and valleys [START_REF] Gottardi | Statistical reanalysis of precipitation fields based on ground network data and weather patterns: Application over French mountains[END_REF]. Instead of fitting the marginal tail parameters separately for each spatial grid cell, a generalized additive regression approach could be implemented to borrow information from nearby sites [START_REF] Gardes | Conditional extremes from heavy-tailed distributions: an application to the estimation of extreme rainfall return levels[END_REF][START_REF] Carreau | Partitioning into hazard subregions for regional peaksover-threshold modeling of heavy precipitation[END_REF]Castro-Camilo et al., 2020, e.g.,). In addition, more sophisticated validation methods for POT stability in large dimensions could be studied. Finally, we note that there are events such as karstic aquifer floods where not only the extreme rainfall but also dry and moderate rainfall periods have to be considered. By extending ideas in [START_REF] Cantet | Using a rainfall stochastic generator to detect trends in extreme rainfall[END_REF] and in [START_REF] Yiou | Anawege: a weather generator based on analogues of atmospheric circulation[END_REF], we plan to implement our method as part of a spatial precipitation generator that simulates complete rainfall series. Rain-flow models will then be fed by simulated series from a precipitation generator, and we will be able to study the impact of the flood by applying risk measures to the outputs of rain-flow models.

Figure 1 :

 1 Figure 1: Empirical return levels at 98% level for each grid cell in our study area from 1997 to 2007. Grey and black contour lines indicate altitude (400 and 800 m respectively).

Figure 2 :

 2 Figure 2: Distribution function G for p 0 = 0.7.

Figure 4 :

 4 Figure 4: Exploratory plots for threshold stability using two cost functionals: spatio-temporal median (top row); spatio-temporal maximum (bottom row). Summary statistics: f u (Y i ) for u = 0 (first panel); same for u = 0.95-quantile (second panel); sd(Y i ) (third panel). QQ-plots (fourth panel) of observed cost functionals at uniform scale 1 -u/ *i with pointwise confidence interval at 95% (dashed lines) .

Figure 5 :

 5 Figure 5: Original precipitation data X(s, t) (left column) and uplifted episodes W (s, t) (right column) based on the spatio-temporal median

  . First row: most extreme episode, shown for t =2000-07-25, 18:00:00; second row: second most extreme episode, shown for t =2003-09-22, 17:00:00. Red dots indicate the center s i of the neighborhood N (s i , t i ) with maximum cost functional. Grey and black contour indicate altitude.

Figure 6 :

 6 Figure 6: Original precipitation data X(s, t) (left column) and uplifted episodes W (s, t) (right column) based on the spatio-temporal maximum

Table 1 :

 1 Last time steps of the 6 most extreme space-time episodes after declustering, using two different cost functionals.

	Episode Spatio-temporal median	(1) s,t	Spatio-temporal maximum	(2) s,t
	1st	2000-07-25 18:00:00		1999-09-14 10:00:00
	2nd	2003-09-22 17:00:00		2005-09-07 01:00:00
	3rd	2002-09-08 07:00:00		1999-09-03 15:00:00
	4th	1999-11-12 15:00:00		1999-08-28 22:00:00
	5th	1999-08-08 14:00:00		2001-07-07 06:00:00
	6th	1997-08-11 05:00:00		2006-10-12 02:00:00

Acknowledgements We thank Météo-France for providing us the data set. We are grateful to the LabEx NUMEV and the French national program LEFE/INSU for financial support.