# Assessing Nonlinear Dynamics and Trends in Precipitation by Ensemble Empirical Mode Decomposition (EEMD) and Fractal Approach in Benin Republic (West Africa)

Authors: Médard Noukpo Agbazo; Gabin Koto N’Gobi; Eric Alamou; Basile Kounouhewa; Abel Afouda
Journal: Complexity
 (2021)

---

## Abstract

Climate dynamics and trends have significant environmental and socioeconomic impacts; however, in the Benin Republic, they are generally studied with diverse statistical methods ignoring the nonstationarity, nonlinearity, and self-similarity characteristics contained in precipitation time series. This can lead to erroneous conclusions and an unclear understanding of climatic dynamics. Based on daily precipitation data observed in the six synoptic stations of Benin Republic, in the period from 1951 to 2010, we have proposed (i) determining the local trends of precipitations, (ii) investigating precipitation nonlinear dynamics, and (iii) assessing climatic shift in the study period by Ensemble Empirical Mode Decomposition (EEMD) and Multifractal Detrended Fluctuation Analysis (MFDFA) method. To overcome the detrending issue in the standard MFDFA method, the EEMD algorithm is embedded into the MFDFA. The study period is subdivided into three subperiods: 1951–1970, 1971–1990, and 1991–2010. Intrinsic Mode Functions (IMFs) are obtained according to the climatic region in which the stations are located. Results show that precipitation variation is significantly governed by the five first IMFs, in which oscillation periods vary from 1 to 25 days. The trend curves decrease at all the synoptic stations, and their slope values vary accordingly to the subperiods. Referring to the values of the multifractal spectrum parameters,α0, and the width of the spectrum w, consistent changes are observed regardless of the subperiods and the concerned stations. The spatial and temporal variability of precipitation indicates that the multifractal properties are good indicators for assessing changes in precipitation dynamics, and the changes in their features could be explained by the global change than by the local climate variation (climatic zones). Despite the observed differences in multifractal spectra properties from the three subperiods, it is not possible to verify the subdivision of 1951–2010 in three subperiods as it is done by previous studies in West Africa. Our findings can be used in the validation of global and regional climate models since a valid model should explain empirically detected scaling properties in observed data.


---

## Body

## 1. Introduction

Climate variability and trends have enormous influences on the environment and social development of locations with a growing human population [1, 2]. Because many global challenges such as food insecurity, water crisis, biodiversity loss, and health issues are tied to the changing climate, understanding climatic patterns is of great significance [3, 4]. According to Pelletier and Trucotte [5], variability indicates the degree of fluctuation and uncertainty of the climate change process. Climate variability has enormous influences on agricultural and economic development [2]. Precipitations are variable widely collected all over the world and the most affected by climate change and climate variability [6, 7]. Thus, they are the main meteorological variables, used for understanding climate dynamics and climate changes [8]. However, it is not always easy to analyze them because precipitations are highly complex and variable process, exhibiting enormous variability over a wide range of time and space scales [9, 10]. Furthermore, precipitation variations are extremely complex [10], and their time series show obvious nonstationary characteristics in the climatic system [6, 11]. In fact, it is not feasible to detect intrinsic dynamic properties of precipitations due to nonstationarity and nonlinearity by using some traditional approaches, such as the power spectrum, correlation analysis, and the nonparametric Mann–Kendall (MK) test (the most popularly used for trend analysis). These methods are, in general, insufficient and unreliable for a complete analysis of meteorological time series like precipitations [9, 11–13].Wu and Huang [14] have proposed the Ensemble Empirical Mode Decomposition (EEMD) method for processing nonlinear and nonstationary signals, widely applicable to hydrology, climatology, meteorology, and other fields [15]. EEMD can be used to correctly (with higher accuracy) represent the trend in meteorological variables like precipitations without ignorance of its nonstationarity and nonlinearity characteristics, compared to the traditional statistical methods mentioned above [14, 15]. In the Benin Republic, precipitation trend assessment is generally studied using the Mann–Kendall test, Pettit’s test, or linear regression statistical tools. Sometimes, the results can lead to erroneous conclusions because they are affected by the presence of nonstationarity and nonlinearity contained in time series. However, the EEMD method is rarely applied to detect precipitation trends.Many case studies indicated that the regional and local climate is a huge and complex and dynamic system that needs strong and appropriate study tools [16–18]. Shen et al. [19] have demonstrated that climate change has self-memory and self-similarity characteristics, which are important features of fractal and multifractal time series [19]. Moreover, it has been recognized that the climate system has chaotic dynamics with high variability [20–22]. Trend analysis could only reflect the overall change of a climatic variable over one period of time and ignored the variability of climatic dynamics to which human health, crop production, and plant growth are sensitive [22, 23]. Also, the climatic system has similar characteristics on different temporal scales (“self-similarity”). Morata et al. [24], applying fractal theory to assess the variability of climate change, show a more comprehensive picture of the variability of climatic dynamics [23, 24]. Multifractals describe the complex dynamic characteristics of systems more carefully and comprehensively and fully characterize their properties both locally and globally [25, 26]. Climate variability can be described also by using the multifractal characteristics of data sets [23, 26]. Indeed, a multifractal approach leads us to derive trend and seasonality within time series that can help to describe climate variability among specific regions and periods [9, 27, 28].To completely understand the complexity of precipitation dynamics, three scientific questions are arising: What is the real trend within a precipitation time series in Benin synoptic stations? Can the multifractal approach reveal shifts in climate? Do the dynamics of the underlying process govern precipitation changes?The purpose of this paper is, firstly, to decompose precipitation time series into a limited number of intrinsic mode functions (IMFs) and a residual component and, secondly, to investigate changes in precipitation dynamics by EEMD method and improved multifractal approach. The remaining parts of the paper are organized as follows: in Section2, study sites, the dataset, and methods are described. Analysis of results and their interpretations are made in Section 3. Finally, the paper is concluded with a summary and outlook for further research in Section 4.


## 2. Materials and Methods

### 2.1. Materials

#### 2.1.1. Study Sites

The study covers all the synoptic stations of Benin, the geographical positions of which are presented in Figure1. The occurrence of the rainy season is related to the seasonal variations of the Intertropical Convergence Zone (ITCZ), in the southern part of the ITF [29, 30], between 100 and 200 kilometers, where the monsoon wedge is deeper [31]. The seasonal variations of the ICTZ, as well as the distribution of the precipitations, depend on the regions. They depend also on the geographic position of the stations (South-North gradient of 2°) and their intrinsic characteristics. For example, the presence of mountain chains in Natitingou allows the development of convective system cells. So, Natitingou and Parakou remain more humid (∼1300 mm/year and 1200 mm/year annual amount of precipitation, respectively) than Kandi (∼700 mm/year) [32]. Benin is characterized from the South to North by three climatic zones in which stations are located [11, 33]: Cotonou, Bohicon, and Save are located in the subequatorial region where March is the hottest month (∼26°C), while August is the coldest month (∼24°C). The daily and annual thermal amplitudes are, respectively, ∼10°C and∼5°C. The relative humidity ranges between 70% and 95% because of the proximity to the Atlantic Ocean. A coastline of around 121 km long (in the Gulf of Guinea) separated the country from the Atlantic Ocean. The subequatorial climate has four seasons: a long rainy season (April to July) followed by a short dry season (July-August to September) and a short rainy season (September to October) followed by a long dry season (November-December to March) in the year. However, the stations of Parakou, Kandi, and Natitingou are located in the Sudanian region in the northern part of the country. The daily mean of air temperatures in Natitingou, Parakou, and Kandi is, respectively, ∼25°C, ∼27°C, and ∼35°C. Parakou and Save (200 km far) are located in the transition area between the two kinds of climatic zone. The “Harmattan” flux blows from the Sahara Desert toward the Atlantic Ocean and covers, each season, the study area with dusty and dry air (known as “Harmattan”) and influences circulations in the northern part of the country than the other locations.Figure 1
Study site location: (a) Benin location in West Africa; (b) synoptic station locations in Benin.


#### 2.1.2. Data Records

Data were provided by the Agency for the Aerial Navigation’s Security in Africa and Madagascar (ASECNA). Daily precipitation data from Benin synoptic stations are used during a period from 1951 to 2010. The descriptive statistics of the precipitation time series collected during the study period are presented in Table1.Table 1
Statistical description of the whole sixty-year daily meteorological time series from Benin’s synoptic stations. .
Synoptic stationsMinMaxSkewnessKurtosisCotonou0194.24.6833.15Bohicon0128.54.0023.60Save0184.24.3130.60Parakou01764.0027.17Natitingou0149.73.2017.487Kandi0157.83.6922.62Co: Cotonou; Bo: Bohicon; Sa: Save; Pa: Parakou; Na: Natitingou; Ka: Kandi. Min and max have units corresponding to the units of the meteorological variable (mm.d-1); skewness and kurtosis are nondimensional.The highest maximum value of precipitation during the study period is observed at Cotonou station. In contrast, the lowest value is noticed at Bohicon station. The skewness and kurtosis parameters of precipitation time series present information about differences in their statistical distributions. Precipitation reveals a positive skewness at all the synoptic stations which indicates that the distributions are right-tailed. The lowest skewness and kurtosis values are observed at Natitingou, informing that the precipitation’s distribution at this station has a more rounded peak and thinner tails, compared to the other stations. However, the highest skewness and kurtosis values are observed at Cotonou, indicating strongly right-tailed distributions with sharp peaks and fat tails.Figure2 represents the occurrence of daily rainfall during the 60 years of observations in Cotonou synoptic station, as a disjoint geometric set composed of elementary segments supported by the time axis. This figure evokes the object obtained after a certain number of iterations in the random process that generates Cantor's dust. The figure also reveals the random and nonlinear nature of the rainfall process in the measurement site.Figure 2
Daily rainfall occurrence during 60 years of observations in the Cotonou synoptic station. Each row corresponds to a year of observations. The occurrence of rain is defined by a precipitation threshold that distinguishes between dry and rainy days. The precipitation threshold chosen in this representation is zero.


### 2.2. Methods

#### 2.2.1. Seasonal Detrending

Natural time series like precipitations generally exhibit a seasonal cycle; thus, there are usually seasonal variations in precipitation data. To make sure that seasonality does not affect the Multifractal Detrended Fluctuation Analysis (MFDFA) method, the seasonal detrending is done. In this study, a deseasonalized time serieszt is obtained from a procedure of the original series xt. The method is performed by adjusting the data with the seasonal mean and standard deviation as in [18, 34, 35]:(1)zt=xt−xtdσd,where xtd is the mean of the observed daily precipitation’s quantity calculated for each calendar date d (obtained by averaging over all the years in the record) and σd is the corresponding standard deviation (also calculated for each calendar date). The deseasonalized precipitation forms a new time series for MFDFA analysis.


#### 2.2.2. Description of MFDFA Method

Based on the Detrended Fluctuation Analysis (DFA), the Multifractal Detrended Fluctuation Analysis was proposed by Kantelhardt et al. [27]. In literature, the MFDFA method is an effective tool to detect the scaling behaviors and multifractal properties of nonlinear and nonstationary time series such as precipitation records in hydrology and hydrometeorology [27, 28, 36]. Owing to its simplicity and robustness, MFDFA becomes very popular for the fractal characterization of daily, monthly, or annual precipitation records worldwide [9, 10, 13, 27, 28, 34, 35, 37–39].According to Kantelhardt et al. [27], Li et al. [39], Krzyszczak et al. [26], Adarsh et al. [40], and Baranowski et al. [13], the different steps involved in MFDFA computational procedure can be described as follows:Step 1: letzt|t=1,2,…,N be a deseasonalized precipitation time series of N equidistant measurements to which the procedure of the MFDFA method is applied. The accumulated deviation of the series (known as “profile”) is calculated as(2)Yk=∑t=1kzt−z,wherez is the mean value of the series, zt and k = 1, 2, …, N, and t = 1, 2, …, N.Step 2: divide the profileYk into Ns = int(N/s) equal-sized nonoverlapping windows with a length s. While considering multiple timescales, sometimes, a small portion of the time series at the end may be left out, as N need not be a multiple of s always. To retain this part of the series, the same procedure is repeated starting from the opposite end resulting in 2 Ns segments.Step 3: calculate the local trend for each of the 2Ns segments by a least-squares fit of the series as(3)F2ν,s=1s∑i=1sYν−1s+i−yνmi2,forν=1,2,…,Ns,1s∑i=1sYN−ν−Nss+i−yνmi2,forν=Ns+1,…,2Ns,whereyνi is the trend function in segment ν. yνi stands for the fitted polynomial trend in the vth segment. Different types of fittings such as linear (i.e., polynomial order m = 1), quadratic (i.e., m = 2), cubic (i.e., m = 3), or higher-order polynomials yνmi can be used to ﬁt the local trend. Determine the qth-order fluctuation function by averaging over all segments from the following equation:(4)Fqs=12Ns∑ν=12NsF2ν,sq/21/q,forq≠0,exp14Ns∑ν=12NslnF2ν,s,forq=0.Step 4: determine the scaling behavior of the fluctuation functions by analyzing thelog-log plots of Fqs versus s for each value of q (q is the order of moment). If the time series is long-range power-law correlated, Fqs increases as(5)Fqs∝shq,where h (q) is the generalized Hurst scaling function [27]. For monofractal time series, h (q) is independent of q, whereas, for a multifractal time series, h (q) depends on q. To avoid divergence of moments in the fat tails of the fluctuation distribution as mentioned by some authors [41, 42], one restricts the order q within the range −5 ≤ q ≤ 5. The density distributions of all the studied meteorological time series had heavy tails. To prevent and avoid potential distortion of the results by the so-called “freezing” phenomenon [9, 13, 28], the range of q has to be limited in the interval −5,5, as in [9, 13, 41].


#### 2.2.3. Multifractal Spectra Analysis (MFS)

According to [9, 13, 27, 43–45], using the Legendre transform, the relationship between the multifractal spectrum fα and the generalized Hurst index hq can be obtained and written as(6)αq=hq+qdhqdq,fα=qαq−hq+1,where αq and fα are the singularity exponent (or Hölder exponent) and the multifractal spectrum, respectively. Additionally, the mass exponent function τ (q) can be calculated by the generalized Hurst exponent hq as(7)τq=qhq−1.The schematic representation of a multifractal spectrum with its most important parametersαmax, αmin, α0, as, and w is shown in Figure 3 (obtained from Baranowski et al. [9], Krzyszczak et al. [26], and Baranowski et al. [13]).Figure 3
Schematic presentation of the main parameters (position of maximumα0, asymmetry as, and width w) of a multifractal spectrum fα. Parameters α0, as, and w are dimensionless (modified figure from Baranowski et al., [13]).The parameterαmin represents the most extreme event and αmax the smoothest event in the studied process. The value of α0 delivers valuable information about the structure of the studied process. It indicates at which value of α multifractal spectrum reaches its maximum. A low value of α0 indicates that the underlying process becomes correlated and loses fine structure, becoming more regular in appearance.The negative values of the asymmetry parameteras (left-skewed shapes of the multifractal spectra) indicate low fractal exponents of small weights. They suggest the dominance of large fluctuations which can be connected with extreme events [9, 13]. A right skewness of the multifractal spectrum (positive value of as) denotes high fractal exponents with large weights, which are characteristic of fine structures and indicate a relative dominance of the small fluctuations. The width of the spectrum w (the difference between αmax and αmin) represents the length of the fractal exponent range in the series. It is strongly connected to the “richness” of the signal structure. The greater the value of w is, the more developed the multifractality is. In contrast, for pure monofractal, the width w of the spectra is equal to 0, but, according to Makowiec and Fuliński [46], if w is less than 0.05, then monofractal behavior of the spectrum should be assumed. Therefore, the values of the multifractal spectrum parameters (α0,w, and as) can be used as quantitative and qualitative indicators for studying the dynamics of meteorological processes [9, 13, 26].Additionally, according to Hou et al. [25], the asymmetric index (R), as presented in equation (8), is a useful parameter of multifractal analysis. It is computed based on the width of the left (ΔαL) and right (ΔαR) hand branches of the multifractal spectrum, respectively. Their values describe the distribution patterns of high and low fluctuations [47]. The values of R range from −1 to 1, quantifying the deviations of the multifractal spectrum curve [25, 48]. Positive values of R suggest a left-hand deviation of the multifractal spectrum, with high local fluctuations; negative values of R suggest a right-hand deviation with local low fluctuations; R-value equal to zero represents a symmetrical multifractal spectrum. The asymmetry index (R) is obtained as(8)R=ΔαL−ΔαRΔαL+ΔαR,where ΔαL=α0−αmin and ΔαR=αmax−α0.


#### 2.2.4. Description of the EEMD Algorithm

The empirical mode decomposition (EMD) method proposed by Huang and Wu [49] is widely used for processing nonlinear and nonstationary signals in hydrology, hydrometeorology, and other fields [15, 49]. EMD method can decompose time series data into a limited number of intrinsic mode functions (IMFs) and a residual component (trend). The IMFs satisfy the following two main requirements [15]: (i) in the whole set of data, the number of local extrema and the numbers of zero-crossings must be equal to or different from 1 at most and (ii) at any point, the mean value of the ‘‘upper envelope’’ (defined by the local maxima) and the ‘‘lower envelope’’ (defined by the local minima) must be zero. According to these properties, some meaningful IMFs can be well defined. Generally, an IMF indicates a simple oscillatory mode, compared to a simple harmonic function. Based on this definition, a shifting process of the original time series can be briefly expressed by Liu et al. [15] in steps as follows:(1)
Identify all extrema (local maxima and minima) points of the given time seriesy (t)(2)
Connect that extrema by spline interpolation, to create an upper envelope of local maxima pointsemaxt, and a lower envelope emint, of all minima points(3)
Compute the meanm (t) of two envelopes by m=emaxt+emint/2(4)
Subtract the mean from the data to establish an IMF candidate;ht=yt−mt(5)
Ifht satisfies the two properties of IMF as indicated by a predefined stopping criterion, ht is indicated as the first IMF and is recorded as IMFit; if ht is not an IMF, yt is replaced by ht, and iterate steps 1–4 until h (t) satisfies the two conditions of IMFs(6)
The residuerit=yt−IMFit is then treated as new data subjected to the same shifting procedure as defined above for the next IMF from rit. In the final step, the shifting process can be stopped when the residue rt has a monotonic trend or has one local maxima and minima point from which no more IMFs can be extracted. At the end of this shifting process, the original signal yt can be reconstructed as the sum of IMFs and residual as(9)yt=∑i=1nIMFit+rnt,where rnt represents the final residue and n is the number of IMFs. The residue rnt is known as the trend of the signal yt [15].Although EMD has obvious advantages in time series data processing, there are also some unavoidable weaknesses [14, 15]. The major weakness is mode-mixing. Mode-mixing means that the same IMF contains different frequency components or the frequency of the same and similar scales is distributed in different IMFs. Thus, mode-mixing not only causes the mixing of various scale vibration modes but can even change the physical meaning of the individual IMF. To overcome these drawbacks of EMD, Ensemble Empirical Mode Decomposition (EEMD) is proposed in Wu and Huang [14] using white noises. The EEMD algorithm is straight forward and can be summarized as follows:(1)
Add the white noise series into the original series(2)
Decompose the signal with the added white noise into the IMFs using EMD(3)
Repeat steps (1) and (2) with a different white noise series each time(4)
Obtain the corresponding IMF components of the decomposition and calculate the means of the ensemble corresponding to the IMFs of the decomposition as the final result [14]For the EEMD method, the first important step is to determine the ensemble times and the amplitude of the added noise. But the way of selecting the best ensemble number and amplitude of added noise is still an open question. The effect of adding white noise should obey the following statistical rule:(10)εn=εN,where εn is the final standard deviation of the error, which is described as the difference between the input signal and the corresponding IMFs. ε denotes the amplitude of the added noise, and N is the number of ensemble numbers. In this study, the ensemble number was set to 100, and the amplitude of added white noise was set to 0.2 times the standard deviation of the corresponding data as in [10, 15].


#### 2.2.5. Period Computation Method and Criterion for Selecting Relevant IMFs

Since each IMF represents different intrinsic modes of oscillation, it is possible to calculate the period for each of them. The average periods are calculated in this study by the time intervals between consecutive zero-crossings on successive waves as in [50].Since the IMFs are supposed to be almost orthogonal components of the original signal, each IMF would have a relatively good correlation with the original signal; this presupposes that the irrelevant components would have a relatively poor correlation with the original signal [51]. For discriminating between relevant and irrelevant IMFs, the following threshold has been proposed by Ayenu-Prah and Attoh-Okine [52].(11)λ=maxμi10∗maxμi−3,where i=1,2,…,n and λ is the threshold. Moreover, μi is Pearson’s correlation coefficient between the ith IMF and the original signal, and n is the total number of IMFs; maxμi is the maximum observed correlation coefficient. The selection criterion for IMFs is given as follows: if μi>λ, then keep the ith IMF, else eliminate the ith IMF, and add it into the residue.


#### 2.2.6. EEMD-MFDFA-Based Method

The selection of the scale range for computing the fluctuation functionFqs and the type of polynomial chosen (detrending) is one of the major issues in applying the MFDFA method [10, 40]. Indeed, firstly, if there is a trend in a real-time series, the functional form of this trend is usually unknown. For real-time series, it is usually unknown whether there is a trend component and, if so, what the functional form of this trend is. A conventional strategy is to assume that local trends are in the form of a polynomial. The discontinuity at the segmentation points using polynomial fitting comes from the new pseudofluctuation errors and estimation errors in the qth-order generalized Hurst exponent hq. Secondly, the qth-order fluctuation function Fqs depends strongly on the selection of the multiple segment sizes s. To overcome the detrending issue in the standard MFDFA method, the EEMD algorithm is embedded into the MFDFA to modify the third step of the algorithms, while keeping all other steps unchanged as in [10, 45]. Thus, in the standard MFDFA method, the third step is reformulated as follows:Step 3: For each segmentν, one obtains the EEMD-based local trend yνi=rνi with the shifting process. Note that the trend rνi should be determined for each segment separately at each time scale. One can then obtain the variance as(12)F2ν,s=1s∑i=1sYν−1s+i−rνi2,forν=1,2,…,Ns,1s∑i=1sYN−ν−Nss+i−rνi2,forν=Ns+1,…,2Ns.We denote this MFDFA with the EEMD-based detrending method as EEMD-MFDFA. The criterion for selecting the range of scale (s) is the highest possible stability of the obtained spectra as in [9, 13, 26].According to Verrier et al. [53], zero values are susceptible to bias multifractal analysis results. To avoid the effect of the zeros on multifractal analysis results, the precipitation data recorded from April to October are selected to take into account the rainy season in all the Benin synoptic stations. To determine whether changes in the dynamics of precipitation time series can be assessed with the MFDFA method, the sixty-year data were divided into three separated subdatasets: the first is from 1951 to 1970, the second is from 1971 to 1990, and the third is from 1991 to 2010. This division was made considering the classification often mentioned by some studies in West Africa [54, 55]. The period 1951–1970 is considered as the wet period, 1971–1990 is the drought period, and 1991–2010 is considered as the end of the drought. In each synoptic station, the comparison of precipitation multifractal spectrum parameters values (α0,w,as) in the subsets is done.


## 2.1. Materials

### 2.1.1. Study Sites

The study covers all the synoptic stations of Benin, the geographical positions of which are presented in Figure1. The occurrence of the rainy season is related to the seasonal variations of the Intertropical Convergence Zone (ITCZ), in the southern part of the ITF [29, 30], between 100 and 200 kilometers, where the monsoon wedge is deeper [31]. The seasonal variations of the ICTZ, as well as the distribution of the precipitations, depend on the regions. They depend also on the geographic position of the stations (South-North gradient of 2°) and their intrinsic characteristics. For example, the presence of mountain chains in Natitingou allows the development of convective system cells. So, Natitingou and Parakou remain more humid (∼1300 mm/year and 1200 mm/year annual amount of precipitation, respectively) than Kandi (∼700 mm/year) [32]. Benin is characterized from the South to North by three climatic zones in which stations are located [11, 33]: Cotonou, Bohicon, and Save are located in the subequatorial region where March is the hottest month (∼26°C), while August is the coldest month (∼24°C). The daily and annual thermal amplitudes are, respectively, ∼10°C and∼5°C. The relative humidity ranges between 70% and 95% because of the proximity to the Atlantic Ocean. A coastline of around 121 km long (in the Gulf of Guinea) separated the country from the Atlantic Ocean. The subequatorial climate has four seasons: a long rainy season (April to July) followed by a short dry season (July-August to September) and a short rainy season (September to October) followed by a long dry season (November-December to March) in the year. However, the stations of Parakou, Kandi, and Natitingou are located in the Sudanian region in the northern part of the country. The daily mean of air temperatures in Natitingou, Parakou, and Kandi is, respectively, ∼25°C, ∼27°C, and ∼35°C. Parakou and Save (200 km far) are located in the transition area between the two kinds of climatic zone. The “Harmattan” flux blows from the Sahara Desert toward the Atlantic Ocean and covers, each season, the study area with dusty and dry air (known as “Harmattan”) and influences circulations in the northern part of the country than the other locations.Figure 1
Study site location: (a) Benin location in West Africa; (b) synoptic station locations in Benin.


### 2.1.2. Data Records

Data were provided by the Agency for the Aerial Navigation’s Security in Africa and Madagascar (ASECNA). Daily precipitation data from Benin synoptic stations are used during a period from 1951 to 2010. The descriptive statistics of the precipitation time series collected during the study period are presented in Table1.Table 1
Statistical description of the whole sixty-year daily meteorological time series from Benin’s synoptic stations. .
Synoptic stationsMinMaxSkewnessKurtosisCotonou0194.24.6833.15Bohicon0128.54.0023.60Save0184.24.3130.60Parakou01764.0027.17Natitingou0149.73.2017.487Kandi0157.83.6922.62Co: Cotonou; Bo: Bohicon; Sa: Save; Pa: Parakou; Na: Natitingou; Ka: Kandi. Min and max have units corresponding to the units of the meteorological variable (mm.d-1); skewness and kurtosis are nondimensional.The highest maximum value of precipitation during the study period is observed at Cotonou station. In contrast, the lowest value is noticed at Bohicon station. The skewness and kurtosis parameters of precipitation time series present information about differences in their statistical distributions. Precipitation reveals a positive skewness at all the synoptic stations which indicates that the distributions are right-tailed. The lowest skewness and kurtosis values are observed at Natitingou, informing that the precipitation’s distribution at this station has a more rounded peak and thinner tails, compared to the other stations. However, the highest skewness and kurtosis values are observed at Cotonou, indicating strongly right-tailed distributions with sharp peaks and fat tails.Figure2 represents the occurrence of daily rainfall during the 60 years of observations in Cotonou synoptic station, as a disjoint geometric set composed of elementary segments supported by the time axis. This figure evokes the object obtained after a certain number of iterations in the random process that generates Cantor's dust. The figure also reveals the random and nonlinear nature of the rainfall process in the measurement site.Figure 2
Daily rainfall occurrence during 60 years of observations in the Cotonou synoptic station. Each row corresponds to a year of observations. The occurrence of rain is defined by a precipitation threshold that distinguishes between dry and rainy days. The precipitation threshold chosen in this representation is zero.


## 2.1.1. Study Sites

The study covers all the synoptic stations of Benin, the geographical positions of which are presented in Figure1. The occurrence of the rainy season is related to the seasonal variations of the Intertropical Convergence Zone (ITCZ), in the southern part of the ITF [29, 30], between 100 and 200 kilometers, where the monsoon wedge is deeper [31]. The seasonal variations of the ICTZ, as well as the distribution of the precipitations, depend on the regions. They depend also on the geographic position of the stations (South-North gradient of 2°) and their intrinsic characteristics. For example, the presence of mountain chains in Natitingou allows the development of convective system cells. So, Natitingou and Parakou remain more humid (∼1300 mm/year and 1200 mm/year annual amount of precipitation, respectively) than Kandi (∼700 mm/year) [32]. Benin is characterized from the South to North by three climatic zones in which stations are located [11, 33]: Cotonou, Bohicon, and Save are located in the subequatorial region where March is the hottest month (∼26°C), while August is the coldest month (∼24°C). The daily and annual thermal amplitudes are, respectively, ∼10°C and∼5°C. The relative humidity ranges between 70% and 95% because of the proximity to the Atlantic Ocean. A coastline of around 121 km long (in the Gulf of Guinea) separated the country from the Atlantic Ocean. The subequatorial climate has four seasons: a long rainy season (April to July) followed by a short dry season (July-August to September) and a short rainy season (September to October) followed by a long dry season (November-December to March) in the year. However, the stations of Parakou, Kandi, and Natitingou are located in the Sudanian region in the northern part of the country. The daily mean of air temperatures in Natitingou, Parakou, and Kandi is, respectively, ∼25°C, ∼27°C, and ∼35°C. Parakou and Save (200 km far) are located in the transition area between the two kinds of climatic zone. The “Harmattan” flux blows from the Sahara Desert toward the Atlantic Ocean and covers, each season, the study area with dusty and dry air (known as “Harmattan”) and influences circulations in the northern part of the country than the other locations.Figure 1
Study site location: (a) Benin location in West Africa; (b) synoptic station locations in Benin.


## 2.1.2. Data Records

Data were provided by the Agency for the Aerial Navigation’s Security in Africa and Madagascar (ASECNA). Daily precipitation data from Benin synoptic stations are used during a period from 1951 to 2010. The descriptive statistics of the precipitation time series collected during the study period are presented in Table1.Table 1
Statistical description of the whole sixty-year daily meteorological time series from Benin’s synoptic stations. .
Synoptic stationsMinMaxSkewnessKurtosisCotonou0194.24.6833.15Bohicon0128.54.0023.60Save0184.24.3130.60Parakou01764.0027.17Natitingou0149.73.2017.487Kandi0157.83.6922.62Co: Cotonou; Bo: Bohicon; Sa: Save; Pa: Parakou; Na: Natitingou; Ka: Kandi. Min and max have units corresponding to the units of the meteorological variable (mm.d-1); skewness and kurtosis are nondimensional.The highest maximum value of precipitation during the study period is observed at Cotonou station. In contrast, the lowest value is noticed at Bohicon station. The skewness and kurtosis parameters of precipitation time series present information about differences in their statistical distributions. Precipitation reveals a positive skewness at all the synoptic stations which indicates that the distributions are right-tailed. The lowest skewness and kurtosis values are observed at Natitingou, informing that the precipitation’s distribution at this station has a more rounded peak and thinner tails, compared to the other stations. However, the highest skewness and kurtosis values are observed at Cotonou, indicating strongly right-tailed distributions with sharp peaks and fat tails.Figure2 represents the occurrence of daily rainfall during the 60 years of observations in Cotonou synoptic station, as a disjoint geometric set composed of elementary segments supported by the time axis. This figure evokes the object obtained after a certain number of iterations in the random process that generates Cantor's dust. The figure also reveals the random and nonlinear nature of the rainfall process in the measurement site.Figure 2
Daily rainfall occurrence during 60 years of observations in the Cotonou synoptic station. Each row corresponds to a year of observations. The occurrence of rain is defined by a precipitation threshold that distinguishes between dry and rainy days. The precipitation threshold chosen in this representation is zero.


## 2.2. Methods

### 2.2.1. Seasonal Detrending

Natural time series like precipitations generally exhibit a seasonal cycle; thus, there are usually seasonal variations in precipitation data. To make sure that seasonality does not affect the Multifractal Detrended Fluctuation Analysis (MFDFA) method, the seasonal detrending is done. In this study, a deseasonalized time serieszt is obtained from a procedure of the original series xt. The method is performed by adjusting the data with the seasonal mean and standard deviation as in [18, 34, 35]:(1)zt=xt−xtdσd,where xtd is the mean of the observed daily precipitation’s quantity calculated for each calendar date d (obtained by averaging over all the years in the record) and σd is the corresponding standard deviation (also calculated for each calendar date). The deseasonalized precipitation forms a new time series for MFDFA analysis.


### 2.2.2. Description of MFDFA Method

Based on the Detrended Fluctuation Analysis (DFA), the Multifractal Detrended Fluctuation Analysis was proposed by Kantelhardt et al. [27]. In literature, the MFDFA method is an effective tool to detect the scaling behaviors and multifractal properties of nonlinear and nonstationary time series such as precipitation records in hydrology and hydrometeorology [27, 28, 36]. Owing to its simplicity and robustness, MFDFA becomes very popular for the fractal characterization of daily, monthly, or annual precipitation records worldwide [9, 10, 13, 27, 28, 34, 35, 37–39].According to Kantelhardt et al. [27], Li et al. [39], Krzyszczak et al. [26], Adarsh et al. [40], and Baranowski et al. [13], the different steps involved in MFDFA computational procedure can be described as follows:Step 1: letzt|t=1,2,…,N be a deseasonalized precipitation time series of N equidistant measurements to which the procedure of the MFDFA method is applied. The accumulated deviation of the series (known as “profile”) is calculated as(2)Yk=∑t=1kzt−z,wherez is the mean value of the series, zt and k = 1, 2, …, N, and t = 1, 2, …, N.Step 2: divide the profileYk into Ns = int(N/s) equal-sized nonoverlapping windows with a length s. While considering multiple timescales, sometimes, a small portion of the time series at the end may be left out, as N need not be a multiple of s always. To retain this part of the series, the same procedure is repeated starting from the opposite end resulting in 2 Ns segments.Step 3: calculate the local trend for each of the 2Ns segments by a least-squares fit of the series as(3)F2ν,s=1s∑i=1sYν−1s+i−yνmi2,forν=1,2,…,Ns,1s∑i=1sYN−ν−Nss+i−yνmi2,forν=Ns+1,…,2Ns,whereyνi is the trend function in segment ν. yνi stands for the fitted polynomial trend in the vth segment. Different types of fittings such as linear (i.e., polynomial order m = 1), quadratic (i.e., m = 2), cubic (i.e., m = 3), or higher-order polynomials yνmi can be used to ﬁt the local trend. Determine the qth-order fluctuation function by averaging over all segments from the following equation:(4)Fqs=12Ns∑ν=12NsF2ν,sq/21/q,forq≠0,exp14Ns∑ν=12NslnF2ν,s,forq=0.Step 4: determine the scaling behavior of the fluctuation functions by analyzing thelog-log plots of Fqs versus s for each value of q (q is the order of moment). If the time series is long-range power-law correlated, Fqs increases as(5)Fqs∝shq,where h (q) is the generalized Hurst scaling function [27]. For monofractal time series, h (q) is independent of q, whereas, for a multifractal time series, h (q) depends on q. To avoid divergence of moments in the fat tails of the fluctuation distribution as mentioned by some authors [41, 42], one restricts the order q within the range −5 ≤ q ≤ 5. The density distributions of all the studied meteorological time series had heavy tails. To prevent and avoid potential distortion of the results by the so-called “freezing” phenomenon [9, 13, 28], the range of q has to be limited in the interval −5,5, as in [9, 13, 41].


### 2.2.3. Multifractal Spectra Analysis (MFS)

According to [9, 13, 27, 43–45], using the Legendre transform, the relationship between the multifractal spectrum fα and the generalized Hurst index hq can be obtained and written as(6)αq=hq+qdhqdq,fα=qαq−hq+1,where αq and fα are the singularity exponent (or Hölder exponent) and the multifractal spectrum, respectively. Additionally, the mass exponent function τ (q) can be calculated by the generalized Hurst exponent hq as(7)τq=qhq−1.The schematic representation of a multifractal spectrum with its most important parametersαmax, αmin, α0, as, and w is shown in Figure 3 (obtained from Baranowski et al. [9], Krzyszczak et al. [26], and Baranowski et al. [13]).Figure 3
Schematic presentation of the main parameters (position of maximumα0, asymmetry as, and width w) of a multifractal spectrum fα. Parameters α0, as, and w are dimensionless (modified figure from Baranowski et al., [13]).The parameterαmin represents the most extreme event and αmax the smoothest event in the studied process. The value of α0 delivers valuable information about the structure of the studied process. It indicates at which value of α multifractal spectrum reaches its maximum. A low value of α0 indicates that the underlying process becomes correlated and loses fine structure, becoming more regular in appearance.The negative values of the asymmetry parameteras (left-skewed shapes of the multifractal spectra) indicate low fractal exponents of small weights. They suggest the dominance of large fluctuations which can be connected with extreme events [9, 13]. A right skewness of the multifractal spectrum (positive value of as) denotes high fractal exponents with large weights, which are characteristic of fine structures and indicate a relative dominance of the small fluctuations. The width of the spectrum w (the difference between αmax and αmin) represents the length of the fractal exponent range in the series. It is strongly connected to the “richness” of the signal structure. The greater the value of w is, the more developed the multifractality is. In contrast, for pure monofractal, the width w of the spectra is equal to 0, but, according to Makowiec and Fuliński [46], if w is less than 0.05, then monofractal behavior of the spectrum should be assumed. Therefore, the values of the multifractal spectrum parameters (α0,w, and as) can be used as quantitative and qualitative indicators for studying the dynamics of meteorological processes [9, 13, 26].Additionally, according to Hou et al. [25], the asymmetric index (R), as presented in equation (8), is a useful parameter of multifractal analysis. It is computed based on the width of the left (ΔαL) and right (ΔαR) hand branches of the multifractal spectrum, respectively. Their values describe the distribution patterns of high and low fluctuations [47]. The values of R range from −1 to 1, quantifying the deviations of the multifractal spectrum curve [25, 48]. Positive values of R suggest a left-hand deviation of the multifractal spectrum, with high local fluctuations; negative values of R suggest a right-hand deviation with local low fluctuations; R-value equal to zero represents a symmetrical multifractal spectrum. The asymmetry index (R) is obtained as(8)R=ΔαL−ΔαRΔαL+ΔαR,where ΔαL=α0−αmin and ΔαR=αmax−α0.


### 2.2.4. Description of the EEMD Algorithm

The empirical mode decomposition (EMD) method proposed by Huang and Wu [49] is widely used for processing nonlinear and nonstationary signals in hydrology, hydrometeorology, and other fields [15, 49]. EMD method can decompose time series data into a limited number of intrinsic mode functions (IMFs) and a residual component (trend). The IMFs satisfy the following two main requirements [15]: (i) in the whole set of data, the number of local extrema and the numbers of zero-crossings must be equal to or different from 1 at most and (ii) at any point, the mean value of the ‘‘upper envelope’’ (defined by the local maxima) and the ‘‘lower envelope’’ (defined by the local minima) must be zero. According to these properties, some meaningful IMFs can be well defined. Generally, an IMF indicates a simple oscillatory mode, compared to a simple harmonic function. Based on this definition, a shifting process of the original time series can be briefly expressed by Liu et al. [15] in steps as follows:(1)
Identify all extrema (local maxima and minima) points of the given time seriesy (t)(2)
Connect that extrema by spline interpolation, to create an upper envelope of local maxima pointsemaxt, and a lower envelope emint, of all minima points(3)
Compute the meanm (t) of two envelopes by m=emaxt+emint/2(4)
Subtract the mean from the data to establish an IMF candidate;ht=yt−mt(5)
Ifht satisfies the two properties of IMF as indicated by a predefined stopping criterion, ht is indicated as the first IMF and is recorded as IMFit; if ht is not an IMF, yt is replaced by ht, and iterate steps 1–4 until h (t) satisfies the two conditions of IMFs(6)
The residuerit=yt−IMFit is then treated as new data subjected to the same shifting procedure as defined above for the next IMF from rit. In the final step, the shifting process can be stopped when the residue rt has a monotonic trend or has one local maxima and minima point from which no more IMFs can be extracted. At the end of this shifting process, the original signal yt can be reconstructed as the sum of IMFs and residual as(9)yt=∑i=1nIMFit+rnt,where rnt represents the final residue and n is the number of IMFs. The residue rnt is known as the trend of the signal yt [15].Although EMD has obvious advantages in time series data processing, there are also some unavoidable weaknesses [14, 15]. The major weakness is mode-mixing. Mode-mixing means that the same IMF contains different frequency components or the frequency of the same and similar scales is distributed in different IMFs. Thus, mode-mixing not only causes the mixing of various scale vibration modes but can even change the physical meaning of the individual IMF. To overcome these drawbacks of EMD, Ensemble Empirical Mode Decomposition (EEMD) is proposed in Wu and Huang [14] using white noises. The EEMD algorithm is straight forward and can be summarized as follows:(1)
Add the white noise series into the original series(2)
Decompose the signal with the added white noise into the IMFs using EMD(3)
Repeat steps (1) and (2) with a different white noise series each time(4)
Obtain the corresponding IMF components of the decomposition and calculate the means of the ensemble corresponding to the IMFs of the decomposition as the final result [14]For the EEMD method, the first important step is to determine the ensemble times and the amplitude of the added noise. But the way of selecting the best ensemble number and amplitude of added noise is still an open question. The effect of adding white noise should obey the following statistical rule:(10)εn=εN,where εn is the final standard deviation of the error, which is described as the difference between the input signal and the corresponding IMFs. ε denotes the amplitude of the added noise, and N is the number of ensemble numbers. In this study, the ensemble number was set to 100, and the amplitude of added white noise was set to 0.2 times the standard deviation of the corresponding data as in [10, 15].


### 2.2.5. Period Computation Method and Criterion for Selecting Relevant IMFs

Since each IMF represents different intrinsic modes of oscillation, it is possible to calculate the period for each of them. The average periods are calculated in this study by the time intervals between consecutive zero-crossings on successive waves as in [50].Since the IMFs are supposed to be almost orthogonal components of the original signal, each IMF would have a relatively good correlation with the original signal; this presupposes that the irrelevant components would have a relatively poor correlation with the original signal [51]. For discriminating between relevant and irrelevant IMFs, the following threshold has been proposed by Ayenu-Prah and Attoh-Okine [52].(11)λ=maxμi10∗maxμi−3,where i=1,2,…,n and λ is the threshold. Moreover, μi is Pearson’s correlation coefficient between the ith IMF and the original signal, and n is the total number of IMFs; maxμi is the maximum observed correlation coefficient. The selection criterion for IMFs is given as follows: if μi>λ, then keep the ith IMF, else eliminate the ith IMF, and add it into the residue.


### 2.2.6. EEMD-MFDFA-Based Method

The selection of the scale range for computing the fluctuation functionFqs and the type of polynomial chosen (detrending) is one of the major issues in applying the MFDFA method [10, 40]. Indeed, firstly, if there is a trend in a real-time series, the functional form of this trend is usually unknown. For real-time series, it is usually unknown whether there is a trend component and, if so, what the functional form of this trend is. A conventional strategy is to assume that local trends are in the form of a polynomial. The discontinuity at the segmentation points using polynomial fitting comes from the new pseudofluctuation errors and estimation errors in the qth-order generalized Hurst exponent hq. Secondly, the qth-order fluctuation function Fqs depends strongly on the selection of the multiple segment sizes s. To overcome the detrending issue in the standard MFDFA method, the EEMD algorithm is embedded into the MFDFA to modify the third step of the algorithms, while keeping all other steps unchanged as in [10, 45]. Thus, in the standard MFDFA method, the third step is reformulated as follows:Step 3: For each segmentν, one obtains the EEMD-based local trend yνi=rνi with the shifting process. Note that the trend rνi should be determined for each segment separately at each time scale. One can then obtain the variance as(12)F2ν,s=1s∑i=1sYν−1s+i−rνi2,forν=1,2,…,Ns,1s∑i=1sYN−ν−Nss+i−rνi2,forν=Ns+1,…,2Ns.We denote this MFDFA with the EEMD-based detrending method as EEMD-MFDFA. The criterion for selecting the range of scale (s) is the highest possible stability of the obtained spectra as in [9, 13, 26].According to Verrier et al. [53], zero values are susceptible to bias multifractal analysis results. To avoid the effect of the zeros on multifractal analysis results, the precipitation data recorded from April to October are selected to take into account the rainy season in all the Benin synoptic stations. To determine whether changes in the dynamics of precipitation time series can be assessed with the MFDFA method, the sixty-year data were divided into three separated subdatasets: the first is from 1951 to 1970, the second is from 1971 to 1990, and the third is from 1991 to 2010. This division was made considering the classification often mentioned by some studies in West Africa [54, 55]. The period 1951–1970 is considered as the wet period, 1971–1990 is the drought period, and 1991–2010 is considered as the end of the drought. In each synoptic station, the comparison of precipitation multifractal spectrum parameters values (α0,w,as) in the subsets is done.


## 2.2.1. Seasonal Detrending

Natural time series like precipitations generally exhibit a seasonal cycle; thus, there are usually seasonal variations in precipitation data. To make sure that seasonality does not affect the Multifractal Detrended Fluctuation Analysis (MFDFA) method, the seasonal detrending is done. In this study, a deseasonalized time serieszt is obtained from a procedure of the original series xt. The method is performed by adjusting the data with the seasonal mean and standard deviation as in [18, 34, 35]:(1)zt=xt−xtdσd,where xtd is the mean of the observed daily precipitation’s quantity calculated for each calendar date d (obtained by averaging over all the years in the record) and σd is the corresponding standard deviation (also calculated for each calendar date). The deseasonalized precipitation forms a new time series for MFDFA analysis.


## 2.2.2. Description of MFDFA Method

Based on the Detrended Fluctuation Analysis (DFA), the Multifractal Detrended Fluctuation Analysis was proposed by Kantelhardt et al. [27]. In literature, the MFDFA method is an effective tool to detect the scaling behaviors and multifractal properties of nonlinear and nonstationary time series such as precipitation records in hydrology and hydrometeorology [27, 28, 36]. Owing to its simplicity and robustness, MFDFA becomes very popular for the fractal characterization of daily, monthly, or annual precipitation records worldwide [9, 10, 13, 27, 28, 34, 35, 37–39].According to Kantelhardt et al. [27], Li et al. [39], Krzyszczak et al. [26], Adarsh et al. [40], and Baranowski et al. [13], the different steps involved in MFDFA computational procedure can be described as follows:Step 1: letzt|t=1,2,…,N be a deseasonalized precipitation time series of N equidistant measurements to which the procedure of the MFDFA method is applied. The accumulated deviation of the series (known as “profile”) is calculated as(2)Yk=∑t=1kzt−z,wherez is the mean value of the series, zt and k = 1, 2, …, N, and t = 1, 2, …, N.Step 2: divide the profileYk into Ns = int(N/s) equal-sized nonoverlapping windows with a length s. While considering multiple timescales, sometimes, a small portion of the time series at the end may be left out, as N need not be a multiple of s always. To retain this part of the series, the same procedure is repeated starting from the opposite end resulting in 2 Ns segments.Step 3: calculate the local trend for each of the 2Ns segments by a least-squares fit of the series as(3)F2ν,s=1s∑i=1sYν−1s+i−yνmi2,forν=1,2,…,Ns,1s∑i=1sYN−ν−Nss+i−yνmi2,forν=Ns+1,…,2Ns,whereyνi is the trend function in segment ν. yνi stands for the fitted polynomial trend in the vth segment. Different types of fittings such as linear (i.e., polynomial order m = 1), quadratic (i.e., m = 2), cubic (i.e., m = 3), or higher-order polynomials yνmi can be used to ﬁt the local trend. Determine the qth-order fluctuation function by averaging over all segments from the following equation:(4)Fqs=12Ns∑ν=12NsF2ν,sq/21/q,forq≠0,exp14Ns∑ν=12NslnF2ν,s,forq=0.Step 4: determine the scaling behavior of the fluctuation functions by analyzing thelog-log plots of Fqs versus s for each value of q (q is the order of moment). If the time series is long-range power-law correlated, Fqs increases as(5)Fqs∝shq,where h (q) is the generalized Hurst scaling function [27]. For monofractal time series, h (q) is independent of q, whereas, for a multifractal time series, h (q) depends on q. To avoid divergence of moments in the fat tails of the fluctuation distribution as mentioned by some authors [41, 42], one restricts the order q within the range −5 ≤ q ≤ 5. The density distributions of all the studied meteorological time series had heavy tails. To prevent and avoid potential distortion of the results by the so-called “freezing” phenomenon [9, 13, 28], the range of q has to be limited in the interval −5,5, as in [9, 13, 41].


## 2.2.3. Multifractal Spectra Analysis (MFS)

According to [9, 13, 27, 43–45], using the Legendre transform, the relationship between the multifractal spectrum fα and the generalized Hurst index hq can be obtained and written as(6)αq=hq+qdhqdq,fα=qαq−hq+1,where αq and fα are the singularity exponent (or Hölder exponent) and the multifractal spectrum, respectively. Additionally, the mass exponent function τ (q) can be calculated by the generalized Hurst exponent hq as(7)τq=qhq−1.The schematic representation of a multifractal spectrum with its most important parametersαmax, αmin, α0, as, and w is shown in Figure 3 (obtained from Baranowski et al. [9], Krzyszczak et al. [26], and Baranowski et al. [13]).Figure 3
Schematic presentation of the main parameters (position of maximumα0, asymmetry as, and width w) of a multifractal spectrum fα. Parameters α0, as, and w are dimensionless (modified figure from Baranowski et al., [13]).The parameterαmin represents the most extreme event and αmax the smoothest event in the studied process. The value of α0 delivers valuable information about the structure of the studied process. It indicates at which value of α multifractal spectrum reaches its maximum. A low value of α0 indicates that the underlying process becomes correlated and loses fine structure, becoming more regular in appearance.The negative values of the asymmetry parameteras (left-skewed shapes of the multifractal spectra) indicate low fractal exponents of small weights. They suggest the dominance of large fluctuations which can be connected with extreme events [9, 13]. A right skewness of the multifractal spectrum (positive value of as) denotes high fractal exponents with large weights, which are characteristic of fine structures and indicate a relative dominance of the small fluctuations. The width of the spectrum w (the difference between αmax and αmin) represents the length of the fractal exponent range in the series. It is strongly connected to the “richness” of the signal structure. The greater the value of w is, the more developed the multifractality is. In contrast, for pure monofractal, the width w of the spectra is equal to 0, but, according to Makowiec and Fuliński [46], if w is less than 0.05, then monofractal behavior of the spectrum should be assumed. Therefore, the values of the multifractal spectrum parameters (α0,w, and as) can be used as quantitative and qualitative indicators for studying the dynamics of meteorological processes [9, 13, 26].Additionally, according to Hou et al. [25], the asymmetric index (R), as presented in equation (8), is a useful parameter of multifractal analysis. It is computed based on the width of the left (ΔαL) and right (ΔαR) hand branches of the multifractal spectrum, respectively. Their values describe the distribution patterns of high and low fluctuations [47]. The values of R range from −1 to 1, quantifying the deviations of the multifractal spectrum curve [25, 48]. Positive values of R suggest a left-hand deviation of the multifractal spectrum, with high local fluctuations; negative values of R suggest a right-hand deviation with local low fluctuations; R-value equal to zero represents a symmetrical multifractal spectrum. The asymmetry index (R) is obtained as(8)R=ΔαL−ΔαRΔαL+ΔαR,where ΔαL=α0−αmin and ΔαR=αmax−α0.


## 2.2.4. Description of the EEMD Algorithm

The empirical mode decomposition (EMD) method proposed by Huang and Wu [49] is widely used for processing nonlinear and nonstationary signals in hydrology, hydrometeorology, and other fields [15, 49]. EMD method can decompose time series data into a limited number of intrinsic mode functions (IMFs) and a residual component (trend). The IMFs satisfy the following two main requirements [15]: (i) in the whole set of data, the number of local extrema and the numbers of zero-crossings must be equal to or different from 1 at most and (ii) at any point, the mean value of the ‘‘upper envelope’’ (defined by the local maxima) and the ‘‘lower envelope’’ (defined by the local minima) must be zero. According to these properties, some meaningful IMFs can be well defined. Generally, an IMF indicates a simple oscillatory mode, compared to a simple harmonic function. Based on this definition, a shifting process of the original time series can be briefly expressed by Liu et al. [15] in steps as follows:(1)
Identify all extrema (local maxima and minima) points of the given time seriesy (t)(2)
Connect that extrema by spline interpolation, to create an upper envelope of local maxima pointsemaxt, and a lower envelope emint, of all minima points(3)
Compute the meanm (t) of two envelopes by m=emaxt+emint/2(4)
Subtract the mean from the data to establish an IMF candidate;ht=yt−mt(5)
Ifht satisfies the two properties of IMF as indicated by a predefined stopping criterion, ht is indicated as the first IMF and is recorded as IMFit; if ht is not an IMF, yt is replaced by ht, and iterate steps 1–4 until h (t) satisfies the two conditions of IMFs(6)
The residuerit=yt−IMFit is then treated as new data subjected to the same shifting procedure as defined above for the next IMF from rit. In the final step, the shifting process can be stopped when the residue rt has a monotonic trend or has one local maxima and minima point from which no more IMFs can be extracted. At the end of this shifting process, the original signal yt can be reconstructed as the sum of IMFs and residual as(9)yt=∑i=1nIMFit+rnt,where rnt represents the final residue and n is the number of IMFs. The residue rnt is known as the trend of the signal yt [15].Although EMD has obvious advantages in time series data processing, there are also some unavoidable weaknesses [14, 15]. The major weakness is mode-mixing. Mode-mixing means that the same IMF contains different frequency components or the frequency of the same and similar scales is distributed in different IMFs. Thus, mode-mixing not only causes the mixing of various scale vibration modes but can even change the physical meaning of the individual IMF. To overcome these drawbacks of EMD, Ensemble Empirical Mode Decomposition (EEMD) is proposed in Wu and Huang [14] using white noises. The EEMD algorithm is straight forward and can be summarized as follows:(1)
Add the white noise series into the original series(2)
Decompose the signal with the added white noise into the IMFs using EMD(3)
Repeat steps (1) and (2) with a different white noise series each time(4)
Obtain the corresponding IMF components of the decomposition and calculate the means of the ensemble corresponding to the IMFs of the decomposition as the final result [14]For the EEMD method, the first important step is to determine the ensemble times and the amplitude of the added noise. But the way of selecting the best ensemble number and amplitude of added noise is still an open question. The effect of adding white noise should obey the following statistical rule:(10)εn=εN,where εn is the final standard deviation of the error, which is described as the difference between the input signal and the corresponding IMFs. ε denotes the amplitude of the added noise, and N is the number of ensemble numbers. In this study, the ensemble number was set to 100, and the amplitude of added white noise was set to 0.2 times the standard deviation of the corresponding data as in [10, 15].


## 2.2.5. Period Computation Method and Criterion for Selecting Relevant IMFs

Since each IMF represents different intrinsic modes of oscillation, it is possible to calculate the period for each of them. The average periods are calculated in this study by the time intervals between consecutive zero-crossings on successive waves as in [50].Since the IMFs are supposed to be almost orthogonal components of the original signal, each IMF would have a relatively good correlation with the original signal; this presupposes that the irrelevant components would have a relatively poor correlation with the original signal [51]. For discriminating between relevant and irrelevant IMFs, the following threshold has been proposed by Ayenu-Prah and Attoh-Okine [52].(11)λ=maxμi10∗maxμi−3,where i=1,2,…,n and λ is the threshold. Moreover, μi is Pearson’s correlation coefficient between the ith IMF and the original signal, and n is the total number of IMFs; maxμi is the maximum observed correlation coefficient. The selection criterion for IMFs is given as follows: if μi>λ, then keep the ith IMF, else eliminate the ith IMF, and add it into the residue.


## 2.2.6. EEMD-MFDFA-Based Method

The selection of the scale range for computing the fluctuation functionFqs and the type of polynomial chosen (detrending) is one of the major issues in applying the MFDFA method [10, 40]. Indeed, firstly, if there is a trend in a real-time series, the functional form of this trend is usually unknown. For real-time series, it is usually unknown whether there is a trend component and, if so, what the functional form of this trend is. A conventional strategy is to assume that local trends are in the form of a polynomial. The discontinuity at the segmentation points using polynomial fitting comes from the new pseudofluctuation errors and estimation errors in the qth-order generalized Hurst exponent hq. Secondly, the qth-order fluctuation function Fqs depends strongly on the selection of the multiple segment sizes s. To overcome the detrending issue in the standard MFDFA method, the EEMD algorithm is embedded into the MFDFA to modify the third step of the algorithms, while keeping all other steps unchanged as in [10, 45]. Thus, in the standard MFDFA method, the third step is reformulated as follows:Step 3: For each segmentν, one obtains the EEMD-based local trend yνi=rνi with the shifting process. Note that the trend rνi should be determined for each segment separately at each time scale. One can then obtain the variance as(12)F2ν,s=1s∑i=1sYν−1s+i−rνi2,forν=1,2,…,Ns,1s∑i=1sYN−ν−Nss+i−rνi2,forν=Ns+1,…,2Ns.We denote this MFDFA with the EEMD-based detrending method as EEMD-MFDFA. The criterion for selecting the range of scale (s) is the highest possible stability of the obtained spectra as in [9, 13, 26].According to Verrier et al. [53], zero values are susceptible to bias multifractal analysis results. To avoid the effect of the zeros on multifractal analysis results, the precipitation data recorded from April to October are selected to take into account the rainy season in all the Benin synoptic stations. To determine whether changes in the dynamics of precipitation time series can be assessed with the MFDFA method, the sixty-year data were divided into three separated subdatasets: the first is from 1951 to 1970, the second is from 1971 to 1990, and the third is from 1991 to 2010. This division was made considering the classification often mentioned by some studies in West Africa [54, 55]. The period 1951–1970 is considered as the wet period, 1971–1990 is the drought period, and 1991–2010 is considered as the end of the drought. In each synoptic station, the comparison of precipitation multifractal spectrum parameters values (α0,w,as) in the subsets is done.


## 3. Results and Discussion

### 3.1. Multiscale Decomposition of Observed Precipitation in Benin Synoptic Stations

Precipitation observed in Benin synoptic stations is decomposed into different intrinsic mode functions (IMFs) and residue representing different scales. Thirteen (13) IMFs (not shown) are obtained in the subequatorial region and twelve (12) IMFs (not shown) in the Sudanian region (except Kandi station). The IMFs with lower numerical numbers correspond to higher-frequency (smaller scale) oscillations, whereas IMFs with higher numerical numbers correspond to lower-frequency oscillations at larger scales.In Figure4, the curve of residue r (t) obtained during 1951–2010 at Cotonou synoptic station is presented. The curve of the residue obtained at the other studied synoptic stations is not shown. These results indicate that, in Benin, the EEMD method can be used for decomposing precipitation data into different IMFs with a residue. Since residue r (t) represents the trend of the original data, this curve displays the precipitation trend in the study period. A decreasing trend occurred in the whole synoptic stations during 1951–2010. This result is not systematically in agreement with those of other studies [54, 55]. Indeed, based on the previous findings in West Africa, contrasting the existence of a general trend in the study period, 1951–1970 is considered as the wet period, 1971–1990 is the drought period, and 1991–2010 is the end of the drought [55]. Considering the EEMD approach, our results show that the subperiod 1991–2010 is not a transition period as mentioned by these authors. The drought is prolonged until 2010. Generally, in West Africa, linear methods (linear regression and Mann–Kendall test) are used in many types of research to study precipitation trends. But, according to Kulkarni and von Storch [56] and Yue et al. [57, 58], the results of linear methods in detecting and interpreting trends in hydrologic time series are affected by the presence of serial correlation in time series. They are, in general, insufficient and unreliable for a complete analysis of meteorological time series (e.g., precipitations). Our findings reveal questionable statements concerning the reliability of the approach used in previous studies.Figure 4
Residues (trend curve) of Cotonou synoptic station based on EEMD in 1951–2010.The dominant components in terms of a good correlation between the IMFs and the original signal are identified. Pearson’s correlation coefficient (PCC) between theith IMF and the original signal at Cotonou synoptic station is shown in Figures 5. The results obtained at the other studied synoptic stations are not shown. It is observed that whatever the synoptic stations are, the first five obtained IMFs (IMF1, IMF2, IMF3, IMF4, and IMF5) reveal high PCC values than the threshold value. Thus, according to the selection criterion, these IMFs are the relevant IMFs in all the stations. These results mean that, in Benin, precipitation is mainly controlled and governed by these five IMFs, influencing significantly the precipitation variations.Figure 5
Comparison of Pearson’s correlation coefficient values between IMFs and original data and the threshold values at Cotonou synoptic station.Table2 presents the average periods calculated for each of the IMFs obtained during the period 1951–2010 at each of the synoptic stations. The values of the period (in the day) are highlighted in bold.Table 2
Periods (in days) of decomposed data series with EEMD.
CotonouBohiconSaveParakouNatitingouKandiIMF12.72.62.62.62.62.6IMF25.35.55.55.45.05.03IMF38.68.58.68.57. 98.1IMF413.113.913.813.313.413.1IMF521.423.222.724.524.021.2IMF635.339.640.242. 443.436.3IMF763.081.078.7101.7103.486.6IMF8107.1142.4156.6210. 2213.6213.4IMF9213.9291.0288220.7226.3214.1IMF10366.4532.8647.4651.3737.9480.3IMF11730.1986.2986.21325. 61221899.4IMF121288.22094.82212.2304330652243.6IMF133301.740624527.5——6759In Table2, it is observed that the decomposed data presents nearly identical periods for the lower IMFs such as IMF1 to IMF5 in all synoptic stations. The periods of the relevant IMFs are between 1 and 25 days. The higher IMFs, which represent the longer periodic oscillations of the original data, have larger differences because there are fewer cycles to average overfluctuations in instantaneous frequencies, as we can find in IMF10,11,12,13.


### 3.2. Multifractal Spectrum Analysis

Figure6 shows the variation of scaling exponent function τ (q) with q for the daily precipitation time series at Cotonou synoptic station during the subperiods 1951–1970, 1971–1990, and 1991–2010. The variation of τ (q) with q obtained for the other studied synoptic stations is not shown. The function τ (q) increased nonlinearly with the increasing values of q, indicating significant multifractal characteristics of the daily precipitation time series in Benin, whatever the synoptic stations and the period considered. Moreover, this nonlinear dependence of τ (q) function reveals the presence of nonlinear interaction between the different scale events and the multifractal nature of the precipitation time series. So, in the synoptic stations, the function τ (q) is nonlinear, illustrating the presence of multiple scaling in precipitation observations during the study periods. The degree of the nonlinearity of the function τ (q) reflects the degree of multifractality [59]. Furthermore, whatever the synoptic stations and the studied periods, different relationships between τ (q) and q for −5 < q < 0 and 0 < q < 5 are found as in Figure 6 for Cotonou synoptic station. Therefore, the scaling exponent function’s behavior is different for q < 0 (for q > 0) and possesses different slopes before and after the value for q = 0 (the values of the slopes are not shown).Figure 6
Plot of scaling exponent with moment order (τ(q) versus (q) of Cotonou synoptic station precipitation time series: 1951–1970 (black color), 1971–1990 (red color), and 1991–2010 (green color).The relationship between singularity spectrumf (α) and singularity exponent α of Cotonou synoptic station' precipitation in the subsets is displayed in Figure 7. Similar results are obtained for the other synoptic stations (not shown). From the figure, one can notice that whatever the synoptic stations and the subperiods considered, all the multifractal spectra f (α) exhibit the shape of a parabolic curve, indicating the multifractal structure of the precipitation time series but with varying degree, since the curves are all different. This conforms to the results obtained by Agbazo et al. [11] in the same study stations. The analysis of the multifractal spectra obtained at all synoptic stations reveals evident differences in the dynamics of the precipitations processes for Benin. The multifractal spectra of the precipitation quantities vary regarding the climate region of the synoptic stations and the considered subperiods. The structures of the multifractal spectrum vary relatively to the synoptic stations and the study period, indicating a different evolution law. This indicates also that the long-term characteristics of precipitation in different subdivisions are not stable. Furthermore, the multifractal spectra at all the synoptic stations in the different periods are not symmetric, and all of the stations have left truncation. On examining the plots, it can be noted that the multifractal spectra have a left truncation and a long right tail, which can be attributed to the sensibility of the time series to the local climate fluctuations with small magnitudes [42]. Therefore, multifractal spectra analysis revealed that the effect of small fluctuations plays a dominant role in daily precipitation time series.Figure 7
Multifractal spectra of precipitation recorded at Cotonou synoptic station, plotted for 1951–1970 (black color), 1971–1990 (red color), and 1991–2010 (green color).The result of the asymmetry index values (R) obtained for different synoptic stations and periods is presented in Figure 8(a)). The findings show that R is systematically negative for all the synoptic stations regardless of the periods. According to Hou et al. [25] and Adarsh et al. [40], the negative values of R from a spectrum with right-hand deviation indicates that extreme events play a prominent role in the temporal structure of daily precipitation. That is the case for our dataset. Whatever the period, the multifractal spectra of precipitation in synoptic stations are right-deviant with local low fluctuations, meaning that the singularity of the low values is larger than the high values. Examining the plots, it can be noted that, for the synoptic stations located in the Sudanian region R1991-2010 < R1971-1990 < R1951-1970, whereas in the subequatorial region the variation on R is random.Figure 8
Spatiotemporal variation of the asymmetry index values (R) and the multifractal spectrum parameters (as,α0, and w) of precipitation data. Co, Bo, Sa, Pa, Na, and Ka mean Cotonou, Bohicon, Save, Parakou, Natitingou, and Kandi synoptic stations, respectively. The spatiotemporal variation of (R) as,α0, and w are shown at panels (a), (b), (c), and (d), respectively.
(a)(b)(c)(d)The multifractal spectrum parameters (as,α0, and w) of precipitation data from the synoptic stations in Benin in the different periods are presented in Figures 8(b), 8(c), and 8(d), respectively. The width of the multifractal spectrum w is used to measure the degree of multifractality, which denoted a nonuniform clustering structure for daily precipitation. The analysis of changes in multifractal properties shows that there is a clear difference between the subperiods in multifractal properties. This result shows a change in the dynamics of precipitation time series, namely, 1951–1970, 1971–1990, and 1991–2010 as shown in the literature (Balme et al. [54] and Nicholson [55]). Climatic modifications are observed between 1951–1970 and 1971–1990, and between 1971–1990 and 1991–2010.In the case of the asymmetry parameteras (Figure 8(b)), more consistent changes are observed. The asymmetry parameters are positive in all the synoptic stations regardless of the study period, indicating that the fine structure of the precipitation is preserved during all the subperiods. For the synoptic stations located in the Sudanian region, there is a clear increase in as values for precipitation series from 1951 to 1970, 1971 to 1990, and 1991 to 2010. For the synoptic stations located in the Sudanian region, one can note systematically that as1951−1970<as1971−1990<as1991−2010. However, no explicit direction of changes is observed in subequatorial stations. In the Sudanian region, as values change from a positive value to a higher positive from 1951 to 1970, 1971 to 1990, and 1991 to 2010. Specifically, at Parakou, Natitingou, and Kandi, as vary from 0.3 to 0.5, 0.4 to 0.9, and 0.5 to 0.8, respectively. In the Sudanian region, by comparing 1951–1970 and 1971–1990, one can note that as1951−1970<as1971−1990; therefore, more extreme events are observed in the dynamics of the precipitation processes during the 1951–1970 periods than in 1971–1990, where as as1971−1990<as1991−2010; thus, the 1971–1990 period had recorded more extreme events than 1991–2010.At Bohicon and Save stations (subequatorial region), the lowest positive value ofas is obtained from 1971 to 1990; thus, comparing the three subperiods, 1971–1990 is the period where more extreme events happened. At Cotonou station (subequatorial region), as changes from higher positive values to lower positive ones, indicating that the fine structure of the signal is preserved, but in the 1991–2010 period, more extreme events happened. The asymmetry as values for the subequatorial region changed differently compared to the Sudanian region. This result reveals evident differences in the dynamics of the precipitation processes in the southern part of Benin. This difference could be explained by the proximity of the subequatorial region to the Atlantic Ocean, which strongly influences the precipitation regime in the subequatorial region than Harmattan fluxes in the Sudanian area.In the case of the multifractal spectrum widthw (also known as the richness of the signal) presented in Figure 8(c)), its values changed during the studied periods. Therefore, multifractal properties could be regarded as good indicators of changes in the dynamics of precipitation.One can notice that the greaterw is, the larger the variety of daily precipitation is. Also, the stronger the degree of multifractality is obtained, the more complex the structure of daily precipitation is. Generally, except for Save and Parakou stations, the greatest value of the multifractal spectrum width is obtained from 1951 to 1970. This result indicates that the distribution of daily precipitation during 1951–1970 at the oversynoptic stations, except Save and Parakou, is relatively uneven, and the multifractal strength of daily precipitation is greater than the other studied periods. At Save and Parakou stations, the greatest value of the multifractal spectrum width is obtained during 1971–1990. Therefore, in this subperiod, the degree of multifractality is stronger and the structure of daily precipitation is more complex.The spatial and temporal variability ofα0 is provided in Figure 8(d). Generally, the highest values of α0 are obtained from 1951 to 1970 at all the stations, indicating that the precipitation process is less correlated and possesses fine structure in the studied period. For all the stations located in the subequatorial region and at Kandi (Sudanian region), the lowest value of α0 is obtained during 1991–2010, indicating that the underlying process becomes correlated and loses the fine structure, becoming more regular in appearance, whereas at Parakou and Natitingou synoptic stations (in the Sudanian region), the lowest value is observed in 1971–1990.Referring to the values ofα0 or the width of the spectrum w, generally consistent changes are observed. According to Baranowski et al. [9], this result implies that the impact on these two features (α0 and w) is more from the global changes in climate dynamics than in the local changes (climatic zones). The results of precipitation multifractal parameters analysis indicate that the spatial and temporal variability of α0 Figure 8(d) differs from the asymmetry as and the width w parameters (Figures 8(b) and 8(c)). Therefore, multifractal properties could be regarded as good indicators to assess changes in the dynamics of precipitation. Despite the observed differences in multifractal spectra properties in the three studied periods, it is not possible to identify similarities in the direction of changes for all the stations at the same time, maybe because of the climate difference.


## 3.1. Multiscale Decomposition of Observed Precipitation in Benin Synoptic Stations

Precipitation observed in Benin synoptic stations is decomposed into different intrinsic mode functions (IMFs) and residue representing different scales. Thirteen (13) IMFs (not shown) are obtained in the subequatorial region and twelve (12) IMFs (not shown) in the Sudanian region (except Kandi station). The IMFs with lower numerical numbers correspond to higher-frequency (smaller scale) oscillations, whereas IMFs with higher numerical numbers correspond to lower-frequency oscillations at larger scales.In Figure4, the curve of residue r (t) obtained during 1951–2010 at Cotonou synoptic station is presented. The curve of the residue obtained at the other studied synoptic stations is not shown. These results indicate that, in Benin, the EEMD method can be used for decomposing precipitation data into different IMFs with a residue. Since residue r (t) represents the trend of the original data, this curve displays the precipitation trend in the study period. A decreasing trend occurred in the whole synoptic stations during 1951–2010. This result is not systematically in agreement with those of other studies [54, 55]. Indeed, based on the previous findings in West Africa, contrasting the existence of a general trend in the study period, 1951–1970 is considered as the wet period, 1971–1990 is the drought period, and 1991–2010 is the end of the drought [55]. Considering the EEMD approach, our results show that the subperiod 1991–2010 is not a transition period as mentioned by these authors. The drought is prolonged until 2010. Generally, in West Africa, linear methods (linear regression and Mann–Kendall test) are used in many types of research to study precipitation trends. But, according to Kulkarni and von Storch [56] and Yue et al. [57, 58], the results of linear methods in detecting and interpreting trends in hydrologic time series are affected by the presence of serial correlation in time series. They are, in general, insufficient and unreliable for a complete analysis of meteorological time series (e.g., precipitations). Our findings reveal questionable statements concerning the reliability of the approach used in previous studies.Figure 4
Residues (trend curve) of Cotonou synoptic station based on EEMD in 1951–2010.The dominant components in terms of a good correlation between the IMFs and the original signal are identified. Pearson’s correlation coefficient (PCC) between theith IMF and the original signal at Cotonou synoptic station is shown in Figures 5. The results obtained at the other studied synoptic stations are not shown. It is observed that whatever the synoptic stations are, the first five obtained IMFs (IMF1, IMF2, IMF3, IMF4, and IMF5) reveal high PCC values than the threshold value. Thus, according to the selection criterion, these IMFs are the relevant IMFs in all the stations. These results mean that, in Benin, precipitation is mainly controlled and governed by these five IMFs, influencing significantly the precipitation variations.Figure 5
Comparison of Pearson’s correlation coefficient values between IMFs and original data and the threshold values at Cotonou synoptic station.Table2 presents the average periods calculated for each of the IMFs obtained during the period 1951–2010 at each of the synoptic stations. The values of the period (in the day) are highlighted in bold.Table 2
Periods (in days) of decomposed data series with EEMD.
CotonouBohiconSaveParakouNatitingouKandiIMF12.72.62.62.62.62.6IMF25.35.55.55.45.05.03IMF38.68.58.68.57. 98.1IMF413.113.913.813.313.413.1IMF521.423.222.724.524.021.2IMF635.339.640.242. 443.436.3IMF763.081.078.7101.7103.486.6IMF8107.1142.4156.6210. 2213.6213.4IMF9213.9291.0288220.7226.3214.1IMF10366.4532.8647.4651.3737.9480.3IMF11730.1986.2986.21325. 61221899.4IMF121288.22094.82212.2304330652243.6IMF133301.740624527.5——6759In Table2, it is observed that the decomposed data presents nearly identical periods for the lower IMFs such as IMF1 to IMF5 in all synoptic stations. The periods of the relevant IMFs are between 1 and 25 days. The higher IMFs, which represent the longer periodic oscillations of the original data, have larger differences because there are fewer cycles to average overfluctuations in instantaneous frequencies, as we can find in IMF10,11,12,13.


## 3.2. Multifractal Spectrum Analysis

Figure6 shows the variation of scaling exponent function τ (q) with q for the daily precipitation time series at Cotonou synoptic station during the subperiods 1951–1970, 1971–1990, and 1991–2010. The variation of τ (q) with q obtained for the other studied synoptic stations is not shown. The function τ (q) increased nonlinearly with the increasing values of q, indicating significant multifractal characteristics of the daily precipitation time series in Benin, whatever the synoptic stations and the period considered. Moreover, this nonlinear dependence of τ (q) function reveals the presence of nonlinear interaction between the different scale events and the multifractal nature of the precipitation time series. So, in the synoptic stations, the function τ (q) is nonlinear, illustrating the presence of multiple scaling in precipitation observations during the study periods. The degree of the nonlinearity of the function τ (q) reflects the degree of multifractality [59]. Furthermore, whatever the synoptic stations and the studied periods, different relationships between τ (q) and q for −5 < q < 0 and 0 < q < 5 are found as in Figure 6 for Cotonou synoptic station. Therefore, the scaling exponent function’s behavior is different for q < 0 (for q > 0) and possesses different slopes before and after the value for q = 0 (the values of the slopes are not shown).Figure 6
Plot of scaling exponent with moment order (τ(q) versus (q) of Cotonou synoptic station precipitation time series: 1951–1970 (black color), 1971–1990 (red color), and 1991–2010 (green color).The relationship between singularity spectrumf (α) and singularity exponent α of Cotonou synoptic station' precipitation in the subsets is displayed in Figure 7. Similar results are obtained for the other synoptic stations (not shown). From the figure, one can notice that whatever the synoptic stations and the subperiods considered, all the multifractal spectra f (α) exhibit the shape of a parabolic curve, indicating the multifractal structure of the precipitation time series but with varying degree, since the curves are all different. This conforms to the results obtained by Agbazo et al. [11] in the same study stations. The analysis of the multifractal spectra obtained at all synoptic stations reveals evident differences in the dynamics of the precipitations processes for Benin. The multifractal spectra of the precipitation quantities vary regarding the climate region of the synoptic stations and the considered subperiods. The structures of the multifractal spectrum vary relatively to the synoptic stations and the study period, indicating a different evolution law. This indicates also that the long-term characteristics of precipitation in different subdivisions are not stable. Furthermore, the multifractal spectra at all the synoptic stations in the different periods are not symmetric, and all of the stations have left truncation. On examining the plots, it can be noted that the multifractal spectra have a left truncation and a long right tail, which can be attributed to the sensibility of the time series to the local climate fluctuations with small magnitudes [42]. Therefore, multifractal spectra analysis revealed that the effect of small fluctuations plays a dominant role in daily precipitation time series.Figure 7
Multifractal spectra of precipitation recorded at Cotonou synoptic station, plotted for 1951–1970 (black color), 1971–1990 (red color), and 1991–2010 (green color).The result of the asymmetry index values (R) obtained for different synoptic stations and periods is presented in Figure 8(a)). The findings show that R is systematically negative for all the synoptic stations regardless of the periods. According to Hou et al. [25] and Adarsh et al. [40], the negative values of R from a spectrum with right-hand deviation indicates that extreme events play a prominent role in the temporal structure of daily precipitation. That is the case for our dataset. Whatever the period, the multifractal spectra of precipitation in synoptic stations are right-deviant with local low fluctuations, meaning that the singularity of the low values is larger than the high values. Examining the plots, it can be noted that, for the synoptic stations located in the Sudanian region R1991-2010 < R1971-1990 < R1951-1970, whereas in the subequatorial region the variation on R is random.Figure 8
Spatiotemporal variation of the asymmetry index values (R) and the multifractal spectrum parameters (as,α0, and w) of precipitation data. Co, Bo, Sa, Pa, Na, and Ka mean Cotonou, Bohicon, Save, Parakou, Natitingou, and Kandi synoptic stations, respectively. The spatiotemporal variation of (R) as,α0, and w are shown at panels (a), (b), (c), and (d), respectively.
(a)(b)(c)(d)The multifractal spectrum parameters (as,α0, and w) of precipitation data from the synoptic stations in Benin in the different periods are presented in Figures 8(b), 8(c), and 8(d), respectively. The width of the multifractal spectrum w is used to measure the degree of multifractality, which denoted a nonuniform clustering structure for daily precipitation. The analysis of changes in multifractal properties shows that there is a clear difference between the subperiods in multifractal properties. This result shows a change in the dynamics of precipitation time series, namely, 1951–1970, 1971–1990, and 1991–2010 as shown in the literature (Balme et al. [54] and Nicholson [55]). Climatic modifications are observed between 1951–1970 and 1971–1990, and between 1971–1990 and 1991–2010.In the case of the asymmetry parameteras (Figure 8(b)), more consistent changes are observed. The asymmetry parameters are positive in all the synoptic stations regardless of the study period, indicating that the fine structure of the precipitation is preserved during all the subperiods. For the synoptic stations located in the Sudanian region, there is a clear increase in as values for precipitation series from 1951 to 1970, 1971 to 1990, and 1991 to 2010. For the synoptic stations located in the Sudanian region, one can note systematically that as1951−1970<as1971−1990<as1991−2010. However, no explicit direction of changes is observed in subequatorial stations. In the Sudanian region, as values change from a positive value to a higher positive from 1951 to 1970, 1971 to 1990, and 1991 to 2010. Specifically, at Parakou, Natitingou, and Kandi, as vary from 0.3 to 0.5, 0.4 to 0.9, and 0.5 to 0.8, respectively. In the Sudanian region, by comparing 1951–1970 and 1971–1990, one can note that as1951−1970<as1971−1990; therefore, more extreme events are observed in the dynamics of the precipitation processes during the 1951–1970 periods than in 1971–1990, where as as1971−1990<as1991−2010; thus, the 1971–1990 period had recorded more extreme events than 1991–2010.At Bohicon and Save stations (subequatorial region), the lowest positive value ofas is obtained from 1971 to 1990; thus, comparing the three subperiods, 1971–1990 is the period where more extreme events happened. At Cotonou station (subequatorial region), as changes from higher positive values to lower positive ones, indicating that the fine structure of the signal is preserved, but in the 1991–2010 period, more extreme events happened. The asymmetry as values for the subequatorial region changed differently compared to the Sudanian region. This result reveals evident differences in the dynamics of the precipitation processes in the southern part of Benin. This difference could be explained by the proximity of the subequatorial region to the Atlantic Ocean, which strongly influences the precipitation regime in the subequatorial region than Harmattan fluxes in the Sudanian area.In the case of the multifractal spectrum widthw (also known as the richness of the signal) presented in Figure 8(c)), its values changed during the studied periods. Therefore, multifractal properties could be regarded as good indicators of changes in the dynamics of precipitation.One can notice that the greaterw is, the larger the variety of daily precipitation is. Also, the stronger the degree of multifractality is obtained, the more complex the structure of daily precipitation is. Generally, except for Save and Parakou stations, the greatest value of the multifractal spectrum width is obtained from 1951 to 1970. This result indicates that the distribution of daily precipitation during 1951–1970 at the oversynoptic stations, except Save and Parakou, is relatively uneven, and the multifractal strength of daily precipitation is greater than the other studied periods. At Save and Parakou stations, the greatest value of the multifractal spectrum width is obtained during 1971–1990. Therefore, in this subperiod, the degree of multifractality is stronger and the structure of daily precipitation is more complex.The spatial and temporal variability ofα0 is provided in Figure 8(d). Generally, the highest values of α0 are obtained from 1951 to 1970 at all the stations, indicating that the precipitation process is less correlated and possesses fine structure in the studied period. For all the stations located in the subequatorial region and at Kandi (Sudanian region), the lowest value of α0 is obtained during 1991–2010, indicating that the underlying process becomes correlated and loses the fine structure, becoming more regular in appearance, whereas at Parakou and Natitingou synoptic stations (in the Sudanian region), the lowest value is observed in 1971–1990.Referring to the values ofα0 or the width of the spectrum w, generally consistent changes are observed. According to Baranowski et al. [9], this result implies that the impact on these two features (α0 and w) is more from the global changes in climate dynamics than in the local changes (climatic zones). The results of precipitation multifractal parameters analysis indicate that the spatial and temporal variability of α0 Figure 8(d) differs from the asymmetry as and the width w parameters (Figures 8(b) and 8(c)). Therefore, multifractal properties could be regarded as good indicators to assess changes in the dynamics of precipitation. Despite the observed differences in multifractal spectra properties in the three studied periods, it is not possible to identify similarities in the direction of changes for all the stations at the same time, maybe because of the climate difference.


## 4. Conclusions and Remarks

In the context of global climate change, it is of great importance to accurately determine nonlinear changes in climate factors for understanding the complex change processes of the climate system. In this study, Ensemble Empirical Mode Decomposition (EEMD) is used for decomposing the precipitation time series data recorded at the synoptic stations subdivided into two main climatic areas during 1951–2010. EEMD decomposed Benin precipitation into a finite and low number of oscillatory modes, depending on the local characteristic timescale. The oscillatory modes are revealed by intrinsic mode function (IMFs) components, embedded in the data. Moreover, the Multifractal Detrended Fluctuation Analysis (MFDFA) with EEMD-based detrending is used to investigate the multifractal properties. The major remarks are summarized as follows:(1)
In the subequatorial region, thirteen intrinsic mode functions (IMFs) are obtained, whereas, in the Sudanian region, except Kandi station, twelve IMFs are obtained. The decrease in the obtained residue indicates the existence of a decrease in the precipitation trend, contrasting the known subdivision of the study area from the literature(2)
In Benin synoptic stations, precipitation is mainly controlled and governed by the first five main IMFs (IMF1, IMF2, IMF3, IMF4, and IMF5) in which periods are between 1 and 25 days, corresponding to the local physical phenomenon(3)
The nonlinearity of the mass exponent functionτ (q) indicates the presence of multiple scaling in precipitation in all the stations regardless of the study period(4)
Multifractal spectra, whatever the synoptic stations and the study periods, have a left truncation and a long right tail attributed to the sensibility of the time series to the local fluctuations, confirming the results in point 2. Therefore, multifractal spectra analysis revealed that the effect of large fluctuations plays a dominant role in the daily precipitation series(5)
The analysis of the multifractal spectrum parameters (α0, w, and as) shows that there is a clear difference between the analyzed periods in multifractal spectrum parameters. The changes in the results for different subperiods reveal that these subperiods have different climatic modifications. These results show that multifractal analysis is a good method for assessing the differences in the dynamics of precipitation processes of the Benin Republic. Despite the observed differences in multifractal spectra properties in the subperiods, it is not possible to identify similarities in the direction of changes for all stationsOur findings can be used for the validation of global and regional climate models because a valid model should explain empirically detected scaling properties in observed data. This can be performed in future studies.



---

*Source: 3689397-2021-01-22.xml*