Space-Time Trend Detection and Dependence Modeling in Extreme Event Approaches by Functional Peaks-Over-Thresholds: Application to Precipitation in Burkina Faso

In this paper, we propose a new method for estimating trends in extreme spatiotemporal processes using both information from marginal distributions and dependence structure. We combine two statistical approaches of an extreme value theory: the temporal and spatial nonstationarities are handled via a tail trend function in the marginal distributions. The spatial dependence structure is modeled by a latent spatial process using generalized ℓ -Pareto processes. This methodology for trend analysis of extreme events is applied to precipitation data from Burkina Faso. We show that a significant increasing trend for the 50 and 100year return levels in some parts of the country. We also show that extreme precipitation is spatially correlated with distance for a radius of approximately 200km.


Introduction
In the framework of climate change, the modeling and accurate prediction of the magnitude and extent of extreme events that occur in space and time of climate variables is a particularly di cult task, since it implies taking into account spatial and temporal nonstationarities. Nowadays, there is a general consensus in the scienti c community that climate change has accelerated in recent decades and that the climate will continue to change in the coming decades, mainly due to natural and anthropogenic changes (IPCC 2007(IPCC , 2018(IPCC , 2019.
is change is manifested in most regions of the world by a resurgence of heavy rainfall, heat waves, and pollution peaks with very signi cant economic and social consequences.
In sub-Saharan Africa, works on this topic have shown an increasing trend in the occurrence of extremes in meteorological parameters [1][2][3][4]. In Burkina Faso, a study conducted in the framework of the National Action Program for Adaptation to climate change (PANA, 2007) showed that precipitation is highly variable in space and time.
ere is a decreasing trend in cumulative and daily rainfall. According to recent publications from the National Meteorology Agency, precipitation has returned to humid periods since the late 1980s and decades of 1990 and 2000. e return of rains is more related to a high frequency of high intensity rainfall events than to an increase in the number of rainy days. For example, we can cite the oods of 1st September 2009 that a ected Ouagadougou and its outskirts with a record of 261.3 mm of rainfall not registered since 1919. e oods have a ected more than 150000 peoples, damaged several bridges, and ooded more than 9300 hectares of crops throughout the country. From a statistical point of view, this raises the question of trend detection in extreme events. e classical extreme value theory extended both to nonstationary and non-independent observations provides a rigorous mathematical framework to deal with this question of trend detection in extremes [5][6][7][8]. is issue was rst generally studied in extreme value literature by parametric pointwise approaches in which an extreme value model (GEV or GPD) is tted to the data at each site in turn, leaving the parameters of the marginal distributions to evolve with time or other significant covariates [9][10][11][12][13]. Although it is relatively simple to construct nonstationary models in the univariate framework, it is more difficult to account for spatial and temporal trends in these univariate models. e spatial and temporal dimension of extreme events was first developed for block maxima in a stationary framework and modeled by max-stable processes [14][15][16][17][18]. ese spatial models have been readapted to handle nonstationarities induced by global warming [19,20]. Although attractive, these models are expensive in computing time when tending toward large dimensions and are not adapted for modeling threshold exceedances. e threshold exceedance approach introduced by [21,22] was first extended to multivariate environments [23] before being generalized to functional data [24][25][26][27] to give birth to the family of generalized ℓ-Pareto processes. is family of processes are good candidates for modeling the spatial dependence structure of exceedance data. However, these approaches do not sufficiently take into account marginal nonstationarities and spatial dependence between margins. In this paper, we propose a new method to capture nonstationarities in marginal distributions, while taking into account the spatial dependence structure. To reach this goal, we combine two statistical approaches to the extreme value theory. First temporal and spatial nonstationarities are controlled via a tail trend function in marginal distributions [5][6][7]. Subsequently, the spatial dependence structure is taken into account by a hidden auxiliary spatial process using generalized ℓ-Pareto processes [25,26]. Finally, we use the model to simulate and predict future extreme events. e article is organized as follows: Section 2 details the methodology implemented and the methods used to estimate the parameters. In Section 3, we present our main results from the data analysis, and in particular, we calculate the nonstationary return levels from the developed approaches. Section 4 concludes the study.

Space-Time Trends
Modeling. Let X � X t (s), s ∈ S, t ∈ T} be a continuous nonstationary space-time stochastic process with sample paths in the family of continuous functions C(S × T) equipped with the uniform norm ‖ .‖ ∞ , where S × T ⊂ R d × R + and C + (S × T) its restriction to nonnegative functions deprived of the null function. In practice X is observed at each stations s 1 , . . . , s m and at given dates t � 1, . . . , n. Let F t,s be the continuous univariate marginal distribution with a common right endpoint x F and Z � Z(s), s ∈ S { } an unobserved latent spatial stochastic process with sample paths in C(S × T) satisfying the proportional tail condition such that where c θ : [0, 1] × S ⟶ (0, ∞) is assumed to be a continuous and positive function depending on a parameter vector θ ∈ Θ ⊂ R called a tail trend function or skedasis function [5,7]. e skedasis function describes the evolution of extreme events jointly in space and time. In the framework of the model (1), Mefleh et al. [8] shows that the empirical point measure converges in distribution in the space of point In this case, the times of exceedances for high threshold x and the value of exceedances are asymptotically independent with distributions respectively equal to the trend density function c θ (u, s), s ∈ S and the Pareto distribution of tail index c.
Moreover, we assume that the continuous marginal distributions F Z of the latent process is in the maximum domain of attraction condition for some constant c ∈ R and appropriate real normalization constants a n > 0, b n ∈ R.
anks to the equation (1) and the convergence of Z exceedances to a GPD distribution we deduce a sample of Z t (s), s ∈ S as mentioned in [5,7] in the following manner: where c, a n , b n , and c θ are parameters to be estimated. Parameter inference techniques will be discussed in Section 2.4.1.

Nonstationary Return Period and Return
Level. e concept of return period and level becomes very ambiguous when we leave the stationary context to the nonstationary framework. For example, in the stationary framework, an extreme event with a 100-year return period is likely to occur on average once every 100 years, but with an annual exceedance probability of 1% in a given year. In this case, it is assumed that the underlying probability distribution remains unchanged over time. More formally, for a stationary random variable Z with distribution function F Z , the m-year return period is expressed as m � 1/F Z (Z > z m ), where z m is the return level associated with the return period. However, in the nonstationary context, the underlying probability distribution is no longer constant but is supposed to evolve with time. In this paper, we have chosen to follow return period based approaches, i.e., expected number of exceedances (ENE) [12] and expected waiting time (EWT) [9,11]. In the ENE approach, the number of times the variable X t (s), s ∈ S exceeds the return level value x m in m years is defined by N m � n x m t�1 1 X t (s)>x m ,s∈S { } under nonstationary context. e return level x m can be defined as the value for which the expected number of events exceeding x m in m years equals to one, i.e., the return level x m is the solution of the following equation: where n x is the number of days in the year and θ t (s) the vector of time-dependent marginal parameters or other covariates. Parey et al. [12] uses the ENE method in a pointwise POT model where the parameters of the distribution of exceedances and the intensity of extreme event occurrences are described as polynomial functions of time. e EWT method was first proposed by Olsen et al. [11], and then derived by Salas and Obeysekera [13] using a geometric distribution with time-varying parameters. Under nonstationary conditions, the distribution describing waiting time Y before the first occurrence of an event exceeding the return level x m is as follows: where variable Y is the day of the first occurrence of an event exceeding the quantile and p t,s � 1 − F t,s (x m |θ t (s)) is daily exceedance probability varying with time step t. y max is the time when the daily exceedance probability p t,s is equal to 1 for an increasing-trend series or is equal to 0 for a decreasing-trend series. Reused and simplified by [9], the EWT approach defines the m-year return level x m , as the value for which the expected waiting time until an exceedance of this level is m years, i.e, x m is the solution of the equation: Using the relationship (1), P(X t (s) > x m ) can be rewritten asc θ (t/n, s)P(Z(s) > u)P(Z(s) > x m |Z(s) > u) for u < x m and we obtain the following results. Proposition 1. Let X t (s), s ∈ S, t ∈ T be a nonstationary stochastic process defined on a region S ⊂ R d and Z(s), s ∈ S { } a latent spatial process satisfying equation (1). Given a return period m and threshold u < x m , the return level x m for all s ∈ S is a solution of the following two equations: (i) Return period as expected number of events (ii) Return period as expected waiting time where F Z(s)|ℓ(Z)>u is a survival of generalized Pareto distribution of ℓ-exceedances at position s ∈ S; c θ is the tail trend function and n x is the number of days in the year.
e nonstationary return levels obtained from (6) and (7) are evaluated by numerical algorithms taking into account information from the extrapolation of the trend function c θ . To derive the nonstationary return level at points where we have not observations, we use spatial marginal model (15) described in Section 2.4.1 to extrapolate the parameters at these points. An alternative approach is to calculate the nonstationary return levels from the following result.

Proposition 2.
Let X t (s), s ∈ S, t ∈ T be a nonstationary spatio-temporal stochastic process defined on a region S ⊂ R d and Z(s), s ∈ S { } a latent spatial process. e nonstationary return level x m (s), s ∈ S of the nonstationary process X is deduced from the return level z m of the latent spatial process Z such that where a n , b n , c θ , and c are the respective estimators of a n , b n , c θ , and c. t m � 1 + n x m/n with n x the number of days in the year and n the size of the sample observed.
is result is a consequence of the (2). In practice, the return level z m of the latent spatial process is first computed at grid points where we have not observations using (15). en, the nonstationary return level x m is deduced using the Proposition 2.

Spatial POT Modeling.
After removing the nonstationarity, thanks to the (2), the modeling is then focused on the evaluation of the extreme spatial dependence structure in Z using functional POT [24][25][26]. In the multivariate and spatial framework a threshold exceedances for a random function Z � Z(s), s ∈ S { } is defined by [26] to be an event of the form ℓ(Z) > u { } for some u ≥ 0, where ℓ: C(S) ⟶ R + is a continuous and homogeneous nonnegative risk function, i.e., there exists α > 0 such that ℓ(λy) � λ α ℓ(y) when y ∈ C + (S) and λ > 0. e risk function ℓ determines the type of extreme events of interest. For example, such a function can be the maximum, minimum, average, or value at a specific point s 0 ∈ S. Under minimal assumptions on the risk function, the conditional distribution of ℓ-exceedance for some threshold u ≥ 0 of the process (Z − b n )/a n ) can be approximated by a generalized ℓ-Pareto process, for n large enough [25].
where W ℓ is a nondegenerate stochastic process over S and belongs to the family of generalized ℓ-Pareto processes International Journal of Mathematics and Mathematical Sciences with tail index c (see appendix B for more details). Specifically, W ℓ is a stochastic process taking values in z ∈ C + (S): } ≥ 0} and defined by the following: where a > 0 and b are continuous functions on S, respectively, scale and location functions and Y ℓ is a stochastic process whose probability measure is completely determined by the limit measure Λ. More details on generalized ℓ-Pareto processes can be found in [24,25].
For the modeling of spatial dependence structure of the latent process, we use Y ℓ whose margins are in the Frechet domain of attraction with tail index c � 1, as the process of reference. A pseudo-polar decomposition of Y ℓ in (11) leads to the following formulation [24,26].
where R is a unit Pareto random variable of index c R representing the intensity of process and Q is stochastic process denoted the angular component with state space S and taking values in S � y ∈ C + (S): ‖y‖ 1 � 1 whose probability measure is characterized by limit measure Λ.

Marginal Model of the Spatial Latent Process.
Marginal parameters a n , b n , and c of latent process can be estimated into fixing b n at a high quantile of ℓ(Z), i.e., b n � q 1− α ℓ(Z) { }. In general, a parametric model may be necessary for a n and b n , as in Engelke et al. [28]. However, some forms of parameterization can lead to problems of parameter identifiability and inference. For these reasons and simplicity purpose, we consider in this work that a n (s j ) � a j > 0 and b n (s j ) � b j ∈ R for any j � 1, . . . , m. e threshold stability of the generalized Pareto distributions does not allow us to identify the b n function without additional assumptions, so we set where u q′ Z(s) { } is an empirical quantile of the order q ′ of the ℓ-exceedances at each location s ∈ S, where q ′ is chosen such that ℓ(b n ) � b n in order to impose the identifiability of the parameters.
us, the tail index c ∈ R and the scale parameters a n (s) > 0 are estimated by maximizing the independent log-likelihood; that is, ℓ in de p c, a n s j � n t�1 m j�1 1 x t s j ≥b n s j · log 1 a n s j 1 + c x t s j − b n s j a n s j e assumptions of parameter identifiability (see appendix B) implies that ℓ(a n ) � a n . e trend function parameter c θ is estimated using the maximum likelihood method under the assumption that the dates of exceeding the marginal thresholds b n are asymptotically independent and identically distributed, from the density function c θ (t/n, s) [8]. Several models are eligible to model the function c θ . However, in this study, we opt for parametric models because of their flexibility in trend detection and in order to make extrapolations of the trend beyond the observed data. To facilitate inference by the maximum likelihood method, we are interested in monotonic log-linear and simple linear trend models of c θ such that where θ(s), s ∈ S is a characteristic parameter of the trend function. θ(s) > 0 suggests that extreme are getting more frequent, while θ(s) < 0 indicates that extremes are getting less frequent.

Spatial Extremal
Dependence. e spatial dependence structure is taken into account in the marginal parameters of spatial process Z by letting these parameters evolve as a function of space covariates. More precisely, we use a spatial model of the marginal parameters as functions of the significant space covariates using a generalized linear regression model [9,12,13]. a n (s) � g a n (v(s)) > 0, where v is a vector of covariates which can contain for example the space coordinates, or any other spatialized data explaining significantly the weather variable studied. Furthermore, after removing the nonstationarity, modeling is focused on the evaluation of the remaining spatial dependence structure in Z. us, we approach the limit distribution of ℓ-exceedances of spatial latent process by a generalized ℓ-Pareto process (9). To estimate the parameters of spatial dependence structure, we have chosen to model its angular component Q (11) by a log-Gaussian process with stationary increments and within the framework of the Brown-Resnick model. In this case, we have access to the limit measure Λ. In order to better capture this dependence structure, we use a flexible parametric semi-variogram belonging to the class of power semi-variograms. e parameter of the dependence structure is estimated using the gradient score method or the censored likelihood method [24,29].

Data Set Analysis.
is study uses time series of daily precipitation measurements from 1957 to 2016 provided by ten synoptic stations extracted from the Burkina Faso climatological database. ese stations have been selected to ensure good spatial uniformity and representativeness of different climatic regimes and data quality. Figure 1 gives the spatial distribution of selected stations in our study area. In order to limit the problems related to seasonal rainfall cycles on each station, we worked from the sub-series corresponding to rainy days. e period from May to October was chosen because it is during this period that the most rainfall is recorded in the country. us, a series of 60 years (1957-2016) by 184 days (May to October) is extracted to constitute the time series of daily rainfall. On these time series, we apply a run declustering procedure with a daily step (r � 1 day) to identify the groups of approximately independent extreme observations within the sample in order to avoid short-term dependencies in the time series.

Marginal Characteristics.
e first line of Figure 2 shows the maps of the parameters a n and b n estimated on our selected stations. Map (a) gives us the variability parameter of the precipitation distribution and shows a very high variability of extreme precipitation. While map (b) gives us the estimated level or threshold for each station using equation (12) and we see that the level is higher in the Sudanian climate and lower in the Sahelian climate.In practice, we deduce from the generalized ℓ-Pareto model fitted to the data, the marginal tail behavior of the survival distribution F Z(s),c(s),σ(s) for any s ∈ S from the (16) given a sufficiently large u q > 0 threshold.
with σ(s) � a n (s) + c(u q (s) − b n (s)), s ∈ S, where a n , b n , and c are the marginal parameters estimators described in Section 2.4.1. We can use it to check the quality of the marginal adjustment of the stationary process Z. Figure 3 shows the qq-plots of the local tail distribution due to two stations by climatic zone. Columns 1, 2, and 3 of Figure 3 represent the respective adjustments for the Sudanian, Sudano-Sahelian, and Sahelian climate. Globally, the fit of the marginal models seems convincing, as most of the observations remain within the confidence intervals. e second row of Figure 2 shows the maps of parameter θ estimated for the log-linear and simple linear trend functions. We can observe a strong spatial variability in the occurrence of extreme precipitation. us, we observe a decreasing trend in the frequency of extreme rainfall in Ouaga, Dori, and Bobo. In contrast, in Boromo, Ouahigouya, and Pô, extreme precipitation events will become more and more frequent. Figure 4 gives us the details of the adjustment of the tail trend function by a log-linear model.

Estimated Spatial Dependence Model.
We use a spatial marginal model (15) to estimate the marginal parameters at the different grid points. is is to prepare the ground for the calculation of nonstationary return levels at locations where we have no observations. a n (s) � a 0 + a 1 Long × Lat(s), An alternative to the spatial model of marginal parameters is to consider small subsets N s 0 ⊂ S, where N s 0 is a small neighborhood around s 0 . In what follows, these regional neighborhoods will be determined by a small number D 0 of nearest stations from the site s 0 , thus we will write N s 0 � N s 0 ,D 0 . In principle, the choice of D 0 should be such that the spatial dependence structure and marginal parameters is approximately stationary within each selected neighborhood N s 0 ,D 0 around s 0 ∈ S. Obviously, the choice of neighborhood is important; the assumed stationary marginal parameters could be a poor approximation for large neighborhoods (i.e., for large D 0 ), while the simulation of the process could be cumbersome for small neighborhoods characterized by a small D 0 number of stations. We obtain the relatively homogeneous, nonoverlapping subregions using the k-means clustering method centered on the reference stations [30,31]. is method is extensively used because it is computationally simple and produces accurate results, compared to other more complex clustering methods. e longitude and latitude of the grid points were used as input variables in the k-means clustering algorithm to form the ten clusters centered on the reference stations.  In addition, to better take into account the dependence structure of precipitation data and capture the possible isotropic, we choose a flexible parametric semi-variogram ] belonging to the power class models such that ](h) � (‖h‖/τ) κ , with τ > 0, h ∈ R 2 , κ ∈ ]0, 2]. We check the adequacy of the fit of the model to the data using the extremogram and the variogram. Figures 5(a) and 5(b) show, respectively, the good quality of the fit of the dependence measures such as the semi-variogram and extremogram. In Figure 5(a), the points in gray represent the calculated empirical extremogram and the blue curve is the fitted empirical extremogram, while the red curve represents our fitted dependence model. e red curve in Figure 5(b) is the variogram model fitted to the data as a function of distance.  It is noted that extreme precipitation is spatially correlated for distances of the order of 200 km. In other words, an extreme precipitation event is likely to affect an area within a 200 km radius. is spatial dependence decreases gradually as the distance increases before stabilizing for example for a value of π(h) � 0.0625. We find the very localized nature of extreme precipitation for long distances.
is is generally a property present in rainfall data [32].

Nonstationary Return Level
Results. After estimating the marginal parameters and those of the dependence structure, we now evaluate the nonstationary return levels. e return levels x 50 and x 100 are first estimated (see Figure 6) using the (6) and (7). en, they are spatially interpolated at points where we have no observations using the spatially extrapolated trend function and dependence structure estimated from the observed data. For a given return period m, we first compute the return levels z m on each point of grid from the latent spatial process Z using the estimated spatial model. e nonstationary return levels x m is deducted using the (8) of Proposition 2. We repeat this operation for different values of the return period m and we deduce on each point of the grid, the associated return level (Figure 7). We can then produce maps of the return levels obtained from the dependence structure and the extrapolation of the trend function for a future period m. us, the extreme precipitation, likely to be observed on average at least once every 50 years (resp. 100 years), will be particularly intense in the Sudanian and Sudano-Sahelian zone and less intense in the Sahelian zone, with a potentially quite strong spatial dependence within a radius of 200 km. e southwest and eastern regions of the country will be most affected by extreme precipitations. e results of the return level for a return period m � 50 and 100 years are displayed in Figure 6.

Discussion and Conclusion
In summary, modeling the spatiotemporal trends of extreme precipitation in the sub-Saharan Africa is a challenge and a necessity. Climate models do not often agree on this issue. Advanced statistical methods remain a fairly powerful alternative to obtain pertinent information on possible spatiotemporal changes in the underlying precipitation distribution. A frequently used approach to dealing with nonstationarity in trend analysis of extreme precipitation is to include covariates, such as space and time, in the parameters of the underlying distribution [9,10,12,13]. ese approaches are built under certain important assumptions that the underlying process that generates the data is locally stationary and that the observations generated can be considered as independent approximations and generally do not take into account spatial dependence.
In this study, we proposed a new flexible methodology for trend detection in extremes of spatiotemporal processes.
is is capable of capturing marginal nonstationarities and the dependence structure between margins. Marginal nonstationarity is taken into account by a trend function called the skedasis function which models the frequency of extreme events jointly in space and time. While the spatial dependence structure is governed by a latent spatial process and modeled using generalized ℓ-Pareto processes. We calculated the nonstationary return levels of precipitation in Burkina Faso and exhibit some regions in which there may be a significant increase or decrease of extreme precipitation. We showed that these extreme rainfall events are spatially correlated around a radius of 200 km. is spatial dependence decreases progressively as the distance increases. In sum, we set up a nonstationary stochastic generator of extreme rainfall in Burkina Faso.
Previous studies conducted in the sub-Saharan African region have shown trends in extreme precipitation similar to those we obtained using Burkina Faso data set. According to the latest publications of the National Meteorological Agency of Burkina Faso, rainfall has started to move back toward a wet period since the end of the 1980s and the decade 1990-2000. is return of rain is more related to an increase in the frequency of heavy rainfall events than to an increase in the number of rainy days. Our results are also in agreement with the results found by [2][3][4] in the region which indicate a similar intensification of rainfall events. Based on a point GEV model with time-dependent parameters, these authors show that intense precipitation events will become more frequently occurring in the region. Nkrumah et al. [33] and Mouhamed et al. [34] observed that the average annual precipitation has a decreasing trend in the same region due to less frequent but more intense precipitation during the last decade.
Our method generates extreme events consistent enough with the data already observed. e return levels computed can be improved by introducing also a nonstationarity in the structure of spatial dependence.
Furthermore, it follows that which gives with σ(s) � a n (s) + c(s)(u(s) − b n (s)). (ii) e result (7) is shown in a similar way.

B. Overview on Generalized ℓ-Pareto Process
Let S be a compact subset of R d denoting the spatial domain of the study region. We note by C(S) the set of continuous real functions on S equipped with the uniform norm ‖.‖ ∞ ; C + (S) its restriction to nonnegative functions deprived of the null function and M C + (S) the class of measures associated with C + (S).
Let Z(s), s ∈ S { } be a latent spatial stochastic process indexed by s ∈ S ⊂ R d with sample paths in C(S) of continuous marginal distribution F Z and common right endpoint z F . As mentioned in [25,27], we assume that the spatial process Z is a general functional regular variation (GRV), i.e., that there exists suitable sequences of continuous functions a n : S ⟶ R + , b n : S ⟶ R and c ∈ R such that and noted by Z ∈ GRV(c, a n , b n , Λ), where Λ is a nonzero measure in M C + (S) and homogeneous of order − 1, with Λ(tA) � t − 1 Λ(A) for any positive real t > 0 and Borel set A ⊂ C + (S). e shape parameter c ∈ R is also called the tail index, which determines the strength of the tail and its support. According to the sign of c, we end up with three types of distributions of extreme values known as Gumbel (c ⟶ 0), Frechet (c > 0), and Weibull (c < 0).
In the multivariate and spatial framework, a threshold exceedances for a random function Z � Z(s), s ∈ S { } is defined by Dombry and Ribatet [26] to be an event of the form ℓ(Z) > u { } for some u ≥ 0, where ℓ: C(S) ⟶ R + is a continuous and homogeneous nonnegative risk function, i.e., there exists α > 0 such that ℓ(λy) � λ α ℓ(y) when y ∈ C + (S) and λ > 0. e risk function ℓ determines the type of extreme events of interest. For example, such a function can be the maximum, minimum, average, or value at a specific point s 0 ∈ S. Moreover, as in [25,28], we assume that there exists a continuous and positive real function A Z such that the sequence a n and the risk function ℓ satisfy the following asymptotic decomposition: lim n⟶∞ sup s∈S a n (s) ℓ a n − A Z (s) � 0, i.e a n (s) ≈ ℓ a n A Z (s), n ⟶ ∞.

(B.2)
Furthermore, the marginal distributions of Z are supposed to belong to a location-scale family, thus ensuring a constant c ∈ R shape parameter over any S, i.e., that there exist a H distribution function such that where A Z : S ⟶ (0; ∞] verifies asymptotic decomposition (B.2) and B Z : S ⟶ R are continuous functions. In particular H must belong to the domain of attraction of a generalized extreme value (GEV) distribution of tail index c ∈ R, i.e., F Z ∈ D(G c ). In other words, there are appropriate normalization real sequences a n > 0, b n ∈ R such that lim n⟶∞ H n (a n z + b n ) � G c (z). e condition of the maxdomain of attraction is equivalent to lim n⟶+∞ n 1 − H a n z + b n � − log G c (z), z > 0. (B.4) Assuming these conditions and general functional regular variation, the appropriate sequence of continuous functions a n (s) and b n (s), satisfy the following: a n (s) � a n A Z (s), b n (s) � B Z (s) + A Z (s)b n , s ∈ S.

(B.5)
Functions A Z (s) and B Z (s) as well as the spatial dependence structure of Z are assumed to belong to a family of parametric functions A Z,θ A : θ A ∈ Θ A ; B Z,θ B : θ B ∈ Θ B and Λ θ Λ : θ Λ ∈ Θ Λ . Under minimal assumptions on the risk function and assumptions (18,19) of Z, the conditional distribution of ℓ-exceedance for some threshold u ≥ 0 of the process (Z − b n )/ℓ(a n ) can be approximated by a generalized ℓ-Pareto process, for n large enough [25]. P ⌊ Z − b n ℓ a n ⌋ ∈ A|ℓ Z − b n ℓ a n ≥ u ⟶ P W ℓ ∈ A , n ⟶ ∞, where ⌊z⌋ � max(z, A Z c − 1 ) if c > 0 and ⌊z⌋ � z if c ≤ 0. W ℓ is a nondegenerate stochastic process over S and belongs to the family of generalized ℓ-Pareto processes with tail index c, zero location, scaling function A Z and limit measure Λ. Specifically, a generalized ℓ-Pareto process W ℓ associated to the limit measure Λ and tail index c ∈ R is a stochastic process taking values in z ∈ C + (S): ℓ (z − b)/ℓ(a) { } ≥ 0 and defined by the following: where a � ℓ(a)A Z > 0 and b are continuous functions on S, respectively, scale and location functions. Y ℓ is a stochastic process whose probability measure is completely determined by the limit measure Λ.