Abstract
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 100 year return levels in some parts of the country. We also show that extreme precipitation is spatially correlated with distance for a radius of approximately 200 km.
1. 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 difficult task, since it implies taking into account spatial and temporal nonstationarities. Nowadays, there is a general consensus in the scientific 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, 2018, 2019). This change is manifested in most regions of the world by a resurgence of heavy rainfall, heat waves, and pollution peaks with very significant 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–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. There 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. The 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 floods of 1st September 2009 that affected Ouagadougou and its outskirts with a record of 261.3 mm of rainfall not registered since 1919. The floods have affected more than 150000 peoples, damaged several bridges, and flooded 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.
The 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–8]. This issue was first generally studied in extreme value literature by parametric pointwise approaches in which an extreme value model (GEV or GPD) is fitted to the data at each site in turn, leaving the parameters of the marginal distributions to evolve with time or other significant covariates [9–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. The spatial and temporal dimension of extreme events was first developed for block maxima in a stationary framework and modeled by max-stable processes [14–18]. These 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.
The threshold exceedance approach introduced by [21, 22] was first extended to multivariate environments [23] before being generalized to functional data [24–27] to give birth to the family of generalized -Pareto processes. This 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–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.
The 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.
2. Methodology
2.1. Space-Time Trends Modeling
Let be a continuous nonstationary space-time stochastic process with sample paths in the family of continuous functions equipped with the uniform norm , where and its restriction to nonnegative functions deprived of the null function. In practice is observed at each stations and at given dates . Let be the continuous univariate marginal distribution with a common right endpoint and an unobserved latent spatial stochastic process with sample paths in satisfying the proportional tail condition such thatwhere is assumed to be a continuous and positive function depending on a parameter vector called a tail trend function or skedasis function [5, 7]. The 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 measure to a Poisson point process with intensity measure on . In this case, the times of exceedances for high threshold and the value of exceedances are asymptotically independent with distributions respectively equal to the trend density function and the Pareto distribution of tail index .
Moreover, we assume that the continuous marginal distributions of the latent process is in the maximum domain of attraction condition for some constant and appropriate real normalization constants , . Thanks to the equation (1) and the convergence of exceedances to a GPD distribution we deduce a sample of as mentioned in [5, 7] in the following manner:where , and are parameters to be estimated. Parameter inference techniques will be discussed in Section 2.4.1.
2.2. Nonstationary Return Period and Return Level
The 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 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 with distribution function , the -year return period is expressed as , where 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 exceeds the return level value in years is defined by under nonstationary context. The return level can be defined as the value for which the expected number of events exceeding in years equals to one, i.e., the return level is the solution of the following equation:where is the number of days in the year and 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.
The 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 before the first occurrence of an event exceeding the return level is as follows:where variable is the day of the first occurrence of an event exceeding the quantile and is daily exceedance probability varying with time step t. is the time when the daily exceedance probability 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 -year return level , as the value for which the expected waiting time until an exceedance of this level is years, i.e, is the solution of the equation:
Using the relationship (1), can be rewritten as for and we obtain the following results.
Proposition 1. Letbe a nonstationary stochastic process defined on a regionanda latent spatial process satisfying equation (1). Given a return periodand threshold, the return levelfor allis a solution of the following two equations:(i)Return period as expected number of events(ii)Return period as expected waiting time where is a survival of generalized Pareto distribution of -exceedances at position ; is the probability of exceedances, is the tail trend function and is the number of days in the year.
The nonstationary return levels obtained from (6) and (7) are evaluated by numerical algorithms taking into account information from the extrapolation of the trend function . 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. Letbe a nonstationary spatio-temporal stochastic process defined on a regionanda latent spatial process. The nonstationary return levelof the nonstationary processis deduced from the return levelof the latent spatial processsuch thatwhere,, andare the respective estimators of,,, and.withthe number of days in the year and n the size of the sample observed.
This result is a consequence of the (2). In practice, the return level of the latent spatial process is first computed at grid points where we have not observations using (15). Then, the nonstationary return level is deduced using the Proposition 2.
2.3. 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 using functional POT [24–26]. In the multivariate and spatial framework a threshold exceedances for a random function is defined by [26] to be an event of the form for some , where is a continuous and homogeneous nonnegative risk function, i.e., there exists such that when and . The 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 . Under minimal assumptions on the risk function, the conditional distribution of -exceedance for some threshold of the process can be approximated by a generalized -Pareto process, for large enough [25].where is a nondegenerate stochastic process over S and belongs to the family of generalized -Pareto processes with tail index (see appendix B for more details). Specifically, is a stochastic process taking values in and defined by the following:where and are continuous functions on , respectively, scale and location functions and 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 whose margins are in the Frechet domain of attraction with tail index , as the process of reference. A pseudo-polar decomposition of in (11) leads to the following formulation [24, 26].where is a unit Pareto random variable of index representing the intensity of process and is stochastic process denoted the angular component with state space and taking values in whose probability measure is characterized by limit measure .
2.4. Statistical Inference
2.4.1. Marginal Model of the Spatial Latent Process
Marginal parameters , , and of latent process can be estimated into fixing at a high quantile of , i.e., . In general, a parametric model may be necessary for and , 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 and for any . The threshold stability of the generalized Pareto distributions does not allow us to identify the function without additional assumptions, so we setwhere is an empirical quantile of the order of the -exceedances at each location , where is chosen such that in order to impose the identifiability of the parameters. Thus, the tail index and the scale parameters are estimated by maximizing the independent log-likelihood; that is,
The assumptions of parameter identifiability (see appendix B) implies that .
The trend function parameter is estimated using the maximum likelihood method under the assumption that the dates of exceeding the marginal thresholds are asymptotically independent and identically distributed, from the density function [8]. Several models are eligible to model the function . 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 such thatwhere is a characteristic parameter of the trend function. suggests that extreme are getting more frequent, while indicates that extremes are getting less frequent.
2.4.2. Spatial Extremal Dependence
The spatial dependence structure is taken into account in the marginal parameters of spatial process 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].where 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 . Thus, 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 (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. The parameter of the dependence structure is estimated using the gradient score method or the censored likelihood method [24, 29].
3. Application to Extreme Precipitation in Burkina Faso
3.1. Data Set Analysis
This study uses time series of daily precipitation measurements from 1957 to 2016 provided by ten synoptic stations extracted from the Burkina Faso climatological database. These 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. The period from May to October was chosen because it is during this period that the most rainfall is recorded in the country. Thus, 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.
3.2. Marginal Characteristics
The first line of Figure 2 shows the maps of the parameters and 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 for any from the (16) given a sufficiently large threshold.with , , where , , and 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 . 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.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
The 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. Thus, 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.
3.3. Estimated Spatial Dependence Model
We use a spatial marginal model (15) to estimate the marginal parameters at the different grid points. This is to prepare the ground for the calculation of nonstationary return levels at locations where we have no observations.
An alternative to the spatial model of marginal parameters is to consider small subsets , where is a small neighborhood around . In what follows, these regional neighborhoods will be determined by a small number of nearest stations from the site , thus we will write . In principle, the choice of should be such that the spatial dependence structure and marginal parameters is approximately stationary within each selected neighborhood around . Obviously, the choice of neighborhood is important; the assumed stationary marginal parameters could be a poor approximation for large neighborhoods (i.e., for large ), while the simulation of the process could be cumbersome for small neighborhoods characterized by a small number of stations. We obtain the relatively homogeneous, nonoverlapping subregions using the k-means clustering method centered on the reference stations [30, 31]. This method is extensively used because it is computationally simple and produces accurate results, compared to other more complex clustering methods. The 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 , with , , . 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. The red curve in Figure 5(b) is the variogram model fitted to the data as a function of distance.
(a)
(b)
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. This spatial dependence decreases gradually as the distance increases before stabilizing for example for a value of . We find the very localized nature of extreme precipitation for long distances. This is generally a property present in rainfall data [32].
3.4. Nonstationary Return Level Results
After estimating the marginal parameters and those of the dependence structure, we now evaluate the nonstationary return levels. The return levels and are first estimated (see Figure 6) using the (6) and (7). Then, 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.
(a)
(b)
For a given return period , we first compute the return levels on each point of grid from the latent spatial process using the estimated spatial model. The nonstationary return levels is deducted using the (8) of Proposition 2. We repeat this operation for different values of the return period 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 . Thus, 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. The southwest and eastern regions of the country will be most affected by extreme precipitations. The results of the return level for a return period and 100 years are displayed in Figure 6.
4. 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]. These 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. This 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. This 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. This 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–4] in the region which indicate a similar intensification of rainfall events. Based on a point 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. The return levels computed can be improved by introducing also a nonstationarity in the structure of spatial dependence.
Appendix
A. Proof of Proposition 1
(i)Proof. Let be a random variable describing the number of exceedances for a period of time ; it comes that And according to the (1), we have as follows: Moreover, given a sufficiently large threshold with , Furthermore, it follows that which gives with .(ii)The result (7) is shown in a similar way.B. Overview on Generalized -Pareto Process
Let be a compact subset of denoting the spatial domain of the study region. We note by the set of continuous real functions on equipped with the uniform norm ; its restriction to nonnegative functions deprived of the null function and the class of measures associated with .
Let be a latent spatial stochastic process indexed by with sample paths in of continuous marginal distribution and common right endpoint . As mentioned in [25, 27], we assume that the spatial process is a general functional regular variation (GRV), i.e., that there exists suitable sequences of continuous functions , and such thatand noted by , where is a nonzero measure in and homogeneous of order , with for any positive real and Borel set . The shape parameter is also called the tail index, which determines the strength of the tail and its support. According to the sign of , we end up with three types of distributions of extreme values known as Gumbel , Frechet , and Weibull .
In the multivariate and spatial framework, a threshold exceedances for a random function is defined by Dombry and Ribatet [26] to be an event of the form for some , where is a continuous and homogeneous nonnegative risk function, i.e., there exists such that when and . The 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 . Moreover, as in [25, 28], we assume that there exists a continuous and positive real function such that the sequence and the risk function satisfy the following asymptotic decomposition:
Furthermore, the marginal distributions of are supposed to belong to a location-scale family, thus ensuring a constant shape parameter over any , i.e., that there exist a distribution function such thatwhere verifies asymptotic decomposition (B.2) and are continuous functions. In particular must belong to the domain of attraction of a generalized extreme value (GEV) distribution of tail index , i.e., . In other words, there are appropriate normalization real sequences , such that . The condition of the max-domain of attraction is equivalent to
Assuming these conditions and general functional regular variation, the appropriate sequence of continuous functions and , satisfy the following:
Functions and as well as the spatial dependence structure of are assumed to belong to a family of parametric functions ; and . Under minimal assumptions on the risk function and assumptions (18, 19) of , the conditional distribution of -exceedance for some threshold of the process can be approximated by a generalized -Pareto process, for large enough [25].where if and if . is a nondegenerate stochastic process over S and belongs to the family of generalized -Pareto processes with tail index , zero location, scaling function and limit measure . Specifically, a generalized -Pareto process associated to the limit measure and tail index is a stochastic process taking values in and defined by the following:where and are continuous functions on , respectively, scale and location functions. is a stochastic process whose probability measure is completely determined by the limit measure .
C. Asymptotic Distribution of -Pareto Process
For a threshold vector and for any regularly varying stochastic process , the density function of the -excess of a -Pareto process is obtained by renormalizing suitably the intensity function found by taking the partial derivatives of the previously defined measure (B.6) by :where
while is the intensity function and the region of exceedances. In order to make inferences and to model the dependence structure of -Pareto processes, we focus on the Brown–Resnick model for which the formulas of and are available. The -dimensional intensity function of the Brown–Resnick model is given as follows:where and the covariance matrix.
To better capture the possible dependence structure, we use anisotropic semi-variograms whose parameters change with time and other covariates. The different parameters are estimated using the gradient scoring rule method [24]. For any , the log-density function is given by the following:where is a differentiable weighting function that vanishes on the boundaries of and is an element in the parameter space.
Data Availability
The rainfall data used to support the findings of this study were supplied by the National Meteorology Agency of Burkina Faso. These data were obtained upon request in the framework of the preparation of B. Sawadogo's thesis that led to this study. But the data policy of Burkina Faso do not allow these data to be freely available since they are confidential. However, we can publish the results supporting our study. Requests for access to these data should be addressed to the National Meteorology Agency of Burkina Faso.
Conflicts of Interest
The authors declare that there are no conflicts of interest for the publication of this paper.
Acknowledgments
The authors thank the National Meteorology Agency of Burkina Faso for providing us with the rainfall data. The authors also thank the Cooperation and Cultural Action Service of the French Embassy in Burkina Faso and UMR 518 MIA-Paris-Saclay/AgroParisTech/INRAe for their financial support.