Spatio-temporal modeling of urban extreme rainfall events at high resolution

Spatio-temp oral mo deling of urban extreme rainfall ev en ts at high resolution Chlo e Serre-Com b e 1,* , Nicolas Mey er 1 , Thomas Opitz 2 , and Gwladys T oulemonde 1 1 IMA G, Univ. Montpellier, CNRS, Inria LEMON team, Mon tp ellier, F rance 2 BioSP , Inrae, A vignon, F rance * Corresp onding author: chloe.serre-combe@umontpellier.fr Abstract Mo deling precipitation and its accumulation o ver time and space is essen tial for flo o d risk assessmen t. W e here analyze rainfall data collected o ver sev eral years through a micro- scale precipitation sensor netw ork in Mon tp ellier, F rance, b y the OMSEV observ atory . A no vel spatio-temp oral sto chastic mo del is prop osed for high-resolution urban rainfall and com bines realistic marginal b eha vior and flexible extremal dep endence structure. Rain- fall in tensities are describ ed b y the Extended Generalized P areto Distribution (EGPD), capturing b oth mo derate and extreme even ts without threshold selection. Based on spa- tial extreme-v alue theory , dep endence during extreme episodes is mo deled b y an r -Pareto pro cess with a non-separable v ariogram including episo de-sp ecific advection, allo wing the displacemen t of rainfall cells to b e represented explicitly . P arameters are estimated by a composite likelihoo d based on joint exceedances, and empirical adv ection velocities are deriv ed from radar reanalysis. The mo del accurately repro duces the spatio-temp oral struc- ture of extreme rainfall observ ed in the Montpellier OMSEV netw ork and enables realistic sto c hastic scenario generation for flo o d risk assessmen t. Keyw ords : Extreme v alue theory , EGPD, high spatio-temp oral resolution, r -Pareto pro- cess, adv ection 1. In tro duction Rainfall mo deling is essential in ev aluating flo o d risks and analyzing the impact of urbaniza- tion, particularly in areas suc h as the Cév ennes moun tain range and the coastal agglomeration of the city of Mon tp ellier in the south of F rance. The Cévennes massif is kno wn for its sp ecific climatological features, exp eriencing some of the most extreme rainfall in F rance, esp ecially during the autumn season. The topography of these moun tains acts as a barrier to relativ ely 1 w arm and h umid air masses from the Mediterranean sea making landfall, whic h causes atmo- spheric instability and orographic uplift. These conditions lead to intense rainfall ev ents that can sometimes exceed 500 mm within a 24 -hour p erio d, as evidenced b y historical o ccurrences suc h as the significant flo o ds in Anduze in 2002 with 680 mm of rain in 24 hours ( F ouc hier et al. 2004 ), which is more than the ann ual av erage rainfall amount in P aris (see Hemp elmann et al. 2021 ). Suc h even ts are characterised b y a large amoun t of rainfall o ccurring ov er short durations and are often highly lo calized, p ossibly leading to sev ere flo o ding, esp ecially flash flo o ds. The city of Montpellier also faces such extreme rainfall ev ents, for instance in September 2014 with a total of 260 mm of rain falling in only four hours (see Brunet et al. 2018 ), or more recently with a total of 100 mm in sev en hours recorded on December 22, 2025 by the Meteo F rance station of Mon tp ellier airp ort ( meteo.data.gouv.fr ). The cit y’s pro ximit y to the Cévennes mountains exacerbates these flo o ding risks, underscoring the necessity for precise rainfall mo deling and simulations for effectiv e flo o d risk assessment and urban planning. The highly lo calized nature of these weather even ts emphasizes the imp ortance of understanding rainfall b ehavior at a high spatio-temp oral resolution. In this context, a rain-gauge netw ork has b een deplo y ed in Mon tp ellier by the OMSEV ob- serv atory , detailed in Section 2.1 . This netw ork pro vides high-resolution rainfall measuremen ts that aim to help in the analysis of urban flo o d pro cesses and their impacts within a h ydrological mo del dev elop ed at Inria Mon tp ellier. Ho w ever, these data present several challenges, including a relatively short observ ation p erio d, a high prop ortion of missing v alues, and the sparsit y of non-zero rainfall observ ations, which complicate statistical mo deling and inference. Addressing these challenges constitutes a cen tral ob jective of this pap er in order to dev elop a sto chastic pre- cipitation generator that can b e used to simulate realistic rainfall fields at high spatio-temp oral resolution, while accoun ting for advection effects driving the spatial displacement of rainfall ev ents. This pap er is organized as follows. Section 2 presents the data used in this study , together with related w ork and the metho dological background for the prop osed sto c hastic precipitation generator. Section 3 in tro duces the complete statistical mo del for extreme rainfall, combining an EGPD for the marginal distributions with the r -P areto pro cess to capture spatio-temp oral dep endence. Estimation pro cedures for the mo del parameters are detailed in Section 4 , includ- ing the incorp oration of advectio n effects. Section 5 applies the prop osed framework to the OMSEV rainfall data to estimate the mo del parameters. Section 6 presents the simulations of extreme rainfall fields based on the fitted mo del ov er the OMSEV net work, ev aluates its p erformance in repro ducing the spatio-temp oral structure of extreme rainfall even ts. Finally , Section 7 concludes the pap er with a discussion of the results and p oten tial future research directions. 2 2. Data and metho dological bac kground In this section, we present the datasets used in this work and provide the related w orks and metho dological background for the prop osed spatio-temp oral stochastic precipitation generator. 2.1 Case study and datasets Our study fo cuses on the OMSEV 1 rain gauge netw ork, lo cated within the V erdanson riv er catc hment, a tributary of the Lez riv er (see Figure 1 ), flowing in to the Mediterranean Sea around 10 km South of Montpellier. This netw ork consists of 20 tipping-buc ket rain gauges installed by the urban observ atory of HydroSciences Montpellier ( Finaud-Guy ot et al. 2023 ). Although tipping-buck et rain gauges pro vide high-resolution measuremen ts, they can lead to an underestimation of rainfall intensities. A prop er calibration w as therefore done b y hydro- logists of each rain gauge in order to obtain the most reliable estimates of rainfall intensit y ( Robin 2021 ). W e exclude three rain gauges installed only recently that do not yet pro vide sufficien tly long data series (CINES, Briv es and Hydrop olis, see Figure 1b ). The dataset cov ers the p erio d from September 2019 to January 2025 and provides rainfall measurements at high temp oral resolution, recorded every min ute leading to discretized low v alues. In order to reduce measuremen t errors, to decrease data volume and to a void o v erly discretized v alues, the data are aggregated in 5 -minute interv als whic h preserv e a high temp oral resolution. It allows to reduce the n umber of missing v alues that are considered as 0 ov er these in terv als when there is at least one non-missing v alue, which are more frequent at the original minute resolution. At this temp oral resolution, approximativ ely 23% of data are missing, due to activ ation dates of eac h rain gauges and measurement errors and are not imputed in the present w ork. With this fine resolution, a kind of discretization is still present for v ery low v alues, and the prop ortion of non-zero v alues is very lo w, represen ting only 1 . 2% of the data. Spatial granularit y is also v ery fine with inter-site distance ranging from 77 to 1531 meters (without remov ed gauges). Data from urban microscale rain gauge netw orks like this are very scarce (e.g. Casas et al. 2010 for Barcelona, Spain) but can provide imp ortan t insights rainfall patterns known to b e highly v ariable and complex ev en at very small spatio-temp oral scales. 1 Observ atoire Montpelliérain et au Sud de l’Eau dans la Ville 3 (a) V erdanson river catc hment (b) Rain gauges Figure 1: Lo cation of the 20 rain gauges of the OMSEV netw ork o ver the V erdanson river catc hment in Mon tp ellier The rain gauges of the OMSEV net work are irregularly distributed. A more regular spatial grid would b e desirable for rainfall field mo deling. In addition, the relatively short temp oral co verage of the OMSEV dataset may not fully capture imp ortant features and the v ariability of extreme rainfall ev ents. T o address these limitations, we enric h these data with the F rench COMEPHORE 2 mosaic from Météo F rance (see T abary et al. 2012 ). COMEPHORE is a reanalysis pro duct combining radar and rain gauge observ ations ov er F rance, pro vided on a regular grid with a spatial resolution of 1 km 2 p er pixel and an hourly temp oral resolution. The dataset spans from 1997 to 2024. W e restrict the analysis to after 2008 due to a change in the reanalysis mo del in 2007 inducing p ossible distributional changes. W e extract the data o ver the same area as the OMSEV net work, corresp onding to 111 pixels within a 5 km radius of the net work cen troid, which is sufficien t to represen t the lo cal spatial v ariabilit y . Only fiv e of these pixels con tain OMSEV rain gauges or sev en pixels when considering also remov ed gauges (see Figure 2 ). 2 COm binaison en vue de la Meilleure Estimation de la Précipitation HOraiRE. The dataset is av ailable on AERIS platform ( h ttps://radarsmf.aeris-data.fr/ ) 4 Figure 2: Lo cation of the 20 rain gauges of the OMSEV netw ork (red dots) ov er the COMEPHORE reanalysis (blue grid). 2.2 State of the art and prop osed metho dological framew ork W e aim to mo del and simulate realistic high-resolution spatio-temp oral rainfall fields ov er the Mon tp ellier area and put fo cus on extreme precipitation even ts relev ant for flo o d risk analysis. Ph ysics-based global and regional climate simulation mo dels provide v aluable information at relativ ely large scales, esp ecially for continen tal to global analyses, with resolutions around 10 to 100 km. Ho wev er, their spatial and temp oral resolutions remain to o coarse to accurately represen t short-duration, high-intensit y , and highly lo calized rainfall ev ents, desp ite the recen t dev elopment of con vection-permitting sim ulations models with kilometer resolution ( Lucas- Pic her et al. 2021 ). Moreov er, their high computational cost prev en ts the generation of large ensem bles required for probabilistic risk assessment. In this con text, sto chastic weather generators (SW Gs) constitute a complemen tary and computationally efficient alternativ e ( Wilks and Wilb y 1999 ). Within this class of mo dels, sto c hastic precipitation generators (SPGs) sp ecifically focus on rainfall and aim to reproduce its main statistical characteristics, including o ccurrence, in tensit y , duration, and spatial structure. Designing SPGs is challenging due to the highly non Gaussian distribution of rainfall, t ypically with heavy tails and a large prop ortion of zeros, and the need to pro vide simulations that are explicit and coheren t in space and time, for using them as inputs of numerical impact mo dels. Early dev elopments and metho dological o v erviews highlight the div ersity of SPG frameworks and their ability to capture rainfall v ariabilit y across a range of spatial and temp oral scales ( Allard et al. 2015 ). In particular, geostatistical sim ulation approaches hav e b een shown to effectiv ely repro duce spatial rainfall structures ( Benoit and Mariethoz 2017 ). Several studies ha ve sho wn that SPGs can b e dev elop ed at sub-kilometer and sub-hourly resolutions, and ha ve highlighted the imp ortance of advection or wind consideration for repro ducing fine-scale 5 space–time rainfall dynamics ( Leblois and Creutin 2013 , Schleiss et al. 2014 , Huser and Davison 2014 ). More recently , Benoit et al. 2018 prop osed a meta-Gaussian framew ork to generate realistic high-resolution rainfall fields, considering advection effects. Rainfall intensities exhibit differen t regimes, ranging from dry p erio ds to mo derate and ex- treme precipitation. Since extreme rainfall plays a k ey role in flo o ding, extreme v alue theory (EVT) provides a natural framework for its statistical mo deling. Classical EVT approaches include blo ck maxima mo dels based on the Generalized Extreme V alue (GEV) distribution and P eaks-Over-Threshold (POT) mo dels based on the Generalized P areto Distribution (GPD). In practice, POT approaches are often preferred as they retain more extreme observ ations and offer greater flexibility , esp ecially for mo deling temp oral patterns. How ev er, their p erformance strongly dep ends on the choice of a threshold, which is kno wn to b e challenging and to sig- nifican tly affect inference results, as highligh ted in reviews b y Coles et al. 2001 and Scarrott and MacDonald 2012 . T o av oid explicit threshold selection and to mo del the en tire range of rainfall intensities, mixture mo dels ha ve b een prop osed, combining a distribution for moderate rainfall with a GPD tail ( Carreau and Bengio 2009 ). These approac hes still require a transition threshold and induce a dep endence b et ween bulk and tail b eha vior ( Na v eau et al. 2016 ). T o o vercome these limitations, the class of Extended Generalized Pareto Distributions (EGPDs) ( P apastathop oulos and T a wn 2013 , Na v eau et al. 2016 ) can be used. It provides a unified mo del for strictly p ositiv e v ariables without threshold selection. EGPDs allo w a flexible representation of b oth mo derate and extreme rainfall and hav e already prov en their efficiency in precipitation mo deling applications ( Haruna et al. 2023 , Carrer and Gaetan 2022 ). While marginal mo deling is essen tial, realistic rainfall simulations also require an adequate represen tation of spatio-temp oral dep endence. In the context of extremes, this dep endence is commonly describ ed using asymptotic mo dels suc h as max-stable and P areto pro cesses ( Da vison et al. 2012 , F erreira and De Haan 2014 ), whic h provide a natural framew ork for spatio-temporal extreme even ts and allo w extrap olation b eyond observ ed extremes. How ev er, max-stable mo d- els are computationally demanding and difficult to handle in high-dimensional environmen tal settings ( Davison et al. 2013 ). T o o vercome these limitations, we rely on the Brown-Resnic k r -P areto pro cess ( De F ondeville and Davison 2018 , Dom bry et al. 2024 ), whic h arises as a limiting mo del for threshold exceedances of sto chastic pro cesses. This mo del preserves the de- p endence structure of the max-stable Bro wn-Resnick pro cess while offering greater flexibility and impro ved computational efficiency compared with max-stable mo dels, making it particu- larly suitable for spatio-temp oral mo deling of extreme rainfall. In the present w ork, w e prop ose a nov el framework for the simulation of high-resolution extreme rainfall fields that mo dels mo derate and extreme marginals with a spatio-temp oral extremal dep endence. The prop osed approac h sp ecifically targets extreme rainfall episo des re- lying on on the Extended Generalized P areto Distribution (EGPD) of Na veau et al. 2016 for the marginal modeling of rainfall in tensities. In con trast to recen t high-resolution sto c hastic precip- itation mo deling studies, our approach fo cuses on repro ducing the spatio-temp oral dep endence of extremes through the estimation of v ariogram parameters using a comp osite likelihoo d based 6 on joint exceedances from identified extreme episo des. A k ey contribution of this work is the explicit incorp oration of even t-sp ecific advection effects in to the dep endence structure via a non-separable spatio-temp oral v ariogram. While similar ideas hav e b een explored for max- stable pro cesses using wind information at coarser resolutions and ov er larger spatial domains ( Huser and Da vison 2014 ), our approac h, dev elop ed on r -P areto pro cesses, fo cuses on threshold exceedances and high-resolution urban rainfall data. The episo de-sp ecific advection term ex- plicitly accounts for the displacemen t of rainstorms o v er time. This framework enables the sim ulation of realistic, high-resolution, mo ving extreme precipitation fields o ver a giv en spatial domain. W e apply the prop osed methodology to rainfall sim ulation on a grid cov ering the OMSEV net w ork in the Montpellier area. 3. Spatio-temp oral mo deling for extreme precipitation Let X = { X s ,t | ( s , t ) ∈ S × T } b e a nonnegativ e spatio-temp oral random field represen ting the rainfall in tensity for all lo cations s in a spatial domain S ⊂ R 2 and all times t in a temp oral domain T ⊂ R + . Let Λ S ⊂ R 2 and Λ T ⊂ R be the sets of spatial and temp oral lags, respectively . In the following, we present the prop osed mo del the construction of the sto chastic pro cess X for extreme rainfall, whic h combines a flexible marginal distribution using EGPD and a spatio- temp oral dep endence structure based on the r -Pareto pro cess. 3.1 Marginal rainfall distributions Rainfall at high temp oral resolution is c haracterized by a large prop ortion of zero v alues, cor- resp onding to dry p erio ds. T o accoun t for this structure, rainfall at each site is mo deled using a mixed distribution with a p oin t mass at zero and a distribution for p ositive v alues. 3.1.1 Mo deling p ositiv e rainfall in tensities Let X s b e the rainfall intensit y at a fixed lo cation s . T o mo del the whole marginal distribution of p ositive rainfall intensities, we use the Extended Generalized P areto Distribution (EGPD). The EGPD is constructed as an extension of the classical Generalized P areto Distribution (GPD), whic h corresp onds to the limit distribution of threshold exceedances. F or a sufficiently large threshold u , the distribution of the exceedances can b e appro ximated b y a GPD with shap e parameter ξ and scale parameter σ u : P ( X s − u > y | X s > u ) ≈ H ξ  y σ u  =         1 + ξ y σ u  − 1 /ξ + , if ξ  = 0 , exp  − y σ u  , if ξ = 0 , y > 0 , where ( a ) + = max( a, 0) . The EGPD is then defined b y considering a p olynomial transformation P ( x ) = x κ , with κ > 0 of the GPD. This is the simplest and most efficient choice among other transformations prop osed by Na veau et al. 2016 . The EGPD preserves the GP b ehavior 7 of the upp er tail while extending the mo del to mo derate rainfall v alues through a smo oth transformation of the GPD distribution. The corresp onding shap e is driven b y the parameter κ . W e define F EGPD as the EGPD cum ulative distribution function with parameters ( ξ , σ, κ ) , indep enden t of the threshold u . Then X s | X s > 0 ∼ F EGPD ( ξ , σ, κ ) . 3.1.2 Ov erall marginal mo del Let p 0 b e the marginal probabilit y of zero rainfall for each lo cation s . Giv en that X s ,t > 0 , p ositiv e rainfall amounts are assumed to follo w an EGPD. F or a fine spatio-temp oral resolution it is reasonable to consider constant parameters ( ξ , σ, κ ) o ver the spatial domain and o ver time. Then the EGPD distribution is the same for all s ∈ S and t ∈ T . The resulting marginal distribution for x ∈ R function is given b y F ( x ) =      p 0 , x = 0 , p 0 + (1 − p 0 ) F EGPD ( x ) , x > 0 . This expression allows the dry-p erio d frequency and the distribution of p ositiv e rainfall intens- ities to b e mo deled separately , while k eeping a coheren t marginal mo del. 3.2 Spatio-temp oral extremal dep endence T o mo del the spatio-temp oral dep endence of extreme rainfall, we consider the r -P areto pro cess in tro duced by De F ondeville and Da vison 2018 . This process provides a natural framework for mo deling the dep endence structure of exceedances ab ov e high thresholds. W e consider a regularly v arying spatio-temp oral pro cess X in the maxim um domain of attraction of a max-stable pro cess, and a nonnegativ e and 1 -homogeneous risk function r , that is, a function satisfying r ( a X ) = a r ( X ) for any a > 0 . As sho wn by Dom bry et al. 2024 , the r -exceedances of X conv erge in distribution to an r -Pareto pro cess Y = { Y s ,t | ( s , t ) ∈ S × T } : u − 1 X s ,t | r ( X ) > u d − → Y s ,t , u → ∞ . T o mo del the dep endence structure of Y , we assume a Brown-Resnic k representation induced b y an underlying Gaussian pro cess W = { W s ,t , ( s , t ) ∈ S × T } with stationary increments. The dep endence structure of the process is then c haracterized through its spatio-temp oral v ariogram γ of the underlying Gaussian pro cess defined for an y ( h , τ ) ∈ Λ S × Λ T as γ ( h , τ ) = 1 2 V ( W s + h , t + τ − W s ,t ) . In this pap er w e consider the risk function r ( X ) = X s 0 ,t 0 where ( s 0 , t 0 ) ∈ S × T corresp onds to a fixed reference spatio-temp oral lo cation. The corresp onding r -Pareto pro cess can b e represented 8 as Y s ,t = R exp( W s ,t − W s 0 ,t 0 − γ ( s − s 0 , t − t 0 )) , where R follo ws a univ ariate Pareto distribution, that is, P ( R > v ) = v − 1 for v ≥ 1 . This approac h is relev ant since it links the v ariogram of the underlying Gaussian pro cess to the dep endence structure of the r -P areto pro cess Y . Th us it allo ws to capture the main features of the extremal dep endence structure of the pro cess X with resp ect to the conditional exceedance at ( s 0 , t 0 ) ∈ S × T given by the risk function. In EVT, extremal dep endence can b e summarized through the spatio-temp oral extremo- gram which corresp onds to the probability that an exceedance at one lo cation and time o ccurs giv en an exceedance at another lo cation and time. F ollo wing Coles et al. 1999 , the extremal dep endence co efficient is defined for spatial lag h ∈ Λ S and temp oral lag τ ∈ Λ T as χ ( h , τ ) = lim q → 1 P  X ∗ s ,t > q   X ∗ s + h , t + τ > q  , where X ∗ s ,t denotes the uniform marginal transformation of X s ,t . The co efficient χ ( h , τ ) ranges from 0 (asymptotic indep endence) to 1 (p erfect extremal dep endence). F or the chosen r -P areto pro cess, the extremal dep endence of X can b e summarized through the r -extremogram χ r ( h , τ ) = lim u →∞ P ( X s 0 + h ,t 0 + τ > u | X s 0 ,t 0 > u ) = lim u →∞ P ( X s 0 + h ,t 0 + τ > u ) , whic h describ es the extremal dep endence structure of X with resp ect to the conditional ex- ceedance. This implies that P ( X s 0 ,t 0 > u ) = 1 for large u . Under the Bro wn-Resnic k mo del represen tation, this co efficient dep ends directly on the v ariogram and can b e expressed as χ r ( h , τ ) = 2  1 − Φ  q 1 2 γ ( h , τ )  , (1) where Φ denotes the standard normal distribution function (see Buhl et al. 2019 ) and ( h , τ ) = ( s − s 0 , t − t 0 ) ∈ Λ S × Λ T reprensen ts the space-time lag with resp ect to the reference lo cation. This form ulation entails that small v alues of γ ( h , τ ) corresp ond to strong extremal dep endence, whereas large v ariogram v alues imply w eak dep endence. 3.3 Non-separable space-time v ariogram with adv ection A common approach to mo del spatio-temp oral dep endence consists in considering the v ariogram as a sum of purely spatial and purely temp oral components. In these so-called separable mo dels, the spatio-temp oral v ariogram is expressed as γ ( h , τ ) = γ S ( h ) + γ T ( τ ) , h ∈ Λ S , τ ∈ Λ T , where γ S and γ T are purely spatial and purely temp oral v ariograms, resp ectiv ely . Separable spatio-temp oral v ariogram mo dels pro vide a simple and reasonable description of 9 the o verall dep endence structure of rainfall data. Ho wev er, they fail to capture the complex in teractions b etw een space and time, esp ecially in the con text of meteorological phenomena and in particular for extreme rainfall even ts. Extreme precipitation is often driven b y dynamical mec hanisms, suc h as advection that induce strong interactions b etw een space and time. In meteorology , it corresp onds to the horizontal transp ort of prop erties (heat, moisture) b y the mo vemen t of air masses, such as wind and clouds. This phenomenon therefore influences the spatio-temp oral b ehavior of precipitations, but cannot b e adequately captured by a separable mo del and requires non-separable spatio-temp oral structures. T o construct a non-separable v ariogram accounting for advection, we start from a spatio- temp oral pro cess with stationary increments, whose spatial and temp oral v ariograms are, for spatial lag h ∈ Λ S and temp oral lag τ ∈ Λ T , γ S ( h ) = 2 β 1 ∥ h ∥ α 1 , γ T ( τ ) = 2 β 2 | τ | α 2 , with β 1 , β 2 > 0 and α 1 , α 2 ∈ (0 , 2] . Without advection, this leads to a separable spatio-temp oral v ariogram. A dvection is incorp orated through a deterministic space-time shift defined by a velocity v ector V ∈ R 2 , yielding the non-separable v ariogram γ ( h , τ ; Θ , V ) = 2 ( β 1 ∥ h − τ V ∥ α 1 + β 2 | τ | α 2 ) , where Θ = ( β 1 , β 2 , α 1 , α 2 ) . The corresp onding r -extremogram follo ws from Equation 1 and is denoted b y χ r ( h , τ ; Θ , V ) . 3.4 Sto c hastic precipitation generator The stochastic precipitation generator is constructed by combining the marginal mo del for rainfall amoun ts with the dep endence mo del for extremes based on the r -Pareto pro cess. The r -Pareto pro cess Y describ es the dep endence structure of extreme even ts on a stand- ardized P areto scale and do es not carry information on their real magnitude. T o rein tro duce the intensit y level of extreme episo des, simulations are rescaled by a high threshold u , leading to the spatio-temp oral pro cess Z = { Z s ,t = u Y s ,t | ( s , t ) ∈ S × T } . By construction, Z s ,t has P areto tails with P ( Z s ,t > z ) = u z , for z > u . T o com bine this dep endence mo del with the rainfall marginal distributions, the sim ulated v alues on the Pareto scale are mapp ed to the unit in terv al using a standardisation function G 10 follo wing P alacios-Ro dríguez et al. 2020 . The function G is defined as G ( x ) =                  0 , x < 0 , p 0 , x = 0 , p 0 + (1 − p 0 ) 2 4 x, 0 < x < 2 1 − p 0 , 1 − 1 x , x > 2 1 − p 0 . This transformation preserv es the tail b ehavior of the r -Pareto mo del and puts a p oin t mass at zero. Applying this transformation G to the pro cess Z , w e obtain the simulated v alues on the uniform scale, U s ,t = G ( Z s ,t ) , s ∈ S , t ∈ T . The sim ulated rainfall pro cess are obtained b y applying the in v erse marginal distribution func- tion, X s ,t = F − 1 ( U s ,t ) , s ∈ S , t ∈ T , where F − 1 ( x ) =      0 , x ≤ p 0 , F − 1 EGPD  x − p 0 1 − p 0 ; ξ , σ, κ  , x > p 0 . This approach allo ws the generation of spatio-temp oral rainfall fields repro ducing both the in tensity and the dep endence structure of extreme precipitation ev en ts. 4. Estimation of mo del parameters In this section, we presen t the estimation pro cedures for the parameters of the prop osed mo del, including the marginal parameters, the episo de selection pro cedure and the estimation of the v ariogram parameters with advection effects using a comp osite likelihoo d based on joint ex- ceedances. 4.1 Marginal parameter estimation The marginal parameters ( p 0 , ξ , σ, κ ) are estimated separately at eac h location s ∈ S and then a veraged ov er all lo cations to obtain a single set of parameters for the entire spatial domain. The probabilit y of zero rainfall p 0 is estimated b y the empirical prop ortion of zero v alues at lo cation s and the EGPD parameters ( ξ , σ , κ ) are estimated b y maximum likelihoo d using only p ositiv e rainfall v alues. 4.2 Selection of extreme episo des Let δ > 0 denote a fixed temp oral duration for extreme episodes. F or eac h conditioning exceedance o ccurring at lo cation s 0 ∈ S and time t 0 ∈ T , we define an extreme episo de from 11 the pro cess X with a duration δ and starting at time t 0 : E = { X s ,t | X s 0 ,t 0 > u, s ∈ S , t ∈ [ t 0 , t 0 + δ [ } . Extreme episo des are selected following a pro cedure describ ed by Palacios-Rodríguez et al. 2020 , whic h ensures a weak dep endence b etw een episodes. First w e iden tify all possible conditioning exceedances ab o v e a high threshold u , i.e. all spatio-temp oral lo cations ( s 0 , t 0 ) ∈ S × T suc h that X s 0 ,t 0 > u . Then, these conditioning exceedances are ordered by time and pro cessed in c hronological order. An episo de with conditioning lo cation ( s 0 , t 0 ) ∈ S × T is retained if for ev ery already selected episo de with conditioning lo cation ( s ′ 0 , t ′ 0 ) ∈ S × T we hav e ∥ s ′ 0 − s 0 ∥ ≥ d min or | t ′ 0 − t 0 | ≥ δ, where d min > 0 is a minim um spatial separation and the temp oral separation is giv en b y the episo de duration δ . In other words, all p oten tial episo des falling within the spatio-temp oral neigh b orho o d of an already selected episo de are discarded from further consideration. The pro cedure stops when either a fixed maxim um n umber of episo des is reac hed or no exceedance ab o v e u remains. The resulting set of selected episo des is denoted by E . 4.3 Comp osite likelihoo d based on exceedances F or each episo de E ∈ E , spatio-temp oral neigh b ors are defined relatively to the conditioning p oin t ( s 0 , t 0 ) ∈ S × T . The corresp onding neighborho o d, i.e. the set of all spatio-temp oral lo cations at given lags from the conditioning p oint, is defined as N E ( h , τ ) = { ( s , t ) ∈ S × T | ∥ s − s 0 ∥ ∈ C h , | t − t 0 | ∈ C τ } , where C h and C τ denote predefined spatial and temp oral lag classes. In practice, if a regular spatio-temp oral grid is considered, these classes are defined as direct spatial and temporal lags sets. F or each lo cation ( s , t ) ∈ N E ( h , τ ) , we define the indicator v ariable of join t exceedances as k E ( s , t ) = 1 { X s ,t >u,X s 0 ,t 0 >u } = 1 { X s ,t >u } where the conditioning exceedance X s 0 ,t 0 > u holds by construction. With the r -Pareto mo del, k E ( s , t ) follows a Bernoulli distribution with success probability χ r ( h , τ ; Θ , V E ) with paramet- ers Θ = ( β 1 , β 2 , α 1 , α 2 ) and episo de-sp ecific velocity vectors V E . The num b er of exceedances at lag ( h , τ ) during an episo de is then given by K E ( h , τ ) = X ( s ,t ) ∈N E ( h ,τ ) k E ( s , t ) . 12 After aggregating o ver all episo des, w e obtain the total num b er of exceedances and a total n umber of trials for lag ( h , τ ) that are giv en b y K tot ( h , τ ) = X E ∈E K E ( h , τ ) , These aggregated exceedance counts pro vide the basis for the co mp osite lik eliho o d used to estimate the v ariogram parameters. Under approximate indep endence b et ween episo des and for large u , Bernoulli lik eliho o d contributions are used to construct the comp osite likelihoo d. The resulting comp osite log-likelihoo d satisfies ℓ C (Θ) ∝ X E ∈E X ( h ,τ ) ∈ Λ S × Λ T X ( s ,t ) ∈N E ( h ,τ ) k E ( s , t ) log χ r ( h , τ ; Θ , V E ) +  1 − k E ( s , t )  log  1 − χ r ( h , τ ; Θ , V E )  . (2) The final inference step consists in maximizing the comp osite log-lik eliho o d ℓ C (Θ) with resp ect to the v ariogram parameters Θ . 4.4 A dv ection and parametric transformation A dvection velocities ma y b e either kno wn and directly sp ecified in the mo del, or estimated from the data when unav ailable. In our case, episo de-sp ecific empirical v elo cit y v ectors are inferred from rainfall barycen ter displacemen ts (see Section 5 ), and their estimation uncertain ty is handled with a parametric transformation acting on their magnitude. F or eac h episo de an empirical v elo city vector V emp is estimated from the data. This empirical velocity vector is not used directly in the mo del, but is mapp ed to an effectiv e v elo city v ector with the transformation A ( · ) . The transformation acts only on the magnitude of the velocity vector, while preserving its direction. F or any non-zero empirical v elo cit y v ector V emp , A ( V emp ) = η 1 ∥ V emp ∥ η 2 V emp ∥ V emp ∥ . The advection asso ciated with episo de E ∈ E is given b y V E = A ( V emp E ) . The apparen t episo de-sp ecific advection V E is therefore fully determined by the parameters η 1 and η 2 . The new extended parameter v ector b ecomes e Θ = ( β 1 , β 2 , α 1 , α 2 , η 1 , η 2 ) for all episo des and is es- timated b y maximizing the comp osite log-likelihoo d ℓ C ( e Θ) giv en in (( 2 )) for a sp ecific episo de E , expressed as χ r ( h , τ ; e Θ , V emp E ) . 4.5 V alidation on sim ulated r -P areto pro cesses The prop osed inference metho dology is v alidated using sim ulated r -Pareto pro cesses with known parameters. The simulation of r -Pareto pro cesses is based on the sp ectral representation of spatio-temp oral Brown-Resnic k pro cesses. Sp ecifically , w e adapt the sim ulation algorithm pro- p osed by Leb er 2015 , which itself relies on exact simulation tec hniques for spatial max-stable 13 pro cesses developed in Dombry et al. 2016 . The r -Pareto pro cesses are simulated on a regular grid of n spatial lo cations ov er a p erio d of δ time steps. A threshold u = 1 is used to define the r -exceedances, and the risk function is c hosen as r ( X ) = X s 0 ,t 0 , where s 0 is randomly selected among the spatial lo cations and t 0 = 0 corresp onds to the first time step of each sim ulated episo de. Differen t sets of simulations are generated using known v ariogram parameters e Θ . Three simulation settings are considered to assess the p erformance of the prop osed infer- ence metho dology under differen t adv ection configurations. F or each setting, 50 indep enden t sim ulations of r -P areto pro cesses are generated, each consisting of 500 replicates of n = 49 spa- tial lo cations ov er δ = 24 or 12 time s teps. Figure 3 illustrates the results obtained when the adv ection v aries across episo des and the parameters ( η 1 , η 2 ) are jointly estimated together with the v ariogram parameters. Despite the higher n umber of parameters, ( β 1 , β 2 , α 1 , α 2 ) remain w ell iden tified. As exp ected, higher v ariabilit y is mainly observ ed for the adv ection parameters ( η 1 , η 2 ) , whic h remain reasonably w ell estimated but with a comp ensatory effect b et w een them, i.e. differen t combinaison can lead to similar velocity vectors. Figure 4 depicts the estima- tion results obtained when the adv ection parameters ( η 1 , η 2 ) are fixed to their true v alues. In this case, the estimation results are comparable to those obtained when ( η 1 , η 2 ) are estimated, with accurate recov ery of the v ariogram parameters. How ev er, fixing the advection parameters can reduce the v ariability of the estimates, esp ecially when dealing with real data where the adv ection estimation ma y b e more uncertain. 0.25 0.50 0.75 β 1 β 2 α 1 α 2 Parameters Estimated values (a) β 1 = 0 . 3 , β 2 = 0 . 6 , α 1 = 0 . 3 , α 2 = 0 . 8 0.0 2.5 5.0 7.5 10.0 η 1 η 2 Parameters Estimated values (b) η 1 = 1 . 6 , η 2 = 5 . 2 Figure 3: Estimations of v ariogram parameters with maxim um likelihoo d optimization on r - P areto simulations with 49 sites and 24 time observ ations. The true parameters are indicated b y red crosses. Here random advection by replicate (episo de) is considered. 14 0.5 1.0 1.5 β 1 β 2 α 1 α 2 Parameters Estimated values (a) β 1 = 0 . 4 , β 2 = 0 . 2 , α 1 = 1 . 5 , α 2 = 1 0.00 0.25 0.50 0.75 1.00 β 1 β 2 α 1 α 2 Parameters Estimated values (b) β 1 = 0 . 4 , β 2 = 0 . 8 , α 1 = 0 . 2 , α 2 = 0 . 7 Figure 4: Estimations of v ariogram parameters with maxim um likelihoo d optimization on r - P areto simulations with 49 sites and 12 time observ ations. The true parameters are indicated b y red crosses. W e use a random advection b y replicate (episo de) with fixed η 1 = 0 . 5 and η 2 = 1 . 6 . 5. Application to rainfall data In this section, the prop osed spatio-temp oral sto c hastic precipitation generator is applied to rainfall data from the OMSEV net work in Montpellier. The marginal distributions are estim- ated, and the use of a common distribution is justified. The dep endence structure of extremes is in vestigated by ev aluating the adequacy of the non-separable v ariogram mo del, describing the selection of extreme episo des, and the adv ection estimation pro cedure is explained. This is follo wed by the inference of the spatio-temp oral extremal dep endence structure using b oth COMEPHORE and OMSEV datasets. Finally , the p erformance of the sto c hastic precipitation generator is assessed by comparing the sim ulated and observ ed rainfall fields at OMSEV sites, with a particular fo cus on the representation of extreme rainfall even ts. A final sto chastic precipitation generator is constructed on a regular grid cov ering the entire area of in terest. 5.1 Estimation of marginal distributions (OMSEV) The EGPD mo del is fitted to rainfall data from the OMSEV net w ork in order to assess its abilit y to represen t the marginal distribution of p ositive rainfall intensities at each site. T o accoun t for the discretization induced b y the finite measurement precision of the rain gauges, a small left-censoring is introduced, chosen according to a lo cal go o dness-of-fit criterion based on the Ro ot Mean Square Error (RMSE), following Haruna et al. 2023 . P ossible censoring thresholds are restricted to multiples of the gauge precision, p ≈ 0 . 2153 mm, i.e. v alues of the form k × p with k ∈ N . The final censoring threshold is set to p for most sites, and to 2 p for three sites (CEFE, Archie, and CNRS). The quan tile-quan tile plots display ed in Figure 5 illustrate an ov erall go o d agreemen t b etw een the fitted EGPD and the empirical rainfall distributions for four sites of the OMSEV netw ork. The same level of agreement is observed at all other sites (not shown). No systematic deviations are observ ed in the upp er tail, suggesting that 15 the extreme rainfall b eha vior is well captured b y the EGPD mo del. The distribution of the estimated EGPD parameters across the OMSEV net w ork is summarized in Figure 6 . The shap e parameter ξ is p ositive at all sites, with v alues ranging approximately b et ween 0 . 25 and 0 . 30 , indicating heavy-tailed rainfall intensit y distributions. Estimated scale parameters σ mostly lie b et w een 0 . 60 and 0 . 75 , reflecting mo derate v ariability in rainfall in tensities across sites. It is consisten t with short-duration rainfall data, which t ypically exhibit low er v ariabilit y compared to longer-duration accum ulations. Finally , the low er-tail parameters κ generally fall b etw een 0 . 20 and 0 . 25 , indicating hea vier mass near small rainfall intensities relatively to the central part of the distribution. This is consistent with the short-duration nature of the data that leads to a larger proportion of small rainfall v alues compared to longer-duration accum ulations. The largest κ outlier corresp onds to the Arc hie site. This site is differen t from the others b ecause of an interruption of measuremen ts in 2023, whic h can explain differences in parameter estimates. Ov erall, these results seem to indicate that the EGPD provides a relev ant representation of the marginal rainfall distributions at all sites in the netw ork. (a) CRBM, left-censoring at 0 . 2153 mm (b) CEFE, left-censoring at 0 . 4306 mm (c) Archie, left-censoring at 0 . 4306 mm (d) UM, left-censoring at 0 . 2153 mm Figure 5: Quantile-quan tile plot of the EGPD fitting with b o otstrap confidence interv als on four rain gauges of the OMSEV net work with site-sp ecific left-censoring. 16 0.0 0.5 1.0 1.5 κ σ ξ P arameter V alue Figure 6: Boxplot of EGPD parameter estimates across the OMSEV netw ork with site-sp ecific left-censoring. The use of common EGPD parameters across all sites is motiv ated by the absence of clear evidence of spatial heterogeneit y in extreme rainfall marginals. Based on these complemen tary diagnostics, a common marginal mo del is assumed for the OMSEV netw ork. Global Moran’s I test is applied to the estimated EGPD parameters to detect p ossible spatial auto correlation. F or all parameters, the associated p -v alues are large, indicating that the null hypothesis of no global spatial auto correlation cannot b e rejected for eac h parameter(see T able 1 ). The spatial homo- geneit y of extreme rainfall is examined b y comparing the distributions of exceedances b et ween sites using the pairwise Anderson–Darling test. The results indicate that, for a large ma jority of site pairs, the null h yp othesis of identical upp er-tail distributions cannot b e rejected based on the Anderson–Darling test. Sp ecifically , for approximately 76% of the 136 pairwise com- parisons, no significant difference is detected in the tail b ehavior ab ov e the q = 0 . 90 quan tile. T aken together, these results suggest that the marginal distributions of extreme rainfall are assumed to b e approximately iden tical across the OMSEV net w ork. T able 1: Moran’s I test for global spatial auto correlation in EGPD parameters. The test is p erformed using the 4 nearest neighbors for each site. Rep orted p -v alues corresp ond to tw o- sided tests assessing the null h yp othesis of no global spatial auto correlation. P arameter Moran’s I p -v alue ξ − 0 . 211 0.252 σ − 0 . 002 0.446 κ − 0 . 098 0.765 The marginal common parameters for the OMSEV netw ork are estimated b y fitting an EGPD across all sites together, with left censoring at the gauge precision p , and by join tly es- timating the probabilit y of zero rainfall. The final parameter estimates are rep orted in T able 2 . These will b e used as common marginal parameters for the final sto c hastic precipitation gener- ator and for the transformation of simulated v alues from the P areto scale to the rainfall scale (see Section 3 ). 17 T able 2: Estimated common marginal parameters for the OMSEV net work. P arameter b p 0 b ξ b σ b κ Estimate 0 . 989 0 . 262 0 . 591 0 . 270 5.2 Empirical evidence of space-time non-separabilit y W e inv estigate the space-time structure of extremal dependence using empirical spatio-temp oral extremal v ariogram estimated from the extremogram in Equation (( 1 )). Under a separable space-time structure, the spatial v ariogram curv es are exp ected to b e prop ortional across tem- p oral lags. Figure 7 sho ws that the empirical v ariogram exhibits clear changes in curves shap e across differen t temp oral lags, indicating an in teraction b etw een spatial and temp oral dep end- ence. This b ehavior provides strong empirical evidence against a separable spatio-temp oral v ariogram for extreme rainfall in the OMSEV data. This justifies the non-separable v ariogram in tro duced in our mo del. Figure 7: Empirical spatio-temp oral v ariogram of the extreme OMSEV data for differen t tem- p oral lags (in 5 minutes) o ver the spatial distance (in meters) and quan tile q = 0 . 95 . 5.3 Episo de configuration Extreme episo des are identified by selecting conditioning p oints exceeding a spatio-temp oral quan tile q , a duration δ and a minimum spatial separation d min b et w een episode. It aims to ensure weak dep endence b et w een selected episo des while retaining a sufficient num b er of episo des for inference. T o select an appropriate quan tile threshold q for defining extreme episo des in the OMSEV data, w e examine the num b er of joint exceedances observ ed at v arious spatial distances and temp oral lags. F or this purp ose, we compute the num b er of joint exceedances ab o ve differen t quan tiles q and w e keep the quantile q = 0 . 95 as a compromise b et ween the representation of 18 extremes and sample size requirements. It allo ws a minimum of 40 spatial join t exceedances p er spatial lag and a minimum of 20 temporal join t exceedances per temp oral lag (see App endix A ). An episo de duration of δ = 1 hour ( i.e. , 12 time steps of 5 minutes) is selected for the OMSEV data. W e try several p ossibilities and the choice of δ = 1 hour represents a reasonable compromise b etw een temp oral interpretabilit y , appro ximate indep endence b etw een episo des, and a sufficient sample size for inference. This duration is consistent with the temp oral res- olution of the COMEPHORE data used for the advection estimation. A similar approac h is p erformed to determine the minim um spatial separation d min b et w een episodes. A v alue of d min = 1200 meters is chosen to balance indep endence and sample size for inference (see App endix A ). Regarding the COMEPHORE data, a similar analyses is employ ed to determine the episo de configuration parameters. The same quantile threshold q = 0 . 95 is adopted to define extreme episo des to ensure consistency with the OMSEV data analysis. Given the hourly time resolution of the COMEPHORE data, episo des are defined o v er δ = 24 hours p erio ds, starting from a conditioning exceedance, corresp onding to a daily episode duration. The minim um spatial separation b etw een episo des is set to d min = 5 km preserving a sufficien t n um b er of episo des for inference and ensuring approximate independence b et w een episo des. A summary of the episo de configuration parameters adopted for b oth the OMSEV and COMEPHORE datasets is pro vided in T able 3 . Under suc h configurations, a total of 384 extreme episo des are considered in the OMSEV dataset and 1130 in the COMEPHORE dataset. T able 3: Summary of the episo de configuration parameters adopted for the OMSEV and COMEPHORE datasets and the resulting num b er of selected episo des. Datatset OMSEV COMEPHORE Quan tile q 0 . 95 0 . 95 Episo de duration δ 60 min utes 24 hours Minim um spatial separation d min 1200 m 5 km Num b er of selected episo des |E | 384 1130 5.4 Estimation of adv ection F or eac h extreme rainfall episo de from the OMSEV dataset, a velocity v ector is required to c haracterize the propagation of rainfall storm. Due to the irregular spatial distribution of rain gauges and the limited spatial extent of the OMSEV netw ork, suc h a vector cannot in general b e reliably estimated from OMSEV data alone. T o o vercome this limitation, w e rely on the COMEPHORE dataset which cov ers a m uch larger spatial domain and pro vides rainfall estimates at higher spatial and temp oral resolution. COMEPHORE therefore offers a more comprehensiv e view of the displacement of rainfall systems. The duration of OMSEV extreme episo des is one hour, whic h matc hes the time resolu- tion of COMEPHORE data. Therefore, for eac h OMSEV episo de, the corresp onding hourly COMEPHORE time is iden tified and used to estimate the advection term. The empirical velo- 19 cit y v ector V emp COM is assigned based on COMEPHORE episo des o ccurring o ver a sp ecific time windo w extending tw o hours b efore and after the OMSEV episo de. Adv ection is estimated b y tracking the displacemen t of rainfall intensit y patterns ov er time. At each timestamp, the rainfall field is summarized b y its barycenter, defined as the center of mass of the rainfall field, where rainfall intensities act as weigh ts. The temp oral evolution of these barycenters provides successiv e displacement v ectors. By a veraging these vectors o ver time, we obtain an empirical adv ection vector for the episo de. An illustrative example of this pro cedure is shown in Figure 8 . (a) t 0 − 1 (b) t 0 (c) t 0 + 1 Figure 8: Extreme rainfall episo de from COMEPHORE on May 24, 2022. Rainfall intensit y is represented by p oint size. The empirical velocity vector (red arrow), estimated o v er the episo de, is sho wn at the barycen ter p osition at each time step. Only the barycenter lo cation c hanges b et w een panels. Es timated adv ection: (2 . 10 , − 1 . 96) km/h. Due to differences in spatial cov erage and ev en t detection, some OMSEV episo des do not corresp ond to an y COMEPHORE episo de. In such cases, the advection is estimated directly from the OMSEV net work, yielding an empirical v ector V emp OMSEV . Although this estimate is m uch less precise b ecause of the limited spatial extent of OMSEV netw ork, it still pro vides partial information on the lo cal displacement of precipitation and constitutes the only av ailable 20 option. The final velocity vector V final is defined as V final =    A ( V emp COM ) , if a corresp onding COMEPHORE episo de exists, A ( V emp OMSEV ) , otherwise. The function A ( · ) denotes the transformation describ ed in Section 4.4 . 5.5 V ariogram estimation The OMSEV v ariogram is estimated using a strategy based on COMEPHORE data. V ariogram parameters, including the advection transformation parameters η 1 and η 2 , are estimated from COMEPHORE episo des using the comp osite lik eliho o d describ ed in Section 4.3 . Empirical v elo cit y vectors are computed as detailed in Section 5.4 , and weigh ted least-squares estimates follo wing Buhl et al. 2019 are used to initialize the optimization. These estimates are then used b oth to initialize the OMSEV v ariogram optimization and to fix the advection parameters. Preliminary analyses rev eal that the range of the adv ection speeds has a non-negligible impact on the estimated v ariogram parameters. T o ensure comparabilit y betw een both datasets empirical adv ection terms (b efore transformation), w e perform the COMEPHORE optimization on a restricted range of adv ection speeds matc hing those observ ed o v er OMSEV data, i.e. ∥ V ∥ ≤ 5 . 6 km/h. This restriction remov es around 3% of the episo des and we keep a total of 1097 episo des for the COMEPHORE estimation. The resulting COMEPHORE parameter estimates are rep orted in T able 4 . W e observe that the temp oral comp onen t b β 2 is larger than the spatial one b β 1 , indicating that temporal v ariabilit y dominates. The small v alues of b α 1 and b α 2 suggest a high degree of irregularity in the spatial and temp oral dep endence structure of extremes. This is consistent with conv ectiv e extreme rainfall, esp ecially giv en the small spatial exp onent b α 1 . The estimated advection parameters b η 1 = 1 . 621 and b η 2 = 5 . 219 indicate a strong non-linear transformation of the empirical v elo cit y v ectors. In particular, the large v alue of b η 2 can lead to unrealistically high transformed advection sp eeds, although this affects only a small n umber of episo des. Indeed, there are 43 episo des (appro ximately 4% of COMEPHORE episodes) for which the transformed adv ection speed exceeds 150 km/h, and 13 for which it exceeds 1000 km/h, corresp onding to outliers. W e use 150 km/h as an upp er b ound for the transformed advection sp eed in the subsequent analysis. This restriction do es not remov e muc h episo des from the OMSEV dataset after transformation on their corresp onding adv ection sp eeds. Only 9 episo des are remov ed (appro ximately 2% ) b efore the v ariogram parameter estimation. The COMEPHORE parameter estimates are used as initial v alues for the OMSEV optimization, with the adv ection parameters η 1 and η 2 fixed. The resulting OMSEV parameter estimates, with mon th-leav e-one-out Jac kknife confidence in terv als, are presen ted in T able 5 . As for COMEPHORE, the temp oral comp onent remains larger than the spatial one and the estimated smo othness parameters α 1 and α 2 indicate strong irregularit y . The adequacy of the fitted mo del is assessed b y comparing empirical and fitted r -extremograms 21 for b oth datasets in Figure 9 . In b oth cases, the fitted r -extremograms closely match the em- pirical ones across a range of spatial and temp oral lags, indicating that the prop osed mo del effectiv ely captures the extremal dep endence structure of b oth COMEPHORE and OMSEV data. T able 4: COMEPHORE results for OMSEV-like adv ection class ( ∥ V ∥ ∈ [0 , 5 . 6] km/h) with |E | = 1097 . b β 1 b β 2 b α 1 b α 2 b η 1 b η 2 0.308 0.602 0.342 0.761 1.621 5.219 T able 5: Parameter estimates for OMSEV advection classes with month-lea v e-one-out Jackknife confidence in terv als. Parameters η 1 and η 2 are fixed to the COMEPHORE estimates. Number of episo des: |E | = 375 . km/h m/5 min b β 1 1 . 090 [0 . 780 , 1 . 401] 0 . 230 [0 . 087 , 0 . 560] b β 2 4 . 628 [3 . 383 , 5 . 873] 0 . 786 [0 . 324 , 1 . 767] b α 1 0 . 225 [0 . 133 , 0 . 317] 0 . 225 [0 . 133 , 0 . 317] b α 2 0 . 713 [0 . 483 , 0 . 944] 0 . 713 [0 . 483 , 0 . 944] (a) OMSEV data (b) COMEPHORE data Figure 9: Empirical against theoretical r -extremogram for each spatio-temp oral lags with es- timated parameters according to pairs num b er within the same spatio-temp oral lag. 6. Sto c hastic rainfall generation This section presen ts the final sto c hastic rainfall generator obtained by com bining the marginal and dep endence models dev elop ed in the previous sections and ev aluates its abilit y to reproduce the marginal and dep endence structure of extreme rainfall episo des observed in the OMSEV dataset. 22 6.1 V alidation of sim ulated rainfall episo des at OMSEV sites Rainfall sim ulations are p erformed by generating episo des of the r -Pareto pro cess given the o c- currence of an extreme ev ent at a conditioning location within the OMSEV rain gauge netw ork. F or eac h observed episo de, w e sim ulate a new episo de by conditioning on the same site using the same threshold and the same empirical velocity v ector, which will b e transformed using the estimated parameters b η 1 and b η 2 . The simulated r -P areto pro cesses are then transformed bac k to the original rainfall scale as describ ed in Section 3.4 . F or the v alidation we consider the marginal distributions fitted at each site separately and for the final sim ulations we use the common marginal distribution fitted across all sites, as justified in Section 5.1 . T o v alidate the marginal b ehavior and dep endence b eha vior of rainfall in tensities, w e con- sider the 384 extreme episodes iden tified in the OMSEV dataset (see Section 5.3 ) and w e generate 100 conditional simulations for eac h episo de. W e compare the distribution of strictly p ositiv e rainfall v alues at eac h site observed during extreme episo des with those pro duced b y the sim ulations. Figure 10 compares the densities of strictly p ositive rainfall v alues b et ween observ ations and sim ulations for t wo sites of the OMSEV netw ork. Without correction (top ro w), sim ulated distributions agree well with observ ations at mo derate and high in tensities but underestimate the frequency of low rainfall v alues (b elo w 1 . 5 mm). This difference is mainly due to the dis- cretization of observ ed rainfall, whereas simulations are contin uous by construction, making comparisons at lo w v alues unreliable. T o account for this effect, a discretization correction is applied to sim ulations: v alues strictly b etw een 0 and the gauge precision p are mapp ed to p , while larger v alues remain unchanged. As shown in the b ottom row, this adjustment sub- stan tially impro ves the agreement at lo w intensities while preserving the upp er-tail b ehavior. Ov erall, the simulated marginal distributions are consisten t with observ ations across the net- w ork (not shown), although sim ulations exhibit slightly more high-intensit y v alues, likely due to the absence of an upp er b ound and the presence of o ccasional extreme simulated outliers. These results indicate that the prop osed mo del adequately repro duces the marginal distribution of rainfall intensities during extreme episo des. 23 (a) IEM (without correction) (b) UM (without correction) (c) IEM (with correction) (d) UM (with correction) Figure 10: Density of strictly positive rainfall v alues for observ ed against simulated episodes for t wo selected sites of the OMSEV net work. T op row: without correction for the discretization effect. Bottom row: with a correction on low simulated v alues to account for discretization in observ ations. F or b etter visibility , only v alues b elow 10 mm are sho wn. T o ev aluate the abilit y of the model to capture the dep endence structure of rainfall extremes, w e compute the conditional triv ariate exceedance probabilities P ( X s 1 > u, X s 2 > u | X s > u ) for an y pairs of sites ( s 1 , s 2 ) ∈ S 2 , giv en an exceedance at a conditioning site s ∈ S . The threshold u is c hosen as the empirical 95% quan tile of the data. These probabilities are es- timated empirically from b oth observ ed and simulated data by coun ting join t exceedances at the three sites ( s 1 , s 2 , s ) and dividing by the n umber of exceedances at the conditioning site s . The results are sho wn in Figure 11 for t wo selected conditioning sites. W e observ e that the largest conditional probabilities are generally well reproduced by the sim ulations, whereas smaller probabilities tend to b e less accurately captured. As a complementary summary indicator of the dep endence structure, w e consider the dis- tribution of cumulativ e rainfall within episo des. This quantit y reflects b oth the num b er of exceedances o ccurring across the netw ork and their asso ciated intensities. Cumulativ e rainfall is defined as the sum of rainfall in tensities ov er all sites and time steps within an episo de, 24 pro viding a global measure of even t intensit y . The observed and sim ulated distributions are compared in Figure 12 . The sim ulated distribution exhibit a slightly larger prop ortion of lo w cum ulative rainfall v alues with a low er median than observ ed data while the bulk of the dis- tribution and the upp er tail is w ell captured. This b eha vior is exp ected from the conditional sim ulation framew ork: eac h episo de is generated conditionally on an exceedance at a single conditioning site, so that some simulated realizations hav e a single exceedance while rainfall at other sites remains close to zero. In con trast, observed episo des more frequently inv olve m ultiple exceedances during the episo de across the netw ork, whic h leads to larger cumulativ e totals ev en for mo derate even ts. Consequently , the sim ulated cum ulative distribution tends to underestimate the low er range of cumulativ e rainfall v alues. (a) P ( X s 1 > u, X s 2 > u | X IEM > u ) (b) P ( X s 1 > u, X s 2 > u | X UM > u ) Figure 11: Conditional triv ariate exceedance probabilities according to an excess at one site from observ ations v ersus sim ulations. 25 Figure 12: Densit y of cumulativ e rainfall o ver all sites within episo des for observed and s im u- lated data. F or b etter visibility , only v alues b elo w 400 mm are shown. The prop osed mo del is able to reproduce b oth the marginal distribution of rainfall intensities during extreme episo des and the main features of the extremal dep endence structure. The agreemen t is b etter for larger probabilities and higher rainfall accumulations. This is consistent with the primary ob jective of the mo del, whic h is to capture the b ehavior of extremes and their extremal dep endence. 6.2 Final sim ulations on a regular grid W e no w use the fitted mo del to sim ulate rainfall fields on a regular grid cov ering the OMSEV domain at fine spatial resolution. Sim ulations are p erformed on a regular grid with 100 m pixels, assuming spatially stationary marginal distributions with parameters given in T able 2 , and using the fitted dep endence mo del describ ed in T able 5 . An example of a sim ulated extreme rainfall episode is shown in Figure 13 . In this example, the conditional exceedance is imposed at a randomly selected grid pixel, and a randomly sampled velocity vector is applied. The resulting rainfall field exhibits a coherent displacement ov er time, consistent with the given advection direction. These simulations demonstrate the ability of the framework to generate con tinuous extreme rainfall episo de scenarios on a regular grid at fine resolution. The contin uous nature of the simulations allo ws for more realistic rainfall patterns that cannot b e directly captured from the discrete rain gauge netw ork alone. 26 Figure 13: Sim ulated rainfall episo de on a regular grid o ver the OMSEV area at the first four time steps. Rain gauge lo cations are shown as grey dots. The selected conditioning grid pixel is highlighted by a white square and the adv ection direction is indicated by a red arrow plotted at the barycen ter of the rainfall field. Here V E = ( − 1 , 2) km/h. A gif version is a v ailable here . 7. Discussion and p ersp ectiv es In this work, w e prop osed a spatio-temp oral sto c hastic framework for the simulation of ex- treme rainfall even ts at fine spatial and temp oral scales that explicitly accounts for advection within the dependence structure. The approac h com bines flexible marginal mo deling based on the EGPD with an adv ection-informed spatio-temp oral v ariogram for the dependence structure. The framew ork w as shown to accurately repro duce b oth the marginal b eha vior of extreme rain- fall and its spatio-temp oral dep endence structure. The resulting sto chastic rainfall generator can b e applied on a regular grid at fine spatial resolution, making it suitable for applications suc h as risk assessment and hydrological impact studies. Although developed and v alidated on the OMSEV dataset, the metho dology can be adapted to other high-resolution rainfall net works. Despite its flexibilit y , the prop osed mo deling framework has some limitations at such a fine spatial and temp oral scale and can b e improv ed. A ma jor source of uncertaint y is the estimation of advection velocities. In this study , velocity vectors were mainly inferred from COMEPHORE data due to the limited spatial cov erage and the irregularity of the OMSEV rain gauge netw ork. Even with COMEPHORE data, the spatial exten t of the study area (appro ximately 5 km around the OMSEV net work) strongly constrains the range of detectable adv ection v elo cities. As a result, sp eeds ab o ve appro ximately 10 km/h cannot b e iden tified and ma y b e in terpreted as null advection. In addiction, differences in ev en ts detection b et w een the 27 OMSEV and COMEPHORE datasets ma y lead to inconsistencies in the advection estimates, as extreme even ts detected by OMSEV data are not alw ays captured by COMEPHORE data within the corresp onding hour. A natural p ersp ective would b e to directly incorp orate real wind data or adv ection information in to the mo del. Suc h information w ould need to b e av ailable at sufficiently high spatial and temp oral resolution to remain consistent with the rainfall data considered here. Ov erall, the prop osed framework provides a coheren t and ph ysically in terpretable approach for the simulation of spatio-temp oral extreme rainfall fields, and constitutes a promising basis for future developmen ts and applications in this area. Co de a v ailabilit y The co de used in this pap er is av ailable on the GitHub rep ository chloesrcb/extreme-rainfall- generator . References D. Allard, P . Ailliot, V. Monbet, and P . Nav eau. Sto chastic we ather gener ators: An overview of we ather typ e mo dels , v olume 156. 2015. L. Benoit and G. Mariethoz. Generating syn thetic rainfall with geostatistical simulations. Wiley Inter disciplinary R eviews: W ater , 4(2):e1199, 2017. L. Benoit, D. Allard, and G. Mariethoz. Stochastic Rainfall Mo deling at Sub-kilometer Scale. 54(6):4108–4130, 2018. ISSN 1944-7973. P . Brunet, C. Bouvier, and L. Nepp el. Retour d’exp érience sur les crues des 6 et 7 o ctobre 2014 à Mon tp ellier-Grab els (Hérault, France) : caractéristiques h ydro-météorologiques et contexte historique de l’épiso de. Physio-Gé o. Gé o gr aphie physique et envir onnement , 12:43–59, 2018. ISSN 1958-573X. S. Buhl, R. A. Davis, C. Klüppelb erg, and C. Steinkohl. Semiparametric estimation for isotropic max-stable space-time pro cesses. Bernoul li , 25:2508–2537, 2019. J. Carreau and Y. Bengio. A h ybrid Pareto model for asymmetric fat-tailed data: the univ ariate case. Extr emes , 12(1), 2009. N. L. Carrer and C. Gaetan. Distributional regression mo dels for Extended Generalized Pareto distributions. , 2022. M. C. Casas, R. Ro dríguez, and Á. Redaño. Analysis of extreme rainfall in Barcelona using a microscale rain gauge netw ork. Mete or olo gic al Applic ations: A journal of for e c asting, pr actic al applic ations, tr aining te chniques and mo del ling , 17(1):117–123, 2010. 28 S. Coles, J. Heffernan, and J. T a wn. Dep endence measures for extreme v alue analyses. Extr emes , 2:339–365, 1999. S. Coles, J. Ba w a, L. T renner, and P . Dorazio. An intr o duction to statistic al mo deling of extr eme values , v olume 208. Springer, 2001. A. C. Davison, S. A. P adoan, and M. Ribatet. Statistical mo deling of spatial extremes. 2012. A. C. Davison, R. Huser, and E. Thibaud. Geostatistics of dep endent and asymptotically indep enden t extremes. Mathematic al Ge oscienc es , 45(5):511–529, 2013. R. De F ondeville and A. C. Da vison. High-dimensional p eaks-ov er-threshold inference. 105(3): 575–592, 2018. ISSN 0006-3444. C. Dombry , S. Engelk e, and M. Oesting. Exact sim ulation of max-stable pro cesses. 103(2): 303–317, 2016. ISSN 0006-3444, 1464-3510. C. Dombry , J. Legrand, and T. Opitz. Pareto pro cesses for threshold exceedances in spatial extremes, 2024. A. F erreira and L. De Haan. The generalized Pareto pro cess; with a view tow ards application and sim ulation. 2014. P . Finaud-Guyot, V. Guinot, P . Marc hand, L. Nepp el, C. Salles, and G. T oulemonde. Rainfall data collected by the HSM urban observ atory (OMSEV), 2023. C. F ouchier, J. La v abre, P . Ro y et, and H. F elix. Inondations de septembre 2002 dans le Sud de la France - Analyse hydrologique et h ydraulique au niv eau des barrages écrêteurs du Vidourle. 2004. A. Haruna, J. Blanc het, and A.-C. F a vre. Mo deling in tensity-duration-frequency curves for the whole range of non-zero precipitation: a comparison of mo dels. W ater R esour c es R ese ar ch , 59(6):e2022WR033362, 2023. N. Hemp elmann, C. Adeline, C. Pagé, and Y. Ro driguez-Dartois. Ville de Paris : Actualisation du diagnostic de vulnérabilités et de robustesses de Paris face aux changemen ts climatiques et à la raréfaction des ressources - Cahier 2 : Les évolutions climatiques à Paris. R amb ol l , 4, 2021. R. Huser and A. C. Da vison. Space–time mo delling of extreme ev ents. Journal of the R oyal Statistic al So ciety Series B: Statistic al Metho dolo gy , 76(2):439–461, 2014. D. Leb er. Comparison of simulation metho ds of Bro wn–Resnic k pro cesses. Master’s thesis, T echnisc he Universität München, Munic h, Germany , 2015. E. Leblois and J.-D. Creutin. Space-time sim ulation of intermitten t rainfall with prescrib ed adv ection field: Adaptation of the turning band metho d. 49(6):3375–3387, 2013. ISSN 1944-7973. 29 P . Lucas-Pic her, D. Argüeso, E. Brisson, Y. T ram blay , P . Berg, A. Lemonsu, S. Kotlarski, and C. Caillaud. Con vection-permitting mo deling with regional climate mo dels: Latest dev elopments and next steps. Wiley Inter disciplinary R eviews: Climate Change , 12(6):e731, 2021. P . Nav eau, R. Huser, P . Rib ereau, and A. Hannart. Mo deling join tly low, mo derate, and hea vy rainfall in tensities without a threshold selection. W ater R esour c es R ese ar ch , 2016. F. P alacios-Ro dríguez, G. T oulemonde, J. Carreau, and T. Opitz. Generalized Pareto pro- cesses for simulating space-time extreme ev en ts: an application to precipitation reanalyses. Sto chastic Envir onmental R ese ar ch and R isk Assessment , 34:2033–2052, 2020. I. Papastathopoulos and J. A. T a wn. Extended generalised Pareto mo dels for tail estimation. 143(1):131–143, 2013. ISSN 0378-3758. R. Robin. Calibration du réseau de pluviographes de l’observ atoire HSM « Eau dans la ville » , 2021. C. Scarrott and A. MacDonald. A review of extreme v alue threshold estimation and uncertaint y quan tification. REVST A T-Statistic al journal , 10(1):33–60, 2012. M. Sc hleiss, S. Chamoun, and A. Berne. Stochastic sim ulation of in termitten t rainfall using the concept of “dry drift”. 50(3):2329–2349, 2014. ISSN 1944-7973. P . T abary , P . Dupuy , G. L’Henaff, C. Gueguen, L. Moulin, and O. Lauran tin. A 10-y ear (1997—2006) reanalysis of Quantitativ e Precipitation Estimation o ver France: metho dology and first results. IAHS-AISH public ation , (351):255–260, 2012. D. S. Wilks and R. L. Wilb y . The w eather generation game: a review of sto chastic w eather mo dels. Pr o gr ess in physic al ge o gr aphy , 23(3):329–357, 1999. 30 A. A dditional figures for episo de configuration (a) Spatial joint exceedances (b) T emp oral joint exceedances Figure 14: Number of joint exceedances according to spatial distance or temp oral lag for the 0 . 95 quan tile. (a) By episo de duration δ . (b) By minimum inter-episode distance d min . Figure 15: Number of selected episo des according to episo de configuration parameters on the OMSEV dataset for a quan tile threshold q = 0 . 95 . Figure (a) shows the effect of the minim um spatial separation d min on the num b er of selected episo des for a fixed episo de duration δ = 60 min utes. Figure (b) shows the effect of the episode duration δ on the num b er of selected episo des for a fixed spatial separation d min = 1200 meters. In b oth figures, the final selected configuration corresp onds to the red dashed v ertical line. 31