Evolution of the East Asian winter land temperature trends during 1961–2018: role of internal variability and external forcing

Detecting the contributions of internal variability and external forcing to the evolution of surface air temperature (SAT) trend at regional scales is a challenge. Based on the observations and large-ensemble simulations of climate models, we estimate the contribution of the internal and forced components to the evolution of East Asian winter land SAT (EAWT) during 1961–2018. Although the external forcing induced EAWT trends show a slow increase, both the total and internally generated EAWT trends exhibit a decrease with the extension of the time period, suggesting a critical role of internal variability in the evolution of the EAWT trends. The internal variability contributes to about 70% of total EAWT trends during 1961–1995. With the extension of the time period, the contribution of internal variability decreases, whereas the contribution of external forcing gradually grows to dominate the EAWT trends. Based on the dynamical adjustment method, we identify that the internal dynamics and forced thermodynamics account for a majority of internal and forced EAWT variations, respectively. We further identify that the multidecadal fluctuation of internal component of autumn Arctic sea ice is a critical precursor of the internal variability, especially the internal dynamically induced EAWT variations, through triggering a meridional stationary Rossby wave response in the following boreal winter. Our findings provide an insight into the understanding of the present and future climate change over East Asia.


Introduction
Despite the increasing emission of greenhouse gases (GHGs) in the past decades, the surface air temperature (SAT) trends exhibit considerable spatial and temporal differences (e.g. IPCC 2013, Gong et al 2019a). Understanding the causes of the observed SAT trends, particularly at regional scales, is a major challenge because the climate system varies on multidecadal time scales in response to a variety of external forcing as well as due to its own internal variability (e.g. Dong and Mcphaden 2017, Perkins-Kirkpatrick et al 2017. East Asia is one of the most populated regions of the world, and climate change exerts substantial impacts on the society and economy in this region (e.g. Wang et al 2017). In sharp contrast to the global warming due to increasing GHGs, many studies have demonstrated that East Asia has experienced frequent damaging snowstorms and severe cold spells in winter which is associated with the recovered intensity of the East Asian winter monsoon (EAWM) during recent two decades (e.g. Zhou et al 2009, Jeong et al 2011, Wu et al 2011, Wang and Chen 2014, Zhou et al 2017, Gong et al 2019b. The abnormal phenomenon implies that the internal variability may modulate the SAT trend over East Asia on the multidecadal time scale. Therefore, quantifying the contributions of internal variability and external forcing to the recent East Asian winter SAT trends and their temporal evolutions during past decades is particularly essential for us to understand the wintertime climate change over East Asia in the past and future. A recent study has indicated that the observed East Asian winter SAT trends indeed can be largely influenced by internal variability (Gong et al 2019b).
During the past 40 years , the internal variability offsets more than 70% of forced warming over northern East Asia (Gong et al 2019b). However, some questions remain unanswered. How have the East Asian winter SAT trends evolved during the past decades? What is the role of internal variability and external forcing in the evolution of East Asian winter SAT trends? Moreover, if the internal variability plays a certain role in the evolution of the East Asian winter SAT, what are the key internal physical processes and precursors? Previous studies have indicated that the Arctic sea ice loss may play a role in the widespread winter cooling over the Eurasian continent during recent two decades (e.g. Honda et al 2009, Wu et al 2011, Cohen et al 2012, 2014, Liu et al 2012, Gao et al 2015. How does the importance of the preceding Arctic sea ice change in the evolution of the East Asian winter SAT trends during past decades?
This work aims to address the above questions in detail. The 'dynamical adjustment' approach is employed to estimate the time-varying contributions of internal and forced dynamic and thermodynamic components to the East Asian winter SAT changes. This approach can linearly separate the effects of dynamical and thermodynamical variability on surface climate change (e.g. Smoliak et al 2015, Deser et al 2016, Lehner et al 2018.

Data
The observational monthly-mean land SAT data are obtained from the Climatic Research Unit (CRU) of the University of East Anglia (CRU_TS version 4.03) with a horizontal resolution of 0.5 • × 0.5 • covering the period from 1901 to 2018 (Harris et al 2014). The CRU data are bilinearly interpolated to a horizontal resolution of 1 • × 1 • prior to the analysis. The observed monthly sea ice concentration data are obtained from the Hadley Centre Sea Ice and Sea Surface Temperature since 1870 to present with a horizontal resolution of 1 • × 1 • (Rayner et al 2003). The observational proxies of monthly sea level pressure (SLP) are derived from the Japanese 55 year (JRA-55) reanalysis, which spans the period from January 1958 to the present with a horizontal resolution of 1.25 • × 1.25 • (Kobayashi et al 2015). Hereinafter, the results based on the above data are referred to as 'observations' .

Model simulation
To estimate the forced response of observations to the external forcing, the large-ensemble simulation of the state-of-the-art climate models is used. Previous studies have reported that the MPI-ESM and CESM1 models are two of the few climate models that simulate the Northern Hemisphere climate change well (e.g. Guo et al 2019, Huang et al 2020). The MPI and CESM models can also reasonably simulate the winter SAT changes over East Asia (figure S1 (available online at stacks.iop.org/ERL/16/024015/mmedia)). Therefore, the outputs of historical and representative concentration pathway (RCP) 8.5 simulations in 100-member ensemble simulations with Max-Planck-Institute Earth System Model (MPI-ESM) and 40-member ensemble simulations with the Community Earth System Model, version 1 (CESM1) are employed (Giorgetta et al 2013, Deser et al 2016 to investigate the relative contributions of external forcing and internal variability to the East Asian winter SAT changes, similar to previous studies (e.g. Deser et al 2016, Guo et al 2019. These experiments follow the design from the fifth phase of Coupled Model Intercomparison Project (CMIP5). Each ensemble member in MPI-ESM and CESM1 is subjected to the same radiative forcing (historical for 1920-2005 and RCP 8.5 emissions scenario for 2006-2100), but starts from slightly different initial conditions. Therefore, the large-ensemble mean SATs in MPI-ESM and CESM1 contains a primarily externally forced response (e.g. Dai et al 2015). Here the external forcing response is a general term including natural and anthropogenic sources. Prior to the analysis, all the models' data are bilinearly interpolated to a common resolution of 1 • × 1 • . It should be noted that the results are similar if a lower resolution is used in the analysis.

Analysis methods
Given the quality of circulation and temperature data (Wu et al 2005), the present analysis focuses on the period from 1961 to 2018. Here, our analysis focuses on boreal winter, and the winter means are constructed by averaging monthly data of December, January, and February. Our convention is that the winter of 1961 refers to the 1960/61 winter. The linear trends are calculated using the least-square regression analysis. To eliminate the noise, all the results are smoothed using 3 year running averages before computing the trends. The significance of the trends is evaluated with a two-tailed Student's t-test.
Note that the forced responses in both global and East Asian mean winter SAT are basically consistent between CESM1 and MPI-ESM ( figure S2). Therefore, the averaged SAT in the 40-ensemble mean of CESM1 and 100-ensemble mean of MPI-ESM is used to provide the best estimate of the forced responses of the observed SAT trends. Similar to Dai et al (2015), the component of regional-scale internal variability in East Asian winter SAT is obtained by removing the variations and changes associated with the area-averaged forced SAT over East Asia (20 • -60 • N, 80 • -150 • E) at each grid point through linear regression.
The dynamical adjustment method based on partial least squares (PLS) regression (Smoliak et al 2015) is employed to further separate the contributions of forced and internal dynamical and thermodynamical components to SAT trends. The thermodynamically induced refers to SAT changes induced by timevarying radiative fluxes or by time-varying fluxes of sensible and latent heat at the Earth's surface, excluding any concomitant changes in the atmospheric circulation and dynamically-induced denotes SAT changes attributable to changes in the atmospheric circulation. Here, we use the PC-wise dynamical adjustment method to separate the dynamical and thermodynamical components of the East Asian winter SAT field (e.g. Smoliak et al 2015, Gong et al 2019b. First, we decompose the SAT in the domain of East Asia (20 • -60 • N, 80 • -150 • E) into EOFs and perform PLS regression upon the corresponding PCs. Then, the dynamically adjusted SAT field is formed by summing the corresponding products of the SAT PCs with their respective EOFs. The leading three SLP patterns in the region of 10 • -90 • E, 40 • -210 • E are employed because the additional variance explained by the fourth predictor is small. Moreover, the results are similar if the whole Northern Hemisphere SLP is used. This PC-wise dynamical adjustment method is interpreted in Smoliak et al (2015) in details. After the dynamical adjustment, the residual component is regarded as the thermodynamical contribution. To obtain the dynamical contribution, we subtract the residual contribution from total variability. Similarly, the averaged dynamical contributions in all the 140 members of CESM1 and MPI-ESM models are used to derive the best estimate of the forced dynamical contribution of the observed SAT trends. The forced thermodynamic contribution is obtained as a residual of the forced components. Similarly, the observed internal dynamical and thermodynamical components are obtained by removing the variations associated with the area-averaged forced dynamical and thermodynamical components from dynamically and thermodynamically induced SAT variations over East Asia at each grid point through linear regression, respectively. Since the different phase shift of internal variability may induce distinct contribution to the SAT trends, this result implies that the EAWT trends are largely modulated by the internal variability. To further illustrate the role of internal variability in the evolution of EAWT trends, the forced components of SAT variations are linearly removed from the total SAT variability over East Asia (see the details in section 2). It shows that the internally generated variability indeed contributes a large part of the EAWT trend, especially during 1961-1995 (figures 1(d)-(f)). The average warming rate of the internally induced EAWT trend is 0.38 K (10 yr) −1 during 1961-1995 and accounts for 70% of the warming of total EAWT in this period ( figure 1(d)). With the extension of the time period, the internal induced SAT trend is gradually reduced and even becomes insignificant when the time period is extended to 1961-2015. This result suggests that the most pronounced warming during 1961-1995 is largely due to the large positive contribution of the internal variability to the EAWT trend.

Results
The internally generated EAWT variability can be further separated into dynamically and thermodynamically induced components based on the dynamical adjustment method (see the details in section 2). The spatial distribution of both internal dynamically and thermodynamically induced EAWT trends shows a gradually reduced warming rate of EAWT with the extension of the time period (figure 2), which is consistent with the evolution of total internal generated EAWT trends (figures 1(d)-(f)). The patterns of internal dynamically induced EAWT trends resemble those of internal induced EAWT trend patterns in both the magnitudes and spatial distributions. In contrast, the magnitudes of internal thermodynamically induced EAWT trends are much weaker than those related to internal dynamics. The internal thermodynamically induced warming only appears in the northeastern and northwestern China before 2005 as shown in figures 2(d) and (e). During 1961-2015, the internal dynamically and thermodynamically induced EAWT trends are both weak and insignificant in the majority of East Asia. These results indicate that both the internal dynamically and thermodynamically induced EAWT trends are gradually reduced with the extension of the time period, and the internal dynamical component accounts for the major proportion of internally generated EAWT variability.
To quantitatively illustrate the time-varying contributions of forced and internal variability to the EAWT trends, the temporal evolutions of the total, forced, and internal induced area-averaged EAWT trends are displayed in figure 3(a). The x-axis indicates the ending year of the calculated trend for the time period starting from 1961. With the extension of the time period, the total EAWT trends decrease rapidly, and the internally induced trends show a very similar feature to the total trends ( figure 3(a)). The consistent evolution of total and internal EAWT trends suggests that the evolution of total EAWT trends is mainly modulated by the internal variability. In contrast, the forced EAWT trends exhibit a slow increase from the period 1961-1995 to 1961-2018. The internal variability contributes about 70% of total EAWT trends during 1961-1995, but it decreases Three-year running averages were applied to local SAT time series before computing the trends. Dots indicate the SAT trends exceeding the 95% confidence level. The confidence level is calculated from all grids.
to 16% when the time period is extended to the whole period 1961-2018. In contrast, the contribution of the forced component is about 30% during 1961-1995 but increases to over 80% with the extension of the time period to 1961-2018. When the time period is 45 years , the internal variability has a contribution equivalent to the forcing induced EAWT trends ( figure 3(c)). These results indicate that the contributions of the internal variability and external forcing to the EAWT trends depend on the time period, which is somewhat consistent with some previous studies (Hawkins and Sutton 2009. With the extension of the time period, the internal variability shows a decreased contribution to the EAWT trends, and the contribution of external forcing gradually grows to dominate the EAWT trends. Figure 3(b) further displays the temporal evolutions of forced and internal dynamically and thermodynamically induced EAWT trends. The evolutions of internal dynamically and forced thermodynamically induced EAWT trends to resemble the evolutions of internaland forced-induced components (figures 3(a) and (b)). The internal dynamical and forced thermodynamical components account for more than 70% and 60% of internal-and forced-induced EAWT trends, respectively, regardless of the length of the time period during 1961-2018 (figures 3(a) and (b)). The contributions of forced dynamics and internal thermodynamics both play a secondary role in the EAWT trends regardless of the length of the time period ( figure 3(d)).
The above results suggest that the reduced total EAWT trends are mainly attributed to the reduced contribution of internal variability, especially the part due to internal dynamics, with the extension of the time period because the internal variability varies on the multidecadal time scale. In view of the importance of internal dynamics to the EAWT trends, especially when the time period is shorter than 45 years, it should be clarified what is the key internal dynamical process responsible for the time-varying contribution to the EAWT trends? Figures 4(a) and (b) show the correlation patterns between the 9 year low-pass components of internal dynamically induced EAWT (Int_Dyn EAWT) and the corresponding detrended SLP and 500 hPa geopotential height (Z500) fields over the Northern Hemisphere, respectively. In the lower-troposphere, widespread significant negative correlations are found over the northern Eurasian continent extending from the Ural-Mountains region eastward to the central Siberia ( figure 4(a)). In the middle-troposphere, there are negative correlations over the Ural-Mountains region and positive correlations over the eastern Siberia eastward to the East Asian coast, revealing a meridional stationary Rossby wave structure ( figure 4(b)). This circulation pattern is similar to the atmospheric response to the autumn Arctic sea ice variations mentioned in previous studies (  autumn Arctic sea ice anomalies can trigger a stationary Rossby wave response, with an anomalous anticyclone/cyclone over the Ural-Mountains to the surface turbulent heat flux anomalies associated to the changes of autumn Arctic sea ice. And the subsequently dynamical feedback (e.g. eddy-mean flow interaction) processes further facilitate the formation of an anomalous meridional circulation in the following winter that influences the East Asian winter temperature anomalies (e.g. Honda et al 2009, Wu et al 2011, Nakamura et al 2015, Luo et al 2018, 2019. Therefore, the internal dynamically induced EAWT variability may be rooted in the internal autumn Arctic sea ice variations. In view of the long-term decrease of the Arctic sea ice during past decades due to the anthropogenic forcing, the detrended Arctic sea ice is used to represent the internal Arctic sea ice variations. The correlation pattern between the Int_Dyn EAWT index and detrended autumn Arctic sea ice on their lowfrequency components obtained through the 9 year low-pass filter is shown in figure 4(c). Widespread significant positive correlations are found from the Eurasian marginal seas eastward to the Beaufort Sea ( figure 4(c)). Here, we define an Arctic Sea Ice index using the area-averaged Arctic sea ice over the key region (SIC, 72 • -82 • N, 60 • E-130 • W) to further investigate the influence of autumn SIC on the internal dynamical process in the following winter as well as associated EAWT variability. The correlation patterns between the 9 year low-pass components of internal autumn SIC and the corresponding detrended SLP and 500 hPa geopotential height fields are quite consistent with those associated with the Int_Dyn_EAWT (figure S3). The correlation coefficients between the internal SIC and Int_Dyn EAWT indices on their low-frequency components are 0.83 ( figure 4(d)). The preceding autumn SIC accounts for 69% of the total variance of Int_Dyn EAWT variations on the decadal time scales. Moreover, the temporal evolutions of the key regional atmospheric circulation, that is, the Siberian High (SH) in the lower troposphere and Ural blocking (UB) in the middle troposphere, are shown in figure 4(d). The SH and UB indices are defined as the area-mean detrended SLP averaged over the region of 45 • -70 • N, 60 • -120 • E and detrended Z500 over the region of 55 • -85 • N, 30 • -90 • E, respectively ( figure 4(d)). It should be noted that the signs of SH and UB indices shown in figure 4(d) are reversed for convenience of comparison with the SIC and Int_Dyn_EAWT indices. It is clear that the temporal evolutions in the Int_Dyn EAWT, SH and UB indices are highly consistent ( figure 4(d)). The correlation coefficients of the lowfrequency components of Int_Dyn EAWT with the corresponding SH and UB indices are −0.97 and −0.96, respectively. The correlation coefficients of the low-frequency components of internal autumn SIC with the corresponding SH and UB indices reach −0.90 and −0.91, respectively. These results suggest that the Int_Dyn_EAWT variations are closely tied to the regional circulation response to the changes of internal autumn sea ice and the stationary Rossby wave response from the Ural-Mountains to the East Asia associated with the internal autumn Arctic sea ice changes is the key dynamical connection between internal autumn Arctic sea ice and following internal East Asian winter temperature anomalies.
It is noted that there are obvious multidecadal fluctuations in the internal autumn SIC and Int_Dyn EAWT, with an increase from 1961 to the early-1990s and then decrease afterward. To further demonstrate the contribution of internal autumn SIC to the internal winter circulation and EAWT trends, we divide the whole period into two sub-periods at the year of 1993 because there are most significant increases of indices before 1993 and then decreases afterward. To reveal the internal SIC induced winter SLP, 250 hPa geopotential height (Z250) and SAT trends, we reconstruct the three-dimensional winter SLP, Z250 and SAT anomalies associated with the internal SIC based on the linear regression. During 1961During -1993, the increase in internal autumn SIC induces widespread negative SLP anomalies extending southward from the Arctic to northern Eurasia and weakens the SH ( figure S4(a)). In the middleupper troposphere, the strong negative geopotential height anomalies are observed from the Arctic region to Ural-Mountains region and positive geopotential height anomalies from the eastern Siberia to the East Asia, revealing a meridional stationary Rossby wave response to the increase of internal autumn SIC. This anomalous anticyclone over the East Asia further weakens the East Asian trough in situ ( figure S5(a)). The resultant weakened EAWM induces the widespread warming of the SAT over East Asia ( figure 5(a)). In contrast, there are widespread positive SLP anomalies extending southward from the Arctic to northern Eurasia, accompanied  (a) 1961-1993 and (c) 1993-2018. (b), (d) same as (a), (c), but for the corresponding internal dynamically induced land SAT trends. Three-year running averages were applied prior to computing the trends. Dots indicate the SAT trends exceeding the 95% confidence level. The confidence level is calculated from all grids. by the strengthening of SH in response to the decrease of internal autumn SIC during 1993SIC during -2018). In the middle-upper troposphere, strong anticyclonic anomalies extend from the Arctic to the Ural Mountains region and cyclonic anomalies from the eastern Siberia to the East Asia, exhibiting a meridional stationary Rossby wave response to the decrease of internal autumn SIC. This anomalous cyclone further deepens the East Asian trough ( figure S5(c)). The resultant strengthened EAWM induces the widespread cooling trends over East Asia except for the Tibet Plateau ( figure 5(c)). The internal autumn SIC induced SLP, Z250, and EAWT trends resemble those associated with the internal dynamics in both the spatial pattern and magnitudes during the two sub-periods (figures S4(b), S4(d), S5(b), S5(d), 5(b) and (d)). It is noted that the cooling trend patterns over the East Asia associated with the reduction of internal autumn SIC and the internal dynamical processes during recent two decades are similar to the recent observed total EAWT cooling trend pattern shown in Cohen et al (2014). This result suggests that the recent widespread cooling over the Eurasian continent may be largely attributed to the multidecadal decrease of the internal autumn SIC. And the internal autumn SIC induced EAWT trends contribute 79% and 90% of total internal dynamically induced EAWT trends during 1961EAWT trends during -1993EAWT trends during and 1993EAWT trends during -2018 These results further confirm that the multidecadal fluctuation of internal autumn SIC variation is a critical precursor significantly influencing the EAWT trends by modulating the winter atmospheric circulation response. In other words, the evolution of EAWT trends depends on the multidecadal changes of internal autumn SIC to a great extent.

Summary and discussions
In this study, we reveal the evolution of EAWT trends and identify the time-varying contributions of internal and forced dynamical and thermodynamical components to the observed EAWT trends during 1961-2018 based on the dynamical adjustment method and large-ensemble simulation from MPI-ESM and CESM1 models. Our results show that although the external forcing induced EAWT trends show a slow increase, the total EAWT trends exhibit an obvious decrease with an extension of the time period during 1961-2018. The consistent evolutions between the total and internal generated EAWT trends suggest a critical role of internal variability in the evolution of the EAWT trends. The internally generated EAWT trends contribute about 70% of total EAWT trends during 1961-1995, but the contribution gradually decreases to 16% when the time period is extended to the whole period 1961-2018. The internal variability shows an equivalent contribution to the forcing induced EAWT trends during the time period of 45 years . Based on the dynamical adjustment, we identify that the internal dynamical and forced thermodynamical components account for the majority of the internal and forced EAWT variations, respectively. The evolution of total EAWT trends is mainly attributed to the time-varying contribution of internal dynamics to the EAWT trends. Our findings provide an insight into the near-term climate projection that the internal variability is actually a critical factor influencing the SAT trends on the decadal time scale.
Our results further indicate that the internal dynamical processes, a major component of internal variability contributing to the evolution of EAWT, are closely associated with the internal autumn Arctic sea ice variation over the Eurasian marginal seas eastward to the Beaufort Sea, through triggering a meridional stationary Rossby wave response in the following boreal winter that influences the East Asian winter temperature anomalies. It is noted that the hemisphere-scale atmospheric circulation response to the internal autumn Arctic sea ice variation resembles the Arctic Oscillation/North Atlantic Oscillation (AO/NAO) pattern (e.g. Thompson and Wallace 1998, Deser 2000, Ambaum et al 2001, Magnusdottir et al 2004, Wang and Chen 2010, Gong et al 2017, He et al 2017, but the Arctic center shifts southward to the northern Eurasian continent (figures S4 and S5). The correlation coefficient between the low-frequency components of internal AO and Int_Dyn EAWT indices obtained using a 9 year low-pass filter is 0.87, weaker than that in the regional circulation mode (e.g. SH and UB). The correlation coefficients of the internal AO with the lowfrequency components of internal SH and UB are −0.91 and −0.87, respectively. Therefore, it implies that the relationship between the internal AO and Int_Dyn EAWT may be attributed to its close connection with the regional circulation patterns over Eurasia.
The internal autumn Arctic sea ice can be regarded as a critical precursor of the following winter internal dynamically induced EAWT variations, but the reasons responsible for the multidecadal evolution of the internal autumn Arctic sea ice are unclear. Previous studies indicated that the preceding and concurrent atmospheric circulation anomalies over the Arctic region can exert a substantial impact on the evolution of autumn Arctic sea ice variations . Therefore, we speculate that the causes of internal induced autumn Arctic sea ice variations may be associated with the changes of internal preceding and concurrent atmospheric circulation, which needs to be further investigated in the future. Meanwhile, it is noted that the multidecadal fluctuation of internal autumn SIC and the corresponding internal circulation and EAWT variability exhibit a quasi-50 year cycle ( figure 4(d)). Our finding may provide a possibility to narrow the uncertainty of near-term or decadal-scale winter climate projection over East Asia based on the quasi-50 year cycle of internal variability. It is noted that the long-term Arctic sea ice loss due to the anthropogenic forcing is expected to be more severe in the future (e.g. Dai et al 2019) which may strengthen or weaken the climate impacts associated with different phase transition of internal sea ice variations. The synthetical climate impacts of forced and internal generated sea ice changes need to be systematically investigated in the future. It is noted that the decline of autumn SIC can induce a negative (positive) SLP (SAT) trend over the Tibetan Plateau (figures S3(c) and S4(c)). We speculate that it may be related to a meridional stationary wave train in the upper troposphere triggered by the sea ice loss (figures S5(c) and (d)), which needs to be further investigated in the future. It should be pointed out that in addition to the internal component of Arctic sea ice, other internal dynamical processes may also contribute to the internal dynamically induced EAWT variations (Wang and Tan 2020), which need to be identified in future studies.