Decadal-to-Multidecadal Variability of Seasonal Land Precipitation in Northern Hemisphere in Observation and CMIP6 Historical Simulations

Based on the centennial-scale observations and CMIP6 historical simulations, this paper employs the ensemble empirical mode decomposition to extract the decadal-to-multidecadal variability of land precipitation (DMVLP) in the northern hemisphere. The spatial distributions of the dominant mode from the empirical orthogonal function are different in four seasons. Regions with the same sign of precipitation anomalies are likely to be teleconnected through oceanic forcing. The temporal evolutions of the leading modes are similar in winter and spring, with an amplitude increasing after the late 1970s, probably related to the overlap of oceanic multidecadal signals. In winter and spring, the Interdecadal Pacific Oscillation (IPO) and the Atlantic Multidecadal Oscillation (AMO) play a joint role. They were in phase before late 1970s and out of phase after then, weakening/strengthening the impacts of the North Pacific and North Atlantic on the DMVLP before/after late 1970s. In summer and autumn, AMO alone plays a part and the amplitude of time series does not vary as in winter and spring. The ability of the coupled models from CMIP6 historical simulations is also evaluated. The good-models average largely captures the spatial structure in four seasons and the associated oceanic signals. The poor-models average is hardly or weakly correlated with observation.


Introduction
In terms of the variability of land precipitation in a long time scale, many studies focused on the trend variability [1][2][3][4][5][6][7][8][9][10]. Yet it was noted that trend analysis had conceptual problems. For instance, the sign of trends and the significance were sensitive to the choice of period. Besides, trend analysis usually assumed a linear trend, unable to depict nonlinear and cyclic variability [11]. In fact, the contribution of the decadal-to-multidecadal variability (DMV) exceeds the trend variability by a large margin [12]. The DMV of land precipitation (DMVLP) is of great importance in understanding the long-term variability and forecasting at decadal time scale, which has attracted considerable attention in recent decades. In previous studies, the DMVLP has been examined mainly on a regional basis, such as North America [13,14], Africa [15], East Asia [16,17], South Asia [18], and so on. Yet, the DMVLP on a hemispheric scale and its mechanism have not been clearly identified. In particular, there were few discussions on the DMVLP in different seasons, while more on the summer monsoon rainfall [19][20][21][22]. Thus, the characteristics of DMVLP on a hemispheric scale in different seasons are an aim of the current work.
Accordant conclusions have not been reached on the mechanisms of DMVLP. Yet, it has been agreed that the oceanic multidecadal signals should play a critical role. In the Pacific, the Interdecadal Pacific Oscillation (IPO) was positively correlated with precipitation over the western and central U.S., eastern Australia and southern Africa [15,23,24]. The positive phase of the Pacific decadal oscillation (PDO)/IPO was related to more precipitation over southern North America and the middle and lower reaches of the Yangtze River, and less precipitation over northern North America, North China, and South China [13,25,26]. The interdecadal change of East Asian summer rainfall in late 1990s was found to be attributed to the interdecadal sea surface temperature (SST) anomalies of North Pacific and West Pacific [17]. In the Atlantic, the Atlantic multidecadal oscillation (AMO) played a part in the northern hemisphere summer monsoon precipitation [20]. AMO drove the variations of North American summer precipitation at decadal time scales, with warm phase corresponding to less rainfall over central U.S. [27][28][29]. AMO also modulated the multidecadal variability of Indian summer monsoon rainfall [18,30], and accounted for the multidecadal precipitation in Siberian warm season to a great extent [31].
Some studies pointed out a joint impact of oceanic multidecadal signals on the DMVLP. Ault et al. found that the relative magnitude of DMVLP over North America was not a linear response to a single decadal mechanism such as PDO/IPO or AMO [14]. Both PDO and AMO had influence on the spatial patterns of precipitation variations in the twentieth century [32]. The Indian summer rainfall was positively correlated with AMO before mid-1990s, but was after then related to Indian Ocean-western Pacific SST anomalies [33]. The role of oceanic multidecadal signals in the DMVLP on a hemispheric scale in different seasons is another aim of the current work.
The Coupled Model Intercomparison Project (CMIP) simulations have been widely used in investigating the long-term variability of precipitation. The trends of global and hemispheric precipitation were evaluated in phase 5 of CMIP (CMIP5) climate models [34][35][36]. The contribution of DMV of precipitation to the total variance was smaller in CMIP5 historical runs than in observation [37]. The ability of CMIP5 models was assessed in predicting the Sahelian monsoon precipitation at decadal time scale [38]. Compared to the observation, all CMIP5 models underestimated the summer precipitation over western Eurasia and overestimated the summer precipitation over eastern Eurasia, which was possibly caused by little cloud cover and strong local coupling of evaporation and precipitation in models [39]. Lyu et al. compared the observed AMO-related precipitation anomalies to that in CMIP5 models [40]. Joshi et al. assessed the impact of IPO on Indian summer monsoon rainfall in CMIP5 models [41]. Nevertheless, the ability of phase 6 of CMIP (CMIP6) simulations in capturing the observed DMVLP, especially on a hemispheric scale, has not been evaluated yet [42].
To sum up, the seasonal DMVLP on a hemispheric scale and the associated oceanic modes are not clear yet. The performance of state-of-the-art climate models in CMIP6 simulations in capturing the observed DMVLP is unraveled. Therefore, the rest of the paper is organized as follows: Section 2 describes data and methods used in the study. Section 3 examines the observed spatial and temporal characteristics of the dominant mode of DMVLP in northern hemisphere (NH) in four seasons. The oceanic signals related to the DMVLP are investigated in Section 4. The ability of the coupled models in CMIP6 historical runs to simulate the observed DMVLP is evaluated in Section 5. Section 6 provides summary and discussion.

Data
This study used the monthly Hadley Centre Sea Ice and Sea Surface Temperature dataset (HadISST) from 1870 to 2016. The spatial resolution was 1° × 1° latitude by longitude [43].
The Global Precipitation Climatology Centre (GPCC) monthly land-surface precipitation Full Data Product (V2018) from 1891 to 2016, based on quality-controlled data from global stations, was employed. The spatial resolution was 1° × 1° latitude by longitude [44]. There were some grids with missing data before 1940; however, obvious differences were not found between 1891-2016 and 1941-2016 in the spatial and temporal features from empirical orthogonal function (EOF) on the seasonal precipitation. Thus, the total length of time was used in this study.
The monthly SST and precipitation datasets of the coupled models are from CMIP6 historical simulations, with temporal coverage from 1870 to 2014. The models chosen in the current work had at least 3 members in a historical ensemble, and then differences among members in a model ensemble could be examined. The models and the number of ensemble members are listed in Table 1. Table 1. The coupled models from CMIP6 historical simulations.

Models
No. of Ensemble Members An identical period of time from 1891 to 2014 was used for all datasets. The seasonal means were calculated and then climatology was removed. As the current work focuses on the DMVLP in NH, winter/spring/summer/autumn was used instead of boreal winter/boreal spring/boreal summer/boreal autumn for brevity.

Methods
The ensemble empirical mode decomposition (EEMD) was employed to extract the DMV of precipitation and SST. EEMD decomposes a time series into some empirically determined intrinsic mode functions from high to low frequency and separates time scales without subjective criterion selection, which could overcome the confusion of modes with different time scale [45]. Then the EOF analysis was applied to obtain the primary mode of DMVLP.
The auto-correlation of DMVLP was much larger than the seasonal precipitation. Thus, the effective degree of freedom (EDOF) was reduced substantially for DMVLP, and was calculated after Bretherton et al. [46]. The t-test was used to determine the statistical significance.

Spatial and Temporal Features of the DMVLP in GPCC
The observed spatial and temporal characteristics of the DMVLP in four seasons were investigated in this section.

Winter
The dominant mode from EOF analysis of the wintertime DMVLP in NH explained 13% of the total variance. Figure 1a shows its spatial distribution. Significant positive anomalies were found over the central North America, northern South America, and regions to the north of Mediterranean Sea; and significant negative anomalies were over southern U.S., central and western Russia, Arabian Peninsula, and southern China.
The time series displays obvious DMV ( Figure 2a). It was in the positive phase during 1920-1955 and 1995-2014, and in the negative phase during 1895-1920 and 1975-1995. Based on the spectrum analysis, the 30-60-years period was the primary and the 10-15-years period had much less power (Figure 3a). It is noteworthy that the amplitude of variability appeared to increase after the late 1970s, in agreement with the time of the 1976-1977 transition [47,48]. It is probably indicative of a larger rainfall intensity in recent decades.

Spring
The leading EOF mode of the springtime DMVLP in NH explains 11% of the total variance. There were significant positive anomalies over the northern South America, western Russia, and Indo-China Peninsula, and significant negative anomalies over the southwestern North America, Arabian Peninsula, Iran Plateau, and southeastern China (Figure 1b). Similar with winter, the time series shows DMV, with 40 years as the primary period ( Figures  2b and 3b). The larger amplitude after the late 1970s was also found.

Summer
The dominant EOF mode of the summertime DMVLP in NH explained 9% of the total variance. Significant positive anomalies were found over the southern North Africa, western Europe, India, eastern Russia, and northern China; and significant negative anomalies appeared in the central and southern North America and western Russia (Figure 1c).
Unlike winter and spring, the amplitude of variability in summer did not increase after the late 1970s (Figure 2c). Zhu et al. found increased amplitude of the variability of the contiguous United States summer rainfall in the early 1990s [22], which was not seen for the summertime land precipitation in NH. The positive and negative phases of the time series in summer did not agree with winter and spring. There were three significant periods, 60 years, 25-30 years, and 15-20 years, in which the 60-years period was the primary (Figure 3c).

Autumn
The primary EOF of the autumn time DMVLP in NH explained 9% of the total variance. Significant positive anomalies were found over Europe, southern North Africa, and India, and significant negative anomalies were over the western North America, southeastern Greenland, and eastern China (Figure 1d).
The time series in autumn had certain resemblance as summer, without evident change of the amplitude of variability (Figure 2d). The main period was 25 years and the secondary period was 35-40 years (Figure 3d). The 70-years period was not significantly different from red noise.

Summary
There was a seasonal distinction for the leading mode of DMVLP in NH. The temporal evolutions in winter and spring were in phase, with the amplitude of variability increasing after the late 1970s, suggesting a larger intensity of rainfall in recent decades. The temporal evolutions in summer and autumn were similar, without apparent variation of amplitude. In any season, regions with the same sign of anomalies were probably teleconnected. The oceanic signal was likely to be a dominating factor.

Oceanic Signals
In order to identify the critical oceanic regions that influence the DMVLP in NH, regressions of SST anomaly (SSTA) onto the first principle component of DMVLP are plotted in Figure 4.
The SSTA patterns were similar in winter and spring (Figure 4a,b). In the Pacific, it was IPO in its cold phase, with a large negative anomaly in the eastern Pacific triangle and a positive anomaly in the western Pacific K-shape regions. These regions were used to define a mega-ENSO index [20]. In the Atlantic, it was AMO in its warm phase, with a large and significant positive anomaly over the high latitude North Atlantic and tropical Atlantic.
The SSTA in the Pacific changes in summer and autumn (Figure 4c,d). In summer, the positive anomaly in the northwestern Pacific moves northward, and the negative anomaly in eastern Pacific triangle became much weaker and insignificant. The signal in the Pacific disappeared in autumn. Yet, the AMO was still in its warm phase in summer and autumn. The positive SSTA over mid-latitude North Atlantic was larger in summer than in the other seasons. The signal from the Indian Ocean was statistically insignificant in all seasons. In a word, the AMO played a role in the DMVLP in NH in all seasons. In winter and spring, there was joint influence of IPO and AMO. The configuration of the oceanic multidecadal signals was key to understanding the larger amplitude of the time series after the late 1970s in winter and spring.
The positive and negative phases of IPO and AMO overlapped. Previous work showed that different combinations of IPO and AMO phases were related to the multidecadal East Asian summer monsoon precipitation and the multidecadal drought frequency over the U.S. [21,49,50]. IPO and AMO were in phase before the late 1970s, inducing opposite signs of SSTA over the North Pacific and North Atlantic. Thus, the impact of SSTA from these two oceans might counteract in part. However, they were out of phase during the period from the late 1970s to the late 1990s, with positive IPO and negative AMO, corresponding to negative SSTA over majority of NH. They were also out of phase during the period from the late 1990s to 2010s, with negative IPO and positive AMO, corresponding to positive SSTA over majority of NH. The influence of SSTA, with consistent signs over the North Pacific and North Atlantic after late 1970s, could be reinforced, contributing to larger intensity of rainfall.
In summer and autumn, AMO was the only multidecadal oceanic signal that affected the leading DMVLP in NH. The amplitude of the time series did not display a similar variation as in winter and spring.

DMVLP in CMIP6 Simulations
The DMVLP in NH was also investigated in different seasons in state-of-the-art coupled models from CMIP6 historical simulations. The spatial and temporal correlations between observation and each member from all models listed in Table 1 were calculated ( Figure 5). Models differences were evident. The performance of a member was not consistent with other members even in the same model. It shows that the largest spatial correlation was 0.48 between observation and the fourth member of CNRM-ESM2 in spring, and the largest temporal correlation was 0.65 between observation and the seventh member of CanESM5 in winter. A model or an ensemble member from a model, which was able to capture the observed spatial pattern to a certain degree, did not necessarily well simulate the observed temporal evolution; and vice versa. Therefore, both spatial and temporal characteristics should be taken into account when evaluating the ability of a model to simulate the observed DMVLP. In this study, good models were found as those with both spatial and temporal correlations larger than 0.2, indicated by red boxes in Figure 5; poor models were those with spatial and temporal correlations smaller than 0.1, irrespective of sign, indicated by blue boxes. Then the spatial averages of good models and poor models are displayed in Figure 6.
It was evident that the good models average captured the spatial structure of the observed significant precipitation anomalies to a great extent in the four seasons. In winter, the observed positive anomalies over the central North America and northern South America, as well as the negative anomalies over southern U.S., central and western Russia, and southern China, were simulated by the good models average (Figure 6a). In spring, the observed positive anomalies over the northern South America and Indo-China Peninsula, and negative anomalies over the southwestern North America, Iran Plateau, and southeastern China were simulated (Figure 6c). In summer, the observed positive anomalies over the southern North Africa, India, eastern Russia, and northern China, and negative anomalies over the central and southern North America were found (Figure 6e). In autumn, the observed positive anomalies over Europe and India, and negative anomalies over the western North America and eastern and central China were captured (Figure 6g). The observed pattern was better correlated with the good models average than any single model (CC = 0.5 in winter and spring, CC = 0.4 in summer and autumn). The poor models average and the observed pattern were hardly correlated in NH in winter and spring (CC<0.1, Figure 6b,d). They were weakly correlated in summer and autumn (CC = 0.2); however, the magnitudes of the precipitation anomalies in poor models averages were weaker than observations (Figure 6f,h).

Regressions of SSTA
The averages of regressions of SSTA from good models and poor models were also examined ( Figure 7). The good models averages captured the observed AMO pattern in winter, spring, and summer, but not in autumn (Figure 7a,c,e,g). They also captured the observed IPO pattern in winter and spring. Unlike the observation, the signal in the Pacific was still very strong in summer and autumn in good models averages, which could partially account for the smaller correlation of the leading mode of DMVLP between the good models average and observation in summer and autumn than in winter and spring. The poor models averages did not capture the oceanic multidecadal modes (Figure 7b,d,f,h), suggesting that the connection between the IPO/AMO and the DMVLP should be rather weak in the poor models.

Summary
Using EEMD to extract the decadal-to-multidecadal variability of land precipitation (DMVLP) from GPCC, this study examined the spatial and temporal characteristics of the dominant mode from EOF analysis on the seasonal DMVLP in the northern hemisphere. The spatial distributions were different in the four seasons. In winter, precipitation was more than the climatology over the central North America, northern South America, and regions to the north of Mediterranean Sea, and was less over southern U.S., central and western Russia, Arabian Peninsula, and southern China. In spring, there were significant positive anomalies over the northern South America, western Russia, and Indo-China Peninsula, and significant negative anomalies over the southwestern North America, Arabian Peninsula, Iran Plateau, and southeastern China. In summer, significant positive anomalies were found over the southern North Africa, western Europe, India, eastern Russia, and northern China; and significant negative anomalies appeared in the central and southern North America and western Russia. In autumn, precipitation was larger than climatology over Europe, southern North Africa, and India, and smaller over the western North America, southeastern Greenland, and eastern China. The temporal evolutions were similar in winter and spring, with amplitude increasing after the late 1970s, which was probably related to the phase shift of IPO and AMO. In summer and autumn, IPO changed and disappeared, leaving AMO as the only signal that affected the leading DMVLP. The amplitude did not display similar variations as in winter and spring.
There were differences among the coupled models from CMIP6 historical simulations in capturing the observed DMVLP in light of both spatial and temporal correlations. The good models average largely captured the spatial structure of the leading EOF mode in the four seasons, and its correlation with the observation was larger than any single model. The associated oceanic signals, IPO and AMO, were simulated in winter and spring by the good models average. However, the IPO persisted in summer and autumn in good models averages, probably responsible for a weaker correlation between the good models average and the observation in summer and autumn. The poor models average was hardly correlated with observation in winter and spring, and weakly correlated in summer and autumn. The IPO and AMO were not connected with the leading DMVLP in poor models.

Discussion
Recent studies showed that the atmosphere-ocean interactions related to multidecadal variability in an ocean basin could favor the variability in another ocean basin (called "trans-basin interactions"), in particular the Atlantic Ocean played an active role in low frequency variability in the Pacific and in the recent global warming hiatus [51,52]. AMO also accounted for the enhanced global monsoon circulation and precipitation via trans-basin interaction [53]. Smith et al. pointed out that historical anthropogenic aerosols might be the key for the negative PDO phase during global warming hiatus and resultant global-scale atmospheric circulation [54]. In the current work, besides the DMVLP that is considered as intrinsically generated variability, we also examined the externally forced variability of land precipitation in observation and CMIP6 models. In observation, the dominant mode of the externally forced land precipitation in winter shows a significant increasing trend and significant positive anomalies over the majority of Eurasian continent (Figure 8a,d). The associated oceanic signal was a global warming pattern (figure not shown).
In response to external all-forcings, the wintertime land precipitation in a CMIP6 model displays a similar increasing trend. The spatial pattern resembled the observation in many regions, especially in mid and high latitude Eurasian continent (Figure 8b,d). In response to anthropogenic aerosol forcing, the wintertime land precipitation from the same model also shows a significant increasing trend, but the spatial pattern was opposite in sign, indicating that the anthropogenic aerosol forcing mainly induces the precipitation to decrease over most of NH (Figure 8c,d).