Abstract
The theoretical and practical understanding of projected changes in rainfall is desirable for planning and adapting to climate change. In this study, finite size Lyapunov exponents (FSLE) are used to study error growth rates of the system at different timescales. This is done to quantify the impact of enhanced anthropogenic greenhouse gas emissions on the predictability of fast and slow varying components of central Indian rainfall (CIR). The CIR time series for this purpose is constructed using the daily gridded high-resolution India Meteorological Department (IMD) dataset and Coupled Model Inter-comparison Project phase 5 (CMIP5) output for historical run and three representative concentration pathways (RCP2.6, RCP4.5, and RCP8.5) from the HadGEM2-ES, IPSL-CM5A-LR, CCSM4, BCC-CSM1.1, and MPI-ESM-LR models. The analyzed CIR dataset reveals a low dimensional chaotic attractor, suggesting that CIR requires a minimum of 5 and maximum of 11 variables to describe the state of the system. FSLE analysis shows a rapid decrease in the Lyapunov exponent with increasing timescales. This analysis suggests a predictability of about 2–3 weeks for fast varying components at short timescale of the CIR and about 5–9 years for slow varying components at long timescales.
INTRODUCTION
The predictability study of rainfall or another atmospheric variable has its own importance. The projection from the Coupled Model Inter-comparison Project phase 5 (CMIP5)-based model outputs for historical as well as future scenarios has enabled us to utilize these data to study the predictability of the system for historical as well as future scenarios. This quantification of the impact of the increase in greenhouse forcing into the atmospheric predictability can be possible with these long-term data. The influence of anticipated changes in greenhouse gases on climate predictability over timescales of several decades to centuries can be explored using global models (Meehl et al. 2009).
The Lyapunov exponent of a dynamical system measures the average rate at which nearby trajectories of a system diverge. The positive exponent provides evidence that the underlying dynamical system is chaotic and there is a limit on the predictability of the dynamical system, i.e., the dynamical system may be predictable over some finite time periods. The nonlinear dynamical method – finite size Lyapunov exponent (FSLE) – is used for quantifying the predictability. It has several advantages over the traditional methods, since nonlinear dynamical methods like largest Lyapunov exponent (LLE) (Rosenstein et al. 1993; Kantz 1994) and Kolmogorov–Sinai entropy quantify the error growth rate exhibited by a dynamical system (Boffetta et al. 2002). The predictability analysis based on LLE holds for infinitesimal perturbations, which are linked with errors in the smallest, fastest-evolving scales of a dynamical system. However, the smallest errors, if finite, soon grow to nonlinear saturation, thereafter not playing any role in the error growth (Boffetta et al. 1998b; Cencini et al. 1999; Siqueira & Kirtman 2014,) of fast scale dynamics. On the other hand, relatively larger errors have slower error growth rates on larger timescales. Thus, LLE cannot estimate the predictability on longer timescales (Toth & Kalnay 1993; Boffetta et al. 1998a, 1998b). This advantage of FSLE over LLE motivated us to use the method for quantifying the predictability of both long and short timescales.
Aurell et al. (1997) introduced the FSLE as a measure of the divergence of nearby trajectories, to quantify predictability of processes developing on different timescales. Boffetta et al. (1998a) further advanced the concept of FSLE and investigated its applicability to detect large-scale (slow varying) dynamical properties of a measured time series characterized by different timescales. They showed that the large-scale error growth rate due to slow varying components of a dynamical system can be estimated by FSLE, even though data points are limited and embedding dimension is moderate. The FSLE has been used to estimate predictability on both fast and slow timescales (Boffetta et al. 1998a; Cencini & Vulpiani 2013; Siqueira & Kirtman 2014). The climate system is supposed to be a system with many characteristic timescales. In a study by Boffetta et al. (1998a), dynamics of a system with two different timescales were analyzed using the FSLE method. It was found that the system provides two plateaus, one for short-scale perturbation, which reveals FSLE of fast varying component λf and the second plateau at large perturbation, which gives FSLE of slow varying component λs. They also investigated the role of short-scale parameterization by replacing the true dynamics of the fast-varying component with the stochastic process and found that even though the fast dynamics is not resolved it does not have any effect on correct computation of FSLE of slow varying components. It is also pointed out in their study that, even with a smaller number of data points and also with a moderate embedding dimension, the dynamics of slow varying components can be practically captured. Siqueira & Kirtman (2014) investigated the predictability of ENSO using FSLE on the fast timescale. On the other hand, Mitchell & Gottwald (2012) suggested that FSLE can well approximate the predictability at the slow timescale. It is suggested in their study that predictability at the slower timescale can be much higher than the local LLE λmax, since short-scale instabilities grow faster but become nonlinearly saturated at a much smaller time than large-scale instabilities. Thus, in the present study, the strength of FSLE over LLE to efficiently extract predictability at both fast and slow timescales motivated us to quantify the predictability of the central Indian rainfall (CIR) on both fast and slow time. The scale at small perturbation is referred to as short timescale and the scale at large perturbation, where the value of λs saturated (plateau) is referred to as long timescale.
Studies about the impact of global warming on the Indian monsoon at different timescales are very limited. Mani et al. (2009) studied the impact of global warming over the medium range prediction. They showed that the medium range prediction of Indian monsoon will be less predictable due to global warming (Mani et al. 2009). They addressed the impact of climate change over medium range prediction but did not study its impact over long timescales. The short-scale predictability of CIR is about 2–3 weeks, but the clear limit on the predictability and role of climate change on long timescales is not well established. Therefore, in the present study, an attempt is made to quantify the predictability of CIR for historical as well as future projections using the theory of nonlinear dynamics on both short and long timescales.
There are several other studies, both on toy model and real-world scenarios, which suggest that the inclusion of forcing increases the predictability of the system (Lorenz model – Palmer 1993; Dwivedi et al. 2007; Singh et al. 2015). Kirtman et al. (2013a) have pointed out that external forcing like aerosol concentration, greenhouse gases, volcanic eruption, etc. provides predictability to the system, ranging from sub-seasonal to decadal timescales. Latif (2013) states that, on longer lead times of several decades, external forcing will cause weakening in Atlantic meridional overturning circulation, which will give rise to predictability. Luo et al. (2011) inferred that the influence of global warming results in an increase in predictability of the system on seasonal to inter-annual timescales.
The climate system is a combination of several timescales and the longer timescale components play an important role in climate change (Bigg et al. 2003). Climatic variables like rainfall also have variability from sub-seasonal to multidecadal timescales. Climate patterns such as the El-Niño Southern Oscillation (ENSO), Pacific Decadal Oscillation (PDO), etc. reveal climate variability (Stocker 2013) and the Indian monsoon is found to be influenced by ENSO and the decadal oscillations like PDO and Atlantic multidecadal oscillation (AMO) (Krishnamurthy & Krishnamurthy 2016, 2017). The ocean surface variables play a vital role in climate variability on shorter timescales, for example, sea-surface temperature (SST) changes leading to El Niño and La Niña events. On the other hand, the slow varying nature of deep ocean currents produces little influence on climate variability at shorter timescales (Kirtman et al. 2013a). However, on longer timescales, these deep ocean currents may play a significant, even a dominant role (Bigg et al. 2003; Marshall & Zanna 2014).
The fast-varying components of the atmosphere effect predictability at short timescales and the slow varying components of the atmosphere effect predictability at long timescales. The study of predictability at slow timescale is important because future emission scenarios have less influence over near-term predictions than long-term projections (Kirtman et al. 2013b). Over longer periods, these effects can accumulate to cause a critical change (Doblas-Reyes et al. 2006). For example, subtle yearly changes in the thickness of ocean ice in polar locales (influencing the surface temperature of the ice and its albedo) may cause long-term cumulative impacts driving climate change. Therefore, slow processes, which do not influence variability on seasonal to inter-annual timescales, may become important in influencing long-term climate variability and hence predictability (Hurrell et al. 2010; Solomon et al. 2011). At long timescales, the predictability of the system is primarily determined by the ocean, because it has large thermal capacity and longer time for adjustment (Slingo & Palmer 2011).
The predictability horizon depends on the initial condition and on the variable chosen (Anderson et al. 2009). The predictability of fast-varying components is generally governed by the internally generated atmospheric variability, i.e. weather, which loses its memory in around a month (Lorenz 1963; Shukla 1981). On the other hand, the predictability of slow-varying components varies at decadal timescales (Tang et al. 2008). However, the absolute limit on the predictability of slow-varying components is not well established. The Hasselmann mechanism suggests that the ‘scale of predictability is restricted to the order of the autocorrelation time of the slower component of the system, usually the ocean’ (NRC 1998). That is, the long-term memory of the ocean offers some hope for initialized decadal-scale prediction and predictability (Li et al. 2017). There are a large number of other studies which support the presence of predictability of atmospheric variations on multiyear timescale (Griffies & Bryan 1997a, 1997b; Collins et al. 2006; Boer & Lambert 2008; Meehl et al. 2009, 2014; DelSole et al. 2013; Ding et al. 2016; Wu et al. 2016).
The rainfall observation data from IMD and the data from IPCC CMIP5 model output of five models have been used. The main aim of this study is to assess the impact of increase in anthropogenic forcing on the predictability of short and long timescales of CIR using FSLE. Our result for the short timescale is in agreement with the result of Mani et al. (2009). The predictability knowledge of rainfall under future emission scenarios will empower preparedness for optimal water resource usage and disaster management. This paper is organized as follows. Below, the study region and datasets used are described, followed by the methodology adopted. Then, the next section shows the results of the calculated embedding dimensions and FSLE for both fast and slow timescales of IMD observation and CMIP5 model data. The final section explains the conclusions drawn from the results.
STUDY REGION AND DATASETS
Time series of the daily central Indian rainfall (CIR) is constructed by area (spatial) averaging the daily gridded (0.25 × 0.25) IMD dataset over the monsoon trough region, (Figure 1) [15N–30N; 70E–90E]. A lowpass filter of 10 days is applied in each case in order to filter high frequency components having frequency up to 10 days. The choice of 10-day filter is made in order to preserve the information with the smallest loss and produce clearer information from the data, so that the high frequency noise component from the data is filtered. On the other hand, without applying the filter, due to the presence of high frequency components in the data, correlation exponent vs embedding dimension (m) is increasing with the increase in m, which is a behavior of random noise. Therfore, choice of the filtering high frequency components from the data is justified. The time series are also created using similar methodologies for Historical and RCP Scenarios of the IPCC AR5 CMIP5 models (Taylor et al. 2012), namely, HadGEM2-ES, IPSL-CM5A-LR, CCSM4, BCC-CSM1.1, and MPI-ESM-LR. Historical simulations are based on solar and volcanic forcing, land use, observed concentrations of greenhouse gases, and reconstructed aerosol emissions. Future projections are based on the three representative concentration pathways (RCPs) (Moss et al. 2010). In the RCP2.6 pathway the peak in radiative forcing reaches 3 W m−2 before the year 2100 and afterwards it declines to 2.6 W m−2 in the year 2300. RCP4.5 represents the pathways for which radiative forcing reaches 4.5 W m−2 in the year 2300. Similarly, RCP8.5 is the pathway for which radiative forcing reaches 8.5 W m−2 by the year 2300. Most of the GCMs are unable to simulate the mean, intra-seasonal variability of summer monsoon rainfall. According to Ashfaq et al. (2017), among the GCMs under study, only the MPI-ESM-LR model is able to simulate both rainfall maxima in July and less than 25% bias in precipitation for the south Asian summer monsoon during the historical time 1970–1999. However, the model selection is done on the basis of availability of long-term data to get the required data for time series reconstruction. Tsonis et al. (1993) have shown that 10(2+0.4D) number of data points are sufficient to correctly estimate the correlation dimension, therefore, the number of data points listed in Table 1 are adequate for faithful computation on correlation dimension. Therefore the data up to the year 2300 is used, which is a combination of RCP up to 2100 and after that, extended RCP up to 2300; for simplicity it is termed RCP throughout this paper. Table 1 lists the number of data points and durations that have been used and Table 2 lists the model description.
Number of data points used in the CIR time series and time duration
Scenario . | Duration . | Years . | Data points . |
---|---|---|---|
Historical | 1860–2005 | 146 years | 53,290 |
RCP2.6 | 2006–2300 | 295 years | 107,675 |
RCP4.5 | 2006–2300 | 295 years | 107,675 |
RCP8.5 | 2006–2300 | 295 years | 107,675 |
IMD observation | 1901–2013 | 113 years | 41,273 |
Scenario . | Duration . | Years . | Data points . |
---|---|---|---|
Historical | 1860–2005 | 146 years | 53,290 |
RCP2.6 | 2006–2300 | 295 years | 107,675 |
RCP4.5 | 2006–2300 | 295 years | 107,675 |
RCP8.5 | 2006–2300 | 295 years | 107,675 |
IMD observation | 1901–2013 | 113 years | 41,273 |
Model descriptions
Modeling center (or group) . | Model . |
---|---|
Met Office Hadley Centre | HadGEM2-ES |
Institute Pierre-Simon Laplace | IPSL-CM5A-LR |
National Center for Atmospheric Research | CCSM4 |
Beijing Climate Center, China Meteorological Administration | BCC-CSM1.1 |
Max Planck Institute for Meteorology | MPI-ESM-LR |
Indian Metrological Department | IMD observation |
Modeling center (or group) . | Model . |
---|---|
Met Office Hadley Centre | HadGEM2-ES |
Institute Pierre-Simon Laplace | IPSL-CM5A-LR |
National Center for Atmospheric Research | CCSM4 |
Beijing Climate Center, China Meteorological Administration | BCC-CSM1.1 |
Max Planck Institute for Meteorology | MPI-ESM-LR |
Indian Metrological Department | IMD observation |
Map of study region (15N–30N; 70E–90E) central India (monsoon trough region).
METHODOLOGY
From a measured time series, phase space reconstruction is done by using the method of time delay as described in the section immediately below. The time delay embedding was done for the total appended time series. Appropriate time delay is estimated as described in the section Mutual information. The method used for determining the embedding dimension (m) that efficiently accommodates the attractor is then described. This is followed by a section describing Determination of finite size Lyapunov exponents. Figure 2 shows the work flow of the methods shown below.
Phase space reconstruction
Mutual information




Mutual information I(t) vs steps (k) for IMD rainfall time series over central India. Inset panel shows delay = 72 (first minima of I(t)).
Mutual information I(t) vs steps (k) for IMD rainfall time series over central India. Inset panel shows delay = 72 (first minima of I(t)).
Embedding dimension
The calculation of embedding dimension m to acquire the minimum number of variables sufficient to model the dynamics of the system is done by using two methods, those of Grassberger & Procaccia (1983) (termed GP1983) and Cao (1997) (termed Cao 1997). Details about both methods are given in the following two subsections, respectively.
Embedding dimension using GP1983
The minimum number of variables required to map the system dynamics is estimated by using correlation dimension. Consider a single variable time series xi, i=1, 2,…, N of length N. The method described earlier is used for phase space reconstruction for different values of m. The correlation sum C(r) in an m-dimensional phase space is defined as the ratio of the number of pairs of embedded data points having distance less than r to the total number of pairs (GP1983).



This mathematical definition cannot be used for experimental data. However, it is often found that C(r) estimated from a reasonably large dataset of size N satisfies a scaling law for a range of r. The scaling exponent
is called the correlation exponent and for a particular embedding dimension m can be calculated from slope of the scaling region of the
vs
graph.
The correlation exponent is computed for different embedding dimensions m. A graph of
vs m is plotted. If the value of
saturates to a constant value
as m is increased, the time series is considered to be chaotic. In this case, the dynamics of the system can be said to reside on a low-dimension attractor of correlation dimension D=
. On the other hand, if
increases without saturation, as m is increased, the analyzed time series is stochastic because a stochastic process embedded in an m dimensional space always fills that space (Tsonis & Elsner 1990). The data size, sampling frequency, and noise present in the data are the main factors that affect estimation accuracy of the correlation dimension.
The smallest integer which is greater than the correlation dimension D gives the smallest number of suitably chosen variables that can map the dynamics. The value of the embedding dimension ms at which the value of the correlation exponent begins to saturate gives the minimum number of dimensions in which the attractor can be embedded completely (Fraedrich 1986). Abarbanel et al. (1990) has pointed out that although Taken's reconstruction theorem requires D>2Ds+1 to ensure faithful representation of the Ds dimensional attractor as seen in D-dimensional embedding space, in practice, D>Ds is usually sufficient.
Minimum embedding dimension using Cao 1997
In the case of the time series from some random process, the E1(m) will not saturate as m is increased, which is due to the limited data samples. Thus, practically, it may be difficult to determine whether E1(m) is changing slowly or has stopped changing. In such a situation, E2(m) will be useful because for random data future values are independent of the past values and E2(m) will be equal to 1 for any m, whereas for deterministic signals there exists some m such that E2(m) ≠ 1.
Finite size Lyapunov exponent (FSLE)






The FSLE λ(δ) were computed from the time series datasets using the method described by Boffetta et al. (1998a, 1998b). The computation of FSLE from measured time series dataset is done by sequentially applying subsequent steps. First, a phase space reconstruction is performed using the method described earlier and the time delay is calculated with the method described in the section Mutual information. Second, a threshold series of perturbation scales is taken. A value δmin ≪ δ0 is chosen so that perturbations of size less than δmin align with the most unstable direction in the phase space. Then, for each reconstructed vector Xn, its nearest neighbors are obtained. Lastly, from these nearest neighbors, only the nearest neighbors (Xm), which have separation
, smaller than δmin, are chosen. The trajectories starting from reconstructed vector Xn and its nearest neighbor Xm are used to calculate the error doubling time
by measuring the growth time of perturbation from δn to δn+1. The growth of perturbation is followed from initial size less than δmin to the maximum value δmax in the threshold series. Finally, The FSLE λ(δn) at perturbation size δn is obtained by averaging the error doubling time τ(δn) for several nearest neighbors Xm and using Equation (13). A measure of the predictability of a system is the error doubling time. It is thus proportional to the inverse of the FSLE (Mitchell & Gottwald 2012).
RESULTS
The nonlinear time series analysis of different datasets – CIR time series for IMD observation, historical, and various RCP scenarios of CMIP5 model outputs over central India – is performed using the method described in the Methodology section. The method described in Equations (5) and (6) (GP1983) is used for calculation of the correlation dimension D to obtain embedding dimension (m) for each of the datasets. In Figure 4(a), the plot of log C(r) vs log(r) is shown for IMD CIR data with the embedding dimension m ranging from 1 to 20 (from left to right). The slope of the straight line within the appropriately chosen scaling region for each m gives the correlation exponent v corresponding to that m. In Figure 4(b), the plot between correlation exponent v and embedding dimension m shows that Ds saturates to a value of 4.46 at about ms = 11. Thus, the CIR data reveal the existence of a low dimensional chaotic attractor having saturation dimension Ds= 4.46. Hence, the next positive integer above Ds, i.e., 5, is the minimum number of variables sufficient to describe the dynamics of the system (Fraedrich 1986). On the other hand, to confirm the deterministic features in data the surrogate data test is performed (Theiler et al. 1992). The four ensembles of surrogate time series were constructed by shuffling randomly the indices of HadGEM2-ES model historical time series data and then the method of correlation dimension on ensemble surrogate time series is applied. There is no saturation found in these cases, as shown in Figure 5. We also performed surrogate data tests on the rest of the data and similar results were found (figure not shown). Thus, the convergence of ν with m is a genuine characteristic of the data and not an artifact of the method used, confirming the deterministic features of the data. The GP1983 method is also chosen for calculation of correlation dimensions of all the other CMIP5 models listed in Table 2. Figure 6 shows the plot of the correlation exponent vs embedding dimension (m) for historical as well as all the RCP scenarios for the HadGEM2-ES model of CMIP5 (figures for the rest of the models are not shown). The fractal dimension of the historical and RCP scenarios for different models of CMIP5 have been calculated and listed in Table 3. From Table 3, it is clear that similar results have been obtained for all IPCC AR5 CMIP5 models. In the present study, D from Table 3 is found to be 5 for all the models. To validate the correct computation of minimum dimension that is sufficient to represent the dynamics of the system the Cao 1997 method is used. Figure 7 shows E1 and E2 quantities and E1 gives the number of minimum variables sufficient to represent the dynamics of the system and is found to be 5 (saturation value, m0+1) for the HadGEM2-ES model. It is similar to embedding dimension m calculated from the GP1983 method (Table 3). The E2 quantity is also calculated, as mentioned in the Cao 1997 method (represented by diamonds in Figure 7) to check whether the datasets are the outcome of the deterministic or stochastic processes. E2 is also found to be less than 1 for each dimension m, which ensures the dataset under study results from a deterministic process. To check whether the Cao 1997 method is able to capture stochasticity in the surrogate dataset, the method is also applied on the surrogate dataset obtained from the method described above for the HadGEM2-ES RCP2.6 scenario. Figure 8 clearly shows that the E2 ≥ 1 as m increases, so the method is able to capture stochasticity of the surrogate dataset and it is also evident from the figure that the E1 also do not saturate, which is expected for the stochastic dataset. The embedding dimension for the rest of the models listed in Table 2 is also calculated using the Cao 1997 method and similar results are evident as in the GP1983 method, i.e., embedding dimension m= 5 and the E2 quantities for these models also suggest the deterministic nature of datasets since E2 ≤ 1 for each case (figure not shown).
Correlation dimensions of IMD and various scenarios of CMIP5 model
Scenario . | Correlation dimension . | |||||||||
---|---|---|---|---|---|---|---|---|---|---|
HadGEM2-ES . | IPSL-CM5A-LR . | CCSM4 . | BCC-CSM1 . | MPI-ESM-LR . | ||||||
Ds . | ms . | Ds . | ms . | Ds . | ms . | Ds . | ms . | Ds . | ms . | |
Historical | 4.29 | 11 | 4.47 | 11 | 4.51 | 11 | 4.32 | 11 | 4.27 | 11 |
RCP2.6 | 4.78 | 11 | 4.58 | 11 | 4.71 | 11 | 4.82 | 11 | 4.75 | 11 |
RCP4.5 | 4.82 | 11 | 4.65 | 11 | 4.81 | 11 | 4.87 | 11 | 4.82 | 11 |
RCP8.5 | 4.54 | 11 | 4.68 | 11 | 4.85 | 11 | 4.85 | 11 | 4.91 | 11 |
IMD | 4.46 | 11 |
Scenario . | Correlation dimension . | |||||||||
---|---|---|---|---|---|---|---|---|---|---|
HadGEM2-ES . | IPSL-CM5A-LR . | CCSM4 . | BCC-CSM1 . | MPI-ESM-LR . | ||||||
Ds . | ms . | Ds . | ms . | Ds . | ms . | Ds . | ms . | Ds . | ms . | |
Historical | 4.29 | 11 | 4.47 | 11 | 4.51 | 11 | 4.32 | 11 | 4.27 | 11 |
RCP2.6 | 4.78 | 11 | 4.58 | 11 | 4.71 | 11 | 4.82 | 11 | 4.75 | 11 |
RCP4.5 | 4.82 | 11 | 4.65 | 11 | 4.81 | 11 | 4.87 | 11 | 4.82 | 11 |
RCP8.5 | 4.54 | 11 | 4.68 | 11 | 4.85 | 11 | 4.85 | 11 | 4.91 | 11 |
IMD | 4.46 | 11 |
(a) C(r) vs r plot of IMD rainfall time series over central India. (b) Correlation exponent vs embedding dimension (m) plot of IMD rainfall time series over central India. The next positive integer above saturation dimension Ds, i.e., D= 5, is the minimum number of variables sufficient to describe the dynamics of the system.
(a) C(r) vs r plot of IMD rainfall time series over central India. (b) Correlation exponent vs embedding dimension (m) plot of IMD rainfall time series over central India. The next positive integer above saturation dimension Ds, i.e., D= 5, is the minimum number of variables sufficient to describe the dynamics of the system.
Correlation exponent vs embedding dimension (m) plot of four ensemble surrogate data formed by shuffling the indices of HadGEM2-ES historical scenario rainfall time series over central India. The value of
increases as m increases for all the four surrogates; on the other hand, for HadGEM2-ES historical time series Ds saturates to a value of 4.29. Thus, the convergence of ν with m confirms the deterministic features of the data.
Correlation exponent vs embedding dimension (m) plot of four ensemble surrogate data formed by shuffling the indices of HadGEM2-ES historical scenario rainfall time series over central India. The value of
increases as m increases for all the four surrogates; on the other hand, for HadGEM2-ES historical time series Ds saturates to a value of 4.29. Thus, the convergence of ν with m confirms the deterministic features of the data.
Correlation exponent () vs embedding dimension (m) plot of HadGEM2-ES for all four scenarios. Figure indicates a minimum of 5 and maximum of 11 variables are sufficient to efficiently accommodate the attractor.
Correlation exponent () vs embedding dimension (m) plot of HadGEM2-ES for all four scenarios. Figure indicates a minimum of 5 and maximum of 11 variables are sufficient to efficiently accommodate the attractor.
E1 and E2 vs dimension plot of HadGEM2-ES model for all four scenarios: (top left panel: historical; top right panel: RCP2.6; bottom left panel: RCP4.5; bottom right panel: RCP8.5). The quantity E1 (represented by star) starts to saturate at m0 = 4 so that m= 5 is the minimum embedding dimension. It is also evident that E2< 1 (represented by diamond) for each case, thus the time series are deterministic.
E1 and E2 vs dimension plot of HadGEM2-ES model for all four scenarios: (top left panel: historical; top right panel: RCP2.6; bottom left panel: RCP4.5; bottom right panel: RCP8.5). The quantity E1 (represented by star) starts to saturate at m0 = 4 so that m= 5 is the minimum embedding dimension. It is also evident that E2< 1 (represented by diamond) for each case, thus the time series are deterministic.
E1 and E2 vs dimension plot for surrogate data of HadGEM2-ES RCP2.6 scenario. The quantity E1 (represented by star) is not saturated as m increases and it is also evident that E2 ≥ 1 (represented by diamond) for each case, thus the method is able to capture stochasticity in the dataset.
E1 and E2 vs dimension plot for surrogate data of HadGEM2-ES RCP2.6 scenario. The quantity E1 (represented by star) is not saturated as m increases and it is also evident that E2 ≥ 1 (represented by diamond) for each case, thus the method is able to capture stochasticity in the dataset.
The calculation of the FSLE λ(δ) is done by using Equation (13) to be on the safe side embedding dimension equal to 2D + 1, which is 11 for all datasets (as mentioned before, D is the next positive integer above saturation dimension Ds), is used as the embedding parameter for the calculation of FSLE. The FSLE can give information about how a finite initial uncertainty δ grows with time. The perturbation threshold δmin = 0.07 is selected. The k-nearest neighbor method is used for the computation of nearest neighbors and the Euclidian distance is used as a distance parameter. The 20 nearest neighbors corresponding to each point in the reconstructed phase space were calculated and the nearest neighbors which have distance less than δmin are chosen for calculation of FSLE (for instance, 23,305 nearest neighbors found below δmin in the case of HadGEM2-ES historical scenario). It is also verified that increasing the number of nearest neighbors greater than ten, corresponding to each point in the reconstructed phase space, does not contribute to an increase in the number of nearest neighbors below δmin, thus it is the optimal choice for the current analysis. The FSLE λ(δ), as a function of δ for the IMD rainfall observed data, is shown in Figure 9. The FSLE λ(δ) for the short scale, i.e., fast-varying component (λf) of IMD rainfall observed data saturates to a value 0.0343 (Table 4). The error doubling time obtained by using Equation (13) can be calculated from the formula , and is around 20.02 days which is in agreement with the weather predictability found by Dwivedi (2012). This quantifies the predictability time for fast variables. The FSLE λ(δ) for large perturbations, that evolve according to slow-varying components of IMD rainfall dynamics, saturate to a value of 3.34 × 10−4 (Table 5), which indicates a predictability time of 5.7 years (Table 7) by using formula
. Our analysis of the rainfall data suggests that the data can be the outcome of two weakly coupled dynamical systems, one having a fast timescale and the other having a slow timescale. The FSLE λ(δ) for the historical and all RCP scenarios from all five models have been calculated by using the same method as described above for IMD rainfall observed data and, for brevity, only one figure for computation of FSLE λ(δ) for the HadGEM2-ES model is shown in Figure 10. (Figures for the calculation of FSLE for the rest of the models are shown in Figures S1 to S4 in the Supplementary material.)
FSLE of IMD and various scenarios of CMIP5 model for small perturbations
Scenarios . | Finite size Lyapunov (λf) exponent for small perturbations . | ||||
---|---|---|---|---|---|
HadGEM2-ES . | IPSL-CM5A-LR . | CCSM4 . | BCC-CSM . | MPI-ESM-LR . | |
Historical | 0.0352 | 0.0461 | 0.0373 | 0.0359 | 0.0324 |
RCP2.6 | 0.0373 | 0.0425 | 0.0411 | 0.0345 | 0.0484 |
RCP4.5 | 0.0360 | 0.0393 | 0.0445 | 0.0351 | 0.0480 |
RCP8.5 | 0.0315 | 0.0338 | 0.0324 | 0.0352 | 0.0398 |
IMD | 0.0343 |
Scenarios . | Finite size Lyapunov (λf) exponent for small perturbations . | ||||
---|---|---|---|---|---|
HadGEM2-ES . | IPSL-CM5A-LR . | CCSM4 . | BCC-CSM . | MPI-ESM-LR . | |
Historical | 0.0352 | 0.0461 | 0.0373 | 0.0359 | 0.0324 |
RCP2.6 | 0.0373 | 0.0425 | 0.0411 | 0.0345 | 0.0484 |
RCP4.5 | 0.0360 | 0.0393 | 0.0445 | 0.0351 | 0.0480 |
RCP8.5 | 0.0315 | 0.0338 | 0.0324 | 0.0352 | 0.0398 |
IMD | 0.0343 |
FSLE of IMD and various scenarios of CMIP5 model for large perturbations
Scenarios . | Finite size Lyapunov (λs) exponent at large perturbations (×10−4) . | ||||
---|---|---|---|---|---|
HadGEM2-ES . | IPSL-CM5A-LR . | CCSM4 . | BCC-CSM1 . | MPI-ESM-LR . | |
Historical | 3.78 | 3.60 | 3.67 | 3.41 | 3.62 |
RCP2.6 | 3.61 | 3.00 | 3.12 | 2.73 | 3.38 |
RCP4.5 | 3.11 | 2.47 | 3.14 | 2.58 | 2.43 |
RCP8.5 | 2.81 | 2.12 | 2.96 | 2.32 | 2.47 |
IMD | 3.34 |
Scenarios . | Finite size Lyapunov (λs) exponent at large perturbations (×10−4) . | ||||
---|---|---|---|---|---|
HadGEM2-ES . | IPSL-CM5A-LR . | CCSM4 . | BCC-CSM1 . | MPI-ESM-LR . | |
Historical | 3.78 | 3.60 | 3.67 | 3.41 | 3.62 |
RCP2.6 | 3.61 | 3.00 | 3.12 | 2.73 | 3.38 |
RCP4.5 | 3.11 | 2.47 | 3.14 | 2.58 | 2.43 |
RCP8.5 | 2.81 | 2.12 | 2.96 | 2.32 | 2.47 |
IMD | 3.34 |
The FSLE, λ(δ), as a function of scale δ for the IMD observation. The λf and λs (plateau) are the values of FSLE at short timescale and long timescale, respectively.
The FSLE, λ(δ), as a function of scale δ for the IMD observation. The λf and λs (plateau) are the values of FSLE at short timescale and long timescale, respectively.
The FSLE, λ(δ), as a function of scale δ for all the four scenarios of HadGEM2-ES. The λf and λs are the values of FSLE at short timescale and long timescale, respectively.
The FSLE, λ(δ), as a function of scale δ for all the four scenarios of HadGEM2-ES. The λf and λs are the values of FSLE at short timescale and long timescale, respectively.
The FSLE of historical and of future RCP projections for HadGEM2-ES, IPSL-CM5A-LR, CCSM4, BCC-CSM1.1, and MPI-ESM-LR models at small perturbations are shown in Table 4. The FSLE results for all the models and IMD observation for larger perturbations are listed in Table 5. The positive values of the FSLE confirm the chaotic nature of CIR time series (Hilborn 2000; Sprott 2003). Tables 6 and 7 show the estimated predictability times for different CIR datasets on fast and slow timescales, respectively. The predictability of the fast-varying component of CIR for all the datasets examined was found to be between 14 and 22 days. On the other hand, the predictability for slow-varying component of CIR for IMD observation and all historical scenarios was found to be around five to six years while the predictability of the CIR under the RCP scenarios was around six to nine years.
Predictability (in days) of the fast-varying components (short timescale), calculated by using formula , where values of λf are listed in Table 4
Scenarios . | Predictability (in days) of the fast-varying component . | ||||
---|---|---|---|---|---|
HadGEM2-ES . | IPSL-CM5A-LR . | CCSM4 . | BCC-CSM1 . | MPI-ESM-LR . | |
Historical | 19.7 | 15.0 | 18.6 | 19.3 | 21.4 |
RCP2.6 | 18.6 | 16.3 | 16.9 | 17.0 | 14.3 |
RCP4.5 | 19.3 | 17.6 | 15.6 | 18.6 | 14.4 |
RCP8.5 | 22.0 | 20.5 | 21.4 | 16.0 | 17.4 |
IMD | 20.2 |
Scenarios . | Predictability (in days) of the fast-varying component . | ||||
---|---|---|---|---|---|
HadGEM2-ES . | IPSL-CM5A-LR . | CCSM4 . | BCC-CSM1 . | MPI-ESM-LR . | |
Historical | 19.7 | 15.0 | 18.6 | 19.3 | 21.4 |
RCP2.6 | 18.6 | 16.3 | 16.9 | 17.0 | 14.3 |
RCP4.5 | 19.3 | 17.6 | 15.6 | 18.6 | 14.4 |
RCP8.5 | 22.0 | 20.5 | 21.4 | 16.0 | 17.4 |
IMD | 20.2 |
Predictability (in years) of the slow varying components (long timescale), calculated by using, , where values of λs are listed in Table 5
Scenarios . | Predictability (in years) of the slow-varying component . | ||||
---|---|---|---|---|---|
HadGEM2-ES . | IPSL-CM5A-LR . | CCSM4 . | BCC-CSM1 . | MPI-ESM-LR . | |
Historical | 5.0 | 5.3 | 5.2 | 5.6 | 5.2 |
RCP2.6 | 5.3 | 6.3 | 6.1 | 7.0 | 5.6 |
RCP4.5 | 6.1 | 7.7 | 6.0 | 7.4 | 7.8 |
RCP8.5 | 6.8 | 9.0 | 6.4 | 8.2 | 7.7 |
IMD | 5.7 |
Scenarios . | Predictability (in years) of the slow-varying component . | ||||
---|---|---|---|---|---|
HadGEM2-ES . | IPSL-CM5A-LR . | CCSM4 . | BCC-CSM1 . | MPI-ESM-LR . | |
Historical | 5.0 | 5.3 | 5.2 | 5.6 | 5.2 |
RCP2.6 | 5.3 | 6.3 | 6.1 | 7.0 | 5.6 |
RCP4.5 | 6.1 | 7.7 | 6.0 | 7.4 | 7.8 |
RCP8.5 | 6.8 | 9.0 | 6.4 | 8.2 | 7.7 |
IMD | 5.7 |
DISCUSSION AND CONCLUSIONS
India is an agroeconomic country. Most of the areas in the central Indian region are rainfed. Therefore, production and yield of Kharif (April to September (summer)) and Rabi season (October to March (winter)) crops (like rice, wheat, etc.) majorly depend on the CIR. Nowadays, groundwater is heavily used in several places in central India for irrigation, but too much exploitation of groundwater resources has caused diminution of this finite resource. The increase in population, food demand, urbanization, and industrialization have caused some serious concerns over fresh water availability (Mall et al. 2007), and fresh water resources are likely to decrease due to climate change (Gosain et al. 2006). Climate change is supposed to increase extreme events, which could potentially alter the rainfall distribution in terms of time and space. The increase in aridity due to climate warming is likely to add to irrigation water demand (Hanasaki et al. 2013; Leng et al. 2015). These factors will cause serious concerns about water-related disasters over central Indian regions. Therefore, revealing the potential impact of climate change on the predictability of CIR has great significance.
The correlation dimensions for the different datasets of the study are determined by the technique described in the section Embedding dimensions. It is evident that a low dimensional chaotic attractor exists for all the datasets and a minimum of 5 and a maximum of 11 variables will be able to simulate the system (Table 3). Similarly, low dimensional attractors from rainfall data have been found by several investigations in the range of 4–9 (Dhanya & Kumar 2010; Jothiprakash & Fathima 2013). The results from Table 3 do not indicate any clear increasing or decreasing trend in the values of correlation dimensions as a result of climate change. In other words, the correlation dimensions obtained for different RCP scenarios of a particular CMIP5 model do not present any clear influence of increase in radiative forcing on the correlation dimension. It is clear from the table that IMD, historical as well as all future scenarios can be considered as residing on a strange attractor of non-integer dimensions (Mandelbrot 1982; Tsonis & Elsner 1989; Hilborn 2000).
The FSLE for small perturbations is substantially in the same range for all the datasets (Table 4). The limit on predictability for the fast-varying component, obtained from the growth rate of small perturbations, is found to be about 2–3 weeks for IMD rainfall, historical scenarios and the RCP scenarios of all the CMIP5 models (Table 6). The predictability estimates for the fast component are similar to those reported in the literature in the range of 2–3 weeks (Goswami & Xavier 2003; Waliser et al. 2003; Webster & Hoyos 2004; Dwivedi et al. 2006; Rai & Krishnamurthy 2011; Dwivedi 2012; Abhilash et al. 2014). Figure 11 (top panel) shows the potential predictability of the fast-varying components of CIR in the historical scenario and IMD observation is greater than the RCPs. This implies a short timescale will be less predictable due to global warming and this result is consistent with the result of Mani et al. (2009), although some models suggest increased predictability of a few days in RCP8.5 for the fast-varying components. A clear reason for this could not be established.
Top panel: Predictability in days for fast-varying component of the CIR. Figure shows a decrease in the predictability in RCP2.6 and 4.5 than in the historical scenario but a few models in RCP8.5 suggest a slight increase in the potential predictability. Bottom panel: Predictability in years for slow-varying component of the CIR. Figure shows an increase in the potential predictability from historical to RCP8.5 due to anthropogenic forcing and longer memory of the oceans.
Top panel: Predictability in days for fast-varying component of the CIR. Figure shows a decrease in the predictability in RCP2.6 and 4.5 than in the historical scenario but a few models in RCP8.5 suggest a slight increase in the potential predictability. Bottom panel: Predictability in years for slow-varying component of the CIR. Figure shows an increase in the potential predictability from historical to RCP8.5 due to anthropogenic forcing and longer memory of the oceans.
The FSLE obtained for the slow-varying component of CIR (Table 5) for all the historical data and the IMD observation data are significantly larger than the FSLE obtained for each of the RCP scenarios. Most of the models suggest a reduction in the value of λs with the increase in radiative forcing. It can be concluded from the results that the predictability of slow-varying component of CIR for RCP scenarios is in the range of six to nine years (Table 7), whereas for historical data and for IMD observational data, the range is between five and six years. Thus, the lower panel of Figure 11 shows that the increase in radiative forcing from 2.6 W m−2 to 8.5 W m−2 will increase the predictability of CIR. Therefore, the CIR under the RCPs will be more predictable than the historical scenarios in long timescales for central India. This increase in the predictability of slow-varying component under RCPs might be because of the increase in the CO2 concentration from RCP2.6 to RCP8.5. This will lead to increased anthropogenic forcing, which could be the possible source of the potential predictability as suggested by Kirtman et al. (2013a), that external forcing like aerosol concentration, greenhouse gases, volcanic eruption, etc. provide predictability to the system, ranging from sub-seasonal to decadal timescales. Using a conceptual forced Lorenz model, Palmer (1993) demonstrated that inclusion of forcing results in increase in the predictability, in the regime residency time of the Lorenz attractor. This means that the dynamical system is confined to a specific region in phase space for a longer time, leading to increased predictability and suggesting that the climate response will be more predictable even in the presence of weak forcing. Dwivedi et al. (2007) and Singh et al. (2015) also used a conceptual forced Lorenz model and found that the forcing causes increase in predictability. This result from conceptual models is supported by studies on real-world data which suggest that external forcing such as aerosols and greenhouse warming may contribute to the long-term predictability (Meehl et al. 2010; Jia & DelSole 2011; Kirtman et al. 2013a). Beyond the monthly timescales, the atmospheric predictability arises from the forcing induced by slowly varying boundary conditions like SST, soil moisture, etc. (Wu et al. 2016).
However, these results should be taken with the caveat that the climate models are still in their evolutionary stage and substantial improvement is still required for better climate change information and improved ability to track the time evolution of the inherent long-scale variability in combination with the forced response.
ACKNOWLEDGEMENTS
The first author is grateful to CSIR, Govt. of India for providing senior research fellowship letter no. 09/001/0399/2016-EMR-I. We acknowledge the World Climate Research Programme's Working Group on Coupled Modeling, which is responsible for CMIP, and we thank the climate modeling groups (listed in Table 1) for producing and making freely available their model output. The authors are grateful to J. Clint Sprott for the discussion on nonlinear dynamical methods. Thanks are also due to IMD, Govt. of India for producing high-resolution gridded rainfall observations.