Using Co-Behavior Analysis to Interrogate the Performance of CMIP5 GCMs over Southern Africa

As established in earlier research, analysis of the combined roles (co-behavior) of multiple climate processes provides useful insights into the drivers of regional climate variability, especially for regions with no singular large-scale circulation control. Here, we extend the previous study in order to examine the performance of eight models from phase 5 of the Coupled Model Intercomparison Project (CMIP5) in representing co-behavior influence on surface expressions over southern Africa. We find that although models broadly simulate observed precipitation responses over southern Africa, they fail to produce statistically strong response signals for an important drought pattern (El Niño co-behaving with positive Antarctic Oscillation during summer) for the region. We also demonstrate that the models show statistically strong temperature response signals to co-behavior that agree well with observed responses over the region. The multimodel ensemble mean although consistent with observations shows a larger spread. By elucidating the performance of models in representing observed co-behavior of climate processes, we are able to evaluate models while establishing important information for understanding of climate variability.


Introduction
Understanding the drivers of regional-scale climate variability and change is key to understanding and interrogating climate projections for regions, which in turn provides the foundation for informing increasingly urgent adaptation decision making. Studies have often focused on individual processes in attempts to understand regional climate variability and change by establishing statistical relationships with climate variables. While this approach has improved our understanding of regional climate to an extent, not much has been done to understand the response to multiscale interactions among climate processes driving regional climate variability and change. As such alternative approaches are required to help interpret the nature of these interactions and their impact on the present and future climate (Daron et al. 2019).
Interactions among climate processes, or co-behavior, can influence regional weather and climate. Co-behavior recognizes that regional climate variability is an outcome of multiple, potentially interacting climate processes operating across varying temporal and spatial scales, and a better understanding of co-behavior may advance our understanding of regional climate variability. Quagraine et al. (2019) described an approach to identifying and analyzing co-behavior in observation-based and climate reanalysis data. For the purposes of improving understanding of regional climate variability there may be a need to extend this approach to global climate models (GCMs).
GCMs are our primary tool to understand past and future changes in our climate (Frame and Stone 2013; Denotes content that is immediately available upon publication as open access. Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/JCLI-D-19-0472.s1. a ORCID: 0000-0002-7887-6040. Knutti and Sedlá cek 2013). The ability of GCMs to reasonably simulate already identified co-behavior modes in observed datasets may permit climate model weighting and the evaluation of climate model simulations by making available a consistent set of GCMs pertaining to any region of interest.
Over time, we have significantly improved how GCMs are able to simulate the observed climate and to some extent make robust future predictions in the presence of uncertainties through improving our understanding of the physics and by refinement to model configuration and numerical algorithms (Miao et al. 2014;Ramirez-Villegas et al. 2013;Su et al. 2013). Experiments like phase 5 of the Coupled Model Intercomparison Project (CMIP5; Taylor et al. 2012) have certainly helped in this direction. Few studies have used CMIP5 models in tackling a wide range of climate issues over southern Africa. The studies include but are not limited to the modeling of the present and future African climate (Dike et al. 2015), exploration of climate processes associated with rainfall bias (Munday and Washington 2018), and analyzing how models reproduce teleconnections between El Niño-Southern Oscillation (ENSO) and southern African rainfall (Dieppois et al. 2015) and how models have identified drivers of rainfall decline over southern Africa in recent past and near-future projections (James and Washington 2013).
However, co-behavior of climate processes in climate models remain weakly explored, yet may offer a valuable source of information to understand regional climate variability as well as aid in assessing the robustness of climate model results. Few recent studies have directly focused on the role of interactions between climate processes on surface manifestations of climate. Muñoz et al. (2015) investigated the relevance of multiple climate drivers across multiple time scales to improve extreme rainfall predictability. These studies used observation and reanalysis datasets. A study by Muñoz et al. (2017) used five members each of two versions of a GCM to develop a framework for understanding the influence of cross-timescale interactions on weather types. Daron et al. (2019) also used related theoretical concepts in highlighting the relationship that exists among climate processes.
In this paper we examine how a set of GCMs from the CMIP5 experiment response to observed co-behavior over southern Africa (108-358S, 58-458E). We should mention that the scope of the current study does not explore the dynamical processes behind these responses. The rest of the paper is organized as follows: section 2 summarizes the datasets and methodology adopted for the study. The results are discussed in section 3 with a summary and conclusions presented in section 4.

a. Data
ERA-Interim data (Dee et al. 2011) at 700-hPa geopotential height with a grid resolution of 0.758 are used for the self-organizing map (SOM; Kohonen 1982Kohonen , 2001) classification of circulation patterns over southern Africa [see Quagraine et al. (2019) for details of the approach].
In the present study, we select only eight models (Table 1) from the set of CMIP5 GCMs based on the availability of climate variables needed for the study and the ability of the models to reasonably represent the southern African climate (see McSweeney et al. 2015).
For observed temperature, we use the Climatic Research Unit (CRU-TS v4.01; Harris et al. 2014) data, and Climate Hazards Infrared Precipitation with Stations (CHIRPS;Funk et al. 2015) for precipitation.

1) CLIMATE INDICES
We calculate three indices that describe and analyze the state of changes occurring in ENSO, Antarctic  Oscillation (AAO), and the intertropical convergence zone (ITCZ) intensity for each of the CMIP5 models used (see Table 1). Here, the Niño-3.4 index (Trenberth 1997;Trenberth and Stepaniak 2001), is computed using sea surface temperature anomalies, and used to characterize ENSO, while the AAO index is constructed by projecting daily 700-hPa height anomalies poleward of 208S onto the leading pattern of the AAO as described by Thompson and Wallace (2000). The tropical rain belt index (TRBI), which captures the intensity of the ITCZ and is based on methodology from (Nikulin et al. 2012), is generated. For higher rainfall intensities within the tropical rain belt, TRBI is positive. These above mentioned processes are chosen because of their importance to the southern African regional climate (see Klutse et al. 2016;Lennard and Hegerl 2014;Meque and Abiodun 2015;Reason and Rouault 2005;Weldon and Reason 2014). We stress here that other processes could have equally been used but we use these processes as they are well understood in literature and conveniently easy to compute their indices for the study.

2) SELF-ORGANIZING MAP
In this study, we identify archetype synoptic circulation patterns over the region by training a 12-node selforganizing map (SOM; Kohonen 1982Kohonen , 2001 with daily 700-hPa ERA-Interim geopotential anomaly fields to produce 12 characteristic circulation patterns over the study period 1980[Fig. 3 in Quagraine et al. (2019, hereafter Fig. S1 in the online supplemental material]. At 700-hPa geopotential height level the important climate systems such as subtropical high and low pressure, easterly and westerly waves around southern Africa are effectively captured (e.g., Bartman et al. 2003).
A SOM is a topologically sensitive clustering technique that uses unsupervised training to cluster training data (Kim and Seo 2016;Lennard and Hegerl 2014). The SOM analysis typically consists of two main processes: the training phase and the mapping phase. The technique as applied in synoptic climatology is particularly effective as it maintains clusters as a continuum that in turn preserves relationships between weather states and presents the output as a readily understood and visualized array of classified patterns (Hewitson and Crane 2002;Kim and Seo 2016;Lee 2017;Rousi et al. 2015;Sheridan and Lee 2012). The classified patterns can aid in developing an understanding of linkages between largescale regional circulation and local weather expressions (Cassano et al. 2015;Gibson et al. 2017;Hewitson and Crane 1996) and have been shown to capture aspects of interannual rainfall variability (Wolski et al. 2018) and rainfall extremes (Lennard and Hegerl 2014). The implementation of the SOM in this study is shown in Fig. 1a.
The choice of the SOM size is a compromise between, on the one side, the need to generalize circulation field structure for analytical purposes, and on the other side, the need to capture the dominant spatial features of that field that have implications to local climate. Mackellar et al. (2010) and Tadross et al. (2005) found the 12-node SOM size captures the generalized synoptic circulation patterns over southern Africa, and we have used this SOM size in earlier analysis (Quagraine et al. 2019). For this analysis, the benefit of reducing variance to representative time series outweighs the temporal and spatial information loss in using this SOM size, thus we use the 12-node SOM (see Wolski et al. 2018). Here, we first use daily 700-hPa anomaly field of the reanalysis to train the 12-node SOM, where each day in the study period is mapped on to one of the 12 circulation patterns identified by the SOM [ Fig. 1a(i)]. Each model is then mapped through the reanalysis SOM space as shown in Fig. 1a(ii). The purpose of mapping the models to the reanalysis SOM space is to allow for comparison of SOM node frequencies between the GCMs and ERA-Interim representation of the observed frequencies and to further provide comparable inputs to the principal component analysis [PCA; Fig. 1b(i)]. The implicit assumption here is that the GCM patterns fall within the data space represented by the reanalysis, that is, that the archetype synoptic patterns of the ERA-Interim data encompasses the range of synoptic states in the GCMs. To test this assumption, we calculate the mean error of each day's pattern when compared to the reference pattern of the SOM node to which it maps. By comparing the error mapping of the GCMs to that of the ERA-Interim data, one can assess whether the GCM patterns have a greater error than the reanalysis, which would indicate that the GCM is representing synoptic patterns that go out of the data space of the ERA-Interim data. In all cases, except some very few nodes in one model, the GCM patterns have errors smaller than the ERA-Interim data mapped to the ERA-Interim trained SOM as shown in Table 2. This is in some ways not unexpected as the GCMs are a reduced complexity system and hence would be expected to have less variance than the ERA-Interim data, which is constrained by observations. Further, it indicates no significant departure by the GCMs from the archetype patterns determined from the ERA-Interim data and thus allows us to proceed to examine the differences in the frequency of occurrence, which is our initial intention. As a next step, we generate frequency of occurrence histograms for the reanalysis and all eight models. A monthly time series matrix of each synoptic type is then generated from the daily frequency of occurrence and then from this a 3-month centered moving average window that serves as a low-pass filter is used to eliminate short-term variability while highlighting long-term variability (see Quagraine et al. 2019). This seasonal grouping aids in analyzing seasonal variability of the SOM.

3) ROTATED PCA
To explore the frequencies of occurrence of the synoptic circulation types in relation to the co-behavior of the conditioning large-scale drivers (represented by climate indices), we use PCA (Jolliffe 2012;Jolliffe and Cadima 2016;Lever et al. 2017), shown in Fig. 1b. In Fig. 1b(i), the input to the PCA is an extended matrix of seasonal frequencies with additional columns for the climate indices for ENSO, AAO and TRBI derived from the ERA-Interim data (hereafter, training PCA is a multivariate technique that produces a set of abstract variables known as principal components (PCs), which are weighted linear combinations (coefficients) of the original variables (Jolliffe and Cadima 2016) and these coefficients are then stored in a loading matrix with the goal of finding the best summary of the data using limited number of PCs (Lever et al. 2017). The data are rotated by the loading matrix (also known as a rotating matrix) in a way that the projection with the greatest variance maps along the first axis (first PC). By eliminating PCs that represent only small amounts of the total original variance, the complexity of a highdimensional data can be simplified while retaining trends and patterns by transforming the data into fewer dimensions, which are summaries of features and has previously been used as a dimension reduction approach in distribution modeling (Lever et al. 2017;Robertson et al. 2001). Hence the output of our PCA, as applied to our extended matrix, identifies dominant modes of independent variability within our data and three PCs were retained based on randomization and assessment of their significance at 90% confidence level using the N-Rule test as a stopping rule (Peres-Neto et al. 2005). Through the PCA loadings, we are able to establish the relationship between the frequency of occurrence of the synoptic circulation types and the conditioning largescale processes. A modified version of scikit-learn Python module that incorporates a function for varimax rotation is used for the PCA (https://scikit-learn.org/; Pedregosa et al. 2011).

4) GCM ANALYSIS
To analyze how the 8 GCMs represent the observed relationships between synoptic variability and largescale climate indices we first regrid the GCMs to the same resolution as the ERA-Interim dataset using nearest-neighbor interpolation. For the period of study, each GCM is a combination of the historical period  and the future scenario (2006-13) for representative concentration pathway (RCP) 8.5, as this was the likely pathway at that period in time. We then combine the time series of the indices calculated earlier with the SOM node frequencies representing the synoptic circulation state frequencies for each GCM, and these data vectors are transformed (using the fit_transform function from scikit-learn module) through the original reanalysis derived PCA loadings to produce PCA score time series for each GCM. This sort of approach has been frequently and effectively applied in climate science and the reader is referred to Maraun and Widmann (2018) for detailed understanding.

5) ASSESSING REGIONAL PRECIPITATION AND TEMPERATURE RESPONSE, AND SIGNIFICANCE
To assess how each GCM represents local responses to co-behavior of large-scale drivers, we explore the response of regional precipitation and temperature under different combined large-scale processes and circulation states. We do so by using the obtained scores of the retained PCs from the models to identify conditions of strong large-scale and/or circulation variations driving each PC. A PC is considered a ''strong'' driver of the regional climate when its score exceeds a threshold of plus or minus one standard deviation. This allows for the evaluation of the co-behavior role in conditioning the regional climate response on seasonal time scales as each PC relates to the indices of large-scale processes. Composite precipitation and temperature anomalies for each of the subperiods where the threshold is exceeded are computed for each grid cell of the eight models.
For the purposes of quantifying the magnitude of association between the spatial distribution patterns of the observed co-behavior modes and that of the models, we use Spearman's rank-order spatial correlation. This method is a nonparametric rank statistic that measures the strength of the association between two variables (Hauke and Kossowski 2011). It is much more robust in the presence of outliers, does not require the assumption of normality and is unbiased by nonlinear relationships when compared to other methods (Endris et al. 2016).
Furthermore, we use the multimodel ensemble mean to ascertain the average model response to each cobehavior mode. The multimodel ensemble mean is calculated by averaging the patterns of each co-behavior mode for the respective models to obtain a single spatial pattern. We then assess the robustness of the patterns by measuring the level of agreement among models based on the condition that at least 80% of the models must agree on the sign of the anomaly with respect to the ensemble mean (Mastrandrea et al. 2011;Pfeifer et al. 2015).

a. Intermodel comparison of SOM seasonal frequencies
We objectively classify the circulation patterns over the southern African region with the SOM. Overall, we identify two distinct patterns with some transitional patterns at 700-hPa geopotential height (Fig. S1). The SOM ( Fig. S1; nodes 4, 8, and 12) captures disturbances in the easterly flow that develops as a result of continuous interactions between the ITCZ and the easterly wave, causing the southward migration of the semipermanent high pressure system and thereby enhancing rainfall conditions in the interior during (austral) summer. Nodes 1, 2, 5, and 9 (Fig. S1) broadly represent a passing midlatitude frontal system that brings rains over the southwestern parts of South Africa. These circulation patterns represent typical (austral) winter synoptic circulation patterns over the region.
In evaluating how the models represent the SOM node seasonal frequencies (Fig. 2), we find most of the observed node frequencies are overestimated across both summer and winter states (e.g., Fig. 2a; nodes 11 and 12, Fig. 2b; nodes 3, 4, and 8) by the models while there are some few cases of underestimation (e.g., Fig. 2a; nodes 4 and 8, Fig. 2b; nodes 9, 10, and 11).
In winter (Fig. 2b), the models mapped more frequently to all nodes although frequencies in nodes 11 and 12 do not often occur. The models overestimate the observed in nodes 1 to 4 and 8, particularly in node 4; CanESM2 (17%), CNRM-CM5 (19%), GFDL-ESM2M (18%), IPSL-CM5A (16%), and MRI-CGCM3 (13%), respectively. MIROC-ESM (14%) strongly overestimates the observed (1%) in node 8. In nodes with some typical winter circulation states ( Fig. 2b; nodes 1 and 2) and also in node 7, the models do simulate fairly the frequency of occurrence. We suggest here that the overestimation (underestimation) by models compared to ERA-Interim may be as a result of models simulating associated circulation types more (less) frequently due to their inability to resolve highly complex large-scale processes (Munday and Washington 2018;Eyring et al. 2019). This may also account for their inability to capture the seasonal variability in some nodes [e.g., Fig. 2, node 4; associated wet summer circulation state]. We note that these variations in frequencies fall within the natural variations of the distribution.

b. Evaluation of co-behavior in GCMs output against observations
We use rotated PCA as a means of examining the interactions between synoptic circulation variability and large-scale climate indices. We find the first principal component (hereafter called PC1) strongly represents the seasonal variation (SOM node frequency variance) with MEI, AAO, and TRBI showing weak correlations and accounts for 30% of explained variance. The second principal component (hereafter called PC2) is largely ENSO response with moderate and weak correlations from TRBI and AAO, respectively. PC2 accounts for 19.8% of the explained variance while the third principal component (hereafter called PC3) is found to be AAO dominant with weak expressions from ENSO and TRBI and accounts for 12.6% of the explained variance. The different phases of each PCs are found to have a peculiar impact on regional precipitation and temperature response in the southern Africa region when the threshold of a ''strong'' driver is attained. For example, typically, an El Niño [positive phase; PC2 . 1 standard deviation (std)] is linked to notable reduction in precipitation, and warming in central parts of the region while a negative phase (PC3 , 21 std) AAO intensifies precipitation in the southwestern parts of the region.
We assess the potential of the influence of co-behavior on regional precipitation and temperature using eight possible combinations by mixing the alternative phases of the identified PCs from the rotated PCA. Designation for the various combinations are shown in Table 3.
Here we present and discuss the precipitation and temperature response for four co-behavior modes (CMs) that showed statistical significance in observed CMs for all models for a general overview of how the models represent CMs and also as a means to reducing redundancy and repetition in discussion. Results for the remaining CMs are shown in the accompanying supplemental material (Figs. S2-S13). We also present later the mean spatial pattern for each CM with respect to all models for precipitation and temperature.
In CM1, wet summer conditions are characterized by a co-behavior of La Niña (PC2 , 21 std), negative AAO (PC3 , 21 std) during summer (PC1 . 1 std). In this mode, above 60% of models are able to fairly represent the general precipitation response and with most showing some wet bias (Fig. 3). However, there are differences in the magnitude and location of the peak intensity across models although statistically significant. For example, peak precipitation intensity is located in the north for CanESM2 (Fig. 3b) whereas the peak is much southward in CNRM-CM5 (Fig. 3c). HadGEM2-ES (Fig. 3e) best represents the response when compared with the observed (CHIRPS; Fig. 3a).
Similarly, during summer (PC1 . 1 std), El Niño (PC2 . 1 std), and positive AAO (PC3 . 1 std) cobehaving characterizes dry summer conditions. We find a mixed representation of precipitation response across the region (CM4; Fig. 4). However, many models fail to produce statistically strong response signals under these forcing conditions, which is concerning as this is a very important drought pattern for the region. For instance, while MIROC-ESM (Fig. 4g) captures peak intensity south of the region, in GFDL-ESM2M (Fig. 4d) peak intensity is farther southwest. Also, CanESM2 (Fig. 4b) and MIROC-ESM (Fig. 4g) show a wet bias relative to the observed (Fig. 4a). Focusing on warm conditions (CM4; Fig. 5) characterized by El Niño (PC2 . 1 std) co-behaving with positive AAO (PC3 . 1 std) during summer (PC1 . 1 std), over 75% of models are able to represent the significant warming found in CRU (Fig. 5a) although a few overestimate its intensity (e.g., CanESM2; Fig. 5b). On the other hand, the models do well to capture the cooling around the southwestern part of the region with the exception being HadGEM2-ES (Fig. 5e).
Similarly, when La Niña (PC2 , 21 std) co-behaves with a positive AAO (PC3 . 1 std) under winter conditions (PC1 , 21 std), about 70% of the models are able to capture the cold conditions across most parts of the region, especially over South Africa and along the Atlantic coast (CM3; Fig. 6). A few models, GFDL-ESM2M (Fig. 6d), IPSL-CM5A-LR (Fig. 6f), and MRI-CGCM3 (Fig. 6i), predict warming in the north/central parts of the region. Figure 7 shows spatial correlations between model composite response spatial patterns and observed composite response spatial patterns across CMs for precipitation. Here, we find no consistent patterns among models in representing CMs as we generally find weak correlations although significant. A possible reason for this to occur is when the sample size for a model for a particular CM is small (Faber and Fonseca 2014;Leppink et al. 2016). In contrast, the models show coherence in representing CMs for temperature (Fig. 8), as they show strong and significant correlations with observed CMs.
Furthermore, we examined the model agreement and ensemble-mean response for each combination with respect to all models for precipitation ( Fig. 9) and temperature (Fig. 10), showing areas where at least 80% of models agree on the sign of response with hatching. For precipitation (Fig. 9), the general circulation states of the multimodel ensemble mean of the GCMs are fairly consistent with the observed CMs (see Figs. S2a-S13a); however, there is not a strong agreement on the sign of the response under certain CMs and there are also differences in intensity and location. For instance, although the ensemble mean (Fig. 9a) under the co-behavior of La Niña (PC2 , 21 std) and negative AAO (PC3 , 21 std) during summers (PC1 . 1 std) captures intense wet conditions over central (southern Zambia, western Zimbabwe, and northern Botswana) and western (Angola and Botswana coast) parts of the region and show strong model agreement on the sign of response, the magnitude is much lower and the location of the peak is not much spread about as noted in the observed data (see Fig. 3a).
Similarly, the intensity of peak dry conditions is much lower and located southward, specifically central to southern parts of South Africa (Fig. 9d) when El Niño (PC2 . 1 std) co-behaves with positive AAO (PC3 . 1 std) during summers (PC1 . 1 std). In the observed, intensity is higher and located to the northeast of the region; around Mozambique spreading toward Zimbabwe, Zambia, and Angola (see Fig. 4a). About 63% (five out of eight) of models represent the peak intensity much southward than the observed (not shown). We suggest that this southward shift of the peak intensity may be as a result of the models exaggerating the extent of the south Indian Ocean convergence zone (SIOCZ), which brings precipitation to the region (Grimm and Reason 2015;Tang et al. 2019).
Focusing on temperature (Fig. 10), there is also generally a good model agreement of sign of response under the CMs with the exception of CM7 when El Niño (PC2 . 1 std) co-behaves with negative AAO (PC3 , 21 std) during summer (PC1 . 1 std) where models do not agree to the sign of anomaly over a small area in the northern part of South Africa and southern Botswana (Fig. 10g). Additionally, when AAO shifts to positive (Fig. 10d), the peak of the warming is not centered but much more spread around the region when compared to the observed (Fig. S12a). The models (Fig. 10b) do not capture the longitudinal shift as noted in the observed (Fig. S9a).

Summary and conclusions
We assesses how CMIP5 GCMs are able to represent already identified co-behavior modes (CMs) in observation datasets and the nature of the variability. Precipitation and temperature datasets for eight GCMs under eight CMs are analyzed in this study.
As a means to establish the relationship between surface expressions and large-scale processes we first identified the synoptic circulation states over the region using SOM. The seasonal frequency of occurrence for each SOM node with respect to each individual model is examined. An evaluation of ERA-Interim and the GCMs show the models are fairly consistent with the reanalysis in representing the seasonal variability although in some nodes there are overestimations and in others there are underestimations.
For the purpose of identifying co-behavior representation in models, we analyze how individual models and their ensemble mean represent co-behavior modes as captured in precipitation and temperature observations. Results show the individual models reasonably simulate precipitation over southern Africa and is fairly consistent with observations. However, the models tend to overestimate precipitation maxima and its location, especially in northern to central parts of the region. For instance, when El Niño (PC2 . 1 std) co-behaves with positive AAO (PC3 . 1 std) during summer (PC1 . 1 std), some models show a wet bias north of the region while precipitation maxima in other models are farther south when compared with the observed in co-behavior mode 4. This mode depicts an important drought pattern for the region and most models fail to produce statistically strong response signals.
We demonstrate that individual models fairly simulate the temperature response and are consistent with the observed over most parts of southern Africa for all co-behavior modes. Most models are able to represent the maxima locations although the intensity is overestimated. The models do well to represent the temperature response when La Niña (PC2 , 21 std) co-behaves with a positive AAO (PC3 . 1 std) under winter conditions (PC1 , 21 std) where about 70% of the models accurately simulate the response. However, in some instances some models FIG. 7. Correlations of models (listed on the left side) vs observed CMs for composite precipitation based on Spearman rank-order spatial correlation expressed as a heatmap. Hatched boxes denote CMs where correlations are not statistically significant at 95% level. show a warm bias over some areas (north/central parts) of cooling in the observed but does well in representing the maxima and its location. Again, over 75% of models are able to adequately simulate the temperature response when El Niño (PC2 . 1 std) co-behaves with positive AAO (PC3 . 1 std) during summer (PC1 . 1 std). Furthermore, we show that there are generally significant weak correlations for precipitation for most model co-behavior modes when compared with observation. However, some individual models show significant positive correlations on specific co-behavior modes. Similarly, for temperature, a greater number of models show significant strong correlations on most co-behavior modes. Thus, we generally find coherent temperature responses in models as opposed to a not so coherent response in precipitation across the models to large-scale co-behaving processes. On the other hand, the multimodel ensemble mean shows a good representation of observed co-behavior modes with a larger spread in precipitation maxima when compared with the individual GCMs and the observation. However, differences in characteristics such as the location and intensity of the maxima to observation as a result of biases in individual models persists (Eyring et al. 2019). This may account for the low number of grid cells where models agree to the sign of the anomaly of the ensemble mean. For temperature response, the multimodel ensemble mean is consistent with observation. There is a strong model agreement on most cobehavior modes with the exception of a co-behavior of El Niño (PC2 . 1 std) and negative AAO (PC3 , 21 std) during summer (PC1 . 1 std) where most models do not agree.
We assert that evaluating the co-behavior of processes in GCMs can provide a rich source of information for how models represent large-scale processes. This may also be useful for climate model evaluation as it is important that models adequately represent the interactions among climate processes over southern Africa or any region without singular circulation control for that matter. As models become more successful in simulating these climate process interactions, climate scientists will have greater confidence in making future projections.