A synthetic view of rainfall intensification in the West African Sahel

The daily rainfall data are extracted from the BADOPLU data base (BAse de DOnnées PLUviomètres; see Le Barbé et al 2002, Lebel and Ali 2009, Panthou et al 2018) for the 1983–2015 period, following the years of lowest rainfall (figure 1(a)). The AMAX samples are extracted from stations having at least 26 valid years over the study period (i.e. ∼80%, see the quality control procedure in the SM of Panthou et al (2018) for the definition of a valid year). The regional AMAX sample gathers all stations meeting the time series completeness criteria located in a $5^\circ \times 30^\circ$ box framing the 10^∘ N–15^∘ N, 20^∘ W–10^∘ E area and corresponding to the West African Sahel definition used in this study (figure S1(a)). Smaller regional samples are also built for investigating the sub-regional variability. The respective influences of the study area and of the study period starting year have been tested with no major influence on the results obtained (see section 3).

In West Africa, daily rainfall at a given point is essentially associated with the passage of well delimited mesoscale convective system (MCS, e.g. Mathon et al 2002), meaning that daily rainfall amounts can be considered as independent or short-term dependent (Ali et al 2005, Gerbaux et al 2009). Moreover, working with AMAX values ensures the independence between block maxima. Under these conditions and according to the Extreme Value Theory (Coles 2001), annual maxima of daily rainfall amounts can be modeled with one type of GEV distributions (equation (S1)). As illustrated in figure 1(b), the behavior of AMAX series in a supposedly homogeneous region may vary significantly from station to station, due to their high sampling dispersion. Regional approaches based on a pooling of the point series allow for reducing the influence of this sampling dispersion, thus providing a more robust inference of the most uncertain components of the extreme value distribution, such as its tail behavior (e.g. Buishand 1991, Blanchet and Davison 2011). This is a central issue when it comes to estimating the probability—return period—of an extremely rare observed event or, reciprocally, the quantile—return level—for an extremely rare probability of occurrence (Papalexiou and Koutsoyiannis 2013). Moreover, regional pooling also improves the trend detection power of statistical methods (e.g. Fischer et al 2013, Fischer and Knutti 2014, Martel et al 2018). Here these two issues—robust statistical inference and trend detection ability—are tightly connected. Building on the work of Panthou et al (2012), we implement a regional GEV (RGEV) model wherein the location (µ) and scale (σ) parameters vary as a linear function of both the latitude and longitude (equation (S2)).

The RNSGEV model is obtained by expressing one or several GEV parameters as a function of time. Generalized linear models are used to model the time-varying GEV parameters (see Panthou et al 2013 and SM) since they are parsimonious as compared to non-linear models and because a linear form seems acceptable in a first approximation in view of the AMAX evolution displayed in figure 1(a). This means that any non-stationary parameter β of the RNSGEV distribution (β = µ, σ, ξ) may be expressed as:

$$\begin{align} \beta(lat, lon, t) = \beta_0 + \beta_1~\times lat + \beta_2~\times lon + \beta_3~\times t. \end{align} \tag{ 1 }$$

Theoretically, there are thus 12 ($4 \times 3$) parameters to be estimated, a far too large number given the available information. Inferring the most appropriate model, based on a given sample, then pertains to identifying which additional parameter provides a significant improvement in terms of maximizing the model's likelihood function, starting from a reference time-stationary regional model M0. Practically, we rather minimize the model's negative log-likelihood (NLLH) function (see SM, section 1.4). The additional flexibility is then tested relying on the semi-parametric bootstrap procedure thoroughly described in Chagnaud et al (2021) (after Katz 2013, see SM, section 1.5). This procedure provides an empirical assessment of the significance of the various compared models while accounting for the spatial dependence between the pooled station samples. The uncertainty of the inferred parameters is quantified as the 5${\textrm{th}}$–95${\textrm{th}}$ percentile range—corresponding to the 90% confidence interval (CI)—obtained from 200 bootstrapped sets of RNSGEV parameters.

A preliminary visual analysis of our Sahelian dataset allows reducing the number of model formulations to be tested. First, the large sampling variance of ξ prevent any robust estimation of either ξ₁, ξ₂ or ξ₃; ξ is thus assumed stationary in space and time (see also e.g. Renard et al 2006, Hanel et al 2009, Chagnaud et al 2021). The RNSGEV model cumulative distribution function (CDF) therefore reads:

$$\begin{align} &GEV_{lat, lon, t}(i)\nonumber\\ &\quad= \textrm{exp} \left\{- \left[ 1 + \xi \left( \frac{i - \mu(lat, lon, t)}{\sigma(lat, lon, t)} \right) \right] ^{-1/ \xi} \right\}. \end{align} \tag{ 2 }$$

Secondly, since neither the Mann–Kendall AMAX trend estimates (figure 1(b)) nor the point-wise non-stationary GEV models indicate any distinct spatial pattern of trends (figure S2), a spatially uniform relative trend over the study region is assumed. Therefore, a multiplicative expression for µ and σ is preferred to the additive one introduced in equation (1). It then remains to rule on the time non-stationarity parameters µ₃ and σ₃. Among the various tested models (SM, section 1.6), the one that proved the most parsimonious while providing the largest likelihood improvement according to the NLLH/bootstrap procedure is the M3 model (p-value $\lt$ 0.01, meaning that the improvement brought by this non-stationary model would happen by chance less than 1% of the time). The location and scale parameters of the M3 model are defined as follows:

$$\begin{align} \begin{cases} \mu(lat, lon, t) = (\mu_0 + \mu_1~\times lat + \mu_2~\times lon)\\\qquad\qquad\quad\quad \times (1 + \alpha_3~\times t) \\ \sigma(lat, lon, t) = (\sigma_0 + \sigma_1~\times lat + \sigma_2~\times lon)\\\qquad\qquad\quad\quad \times (1 + \alpha_3~\times t) \end{cases} \end{align} \tag{ 3 }$$

where $\mu_3 = \sigma_3 = \alpha_3$. This model formulation is characterized by a similar time evolution of µ and σ (in relative terms). This evolution is estimated at 5%/decade ([2.7:7.5] 90% CI) over the considered period. Therefore, the statistical distribution is both shifted toward larger values and widened, as shown in figure 2. Note that the M1 (time-varying location parameter only) and M4 (independently time-varying location and scale parameters) models are also significant when compared to the time-stationary regional model (M0), but perform slightly worse than M3 (figure S3). Also worth mentioning is the fact that the trend undergoes only small variations when periods starting as far back as 1970 are considered (figure S4), demonstrating its significance even though, of course, its value is dependent on the year chosen to start the exploration (the linear formulation is not very flexible in this regard but again, provides the most parsimonious way to model this trend).

Zoom In Zoom Out Reset image size

The T-year return level, defined as the daily rainfall amount having a probability $1/T$ of being exceeded for a given year, is computed as follows at any point in space ($lat, lon$) and time (t):

$$\begin{align} i_T(lat, lon, t) &= \mu(lat, lon, t) + \frac{\sigma(lat, lon, t)}{\xi}\nonumber\\&\quad\times \left \{\left[-\log\left(1 - \frac{1}{T}\right)\right] ^{-\xi} -1~\right \} \end{align} \tag{ 4 }$$

where the location and scale parameters are computed according to equation (3) using the values shown in the inset box on figure 2. Comparing the value of a given return level at both ends (1983 and 2015) of the study period provides an estimate of the change in severity of an event of a given probability of occurrence (calculation details are provided in section 7 of the SM). This is illustrated in figure 2 for the 10-year (circle) and 50-year (triangle) return levels, with the 1983 distribution plotted in blue and the 2015 distribution in green. Also shown are the values of these two return levels under the stationarity assumption (black markers). It is noteworthy that the stationary CI (grey shading) is almost entirely distinct from the non-stationary 2015 CI (green shading), implying that computing return levels under the stationarity assumption likely leads to a significant (at the 90% level) under-estimation of the present-day value, and thus to the undersizing of infrastructures. Note that the formulation of the M3 model involves a similar relative change for all return levels and for all locations (a 16% increase, [8.8:24.7] 90% CI), a point that will be further discussed in section 3.4.

Alternatively to changes in return levels, the non-stationarity can also be expressed as changes in return periods i.e. in the probability of occurrence of a given rainfall amount, as shown in figure 3: while 1983 return periods (horizontal axis) range from 2 to 1000 years, their 2015 counterparts (vertical axis) range from 1.4 to 360. Hence, the change in probability, defined as the $T_{1983}/T_{2015}$ ratio, ranges from 1.4 ([1.2:1.5] 90% CI) for the 2-year return period to 2.8 ([1.8:4.6]) for the largest 1000-year return period.

Zoom In Zoom Out Reset image size

The potential to derive locally-relevant climate change metrics from the regional diagnostic is illustrated in figure 4 with the time evolution of various return periods corresponding to specific daily rainfall amounts at the Niamey station (13.5^∘ N, 2.2^∘ E). Take for instance a daily rainfall event of magnitude 82 mm d⁻¹ (orange line): it is a 10-year rainfall in 1983 and it has become a 5.2-year rainfall in 2015; this rainfall amount is nearly twice as frequent at the end than it is at the start of the study period.

Zoom In Zoom Out Reset image size

The ability of looking at the effects of non-stationarity either in terms of return levels or in terms of return periods is a clear benefit of the statistical modeling framework and an important step forward from previous observationally-based studies (Panthou et al 2014, Taylor et al 2017). The metric to be used, either a change in magnitude for an event with a given probability or a change in probability of a fixed magnitude event will depend on the application and the related risk management policy.

The regional approach presented above allows for a robust detection of the non-stationarity of the annual rainfall maxima in the Sahel, with both a significant shift toward larger values and an increase of their interannual variability, a conclusion that could not be drawn from analyzing separately the available point series. The significance of this non-stationarity is quantified through the NLLH/bootstrap procedure and its effect is straightforwardly derived in terms of increasing return levels and decreasing return periods.

These achievements leave room for debating whether alternate models would more faithfully depict the rainfall regime evolution in this region. The discussion comes down to asking whether we could use a more complex model—and for which benefits—having in mind that, among the four models presented in the SM, three (M1, M2, M3) have the same number (8) of parameters and one (M4) has 9. Two of the 8-parameter models (M1 and M3) were found to be significantly better than the stationary model, M3 being the most significant. M4 is also significant, but less than M3.

M4 is less significant but its additional free parameter gives it a greater flexibility, allowing for an independent time-evolution of µ (4.9%/decade) and σ (3.1%/decade). In M3, µ and σ are assumed to present the same relative evolution, with the effect that return levels corresponding to any probability of occurrence undergo the same relative change. As further investigated in section 4, a different evolution of µ and σ might correspond to a more realistic account of the underlying changes of the rainfall regime. Choosing between the most significant model and a less significant but possibly more realistic alternative is not a matter easily solved, especially since the M4 model provides return level estimates close to those given by M3, yet with a slightly smaller rate of increase (∼4%/decade). That is a fairly open question, keeping in mind that it is well possible that with longer time series M4 could prove as significant as M3. This, however, is just a conjecture at that point.

A still more complex model would be to assume that the shape parameter of the RGEV distribution, ξ, could also vary in time. Given the high sampling variance of ξ, detecting a trend of its value is out of reach with the dataset used here and we thus kept to the conservative hypothesis that it is constant. Exploring whether the increase of return levels could be magnified for large return periods—typically beyond T = 100 years—remains a challenge, and their evolution as provided by the M3 model should be considered very cautiously.

Another issue relates to the spatial behavior of the trend, assumed here to be uniform in relative terms. This is probably not strictly true, as various mechanisms may modulate the overall regional trend. For instance, two types of atmospheric disturbances called African Easterly Waves (AEWs) support the development and propagation of sahelian storms: those originating to the north of the African Easterly Jet and those to the south, accounting for 71% and 29%, respectively, of the total number of AEWs (Chen 2006). Since these two types of AEWs have different genesis mechanisms, one may expect a distinct response to any large-scale forcing, thus potentially influencing the interannual variability as well as a long-term trend in different ways. Here, however, there is no detectable latitudinal gradient affecting the evolution of the regime of extreme rainfall (see SM, section 8.2). Likewise, the trends computed on smaller areas display a striking zonal similarity (SM, section 8.1). As a consequence of this spatial homogeneity, the signal is more robustly detected when considering the Sahel as a whole. Moreover, the derived metrics are associated with narrower CI as compared to their sub-regional counterparts (figure S7 and table 1 in the SM), evidencing a more robust ability for detecting trends when working on large regions, provided that they are climatologically homogeneous (e.g. Donat et al 2016). Whether the absence of spatial gradient is due to a limited detection ability or to the dominance of a mechanism acting on the whole region remains elusive at this point. Again, it may happen that with longer time series, a finer description of the spatial pattern of the trend could emerge.

All in all, the regional vision provided here may be considered as a first order representation of an undoubtedly more complex reality. It is only by developing and maintaining regional- and local-scale in-situ rainfall monitoring systems that we will obtain the necessary information for a more detailed characterization of this reality, needed to support locally-relevant decision-making in various fields pertaining to agriculture, flood warning or water resources management.