Attributing correlation skill of dynamical GCM precipitation forecasts to statistical ENSO teleconnection using a set-theory-based approach

Climate teleconnections are essential for the verification of valuable precipitation forecasts generated by global climate models (GCMs). This paper develops a novel approach to attributing correlation skill of dynamical GCM forecasts to statistical El Niño–Southern Oscillation (ENSO) teleconnection by using the coefficient of determination (R²). Specifically, observed precipitation is respectively regressed against GCM forecasts, Niño3.4 and both of them, and then the intersection operation is implemented to quantify the overlapping R² for GCM forecasts and Niño3.4. The significance of overlapping R² and the sign of ENSO teleconnection facilitate three cases of attribution, i.e., significantly positive anomaly correlation attributable to positive ENSO teleconnection, attributable to negative ENSO teleconnection and not attributable to ENSO teleconnection. A case study is devised for the Climate Forecast System version 2 (CFSv2) seasonal forecasts of global precipitation. For grid cells around the world, the ratio of significantly positive anomaly correlation attributable to positive (negative) ENSO teleconnection is respectively 10.8 % (11.7 %) in December–January–February (DJF), 7.1 % (7.3 %) in March–April–May (MAM), 6.3 % (7.4 %) in June–July–August (JJA) and 7.0 % (14.3 %) in September–October–November (SON). The results not only confirm the prominent contributions of ENSO teleconnection to GCM forecasts, but also present spatial plots of regions where significantly positive anomaly correlation is subject to positive ENSO teleconnection, negative ENSO teleconnection and teleconnections other than ENSO. Overall, the proposed attribution approach can serve as an effective tool to investigate the sources of predictability for GCM seasonal forecasts of global precipitation.

1 Introduction

Precipitation is one of the most important hydrological variables and is an integral part of the global water cycle (Huffman et al., 2007; Ushio et al., 2009; Scofield and Kuligowski, 2003; Schneider et al., 2016; Beck et al., 2019). It plays a key role in driving hydrological processes at catchment, regional and continental scales (e.g., Yuan et al., 2014; Greuell et al., 2018; Robertson et al., 2013; Wu et al., 2014; Lakshmi and Satyanarayana, 2019). Despite the importance, the forecasting of precipitation remains a formidable task due to complex interactions of ocean, atmosphere and land surface processes (Doblas-Reyes et al., 2013; Zhao et al., 2019; Vano et al., 2014; Johnson et al., 2019; Tesfa et al., 2020). Comparing multiple sets of global temperature and precipitation forecasts from the North American Multi-Model Ensemble (NMME) experiment, Becker et al. (2020) highlighted that there have been substantial improvements in temperature forecasts over both land and ocean during the past decades and that there is still plenty of room for improvement of global precipitation forecasts.

Global climate models (GCMs) generate valuable forecasts of worldwide precipitation for hydrological modeling and water management (Kirtman et al., 2014; Doblas-Reyes et al., 2013; Schepen et al., 2020). Nowadays, GCMs have been employed by major climate centers around the world to produce operational climate outlooks (Demargne et al., 2014; Delworth et al., 2020). For example, the Climate Forecast System version 2 (CFSv2) of the US National Centers for Environmental Prediction (NCEP) has been implemented for coupled ocean–atmosphere forecasting since 2011 (Saha et al., 2010); the European Centre for Medium-Range Weather Forecasts (ECMWF) System 5 model became operational in 2017 (Johnson et al., 2019); and the Seamless System for Prediction and Earth System Research (SPEAR) became the next-generation modeling system at the Geophysical Fluid Dynamics Laboratory (GFDL) in 2020 (Delworth et al., 2020). In the meantime, GCM forecasts have been increasingly incorporated into forecasting systems of streamflow, crop yield and soil water, and they are shown to create enormous socioeconomic benefits (e.g., Vano et al., 2014; Peng et al., 2018; Wang et al., 2019).

Climate teleconnections, which are widely used in conventional statistical hydrological forecasting (Wang et al., 2020; Steinschneider and Lall, 2016; Mendoza et al., 2017; Lima and Lall, 2010; Mortensen et al., 2018), are an essential part of assessing the skill of GCM forecasts (Neelin and Langenbrunner, 2013). That is, a number of teleconnection patterns are usually investigated upon the issuance of a new set of GCM forecasts (Delworth et al., 2020; Jia et al., 2015; Kim et al., 2012). For example, El Niño–Southern Oscillation (ENSO) and Madden–Julian oscillation (MJO) have been investigated for CFSv2 forecasts (Saha et al., 2014). North Atlantic Oscillation (NAO), Arctic Oscillation (AO) and Pacific North American (PNA) teleconnections have been assessed for the Global Earth Observing System (GEOS) modeling and data assimilation system (Molod et al., 2020). ENSO and Pacific decadal oscillation (PDO) have been examined for GFDL-SPEAR forecasts (Delworth et al., 2020). It is generally found that skillful GCM forecasts are owing to effective formulations of teleconnection patterns (e.g., Merryfield et al., 2013; Jia et al., 2015; Delworth et al., 2020; Saha et al., 2014; Molteni et al., 2011).

There have been in-depth investigations of ENSO for GCM forecasts as it is one of the most prominent modes of climate variability (Fu et al., 1997; Wang et al., 2003; Feng and Hao, 2021). For example, the influence of ENSO on the East Asian–western Pacific climate was studied for CFSv2 forecasts (Yang and Jiang, 2014) and also for climate projections in the Coupled Model Intercomparison Project Phase 5 (Gong et al., 2015; Kim et al., 2016; Kim and Kug, 2018). The relationship between forecast skill and the state of sea surface temperatures (SSTs) was evaluated for the seasonal outlook of precipitation over the United States, and the skill was found to be dominantly attributed to ENSO in late autumn to late spring (Quan et al., 2006; Pegion and Kumar, 2013; Shin et al., 2019). Understanding the sources for predictability is important for physically based validations of GCM forecasts (Neelin and Langenbrunner, 2013; Manzanas et al., 2014; Shin et al., 2019). Nevertheless, previous studies tended to pay attention to regions subject to prominent ENSO influences (Vashisht et al., 2021; Kim et al., 2016; Rivera and Arnould, 2020). At the global scale, the attribution of GCM precipitation forecasts to ENSO teleconnection is yet to be conducted.

This paper is devoted to attributing correlation skill of dynamical CFSv2 forecasts (e.g., Saha et al., 2010; Yuan et al., 2011; Jia et al., 2015; Becker et al., 2020; Zhao et al., 2020a) to statistical ENSO teleconnection at the global scale. A novel approach based on set theory is devised to facilitate the attribution. With the coefficient of determination (R²) characterizing the ratio of explained variance and the set operations illustrating the overlapping R², the significance test by grid cell is conducted by random permutations in order to identify where correlation skill is attributable to ENSO teleconnection. As will be demonstrated through the case study of CFSv2 forecasts, three cases are effectively revealed: (1) significantly positive anomaly correlation attributable to positive ENSO teleconnection, (2) significantly positive anomaly correlation attributable to negative ENSO teleconnection and (3) significantly positive anomaly correlation not attributable to ENSO teleconnection.

2 Data description

GCM forecasts comprise a typical high-dimensional dataset (Becker et al., 2020; Chen and Kumar, 2016; Kirtman et al., 2014; Saha et al., 2014; Zhao et al., 2020a). For the CFSv2 forecasts investigated in this paper, there are five dimensions: (1) forecast start time s, (2) lead time l, (3) ensemble size n, (4) latitude y and (5) longitude x. s represents the number of months since the benchmark time that is January 1982; l is the number of months ahead, which ranges from 0 to 9 month for the CFSv2 forecasts; n=1, …, 24; i.e., the total number of ensemble members is 24; y ranges from −90 to 90, while x is from 0 to 359, with a horizontal resolution of 1.0^∘ latitude by 1.0^∘ longitude. The set of forecasts is denoted by

$$\begin{array}{}\text{(1)}& F=\left[f_{s,l,n,y,x}\right],\end{array}$$

in which f represents forecast values specified by the five dimensions, and F is the dataset of forecasts.

There are three dimensions for the dataset of observed precipitation corresponding to forecasts (Infanti and Kirtman, 2015; Xie et al., 2007; Schneider et al., 2016). They are target time t, which is equal to the sum of start time s and lead time l to align observations with forecasts; latitude y; and longitude x. The set of observed precipitation corresponding to the forecasts is denoted as

$$\begin{array}{}\text{(2)}& O=\left[o_{t,y,x}\right](t=s+l).\end{array}$$

The NOAA Climate Prediction Center (CPC) global daily Unified Rain-gauge Database (CPC-URD), which has been widely used in the analysis of regional and global precipitation (Chen et al., 2008), is used as the referenced observed precipitation.

Correlation skill, which is in the form of Pearson's correlation coefficient, is calculated so as to relate CFSv2 precipitation forecasts to CPC-URD observations:

$$\begin{array}{}\text{(3)}& r(o,f)={\displaystyle \frac{\sum _k\left(o_k-\stackrel{\mathrm{‾}}{o}\right)\left(f_k-\stackrel{\mathrm{‾}}{f}\right)}{\sqrt{\sum _k{\left(o_k-\stackrel{\mathrm{‾}}{o}\right)}^{\mathrm{2}}}\sqrt{\sum _k{\left(f_k-\stackrel{\mathrm{‾}}{f}\right)}^{\mathrm{2}}}}}.\end{array}$$

In Eq. (3), k represents the target year (k=1982, 1983, …, 2010), where the other dimensions of o and f are omitted for the sake of simplicity. The value of r(o,f) measures how well CFSv2 forecasts correspond to observed precipitation. A significantly positive r(o,f) implies that large (small) forecasts are indicative of large (small) observations, whereas a neutral or negative r(o,f) indicates non-skillful forecasts (Zhao et al., 2020a, b; Yuan et al., 2011).

It is noted that forecasts/observations, which are monthly, are aggregated by season. The aggregation is meant to facilitate the analysis at the seasonal timescale. Attention is paid to the latest forecasts. That is, seasonal precipitation forecasts generated at the beginning of the season are investigated. For example, the December–January–February (DJF) forecasts are generated at the beginning of December. Similarly, seasonal forecasts for March–April–May (MAM), June–July–August (JJA) and September–October–November (SON) are respectively produced at the start of March, June and September.

The concurrent correlation between Niño3.4 and CPC-URD observations is employed to represent ENSO teleconnection (Kim and Kug, 2018; Cai et al., 2009; Steptoe et al., 2018):

$$\begin{array}{}\text{(4)}& r(o,\text{Niño}\mathrm{3.4})={\displaystyle \frac{\sum _k\left(o_k-\stackrel{\mathrm{‾}}{o}\right)\left(\text{Niño}{\mathrm{3.4}}_k-\stackrel{\mathrm{‾}}{\text{Niño}\mathrm{3.4}}\right)}{\sqrt{\sum _k{\left(o_k-\stackrel{\mathrm{‾}}{o}\right)}^{\mathrm{2}}}\sqrt{\sum _k{\left(\text{Niño}{\mathrm{3.4}}_k-\stackrel{\mathrm{‾}}{\text{Niño}\mathrm{3.4}}\right)}^{\mathrm{2}}}}}.\end{array}$$

Niño3.4 is a commonly used index of ENSO (Mason and Goddard, 2001; Vashisht et al., 2021; Chen and Kumar, 2020). The sign of r(o, Niño3.4) indicates the effects of ENSO. A positive r(o, Niño3.4) means that high (low) values of Niño3.4 correspond to large (small) values of observed precipitation; i.e., El Niño events are associated with wet conditions, whereas La Niña events relate to dry conditions. By contrast, a negative r(o, Niño3.4) indicates that high (low) values of Niño3.4 coincide with below-normal (above-normal) precipitation.

3 Methods

3.1 Mathematical formulation

The approach to attributing correlation skill of GCM seasonal forecasts to ENSO teleconnection is built upon the coefficient of determination, i.e., R² (Koster et al., 2010). Mathematically, R² is equivalent to the squared value of Pearson's correlation coefficient r (Krause et al., 2005):

$$\begin{array}{}\text{(5)}& R^{\mathrm{2}}(Y\sim X)=r^{\mathrm{2}}(Y,X).\end{array}$$

There is a difference in the meanings of r in relating observed precipitation to forecasts and ENSO. As for forecasts, r tends to be positive; i.e., high (low) values of forecasts can be indicative of high (low) values of observations (Zhao et al., 2020a, b; Yuan et al., 2011). However, ENSO teleconnection can be either positive or negative. For example, in DJF, positive r(o, Niño3.4) tends to be dominant over southern North America, and negative r(o, Niño3.4) is generally prevalent over northern South America (Mason and Goddard, 2001).

Both positive and negative correlations contribute to R². R² is respectively derived for forecasts and ENSO in relating them to observed precipitation. Specifically, simple linear regression models are set up to regress observed precipitation respectively against GCM forecasts and Niño3.4:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}& o_k={\mathit{\alpha }}_{\mathrm{1}}+{\mathit{\beta }}_{\mathrm{1}}{f}_k+{\mathit{ϵ}}_{\mathrm{1},k}\Rightarrow R^{\mathrm{2}}(o\sim f)=\mathrm{1}-{\displaystyle \frac{\sum _k{\mathit{ϵ}}_{\mathrm{1},k}^{\mathrm{2}}}{\sum _k{\left(o_k-\stackrel{\mathrm{‾}}{o}\right)}^{\mathrm{2}}}},\\ & o_k={\mathit{\alpha }}_{\mathrm{2}}+{\mathit{\beta }}_{\mathrm{2}}\text{Niño}{\mathrm{3.4}}_k+{\mathit{ϵ}}_{\mathrm{2},k}\Rightarrow R^{\mathrm{2}}(o\sim \text{Niño}\mathrm{3.4})=\mathrm{1}-{\displaystyle \frac{\sum _k{\mathit{ϵ}}_{\mathrm{2},k}^{\mathrm{2}}}{\sum _k{\left(o_k-\stackrel{\mathrm{‾}}{o}\right)}^{\mathrm{2}}}},\end{array}\end{array}$$

where β₁ and β₂ are the regression coefficients; ϵ₁ and  ϵ₂ are the residuals; and k represents the target year. Further, through bivariate linear regression, the variance explained by the union of forecasts and Niño3.4 is calculated:

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{o}_k& ={\mathit{\alpha }}_{\mathrm{3}}+{\mathit{\beta }}_{\mathrm{3},\mathrm{1}}{f}_k+{\mathit{\beta }}_{\mathrm{3},\mathrm{2}}\text{Niño}{\mathrm{3.4}}_k+{\mathit{ϵ}}_{\mathrm{3},k}\\ & \Rightarrow R^{\mathrm{2}}(o\sim f\cup \text{Niño}\mathrm{3.4})=\mathrm{1}-{\displaystyle \frac{\sum _k{\mathit{ϵ}}_{\mathrm{3},k}^{\mathrm{2}}}{\sum _k{\left(o_k-\stackrel{\mathrm{‾}}{o}\right)}^{\mathrm{2}}}},\end{array}\end{array}$$

in which the union operator is introduced to represent the joint effect. If GCM forecasts were independent from Niño3.4, then $R^{\mathrm{2}}(o\sim f\cup \text{Niño}\mathrm{3.4})$ could conceptually be obtained by simply adding up R²(o∼f) and R²(o∼Niño3.4). On the other hand, if GCM forecasts were dependent on Niño3.4, then there would be some overlaps for R²(o∼f) and R²(o∼Niño3.4). As a result, $R^{\mathrm{2}}(o\sim f\cup \text{Niño3.4})$ would not be as large as the sum of R²(o∼f) and R²(o∼Niño3.4).

Figure 1 Flowchart of the three steps that attribute significantly positive anomaly correlation of dynamical GCM forecasts to statistical ENSO teleconnection.

According to set theory, the intersection between R²(o∼f) and R²(o∼Niño3.4) is derived by subtracting $R^{\mathrm{2}}(o\sim f\cup \text{Niño}\mathrm{3.4})$ from the sum of R²(o∼f) and R²(o∼Niño3.4):

$$\begin{array}{}\text{(8)}& \begin{array}{rl}{R}^{\mathrm{2}}(o\sim f\cap \text{Niño}\mathrm{3.4})& =R^{\mathrm{2}}(o\sim f)+R^{\mathrm{2}}(o\sim \text{Niño}\mathrm{3.4})\\ & -R^{\mathrm{2}}(o\sim f\cup \text{Niño}\mathrm{3.4}),\end{array}\end{array}$$

in which the intersection operator is introduced to formulate the overlapping R² between forecasts and Niño3.4. Specifically, the value of $R^{\mathrm{2}}(o\sim f\cap \text{Niño}\mathrm{3.4})$ quantifies the overlapping part of the explained variance of observed precipitation accounted for by both GCM forecasts and Niño3.4.

3.2 Attribution of correlation skill

There are three steps to attributing correlation skill of GCM forecasts to ENSO teleconnection. As shown in Fig. 1, the first step is the implementation of three linear regression models to derive R²(o∼f), R²(o∼Niño3.4) and $R^{\mathrm{2}}(o\sim f\cup \text{Niño}\mathrm{3.4})$ so as to derive $R^{\mathrm{2}}(o\sim f\cap \text{Niño}\mathrm{3.4})$. In the absence of a theoretical distribution function for the overlapping R², the significance is tested by randomly permutating GCM forecasts and Niño3.4 1000 times under the null hypothesis that observed precipitation was independent from either GCM forecasts or Niño3.4. In this way, $R^{\mathrm{2}}(o\sim f\cap \text{Niño}\mathrm{3.4})$ is tested to examine whether the overlapping R² between f and Niño3.4 is significant.

Secondly, the significance of correlation coefficient is calculated for GCM forecasts in relating them to global precipitation. Specifically, r(o,f) is obtained for each grid cell over global land. The two-tailed significance test is implemented, and r(o,f) is therefore identified to be significantly positive, neutral or significantly negative. In this paper, the significance levels in the first and second steps are set to be 0.10.

The third step focuses on significantly positive r(o,f) that indicates informative forecasts (Becker et al., 2020; Zhao et al., 2020a, b). There are two criteria: (1) the significance of overlapping R² and (2) the sign of ENSO teleconnection. Overall, three cases are obtained: (1) significantly positive anomaly correlation attributable to positive ENSO teleconnection, (2) significantly positive anomaly correlation attributable to negative ENSO teleconnection and (3) significantly positive anomaly correlation not attributable to ENSO teleconnection.

Figure 2 Box plots of coefficient of determination (R²) for (a) R²(y∼x₁), (b) R²(y∼x₂), (c) $R^{\mathrm{2}}(y\sim x_{\mathrm{1}}\cup x_{\mathrm{2}})$ and (d) $R^{\mathrm{2}}(y\sim x_{\mathrm{1}}\cap x_{\mathrm{2}})$ in the Monte Carlo experiments.

3.3 An illustrative example

An example based on synthetic data is devised to illustrate how the overlapping R² for the response variable y is influenced by the association between two explanatory variables x₁ and x₂. Samples of x₁, x₂ and y are randomly drawn from a trivariate normal distribution

$$\begin{array}{}\text{(9)}& {\left[x_{\mathrm{1}},x_{\mathrm{2}},y\right]}^T\sim N(\mathit{\mu },\mathbf{\Sigma }),\end{array}$$

where μ and Σ are the mean vector and covariance matrix, respectively:

$$\begin{array}{}\text{(10)}& {\displaystyle }& {\displaystyle }{\mathit{\mu }}^T=[\mathrm{0},\mathrm{0},\mathrm{0}],\text{(11)}& {\displaystyle }& {\displaystyle }\mathbf{\Sigma }=\left[\begin{array}{lll}\mathrm{1}& r\left(x_{\mathrm{1}},x_{\mathrm{2}}\right)& \mathrm{0.5}\\ r\left(x_{\mathrm{1}},x_{\mathrm{2}}\right)& \mathrm{1}& \mathrm{0.5}\\ \mathrm{0.5}& \mathrm{0.5}& \mathrm{1}\end{array}\right].\end{array}$$

In the example, the correlations of y with x₁, x₂ are fixed to be 0.5 respectively. As a result, the focus is on r(x₁,x₂) that determines the intersection between x₁ and x₂. The value of r(x₁,x₂) is set to be 0.0, 0.1, 0.2, 0.3, 0.4 and 0.5. For each pre-specified r(x₁,x₂), 1000 samples of x₁, x₂ and y are drawn to facilitate linear regression models to derive R²(y∼x₁), R²(y∼x₂), $R^{\mathrm{2}}(y\sim x_{\mathrm{1}}\cup x_{\mathrm{2}})$ and $R^{\mathrm{2}}(y\sim x_{\mathrm{1}}\cap x_{\mathrm{2}})$. For R², the median and interquartile ranges are estimated through 1000 Monte Carlo experiments.

Figure 3 Venn diagrams for the variance of y explained by x₁ and x₂ given six levels of correlation between x₁ and x₂.

Figure 2 illustrates the influences of r(x₁,x₂) on R²(y∼x₁), R²(y∼x₂), $R^{\mathrm{2}}(y\sim x_{\mathrm{1}}\cup x_{\mathrm{2}})$ and $R^{\mathrm{2}}(y\sim x_{\mathrm{1}}\cap x_{\mathrm{2}})$. Figure 2a and b show that the median values of R²(y∼x₁) and R²(y∼x₂) approximate the squared value of pre-specified r(y,x₁) and r(y,x₂). That is, they remain around 0.25, i.e., 0.5², as r(x₁,x₂) increases from 0.0 to 0.5. The indication is that the change in correlation between forecasts and Niño3.4 does not influence the amount of information that they respectively provide for observed precipitation. By contrast, Fig. 2c shows that $R^{\mathrm{2}}(y\sim x_{\mathrm{1}}\cup x_{\mathrm{2}})$, which represents the ratio of variance explained by the union of x₁ and x₂, decreases with the increase of r(x₁,x₂). This phenomenon coincides with the increase of $R^{\mathrm{2}}(y\sim x_{\mathrm{1}}\cap x_{\mathrm{2}})$ with r(x₁,x₂).

Figure 4 Correlation coefficients between observed precipitation in DJF with (a) 0-month lead CFSv2 forecasts generated at the beginning of December and (b) concurrent Niño3.4. Significance tests of correlation for (c) CFSv2 forecasts and (d) Niño3.4. Coefficient of determination (R²) for the regression of observed precipitation against (e) CFSv2 forecasts, (f) Niño3.4, (g) the union of CFSv2 forecasts and Niño3.4 and (h) the intersection of CFSv2 forecasts and Niño3.4.

Figure 3 further shows the influence of r(x₁,x₂) using a Venn diagram that illustrates the extent to which R²(y∼x₁) and R²(y∼x₂) intersect. The intersection is represented by the overlapping area between the two circles. From this figure, it can be seen that it is the correlation between x₁ and x₂ that leads to the decrease of $R^{\mathrm{2}}(y\sim x_{\mathrm{1}}\cup x_{\mathrm{2}})$ and the increase of $R^{\mathrm{2}}(y\sim x_{\mathrm{1}}\cap x_{\mathrm{2}})$. For global precipitation forecasting, the intersection reflects the overlapping R² for GCM forecasts and ENSO teleconnection. Figure 3 suggests that the correlation between GCM forecasts and Niño3.4 would lead to a decrease of total information and an increase of overlapping information.

Figure 5 (a) Spatial distribution of grid cells with significant overlapping R² for CFSv2 forecasts generated in December and Niño3.4 in DJF. Scatter plots of anomaly correlation against ENSO teleconnection for grid cells (b) with significant overlapping R² and (c) with non-significant overlapping R². Grey points in panel (c) are marked in white in panel (a).

4 Results

4.1 Correlation skill and ENSO teleconnection in DJF

Global maps of correlation skill and ENSO teleconnection in DJF, which is the peak season of ENSO, are shown in the upper part of Fig. 4. In Fig. 4a, correlation skill is observed to be largely positive, indicating that CFSv2 forecasts are skillful in general (Saha et al., 2010). In Fig. 4b, ENSO teleconnection exhibits both positive and negative values. That is, observed precipitation around the world can be positively or negatively correlated with Niño3.4 (Mason and Goddard, 2001; Vashisht et al., 2021; Chen and Kumar, 2020). The two-tailed significance test is applied to anomaly correlation and ENSO teleconnection at each grid cell. Figure 4c illustrates that correlation skills of CFSv2 forecasts are significantly positive over extensive areas around the globe. Moreover, Fig. 4c is observed to correspond to Fig. 4d to some extent – significantly positive correlation appears over southern North America and East Africa in both Fig. 4c and d. In addition, significantly positive anomaly correlation (Fig. 4c) corresponds to significantly negative ENSO teleconnection (Fig. 4d) in northern South America, southern Africa and southeast Asia.

The results of simple linear regressions that lay the foundation for the attribution analysis are shown in Fig. 4e to h. Figure 4e and f, which are respectively for CFSv2 forecasts and Niño3.4, respectively conform to Fig. 4a and b. This outcome is due to R² being mathematically equal to the squared value of the correlation coefficient (Krause et al., 2005). The union in Fig. 4g exhibits a higher value of R² than that in either Fig. 4e or f. The subtraction of the union from the sum facilitates the intersection. As illustrated in Fig. 4h, dark blue grid cells are seen to distribute in southern North America, northern South America, East Africa, southern Africa and southeast Asia. Over these regions, both GCM forecasts and the Niño3.4 index can explain a considerable ratio of the variance of observed precipitation (Fig. 4e and f). More importantly, their explained variances intersect (Fig. 4h).

Figure 6 (a) Spatial distribution of grid cells for the three cases attributing anomaly correlation of CFSv2 forecasts to ENSO teleconnection. Scatter plots of anomaly correlation against ENSO teleconnection for grid cells (b) with significant overlapping R² and (c) without significant overlapping R². Grey points in panels (b) and (c) are marked in white in panel (a). Four grid cells are selected and labeled with capital letters from panel (a) for the subsequent analysis of the attribution of CFSv2 forecasts to ENSO teleconnection.

4.2 Attribution at the global scale

The significance of the overlapping R² for CFSv2 forecasts and Niño3.4 is shown in Fig. 5. In Figure 5a, grid cells with significant overlapping R² are marked in orange. The corresponding anomaly correlation and ENSO teleconnection, which are respectively obtained from Fig. 4a and b, are illustrated using the scatter plot in Fig. 5b. The scatter points tend to fall towards the upper right and left corners of the plot. The implication is that both anomaly correlation and ENSO teleconnection ought to be large enough to facilitate a significant intersection. Largely owing to overlapping R², anomaly correlation is observed to increase with the increase of positive ENSO teleconnection and also with the decrease of negative ENSO teleconnection. For Fig. 5b, there are notably some outliers that suggest ENSO teleconnection could contribute to negative anomaly correlation. CFSv2 forecasts are generally wrong in these cases, and the potentially problematic grid cells are marked in black in Fig. 5a. Further, the scatter plot in Fig. 5c is for grid cells where overlapping R² is non-significant. For a fair number of grid cells, anomaly correlation can rise above 0.50, but ENSO teleconnection remains at nearly 0.00. The implication is that the corresponding anomaly correlation is not relevant to ENSO teleconnection.

Figure 7 Sunburst diagram of the attribution of significantly positive anomaly correlation to ENSO teleconnection for CFSv2 forecasts in DJF.

Figure 6 presents the three cases of attribution of significantly positive anomaly correlation. The red color marks grid cells where significantly positive anomaly correlation is attributable to positive ENSO teleconnection. Some corresponding grid cells are observed in regions of known positive ENSO teleconnections, such as southern North America (Strazzo et al., 2019) and equatorial eastern Africa (Vashisht et al., 2021); and some are in less-investigated regions, such as parts of central, southern and eastern Asia. The green color indicates grid cells where significantly positive anomaly correlation is attributable to negative ENSO teleconnection. They appear in northern South America and southern Africa, where negative ENSO teleconnection is known to exist (Cai et al., 2020; Howard et al., 2019) and also in parts of Far East and Alaska. There are also grey areas where significantly positive anomaly correlation is not attributable to ENSO teleconnection. The corresponding grid cells are generally located in Europe, North Asia, northwestern Africa and southern Australia. Therein, skillful forecasts can relate to teleconnections other than ENSO, such as AO and NAO (Minami and Takaya, 2020).

Figure 8 Scatter plots of the relationships between normalized observation, forecast and Niño3.4 and Venn diagrams of the union and intersection operations for four selected grid cells.

Figure 7 presents a sunburst diagram that quantifies the percentages of significantly positive anomaly correlation and its attribution results. As shown by the central cycle, anomaly correlation is identified to be significantly positive, neutral or significantly negative (Zhao et al., 2020a, b). The grey segment suggests that more than half of grid cells around the globe are of neutral anomaly correlation, indicating that GCM precipitation forecasts still have plenty of room for improvement (Delworth et al., 2020; Jia et al., 2015; Kim et al., 2012). The pink segment indicates that 39.4 % of the grid cells exhibit significantly positive anomaly correlation. There are three cases of attribution for significantly positive anomaly correlation. The results are shown by the extended segments, of which the color scheme is the same as that of Fig. 6. It can be seen that significantly positive anomaly correlation is attributable to positive (negative) teleconnections for 10.8 % (11.7 %) of grid cells around the globe. Furthermore, significantly positive anomaly correlation is not attributable to ENSO teleconnection for 16.9 % of grid cells around the globe.

4.3 Attribution for selected grid cells

Four grid cells are selected from Figure 6a to showcase the attribution of CFSv2 forecasts to ENSO teleconnection. As is illustrated in the Methods section, there are three variables under investigation. They are observed precipitation, forecast precipitation and Niño3.4. Their values are normalized, i.e., subtracting the mean and dividing by the standard deviation, to facilitate the intercomparison of results across the four grid cells. In Fig. 8, the first column plots normalized observation against normalized forecast; the second column plots normalized observation against normalized Niño3.4; and the third column plots normalized forecasts against normalized Niño3.4. Furthermore, the last column of Fig. 8 presents the Venn diagrams illustrating the set operation of intersection.

Grid cell A shown in the first row of Fig. 8 presents the case of significantly positive anomaly correlation attributable to positive ENSO teleconnection. The coordinates of the grid cell are 38^∘ N, 77^∘ W. In southern North America, DJF precipitation is known to be modulated by the ENSO-induced Pacific North American (PNA) pattern (Jong et al., 2021). Specifically, in southern North America, PNA tends to cause an enhanced DJF Pacific jet stream that extends further east than normal during El Niño events, and there are nearly reversed patterns during La Niña events. Since the jet stream determines the paths of DJF storms, the PNA pattern enables ENSO to affect precipitation in southern North America. It can be seen that forecasts exhibit a high correlation with Niño3.4. This result suggests that CFSv2 forecasts can reasonably represent the influence of ENSO. As a result, R² explained by forecasts and Niño3.4 largely overlap in grid cell A.

Grid cell B shown in the second row of Fig. 8 is for the case of significantly positive anomaly correlation attributable to negative ENSO teleconnection. Its coordinates are 5^∘ N, 60^∘ W. In DJF, there is a negative ENSO teleconnection over northern South America; it is owing to ENSO-related SSTs driving changes in the climatological Walker circulation that promotes anomalous descending (ascending) motion and contributes to negative (positive) precipitation anomalies in El Niño (La Niña) events (Kayano and Andreoli, 2006). The high correlation between forecasts and Niño3.4 highlights the effectiveness of CFSv2 in capturing the negative ENSO teleconnection. There is a considerable intersection between the variability explained by forecasts and Niño3.4 in grid cell B. It is noted that some similarities in the Venn diagrams are observed for grid cells A and B, at which ENSO teleconnection is respectively positive and negative. Specifically, r(o, Niño3.4) is respectively 0.46 and −0.48 at grid cells A and B. R² that is the focus of the Venn diagrams mathematically represents the squared value of the correlation coefficient. Therefore, the similar Venn diagrams highlight that both positive and negative ENSO teleconnection can contribute to correlation skill.

Grid cell C displayed in the third row of Fig. 8 is for the case of significantly positive correlation not attributable to ENSO teleconnection. Its coordinates are 49^∘ N, 2^∘ E. The remote influence of ENSO on Europe has pathways through the North Atlantic or Arctic regions, including the tropospheric and stratospheric bridges (Butler et al., 2014). However, the amplitude of ENSO impacts is weak and generally not significant in the European region (Butler et al., 2014). In addition, DJF precipitation in Europe is known to be modulated by the NAO (Greuell et al., 2018), the Eurasian snow cover extent and the Quasi-biennial Oscillation (Butler et al., 2014). It can be seen that while observed precipitation shows a neutral correlation with Niño3.4, CFSv2 forecasts explain a substantial fraction of observed precipitation variability. This result indicates the capability of CFSv2 in capturing teleconnection patterns other than ENSO.

Grid cell D shown in the last row of Fig. 8 represents the case of neutral anomaly correlation. The coordinates are 40^∘ N, 116^∘ E. Over East Asia, precipitation is known to be influenced by the response of Rossby waves to ENSO (Yang et al., 2018). Also, winter monsoon activities in East Asia are profoundly influenced by wind–SST–evaporation feedbacks over tropical central Pacific to northwestern Pacific (Kim and Kug, 2018). It can be observed that there is a moderate but not significant correlation between observed precipitation and Niño3.4. However, CFSv2 forecasts exhibit a neutral anomaly correlation, suggesting that the information of ENSO teleconnection is not represented in the forecasts.

4.4 Extended analysis of the other seasons

The attribution analysis is further extended to the other seasons, i.e., MAM, JJA and SON. Global maps of the three cases of attribution are illustrated by season in Fig. 9 and also in Figs. S1 to S6 in the Supplement. Overall, the results of attribution vary considerably across the four seasons. This is generally owing to the fact that ENSO teleconnection varies by season (Kim and Kug, 2018; Steptoe et al., 2018; Wang et al., 2019) and that GCMs formulate not only ENSO, but also other teleconnections (Saha et al., 2014; Jia et al., 2015; Delworth et al., 2020). Overall, the percentage of significantly positive anomaly correlation is 27.0 %, 24.0 % and 34.6 % respectively in MAM, JJA and SON, which tends to be smaller than that in DJF. This result can be due to the seasonal cycle of ENSO; i.e., ENSO forcing tends to be the strongest in DJF, and it translates into weaker precipitation variability in the other seasons (Yang et al., 2018).

Figure 9 Spatial maps (a–c) and sunburst diagrams (d–f) of the three cases of attribution of significantly positive anomaly correlation in MAM (a, b), JJA (c, d) and SON (e, f).

The percentage of significantly positive anomaly correlation attributable to positive ENSO teleconnection is respectively 7.1 %, 6.3 % and 7.0 % in MAM, JJA and SON. Representative regions for this case are the western United States in MAM (Pegion and Kumar, 2013), parts of South America in JJA (Cai et al., 2020) and the Middle East in SON (Mariotti, 2007). In total, 7.3 %, 7.4 % and 14.3 % of grid cells are with significantly positive anomaly correlation attributable to negative ENSO teleconnection respectively in MAM, JJA and SON. One representative region is southeast Asia, where precipitation is strongly correlated with ENSO in MAM, JJA and SON (Jiang and Li, 2017); also, in Australia, precipitation in SON is found to be substantially influenced by the extratropical teleconnection pathway of ENSO (Cai and Weller, 2013). Furthermore, 12.6 %, 10.3 % and 13.3 % of grid cells are with significantly positive anomaly correlation not attributable to ENSO teleconnection. This result calls for the investigation of other teleconnection patterns for GCM seasonal precipitation forecasts.

5 Discussion

Correlation skill between forecast and observed precipitation is one of the most important indicators of the usefulness of GCM forecasts (Yuan et al., 2011; Becker et al., 2020; Vano et al., 2014; Johnson et al., 2019; Zhao et al., 2020b). To facilitate forecast applications, correlation skill is conventionally calculated from data and then presented using spatial plotting (Zhao et al., 2020a, b). Focusing on the relationship between correlation skill and ENSO teleconnection, the present paper highlights that significantly positive anomaly correlation, which is always advantageous for practical applications of GCM forecasts (Vano et al., 2014; Yuan et al., 2014; Peng et al., 2018), can be attributed to positive (negative) ENSO teleconnection or not to ENSO. In DJF, significantly positive anomaly correlation for CFSv2 forecasts is attributable to positive ENSO teleconnection in southern North America and East Africa, and it is attributable to negative ENSO teleconnection in northern South America and southern Africa. Moreover, significantly positive anomaly correlation in Europe can be attributable to teleconnections other than ENSO. The different cases of attribution also exist for MAM, JJA and SON, but their spatial extents vary considerably (Figs. S4–S6). These results conform to previous findings that regions exhibiting positive (negative) ENSO teleconnection change substantially by season (Mason and Goddard, 2001; Vashisht et al., 2021; Chen and Kumar, 2020) and that performances of GCM forecasts vary by season (Vano et al., 2014; Johnson et al., 2019; Zhao et al., 2020b).

Lagged ENSO teleconnection and correlation skill by lead time are important issues for seasonal forecasting (Schepen et al., 2012; Peng et al., 2014; Steinschneider and Lall, 2016). A numerical experiment of attribution is performed for Niño3.4 indices at different time lags. Specifically, for precipitation in DJF, the concurrent Niño3.4 in the analysis is replaced by Niño3.4 in November, October and September so as to investigate the overlapping R² for 1-, 2- and 3-month lag ENSO teleconnection. Figures S7 to S10 show that the results tend to be similar at the three lags; the similarities are generally owing to the temporal persistency of Niño3.4 (Yang et al., 2018). Furthermore, another experiment is devised for GCM forecasts at the lead time of 1, 2 and 3 months. That is, for precipitation in DJF, forecasts generated in November, October and September are used to replace forecasts generated in December in the analysis (Figs. S11 to S14). It can be observed that the case of significantly positive anomaly correlation attributable to positive ENSO teleconnection remains for southern North America and that the case of significantly positive anomaly correlation attributable to negative ENSO teleconnection remains for northern South America and southern Africa.

The capability to formulate climate teleconnections is an essential part in the evaluation of GCM forecasts (Jong et al., 2021; Molod et al., 2020). Adding to previous studies that investigated GCM forecasts for regions subject to prominent ENSO influences (Jha et al., 2016; Manzanas et al., 2014; Pegion and Kumar, 2013), this paper presents an investigation of ENSO teleconnection at the global scale. For grid cells around the world, the ratio of significantly positive anomaly correlation attributable to positive (negative) ENSO teleconnection is respectively 10.8 % (11.7 %) in DJF, 7.1 % (7.3 %) in MAM, 6.3 % (7.4 %) in JJA and 7.0 % (14.3 %) in SON. Furthermore, the ratio of significantly positive anomaly correlation not attributable to ENSO teleconnection, which suggests that other climate teleconnections are at play in determining the skill of GCM forecasts, is respectively 16.9 %, 12.6 %, 10.3 % and 13.3 % in DJF, MAM, JJA and SON. Overall, the spatial plots and the attribution results can serve as a reference for further investigations of the effects of ENSO teleconnections and other climate patterns on the predictive performance of GCM forecasts.

6 Conclusions

Climate teleconnections, in particular ENSO, have been extensively used in conventional statistical hydrological forecasting. This paper is devoted to investigating the relationship between statistical ENSO teleconnection and correlation skill of dynamical CFSv2 forecasts. A novel mathematical approach is built upon set theory by making use of the coefficient of determination (R²) that measures the ratio of explained variance to total variance. Specifically, taking advantage of simple linear regression, the ratios of variance explained by GCM forecasts, Niño3.4 and their union are respectively obtained; then, the overlapping R² for GCM forecasts and Niño3.4 is derived based on the intersection operation. Based on the significance of overlapping R² and the sign of ENSO teleconnection, three cases of attribution are derived. They are significantly positive anomaly correlation attributable to positive ENSO teleconnection, attributable to negative ENSO teleconnection and not attributable to ENSO teleconnection. The effectiveness of the developed approach is demonstrated through the case study of CFSv2 seasonal forecasts of global precipitation. Spatial plots of the attribution are illustrated by season. The spatial patterns of forecast skill attributed to different types of ENSO teleconnections confirm previous studies associating seasonal precipitation variability with ENSO and highlight the capability of CFSv2 in capturing the pathways of ENSO teleconnections. The attribution method proposed in this paper can lay a basis for future evaluations of other teleconnections and investigations of predictability sources for GCM seasonal precipitation forecasts.