Improving the predictability of compound dry and hot extremes through complexity science

The definition of dry and hot extremes is based on the 25th and 75th percentile of monthly total precipitation and the average temperature, respectively, at each grid point for a given month. The predictand, CDHE, is defined as the simultaneous occurrence of dry extreme and hot extreme in the same month (see figure S2). The percentiles are calculated using the climate normal from 1991 to 2020, as advised by the World Meteorological Organization (WMO 2017). Predictor events are also defined based on the percentile approach. The index is positive if the 3 month average exceeds the 75th percentile and negative if it falls below the 25th percentile (Hao et al 2020). In our study, we attempted to enhance the severity of these conditions by considering the 10th and 90th percentiles. However, we found that the sample size was insufficient for statistical modeling and previous studies have reported a similar finding (Feng and Hao 2021, Hao et al 2021), highlighting the challenges of sample size restrictions.

ECA proposed by Donges et al (2016) calculates the interdependence between two event series. For a grid point on Land $\left( L \right)$, the time series of Temperature $\left( T \right)$ and Precipitation $\left( P \right)$ are converted into binary series. An event is represented by one, while the absence of an event is indicated by zero at each timestep. Similarly, for each oceanic index $\left( I \right)$ on Sea $\left( S \right)$, the time series is converted into binary series. The time series of $T$, $P$ and $I$ are denoted by ${L^T}\left( {{t_i}} \right)$, ${L^P}\left( {{t_j}} \right)$ and ${S^I}\left( {{t_k}} \right)$, respectively. The observations are taken for each ${t_i},{t_j},{t_k} \in \left[ {1,{ }T} \right]$, where $T$ represents the last timestep of the dataset. The event time series ${Y^T}\left( {{t_i}} \right)$, ${Y^P}\left( {{t_j}} \right)$ and ${Y^I}\left( {{t_k}} \right)$ are defined as:

$$\begin{equation}{Y^T}\left( {{t_i}} \right) = \left\{ {\begin{array}{*{20}{c}} {1,}&{{L^T}\left( {{t_i}} \right) > a} \\ {0,}&{{\text{else }}} \end{array}} \right.\end{equation} \tag{ 1 }$$

$$\begin{equation}{Y^P}\left( {{t_j}} \right) = \left\{ {\begin{array}{*{20}{c}} {1,}&{{L^P}\left( {{t_j}} \right) < b} \\ {0,}&{{\text{else}}} \end{array}} \right.\end{equation} \tag{ 2 }$$

$$\begin{equation}{Y^I}\left( {{t_k}} \right) = \left\{ {\begin{array}{*{20}{c}} {1,}&{{S^P}\left( {{t_k}} \right) < c} \\ {0,}&{{\text{else}}} \end{array}} \right..\end{equation} \tag{ 3 }$$

The percentile cutoff values used for T and P are specific to each grid point, represented by a and b, respectively. Likewise, the percentile cutoff for the Index (I) is specific to each Index, denoted by c (see figure S3). The conditional precursor coincidence rate $\left( {C{P_{{\text{cr}}}}} \right)$ quantifies the number of hot extremes coinciding with low precipitation, conditioned by SST anomaly. We conducted the analysis separately for positive and negative anomalies. The $C{P_{{\text{cr}}}}$ is calculated as follows (Siegmund et al 2017)

$$\begin{align}C{P_{{\text{cr}}}}& = \frac{1}{{{N_T}}}\sum\limits_{i = 1}^{{N_T}} H\left( \sum\limits_{j = 1}^{{N_P}} H\left( {\sum\limits_{k = 1}^{{N_I}} {I_{\left[ {0{ }\Delta T} \right]}}\left( {t_j^P - t_k^I} \right)} \right)\right.\nonumber\\ &\left.\quad\times {I_{\left[ {0{ }\Delta T} \right]}}\left( {t_i^T - t_j^P} \right) \right)\end{align} \tag{ 4 }$$

where ${N_T}$, ${N_P}$ and ${N_I}$ are the number of events in the T, P and I event series, respectively, occurring at times $t_i^T$, $t_j^P$ and $t_K^I$. ${{\Delta }}T$ is the temporal tolerance between the timing of the events. We use an instantaneous coincidence $\left( {{{\Delta }}T = 0} \right)$ for CDHE. The output of the indicator function $I()$ and the Heaviside function $H()$ are defined as follows:

$$\begin{equation}{I_{\left[ {0\;\Delta T} \right]}} = \left\{ {\begin{array}{*{20}{l}} {1,}&{t_j^Pt_k^I = 0} \\ {0,}&{{\text{else}}} \end{array}} \right.\end{equation} \tag{ 5 }$$

$$\begin{equation}H = \left\{ {\begin{array}{*{20}{l}} {1,{\text{ }}\sum\limits_{k = 1}^{{N_P}} {{I_{\left[ {0{\text{ }}\Delta T} \right]}}} \left( {t_j^P - t_k^I} \right) > 0} \\ {0,{\text{ else}}} \end{array}} \right..\end{equation} \tag{ 6 }$$

The value of $C{P_{{\text{cr}}}}$ ranges from zero to one. For example, considering positive anomalies, a value of zero implies that no hot extreme coincided with low precipitation preconditioned by a positive SST anomaly. Conversely, one indicates that all the hot extremes coincided with low precipitation preconditioned by positive SST anomalies. An ensemble of 100 shuffle surrogate event series is generated to assess statistical significance. A significance level of 0.05 is used to reject the null hypothesis that the observed coincidence is purely random (Manoj et al 2023, Singh Raghuvanshi and Agarwal 2023). For more details on ECA, the readers can refer Donges et al (2016).

A network is a set of nodes connected by edges. Mathematically, a network with multiple variables is defined as $G = \left\{ {{X_i},{L_j},E} \right\}$ where ${X_i} = \left\{ {{X_1},{X_2},{X_3}, \ldots .{X_{p - 1}},{X_p}} \right\}$ denotes the set of $i = 1,2, \ldots ,p$ predictor nodes, ${L_j} = \left\{ {{L_1},{L_2},{L_3}, \ldots .{L_{q - 1}},{L_q}} \right\}$ is the set of $j = 1,2, \ldots ,q$ predictand nodes, and $E$ is the set of edges. Two stages are involved in constructing a network: (1) development and (2) pruning. During the development stage, the ${\left[ {C{P_{{\text{cr}}}}} \right]_{{X_i},{L_j}}}$ is calculated between each predictor node $\left( {{X_i}} \right)$ and predictand node $\left( {{L_j}} \right)$. The ${\left[ {C{P_{{\text{cr}}}}} \right]_{{X_i},{L_j}}}$ is done separately for positive and negative anomalies of predictors (see figure 1). Next, $C{P_{{\text{cr}}}} = [ {{{\left[ {C{P_{{\text{cr}}}}} \right]}_{{X_i},{L_j}}}^{{\text{ + ve}}}{ }{{\left[ {C{P_{{\text{cr}}}}} \right]}_{{X_i},{L_j}}}^{ - {\text{ve}}}} ]$ is calculated. In pruning step, an edge-set threshold $\left( \tau \right)$ is applied to filter out weaker connections. The output of the pruning process is represented by an adjacency matrix as (Gadhawe et al 2021)

$$\begin{align} \begin{aligned} {A_{{X_i},{L_j}}}^{{\text{ + ve}}} = \left\{ {\begin{array}{*{20}{l}} {1,{\text{ if }}{{\left[ {C{P_{{\text{cr}}}}} \right]}_{{X_i},{L_j}}}^{{\text{ + ve}}} > \tau } \\ {0,{\text{ otherwise}}} \end{array}} \right.\\ {A_{{X_i},{L_j}}}^{{\text{ - ve}}} = \left\{ {\begin{array}{*{20}{l}} {1,{\text{ if }}{{\left[ {C{P_{{\text{cr}}}}} \right]}_{{X_i},{L_j}}}^{{\text{ - ve}}} > \tau } \\ {0,{\text{ otherwise}}} \end{array}} \right.. \end{aligned} \end{align} \tag{ 7 }$$

Zoom In Zoom Out Reset image size

$\tau$ is defined as the 75th percentile of $C{P_{{\text{cr}}}}$ and this selection is objective rather than subjective. The 75th percentile retains the top 25% connections in the network and creates a well-balanced network structure. Using adjacency matrix, we computed degree $\left( k \right)$ to measure the number of predictors linked to a given predictand node in the network. ${k_j}$ of each predictand node is defined as (Gadhawe et al 2021)

$$\begin{align} \begin{aligned} {k_j}^{ + {\text{ve}}} & = \sum\limits_{i = 1}^n {A_{{X_i},{L_j}}}^{ + {\text{ve}}}\\ {k_j}^{ - {\text{ve}}} & = \sum\limits_{i = 1}^n {A_{{X_i},{L_j}}}^{ - {\text{ve}}}. \end{aligned} \end{align} \tag{ 8 }$$

The range of ${k_j}$ is 0 to the total number of predictors $\left( {12} \right)$. A higher degree indicates that the predictand node is connected to multiple predictors, while a low value suggests that the predictand node has no predictor. The steps involved in the identification of predictors and development of a prediction is illustrated in figure 1.

The obtained predictors for each predictand node are tested for their skill in predicting the occurrence of CDHE using a logistic regression model. CDHE occurrence is a binary variable with two categories: ($Y = 1$ for occurrence and $Y = 0$ for non-occurrence). Two models are used in this study as follows (Agresti 2007):

$$\begin{equation}P\left( {CDHE} \right) = \frac{1}{{\left[ {1 + {e^{\left( { - \left( {{\beta _0} + {\beta _1}{X_1}} \right)} \right)}}} \right]}}\end{equation} \tag{ 9 }$$

$$\begin{equation}P\left( {CDHE} \right) = \frac{1}{{\left[ {1 + {e^{\left( { - \left( {{\beta _0} + {\beta _1}{X_1} + {\beta _2}{X_2}} \right)} \right)}}} \right]}}.\end{equation} \tag{ 10 }$$

The right-hand side term is the logistic function that models the probability of the response variable $\left( {{^{^{\prime}}}CDHE^{\prime}} \right)$ taking the value $1$ given the predictor variables $\left( {{X_1},{X_2}} \right)$. ${\beta _0}$ is the intercept term, and ${\beta _1},{\beta _2}$ are the regression coefficients. Equation (10) is used if the predictand has two predictors, while equation (9) is used for single predictor. The model estimates the regression coefficients by maximizing the likelihood and we use Bayesian Information Criterion (BIC) to prevent overfitting to obtain a parsimonious model (Ben-Ari et al 2018). The significance of the regression coefficients is tested using the Wald test at the 0.05 significance level. If $P\left( {CDHE} \right) > 0.5$, then there is a high probability that the CDHE may occur. Otherwise, the model predicts that CDHE will take the value 0, indicating the non-occurrence. This model is capable of modeling a binary variable and has been widely used in predicting CDHEs (Hao et al 2019, Mukherjee et al 2020). For more details on the model, the readers can refer to Agresti (2007).

10-fold cross-validation is chosen to assess the prediction skill (Gharib and Davies 2021). This scheme is recommended by Rodriguez et al (2010), which concluded that 5-fold or 10-fold is more effective than the leave-one-out approach. Importantly, the uncertainty bandwidth analysis (figure S4), shows that 10-fold results in reduced uncertainty bandwidth compared to 5-fold. To begin with, all samples are randomly shuffled and divided into ten groups. Subsequently, nine groups are used for training, while the hold-out group is used for testing. The process is repeated to ensure that each group is tested at least once. At last, evaluation scores obtained from all ten rounds are summarized and averaged (Yuan et al 2020). The skill is determined using the Brier Skill Score (BSS), and the accuracy of the forecasts are quantified using the using the Critical Succes Index (CSI), False alarm ratio (FAR). The procedure to calculate BSS, CSI, FAR is given in text S1. For more details on the skill assessment procedures the readers can refer Wilks (2011).

The study focuses on the Indian Summer Monsoon (ISM) period and we utilized the Climatic Research Unit (CRU) gridded monthly precipitation and temperature datasets for the Indian mainland (covering 68° E–98° E longitude and 8° N–37° N latitude) from 1951 to 2020 (Harris et al 2020). CRU has been established as reliable in capturing the spatiotemporal variability of precipitation and temperature over India (Bejagam and Sharma 2022). We examined the oceanic indices from three distinct ocean basins (figure 2) based on previous studies since they have a significant association with the Indian monsoon. These include PDO, Niño 3.4, Niño 1.2, and TNI representing the Pacific Ocean (Hassan and Nayak 2020, Sharma et al 2020, Gupta and Jain 2021, Nagaraj and Srivastav 2023); AMO, NAT, SAT, and TASI representing the Atlantic Ocean (Bhatla et al 2016, Malik et al 2017, Borah et al 2020, Sooraj et al 2023); and SWIO, WTIO, SETIO, and DMI representing the Indian Ocean (Sardana et al 2022, Bajrang et al 2023). The monthly time series of these indices were obtained from Extended Reconstructed Sea Surface Temperature version.5 dataset (Huang et al 2017) (see data availability statement).

Zoom In Zoom Out Reset image size

In this study, we aim to do prediction on a seasonal scale (i.e. 2 months in advance). Accordingly, we used the 3 month average SST data of JFM (January–February–March), FMA (February–March–April), MAM (March–April–May), and AMJ (April–May–June) to predict the occurrence of CDHEs in June, July, August, and September, respectively. All the time series were linearly detrended and this is a general practice (Hari et al 2022, Mamalakis et al 2022) to avoid the influence of long-term trends. For discussion, we divided the Indian mainland into ten sub-regions (figure S5) based on previous studies (Guntu et al 2020, 2023, Guntu and Agarwal 2021).

We calculated the $C{P_{{\text{cr}}}}$ considering all the 70 year dataset of oceanic indices and CDHE occurrences. However, in the development of prediction model, we employed 10-fold cross-validation instead of 70 year dataset. Findings with statistical significance at a level of 0.05 are only discussed. In June, positive anomalies in the Niño 3.4 region (figure 3(a)) coincided with the occurrence of CDHEs in Eastern, Western, and North-western India. Similarly, positive anomalies in the other Pacific Ocean indices (figure S7) and DMI (figure 3(i)) were associated with CDHE occurrences in Eastern India. While a negative AMO anomaly coincided with 40% of CDHEs in North-central India (figure 4(e)). Moving to July, a positive Niño 1.2 anomaly (figure S7(f)), coincided with the occurrence of CDHEs in the North-central India and South-central India. Furthermore, positive SETIO anomaly (figure S9(j)) and negative AMO anomaly (figure 4(f)) were linked to CDHEs in South-central India. In the important August month, a positive anomaly in the Pacific Ocean (figures 3(c) and S7) coincided with 40% of CDHEs in Western India and North-western India. Additionally, a positive AMO anomaly (figure 3(g)) was associated with CDHE occurrences over North-eastern India, while positive anomalies in other Atlantic indices (figure S8) were linked to CDHE events in North-western India. Lastly, in September, a positive anomaly in the Atlantic Ocean (figures 3(h) and S8) coincided with the CDHE events in North-western India. A positive DMI anomaly (figure 3(l)) was linked to 40% of CDHEs in North-central India, while negative DMI (figure 4(l)) was associated with CDHEs over North-western India.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our results are in line with the driving mechanism identified by the previous studies. For instance, Niño 3.4 is widely recognized for its global climatic influence, including in India (Gupta and Jain 2021), where positive anomalies indicate warmer SST in the central and eastern equatorial Pacific, triggering atmospheric changes leading to CDHE events. Krishnamurthy and Krishnamurthy (2014), established a robust correlation between the warm phase of the Pacific Ocean and the droughts in Eastern India. The warm temperature triggers changes in the atmospheric circulation patterns, weakening the walker circulation due to the high-pressure system over Darwin (Pandey et al 2020). Similarly, a warm phase in the Niño 3.4 and PDO regions has also been linked to droughts in Central India (Nagaraj and Srivastav 2023) and North-western India (Singh et al 2022). Further, a study by van Oldenborgh et al (2021) highlighted that a warm phase in both the Pacific and Indian Ocean leads to arid conditions in Central India. Additionally, Ratna et al (2021) found a strong correlation between the positive phase of the DMI and droughts in Central India. The warm phase in DMI creates a high-pressure system over WTIO, diverting trade winds towards the SETIO region instead of the Indian Mainland (Shah and Mishra 2020). However, Sardana et al (2022) reported that cold phase of DMI can result in dry conditions in North-western India during the monsoon withdrawal period. Furthermore, previous studies (Borah et al 2020, Singh et al 2021, Sooraj et al 2023) consistently revealed a strong association between the cold phase of the Atlantic Ocean and droughts in Central India. The physical mechanism for this teleconnection has been explained by Goswami et al (2006), where a cold phase of NAO changes wind patterns, leading to a decrease in the meridional gradient of tropospheric temperature, resulting in below normal precipitation.

Complexity science is used in this analysis to identify the important predictors contributing to CDHE. The degree (figure 5) illustrates the number of SST indices linked to each predictand node and figures S13–S16 depicts the network during JJAS. In general, a high degree value for a node indicates that the occurrence of CDHE is influenced by SST variability in multiple oceans. Identifying the most significant predictor linked to CDHE occurrences becomes crucial to reduce the influence of multiple predictors. Conversely, grid points with zero degree do not have any predictors. In June, CDHEs in Eastern India and South-central India were linked to positive (figure 5(a)) and negative (figure 5(b)) anomalies of multiple indices. During July (figure 5(c)), North-central India stood out with a high degree, indicating that CDHEs in this region are connected to positive anomalies of multiple climate indices. Meanwhile, CDHEs in North-western India and North-eastern India are linked to negative anomalies (figure 5(d)). In August, CDHEs in North-western India showed a complex relationship characterized by positive and negative anomalies (figures 5(e) and (f)). In contrast, CDHEs in Eastern India, North-eastern India were associated with positive anomalies (figure 5(e)), while Southern India were primarily linked with negative anomalies only (figure 5(f)). CDHEs during September over Western India, North-western India, and South-central India were driven by both positive and negative anomalies (figures 5(g) and (h)).

Zoom In Zoom Out Reset image size

However, it is important to note that one of the assumptions of the logistic regression model is the independence of predictors. Therefore, using multiple predictors can lead to overfitting of the prediction model (Mamalakis et al 2022). Overfitting results in reduced generalization performance, especially when predictors are interdependent (Meyer et al 2019). To address this concern, it becomes crucial to prioritize the identification of the two most significant predictors that exhibit the strongest association with the occurrence of CDHE. Considering only the top two predictors, we can mitigate the risk of overfitting and ensure the development of robust prediction model. Figure 6 shows the top two predictors. Notably, the Pacific Ocean indices emerge as the most important and second most important predictors for the southern and eastern parts of the country during June and July. In July and September, Atlantic Ocean indices plays an important role in driving CDHEs in North-western India and South-central India. In August, it remains the primary driver for South-central India. Further, indices of Indian Ocean is the primary driver for Eastern India during August, while they are secondary for Central India from July to August.

Zoom In Zoom Out Reset image size

The skill of the proposed predictors is examined using the 10-fold cross validation (see section 2.4). For a grid point with two predictors, the model is chosen based on the low BIC value and considered as the best model. By doing so, the multi-collinearity between the two predictors can be mitigated. A similar procedure has been used in Mukherjee et al (2020) and Mamalakis et al (2022). The statistically significant regression coefficients is shown in figure S17. Notably, it validates the findings in section 3.1, consistently identifying positive and negative relationships between the CDHEs and indices. The prediction skill (BSS) is shown in figure S18. The skill (BSS > 0) illustrates that the proposed predictors can perform better compared to reference prediction. Next, to test the accuracy of the model, the probability of CDHE occurrence was estimated for all the years.

The CSI was computed to gain insight into the accuracy of the proposed predictors compared to stand-alone Niño 3.4 prediction in predicting hits (observed as CDHE). The top panel (figure 7) focuses on the CSI based on Niño 3.4 as a predictor for each month, while the bottom panel displays the CSI based on the proposed predictors. While Niño 3.4 exhibits few grid points with high CSI, the proposed predictors demonstrate higher values for more grid points. For instance, Niño 3.4 prediction in June (figure 7(a)), $CSI > 0.5$ were observed in South-central India only. But utilizing the proposed predictors, higher values are also seen in Eastern India, and North-eastern India. Similarly, in July (figure 7(b)), higher CSI is observed in Central India based on Niño 3.4. Meanwhile with the proposed predictors, higher values are also observed in South-central India. The difference is also prominent in August (figure 7(c)) and September (figure 7(d)). Further, FAR (figure S19) also reveal the accuracy of the negative predictions among all the non-CDHE events. Additionally, the range of skill scores are consistent with the leave-one-out cross validation scheme (figure S20). This comparison highlights the accuracy of the proposed indicators and the ability to provide better prediction across different regions compared to stand-alone Niño 3.4 prediction.

Zoom In Zoom Out Reset image size