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

International Chair in Mathematical Physics and Applications (ICMPA-UNESCO Chair), Université d’Abomey-Calavi, Abomey-Calavi BP: 526 UAC, Benin Laboratoire de Physique du Rayonnement (LPR), Université d’Abomey-Calavi, Abomey-Calavi BP: 526 UAC, Benin Laboratoire de Géoscience, de l’Environnement et Applications (LaGEA), Université Nationale des Sciences Technologies, Ingénierie et Mathématiques, Abomey, Benin Laboratoire d’Hydrologie Appliquée (LHA), Université d’Abomey-Calavi, Abomey-Calavi BP: 526 UAC, Benin


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]. us, 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). ese methods are, in general, insufficient and unreliable for a complete analysis of meteorological time series like precipitations [9,[11][12][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][17][18]. Shen et al. [19] have demonstrated that climate change has selfmemory 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][21][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? e 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. e remaining parts of the paper are organized as follows: in Section 2, 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.

Study Sites.
e study covers all the synoptic stations of Benin, the geographical positions of which are presented in Figure 1. e 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]. e seasonal variations of the ICTZ, as well as the distribution of the precipitations, depend on the regions. ey 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). e daily and annual thermal amplitudes are, respectively, ∼10°C and∼5°C. e 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. e 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. e 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. e "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.  Table 1.

Complexity 3
synoptic stations which indicates that the distributions are right-tailed. e 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. Figure 2 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. is figure evokes the object obtained after a certain number of iterations in the random process that generates Cantor's dust. e figure also reveals the random and nonlinear nature of the rainfall process in the measurement site.

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 series z(t) is obtained from a procedure of the original series x(t). e method is performed by adjusting the data with the seasonal mean and standard deviation as in [18,34,35]:  [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][38][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: let z t |t � 1, 2, . . . , N be a deseasonalized precipitation time series of N equidistant measurements to which the procedure of the MFDFA method is applied.
e accumulated deviation of the series (known as "profile") is calculated as where 〈z〉 is the mean value of the series, z t and k � 1, 2, . . ., N, and t � 1, 2, . . ., N.
Step 2: divide the profile Y(k) into Ns � int(N/s) equalsized 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 N s segments.
Step 3: calculate the local trend for each of the 2 N s segments by a least-squares fit of the series as where y ] (i) is the trend function in segment ]. y ] (i) stands for the fitted polynomial trend in the v th 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 (m) ] (i) can be used to fit the local trend. Determine the q th -order Step 4: determine the scaling behavior of the fluctuation functions by analyzing the log-log plots of F q (s) versus s for each value of q (q is the order of moment). If the time series is long-range power-law correlated, F q (s) increases as 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. e 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].

Multifractal Spectra Analysis (MFS).
According to [9,13,27,[43][44][45], using the Legendre transform, the relationship between the multifractal spectrum f(α) and the generalized Hurst index h(q) can be obtained and written as where α (q) and f(α) are the singularity exponent (or Hölder exponent) and the multifractal spectrum ,   50  100  150  200  250  300  350  1951  1952  1953  1954  1955  1956  1957  1958  1959  1960  1961  1962  1963  1964  1965  1966  1967  1968  1969  1970  1971  1972  1973  1974  1975  1976  1977  1978  1979  1980  1981  1982  1983  1984  1985  1986  1987  1988  1989  1990  1991  1992  1993  1994  1995  1996  1997  1998  1999  2000  2001  2002  2003  2004  2005  2006  2007  2008  2009   Each row corresponds to a year of observations. e occurrence of rain is defined by a precipitation threshold that distinguishes between dry and rainy days. e precipitation threshold chosen in this representation is zero. Complexity 5 respectively. Additionally, the mass exponent function τ (q) can be calculated by the generalized Hurst exponent h(q) as e schematic representation of a multifractal spectrum with its most important parameters α max , α min , α 0 , a s , and w is shown in Figure 3 (obtained from Baranowski et al. [9], Krzyszczak et al. [26], and Baranowski et al. [13]). e parameter α min represents the most extreme event and α max the smoothest event in the studied process. e 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.
e negative values of the asymmetry parameter a s (leftskewed shapes of the multifractal spectra) indicate low fractal exponents of small weights. ey 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 a s ) denotes high fractal exponents with large weights, which are characteristic of fine structures and indicate a relative dominance of the small fluctuations. e 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. e 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. erefore, the values of the multifractal spectrum parameters (α 0 , w, and a s ) 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. eir values describe the distribution patterns of high and low fluctuations [47]. e 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. e asymmetry index (R) is obtained as where Δα L � α 0 − α min and Δα R � α max − α 0 .

Description of the EEMD Algorithm.
e 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). e 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 series y (t) (2) Connect that extrema by spline interpolation, to create an upper envelope of local maxima points e max (t), and a lower envelope e min (t), of all minima points (3) Compute the mean m (t) of two envelopes by m � ((e max (t) + e min (t))/2) (4) Subtract the mean from the data to establish an IMF candidate; h(t) � y(t) − m(t) (5) If h(t) satisfies the two properties of IMF as indicated by a predefined stopping criterion, h(t) is indicated as the first IMF and is recorded as replaced by h(t), and iterate steps 1-4 until h (t) satisfies the two conditions of IMFs (6) e residue r i (t) � y(t) − IMF i (t) is then treated as new data subjected to the same shifting procedure as defined above for the next IMF from r i (t). In the final step, the shifting process can be stopped when the residue r(t) 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 y(t) can be reconstructed as the sum of IMFs and residual as where r n (t) represents the final residue and n is the number of IMFs. e residue r n (t) is known as the trend of the signal y(t) [15]. Although EMD has obvious advantages in time series data processing, there are also some unavoidable weaknesses [14,15]. e 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. us, 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. e EEMD algorithm is straight forward and can be summarized as follows: (1) Add the white noise series into the original series 6 Complexity (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. e effect of adding white noise should obey the following statistical rule: 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].

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. e 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].
where i � 1, 2, . . . , n and λ is the threshold. Moreover, μ i is Pearson's correlation coefficient between the i th IMF and the original signal, and n is the total number of IMFs; max(μ i ) is the maximum observed correlation coefficient. e selection criterion for IMFs is given as follows: if μ i > λ, then keep the i th IMF, else eliminate the i th IMF, and add it into the residue.

EEMD-MFDFA-Based Method.
e selection of the scale range for computing the fluctuation function F q (s) 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 realtime 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. e discontinuity at the segmentation points using polynomial fitting comes from the new pseudofluctuation errors and estimation errors in the q th -order generalized Hurst exponent h(q). Secondly, the q th -order fluctuation function F q (s) 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].
us, in the standard MFDFA method, the third step is reformulated as follows: Step 3: For each segment ], one obtains the EEMDbased 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

Complexity 7
We denote this MFDFA with the EEMD-based detrending method as EEMD-MFDFA. e 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. is division was made considering the classification often mentioned by some studies in West Africa [54,55]. e 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, a s ) in the subsets is done.

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. irteen (13) IMFs (not shown) are obtained in the subequatorial region and twelve (12) IMFs (not shown) in the Sudanian region (except Kandi station). e 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 Figure 4, the curve of residue r (t) obtained during 1951-2010 at Cotonou synoptic station is presented. e curve of the residue obtained at the other studied synoptic stations is not shown. ese 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.
is 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. e 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. ey 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. e dominant components in terms of a good correlation between the IMFs and the original signal are identified. Pearson's correlation coefficient (PCC) between the i th IMF and the original signal at Cotonou synoptic station is shown in Figures 5. e 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 (IMF 1 , IMF 2 , IMF 3 , IMF 4 , and IMF 5 ) reveal high PCC values than the threshold value. us, according to the selection criterion, these IMFs are the relevant IMFs in all the stations.
ese results mean that, in Benin, precipitation is mainly controlled and governed by these five IMFs, influencing significantly the precipitation variations. Table 2 presents the average periods calculated for each of the IMFs obtained during the period 1951-2010 at each of the synoptic stations. e values of the period (in the day) are highlighted in bold.
In Table 2, it is observed that the decomposed data presents nearly identical periods for the lower IMFs such as IMF 1 to IMF 5 in all synoptic stations. e periods of the relevant IMFs are between 1 and 25 days. e 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 IMF 10,11,12,13 . Figure 6 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. e variation of τ (q) with q obtained for the other studied synoptic stations is not shown. e 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. e 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. erefore, 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). e relationship between singularity spectrum f (α) and singularity exponent α of Cotonou synoptic station' precipitation in the subsets is displayed in Figure 7. Similar     the multifractal spectra obtained at all synoptic stations reveals evident differences in the dynamics of the precipitations processes for Benin. e multifractal spectra of the precipitation quantities vary regarding the climate region of the synoptic stations and the considered subperiods. e structures of the multifractal spectrum vary relatively to the synoptic stations and the study period, indicating a different evolution law. is 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]. erefore, multifractal spectra analysis revealed that the effect of small fluctuations plays a dominant role in daily precipitation time series. e result of the asymmetry index values (R) obtained for different synoptic stations and periods is presented in Figure 8(a)). e 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. at 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 R 1991-2010 < R 1971-1990 < R 1951-1970 , whereas in the subequatorial region the variation on R is random. e multifractal spectrum parameters (a s , α 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. e width of the multifractal spectrum w is used to measure the degree of multifractality, which denoted a nonuniform clustering structure for daily precipitation. e analysis of changes in multifractal properties shows that there is a clear difference between the subperiods in multifractal properties. is 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 parameter a s (Figure 8(b)), more consistent changes are observed. e 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 a s 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 a s [1951− 1970] < a s [1971− 1990] < a s [1991− 2010] . However, no explicit direction of changes is observed in subequatorial stations. In the Sudanian region, a s 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, a s 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 a s [1951− 1970] < a s [1971− 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 a s [1971− 1990] < a s [1991− 2010] ; thus, the 1971-1990 period had recorded more extreme events than 1991-2010.

Multifractal Spectrum Analysis.
At Bohicon and Save stations (subequatorial region), the lowest positive value of a s 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), a s 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. e asymmetry a s values for the subequatorial region changed differently compared to the Sudanian region. is result reveals evident differences in the dynamics of the precipitation processes in the southern part of Benin.
is 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 width w (also known as the richness of the signal) presented in Figure 8(c)), its values changed during the studied periods. erefore, multifractal properties could be regarded as good indicators of changes in the dynamics of precipitation.
One can notice that the greater w 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. is 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. erefore, in this subperiod, the degree of multifractality is stronger and the structure of daily precipitation is more complex. e 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). e results of precipitation multifractal parameters analysis indicate that the spatial and temporal variability of α 0 Figure 8(d) differs from the asymmetry a s and the width w parameters (Figures 8(b) and 8(c)). erefore, 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.  Figure 8: Spatiotemporal variation of the asymmetry index values (R) and the multifractal spectrum parameters (a s , α 0 , and w) of precipitation data. Co, Bo, Sa, Pa, Na, and Ka mean Cotonou, Bohicon, Save, Parakou, Natitingou, and Kandi synoptic stations, respectively. e spatiotemporal variation of (R) a s , α 0 , and w are shown at panels (a), (b), (c), and (d), respectively.

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. e 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. e 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. e 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, IMF 3 , IMF 4 , and IMF 5 ) in which periods are between 1 and 25 days, corresponding to the local physical phenomenon (3) e 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. erefore, multifractal spectra analysis revealed that the effect of large fluctuations plays a dominant role in the daily precipitation series (5) e analysis of the multifractal spectrum parameters (α 0 , w, and a s ) shows that there is a clear difference between the analyzed periods in multifractal spectrum parameters. e changes in the results for different subperiods reveal that these subperiods have different climatic modifications. ese 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 stations Our 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. is can be performed in future studies.

Data Availability
e precipitation data used in this study are supplied by the local service of ASCENA in Cotonou. e data are not available online in any database so that we cannot provide a link to reach them. ey are provided when researchers address requests to ASCENA (http://www.asecna.aero).

Conflicts of Interest
On behalf of all authors, the corresponding author states that there are no conflicts of interest.