Changes in Cold Surge Occurrence over East Asia in the Future : Role of Thermal Structure

The occurrence of wintertime cold surges (CSs) over East Asia is largely controlled by the surface air temperature (SAT) distribution at high latitudes and thermal advection in the lower troposphere. The thermodynamic background state over northeastern Asia is associated with the strength of the East Asian winter monsoon and the variation of Arctic Oscillation. This study assesses the importance of the SAT structure with thermal advection in determining the frequency of CS occurrences over East Asia through the analysis of nine atmosphere–ocean coupled global climate models participating in the Coupled Model Intercomparison Project Phase 5. The historical simulations can reproduce the observed typical characteristics of CS development. On the basis of this model performance, ensemble-averaged future simulations under the representative concentration pathway 8.5 project a reduction in CS frequency by 1.1 yr−1 in the late 21st century (2065–2095) compared to the present-day period (1975–2005). The major reason for less frequent CSs in the future is the weakened cold advection, caused by notable SAT warming over the northern part of East Asia. These results suggest that changes in the meridional SAT structure and the associated changes in thermal advection would play a more substantial role than local warming in determining future changes in the frequency of CS occurrences over East Asia.


Introduction
Cold surges (CSs), characterized by a steep temperature drop within one or two days, are energetic systems associated with the East Asian winter monsoon (EAWM).Northerly winds bring cold air from high-latitude regions into East Asia, leading to an abrupt drop in local temperature.These winds often penetrate even further south, altering convective activity off the South China Sea and adjacent regions [1,2].Strong CSs are often accompanied by severe weather conditions such as freezing precipitation and heavy snowfall along coastal regions in East Asia (when sufficient humidity is available) [3][4][5].Consequently, the CSs, which occur about ten times per winter [6], have severe impacts on socio-economical human activities in East Asian countries.
The semi-permanent Siberian High develops in association with radiative cooling at the surface and large-scale subsidence motions during boreal winter [7,8].On the eastern flank of the Siberian High, climatological northerlies exist with a tight meridional temperature gradient near the surface.The interaction between lower-and upper-level atmospheric perturbations known as the vertical westward-tilted structure is the primary mechanism leading to amplification of the Siberian High on intraseasonal time scales.Perturbations of atmospheric circulation often provide strong cold advection by bringing pre-existing colder air from northern Eurasia to extend the Siberian High.Along with the development of the Siberian High, the mean state of thermal advection by anomalous winds plays a critical role in cold air outbreaks into East Asia [9,10].The intensity of cold advection is mainly determined by the meridional distribution of mean state temperature.Changes in the meridional distribution of mean state temperature strongly affect both the Siberian High and the occurrence of CSs.
Previous studies have shown that the frequency of CSs is remotely influenced by large-scale climate variabilities such as the Arctic Oscillation (AO) [11], Madden-Julian Oscillation (MJO) [12], and El Niño-Southern Oscillation (ENSO).During the negative phase of the AO, the large-scale circulation brings cold air to East Asia from northern Eurasia, creating an environment that favors strong cold advection [13,14].Jeong and Ho [13] and Woo et al. [15] reported that more frequent, stronger, and longer-lasting CSs occur in decades of a negative phase of AO.When the MJO convective center is located in the Indian Ocean (i.e., MJO phase 2-3), the large-scale extratropical circulation is changed in a way that can amplify CSs [16].The suppression of convective activity in the maritime continent associated with the weak EAWM in El Niño years provides unfavorable conditions for CS occurrences [6,17].These large-scale climate variabilities affect CS occurrence by changing the mean state temperature distribution.In addition, recent studies have reported that sea ice loss in the Barents-Kara sea can modulate atmospheric circulation in Eurasia by direct Rossby wave propagation [18][19][20] or by indirect modulation via stratospheric circulation [21].These remote modulations influence CS frequency through their effects on the cold advection associated with the meridional temperature gradient.
While global surface warming inevitably reduces the number of cold days [22], its effect on the frequency of cold extremes at regional scales remains controversial [23,24].Under global warming, a huge amount of sea ice is anticipated to melt, which would lead to widespread effects at mid-latitudes [18, [25][26][27].For example, Kug et al. [28] found that the Arctic warming associated with reduced sea ice could lead to cold extremes in the extratropics by inducing upper tropospheric wave trains.Arctic amplification, a larger warming in the lower troposphere over Arctic regions than in other regions, may play an important role in the occurrence of CSs by affecting mid-latitude synoptic temperature variability.Francis and Vavrus [25] suggested that Arctic amplification reduces the meridional temperature gradient, resulting in a weaker jet stream.Under weaker jet conditions, mid-latitude synoptic disturbances tend to be more active, and more blocking events occur, which favors both warm and cold extremes [29].At the same time, it is possible that the weaker meridional temperature gradient resulting from Arctic warming could reduce temperature variability locally through a weakening of cold advection from the north [30].
Previous studies have pointed out the local mean surface air temperature (SAT) warming could not be the possible reason for changes of CS occurrences over East Asia in future climate [23].Here, we explore what is an essential factor for changes of CS occurrences over East Asia in various time scale on the basis of both observation and climate-projected simulation of models participated in the Intergovernmental Panel on Climate Change, Fifth Assessment Report.This study provides an insight into how the future global warming influences CSs occurrences over East Asia.
This paper is organized as follows.Section 2 describes the data used in this study, the definition of CSs adopted here, and the temperature advection analysis.In Section 3, we compare the observed vertical and horizontal structure of CSs to the structure simulated by atmosphere-ocean coupled global climate models (CGCMs), and then examine how CS frequency is expected to change in the future and the mechanisms responsible for those changes.Section 4 provides discussion and summary.

Data
To identify CS occurrences in the observational record, the daily mean SAT and mean sea level pressure (SLP) are obtained from the European Center for Medium-Range Weather Forecasts Re-Analysis Interim [31], which has a horizontal resolution of 1.5 • × 1.5  1.These nine CGCMs are selected because they successfully reproduce the observed climatological center of the Siberian High and provide the variables required to analyze the structure and characteristics of CSs (see Figure S1 for details).We use model simulation data from the historical run and the 8.5 W m −2 representative concentration pathway (RCP8.5)emission scenario [32].The historical run is an experiment with all forcing in the 20th century.The RCP8.5 emission scenario is the pathway with the highest greenhouse gas emissions.We select thirty-year periods from both the historical and the RCP8.5 runs: 1975-2005 and 2065-2095, respectively.All CMIP5 data are re-gridded to 3.75 • × 3.75 • (the coarsest resolution among the nine CGCMs) for consistency.The monthly AO index is calculated in each CGCM by projecting the monthly SLP anomalies poleward of 20 • N onto the first leading mode of their empirical orthogonal functions.

Definition of CS Occurrences
The synoptic criteria adopted in this study to detect wintertime CSs over East Asia are similar to those used in previous studies [15,16,33,34].CSs are defined as follows.First, days associated with an expansion of the Siberian High are identified when the daily SLP exceeds 1030 hPa in an anticyclone center.An anticyclone center is defined as the grid point where the SLP is larger than in the eight surrounding grid points.The domain of the Siberian High is the rectangular region defined by 35 • N-55 • N and 90 • E-115 • E (Figure 1a).Second, the SAT drop should exceed 1.5 standard deviations (σ) of the daily winter SAT anomaly within one or two days over northeastern China or Korea (domains defined by 40 1a).These two 5-degree domains can represent East Asia considering that the amplification and expansion of Siberian High have influenced not only SAT domains, but also other East Asian regions.Third, the daily SAT must be below the climatological mean for that day.This condition excludes relatively warm CSs (i.e., the temperature drops abruptly but remains warmer than the climatological average).To separate consecutive CSs, the termination of the CS is defined the first day when the daily SAT anomaly rises above −0.5σafter the CS occurrence.For explaining detection of CS occurrences, the time series of SAT over Korea and mean SLP averaged over Siberian High domain for the period of January in 2009 are shown in Figure 1b,c, respectively.In this period, two CS occurrences are identified, on 10 and 23 January 2009 (Figure 1b).Before the CS occurrence dates, there are expansions of Siberian High (Figure 1c), which would satisfy the first criterion, daily SLP exceeding 1030 hPa.Second, the large SAT drops are identified in CS occurrence dates in which SAT anomaly is below daily climatology as third criterion.Termination dates are identified before recovering climatology of SAT, which separates two CSs.In addition, the reference date of all CS composites is the CS occurrence date.In the CMIP5 model results, we calculate the present-day and future SAT anomaly as deviations from the daily climatology for 1975-2005 and 2065-2095, respectively.To allow comparison with the 30-year future run, σ is calculated over only 30 years of the historical run .Thus, we use the same SAT drop criteria to define CSs in the historical and future runs.

Thermodynamic Interpretation
We analyze thermal advection terms in the temperature tendency equation at 850 hPa when CSs occur.We note that cold advection is essential for amplification of the Siberian High [9] and the occurrence of CSs.Assuming an isentropic surface, the anomalous temperature tendency at 850 hPa is decomposed as follows: where T is the temperature, V is the horizontal wind vector, VA is the vertical advection, and ∇ is the horizontal gradient operator.Each variable is decomposed into the overbar and prime, which denote daily climatological mean and anomaly, respectively.In this study, we estimate the climatological temperature advection by the anomalous wind ( • ∇ ), the anomalous temperature advection by the climatological wind ( • ∇ ), and the anomalous temperature advection by the anomalous wind ( • ∇ ).We neglect the vertical advection (VA) because its contribution is small.

Simulation of the CS Structure
The general characteristic of CS structure is known for vertically westward tilted geopotential height and accompanies cold advection over East Asia resulting from the expansion of the Siberian High [9,15,33].We compare the spatial patterns and vertical structures of observed CSs to those simulated in the CGCMs prior to investigating future changes in CSs and the contribution of thermal advection to those changes.

Thermodynamic Interpretation
We analyze thermal advection terms in the temperature tendency equation at 850 hPa when CSs occur.We note that cold advection is essential for amplification of the Siberian High [9] and the occurrence of CSs.Assuming an isentropic surface, the anomalous temperature tendency at 850 hPa is decomposed as follows: where T is the temperature, V is the horizontal wind vector, VA is the vertical advection, and ∇ is the horizontal gradient operator.Each variable is decomposed into the overbar and prime, which denote daily climatological mean and anomaly, respectively.In this study, we estimate the climatological temperature advection by the anomalous wind (V •∇T), the anomalous temperature advection by the climatological wind (V•∇T ), and the anomalous temperature advection by the anomalous wind (V •∇T ).We neglect the vertical advection (VA) because its contribution is small.

Simulation of the CS Structure
The general characteristic of CS structure is known for vertically westward tilted geopotential height and accompanies cold advection over East Asia resulting from the expansion of the Siberian High [9,15,33].We compare the spatial patterns and vertical structures of observed CSs to those simulated in the CGCMs prior to investigating future changes in CSs and the contribution of thermal advection to those changes.
The geopotential height patterns at 250 and 850 hPa are shown during the dates of CS occurrence in observation (Figure 2a).The upper-tropospheric ridge-trough-ridge geopotential height pattern stretches from the Ural Mountains (50 The ridge-trough-ridge is arranged in a northwest-southeast direction in the upper troposphere and the anticyclonic-cyclonic couplet stretches from eastern flank of the Siberian High to East Asia in the lower troposphere.These upper-and lower-tropospheric synoptic patterns proceed eastward and imply a westward-tilted vertical structure in which the lower-level circulation leads the upper-level circulation by about a quarter of a wavelength.The westward-tilted vertical structures in geopotential height associated with CSs in the ensemble mean and all nine CGCMs are similar to the observed structures [35].However, the magnitude and location of the geopotential heights are slightly different in each model (Figure 2c-k).In ACCESS1-3, BNU-ESM, CanESM2, CMCC-CESM, and GFDL-ESM2G, the magnitude of the upper-level ridge over the Ural Mountains is underestimated (Figure 2c,f,h).The strong lower-level anticyclone over the Lake Baikal is simulated in CMCC-CMS, MPI-ESM-LR, and MPI-ESM-MR (Figure 2g,i,j).In BNU-ESM, the direction of the lower-level anticyclone-cyclone couplet is more northward than the observation (Figure 2d).The geopotential height patterns at 250 and 850 hPa are shown during the dates of CS occurrence in observation (Figure 2a).The upper-tropospheric ridge-trough-ridge geopotential height pattern stretches from the Ural Mountains (50° E-90° E, 50° N-75° N) to the Korea-Japan region (120° E-140° E, 30° N-45° N).The ridge-trough-ridge is arranged in a northwest-southeast direction in the upper troposphere and the anticyclonic-cyclonic couplet stretches from eastern flank of the Siberian High to East Asia in the lower troposphere.These upper-and lower-tropospheric synoptic patterns proceed eastward and imply a westward-tilted vertical structure in which the lower-level circulation leads the upper-level circulation by about a quarter of a wavelength.The westward-tilted vertical structures in geopotential height associated with CSs in the ensemble mean and all nine CGCMs are similar to the observed structures [35].However, the magnitude and location of the geopotential heights are slightly different in each model (Figure 2c-k).In ACCESS1-3, BNU-ESM, CanESM2, CMCC-CESM, and GFDL-ESM2G, the magnitude of the upper-level ridge over the Ural Mountains is underestimated (Figure 2c,f,h).The strong lower-level anticyclone over the Lake Baikal is simulated in CMCC-CMS, MPI-ESM-LR, and MPI-ESM-MR (Figure 2g,i,j).In BNU-ESM, the direction of the lower-level anticyclone-cyclone couplet is more northward than the observation (Figure 2d).3a).The expansion of the positive SLP anomaly to East Asia is associated with intensification of the Siberian High (not shown).The anomalous northerly located between positive and negative SLP anomalies induces strong cold advection into East Asia across the climatological meridional SAT gradient.The magnitude and location of the SLP and SAT anomaly maxima in the ensemble mean results are in good agreement with those in the observations (Figure 3b).In the ACCESS1-3, the positive SLP anomalies are weaker than the negative SLP anomalies, in contrast to the observations and the results from the other CGCMs (Figure 3c).BNU-ESM, CanESM2, and GFDL-ESM2G simulate relatively wide and intense cold anomalies (Figure 3d,e,h).Overall, the historical simulations capture the westward tilt of the geopotential height patterns and cold SAT anomalies in eastern flank of positive SLP anomalies at 850 hPa, indicating that the CGCMs are able to simulate the general structure of CSs.
Atmosphere 2018, 9, x FOR PEER REVIEW 6 of 16 negative SLP anomaly in the western North Pacific are well shown (Figure 3a).The expansion of the positive SLP anomaly to East Asia is associated with intensification of the Siberian High (not shown).
The anomalous northerly located between positive and negative SLP anomalies induces strong cold advection into East Asia across the climatological meridional SAT gradient.The magnitude and location of the SLP and SAT anomaly maxima in the ensemble mean results are in good agreement with those in the observations (Figure 3b).In the ACCESS1-3, the positive SLP anomalies are weaker than the negative SLP anomalies, in contrast to the observations and the results from the other CGCMs (Figure 3c).BNU-ESM, CanESM2, and GFDL-ESM2G simulate relatively wide and intense cold anomalies (Figure 3d,e,h).Overall, the historical simulations capture the westward tilt of the geopotential height patterns and cold SAT anomalies in eastern flank of positive SLP anomalies at 850 hPa, indicating that the CGCMs are able to simulate the general structure of CSs.

Thermodynamic Environment in Association with CS Occurrence
Takaya and Nakamura [36] suggested that a pre-existing cold anomaly over the northern part of the Siberian High is important to the intensification of the anticyclone that is necessary before a cold air outbreak can occur over East Asia.After the intensification of the Siberian High, low-level thermal advection plays a crucial role in the occurrence of CSs over East Asia.Thus, we examine the thermal advection terms at 850 hPa during CS occurrences.We divide the temperature gradient and wind vector terms into a daily climatology term and an anomaly (deviation from the climatology) term to identify the effects of climatological thermal advection on CS occurrences.

Thermodynamic Environment in Association with CS Occurrence
Takaya and Nakamura [36] suggested that a pre-existing cold anomaly over the northern part of the Siberian High is important to the intensification of the anticyclone that is necessary before a cold air outbreak can occur over East Asia.After the intensification of the Siberian High, low-level thermal advection plays a crucial role in the occurrence of CSs over East Asia.Thus, we examine the thermal advection terms at 850 hPa during CS occurrences.We divide the temperature gradient and wind vector terms into a daily climatology term and an anomaly (deviation from the climatology) term to identify the effects of climatological thermal advection on CS occurrences.

Relative Importance of the Temperature Advection Terms
Figure 4 shows the temperature anomaly (Figure 4a,e) and the temperature advection terms in dates of CS occurrences (Figure 4b,d,f,h).The temperature anomaly patterns show that the cold anomaly approaches East Asia, leading to large temperature drops (Figure 4a).The role of each temperature advection term is explored based on its spatial pattern and spatial correlation coefficient with temperature anomaly pattern.When the negative (positive) area is located over East Asia, this temperature advection term contributes cold (warm) advection over East Asia.The climatological temperature advection by anomalous wind pattern (Figure 4b) induces cold advection over East Asia.This advection contributes to sustain and amplify cold anomaly located over East Asia.Spatial correlation coefficient between temperature anomaly and the climatological temperature advection by anomalous wind is 0.67, which means the advection generally amplify temperature anomaly pattern.Anomalous temperature advection by climatological wind induces cold (warm) advection on the eastern flank of cold (warm) anomaly (Figure 4c).It means climatological north-westerly in winter moves the cold and warm anomaly eastward which does not induce cold advection over East Asia.The anomalous temperature advection by anomalous wind induces slightly weak cold (warm) advection into the southern (northern) part of East Asia (Figure 4d).The cooling (warming) in mid-latitude (high-latitude) regions may imply meridional energy exchanges.Each temperature advection term in the ensemble mean shows a similar spatial pattern to the observations and their spatial correlation coefficients are also similar with the observations (Figure 4e,h).Therefore, the climatological temperature gradient plays a critical role for the intensity of the cold advection that eventually leads to the development of CSs by carrying cold air into East Asia.The temperature distribution in the lower troposphere would change CS occurrence through the variation of cold advection.

Relative Importance of the Temperature Advection Terms
Figure 4 shows the temperature anomaly (Figure 4a,e) and the temperature advection terms in dates of CS occurrences (Figure 4b,d,f,h).The temperature anomaly patterns show that the cold anomaly approaches East Asia, leading to large temperature drops (Figure 4a).The role of each temperature advection term is explored based on its spatial pattern and spatial correlation coefficient with temperature anomaly pattern.When the negative (positive) area is located over East Asia, this temperature advection term contributes cold (warm) advection over East Asia.The climatological temperature advection by anomalous wind pattern (Figure 4b) induces cold advection over East Asia.This advection contributes to sustain and amplify cold anomaly located over East Asia.Spatial correlation coefficient between temperature anomaly and the climatological temperature advection by anomalous wind is 0.67, which means the advection generally amplify temperature anomaly pattern.Anomalous temperature advection by climatological wind induces cold (warm) advection on the eastern flank of cold (warm) anomaly (Figure 4c).It means climatological north-westerly in winter moves the cold and warm anomaly eastward which does not induce cold advection over East Asia.The anomalous temperature advection by anomalous wind induces slightly weak cold (warm) advection into the southern (northern) part of East Asia (Figure 4d).The cooling (warming) in midlatitude (high-latitude) regions may imply meridional energy exchanges.Each temperature advection term in the ensemble mean shows a similar spatial pattern to the observations and their spatial correlation coefficients are also similar with the observations (Figure 4e,h).Therefore, the climatological temperature gradient plays a critical role for the intensity of the cold advection that eventually leads to the development of CSs by carrying cold air into East Asia.The temperature distribution in the lower troposphere would change CS occurrence through the variation of cold advection.

Importance of Meridional SAT Gradient and Local SAT Anomaly
The importance of the meridional SAT gradient and the local SAT anomaly for CS occurrences is assessed by a sensitivity test in which the definition of a CS occurrence is varied.At mid-latitudes, the SAT drop is mainly caused by cold advection, and its intensity is related to the meridional SAT gradient.As described in Section 2, in this study, a CS occurrence is defined as an SAT drop of 1.5σ coincident with a negative SAT anomaly.For the sensitivity test, we check how the frequency of CS occurrences changes as we vary the SAT drop criterion and the daily SAT anomaly criterion (Figure 5).We find that in both the observational data and the CGCMs, the change in the frequency of CS occurrences is more sensitive to changes in the SAT drop criterion than to changes in the anomaly criterion.Changing the daily SAT anomaly criterion can be interpreted as seasonal mean warming/cooling in wintertime.For example, if we change the daily SAT anomaly criterion for a CS to be below −1 • C, ideally, that can be regarded as the region experiencing 1 • C of warming in winter.Thus, the finding that the CS occurrence is less sensitive to change in the daily SAT anomaly criterion than to changes in the SAT drop criterion emphasizes that changes in cold advection would have greater effects on CS occurrence than local mean SAT warming in a future climate.

Importance of Meridional SAT Gradient and Local SAT Anomaly
The importance of the meridional SAT gradient and the local SAT anomaly for CS occurrences is assessed by a sensitivity test in which the definition of a CS occurrence is varied.At mid-latitudes, the SAT drop is mainly caused by cold advection, and its intensity is related to the meridional SAT gradient.As described in Section 2, in this study, a CS occurrence is defined as an SAT drop of 1.5σ coincident with a negative SAT anomaly.For the sensitivity test, we check how the frequency of CS occurrences changes as we vary the SAT drop criterion and the daily SAT anomaly criterion (Figure 5).We find that in both the observational data and the CGCMs, the change in the frequency of CS occurrences is more sensitive to changes in the SAT drop criterion than to changes in the anomaly criterion.Changing the daily SAT anomaly criterion can be interpreted as seasonal mean warming/cooling in wintertime.For example, if we change the daily SAT anomaly criterion for a CS to be below −1 °C, ideally, that can be regarded as the region experiencing 1 °C of warming in winter.Thus, the finding that the CS occurrence is less sensitive to change in the daily SAT anomaly criterion than to changes in the SAT drop criterion emphasizes that changes in cold advection would have greater effects on CS occurrence than local mean SAT warming in a future climate.

Influence of the SAT Gradient in Interannual Time Scale
On the interannual time scale, the thermodynamic environment (i.e., mean state of SAT) affects the frequency of CSs. Figure 6 shows the spatial patterns of the interannual correlation coefficients between the mean number of CS occurrences over East Asia and the mean SAT in each winter for both the observations and the ensemble means of the historical and future CGCM runs.In the observations, negative and positive correlation coefficients are located over northern and southern regions, respectively (Figure 6a).The area of negative correlation coefficients spreads zonally from the northern regions of Mongolia to East Asia.If the northern surface air is cold and the meridional SAT gradient steepens, the conditions are favorable for CS genesis.While the regions with negative correlation coefficient are rather different in the historical run and the RCP8.5, both display negative correlation coefficients over the northern part of East Asia (Figure 6b,c).In the historical run, the center of negative correlation coefficient is located over the northeastern regions of the analysis domain, while in the RCP8.5 run, the center is located over the northern regions.The difference in the location of the negative correlation coefficients may be due to the difference between blocking and wave-train CSs [33], which were not considered in this study and there are limitations to simulate the positive correlation coefficient over the southern regions in CGCMs.The meridional SAT gradient over eastern Eurasia has a crucial role for the frequency of CSs in interannual time scale, which could support the importance of thermal structure for predicting the number of CS occurrences in future climate.

Influence of the SAT Gradient in Interannual Time Scale
On the interannual time scale, the thermodynamic environment (i.e., mean state of SAT) affects the frequency of CSs. Figure 6 shows the spatial patterns of the interannual correlation coefficients between the mean number of CS occurrences over East Asia and the mean SAT in each winter for both the observations and the ensemble means of the historical and future CGCM runs.In the observations, negative and positive correlation coefficients are located over northern and southern regions, respectively (Figure 6a).The area of negative correlation coefficients spreads zonally from the northern regions of Mongolia to East Asia.If the northern surface air is cold and the meridional SAT gradient steepens, the conditions are favorable for CS genesis.While the regions with negative correlation coefficient are rather different in the historical run and the RCP8.5, both display negative correlation coefficients over the northern part of East Asia (Figure 6b,c).In the historical run, the center of negative correlation coefficient is located over the northeastern regions of the analysis domain, while in the RCP8.5 run, the center is located over the northern regions.The difference in the location of the negative correlation coefficients may be due to the difference between blocking and wave-train CSs [33], which were not considered in this study and there are limitations to simulate the positive correlation coefficient over the southern regions in CGCMs.The meridional SAT gradient over eastern Eurasia has a crucial role for the frequency of CSs in interannual time scale, which could support the importance of thermal structure for predicting the number of CS occurrences in future climate.

Decadal Variability of CS Occurrence
The importance of cold advection associated with SAT structure to CS occurrence can also be seen on decadal time scales.We divide the observational analysis period of 38 years into three periods: 1979-1987, 1988-1994, and 1995-2016 (Figure 7).The number of CS occurrences in 1979-1987 and 1995-2016 are 10.8 yr −1 and 10.5 yr −1 , respectively.However, fewer CSs occur (8.8 yr −1 ) in 1988-1994.During the period of 1988-1994, the SAT warming over Northern Eurasia decreases the meridional SAT gradient, especially since the center of the warming occurs in the northern part of the Siberian High (Figure 7a).The reduced meridional SAT gradient from the center of Siberian High to East Asia would lead to relatively warmer air by northerly in-weather systems in the region, providing unfavorable conditions for CS occurrence.The decadal changes in CS frequency described in Figure 7b are consistent with previous studies [15,37].In sum, CS occurrence is closely associated with low-level thermodynamic environments, especially the meridional SAT gradient over East Asia on decadal time scale.

Decadal Variability of CS Occurrence
The importance of cold advection associated with SAT structure to CS occurrence can also be seen on decadal time scales.We divide the observational analysis period of 38 years into three periods: 1979-1987, 1988-1994, and 1995-2016 (Figure 7).The number of CS occurrences in 1979-1987 and 1995-2016 are 10.8 yr −1 and 10.5 yr −1 , respectively.However, fewer CSs occur (8.8 yr −1 ) in 1988-1994.During the period of 1988-1994, the SAT warming over Northern Eurasia decreases the meridional SAT gradient, especially since the center of the warming occurs in the northern part of the Siberian High (Figure 7a).The reduced meridional SAT gradient from the center of Siberian High to East Asia would lead to relatively warmer air by northerly in-weather systems in the region, providing unfavorable conditions for CS occurrence.The decadal changes in CS frequency described in Figure 7b are consistent with previous studies [15,37].In sum, CS occurrence is closely associated with low-level thermodynamic environments, especially the meridional SAT gradient over East Asia on decadal time scale.

Projection of Future CS Occurrence Frequency under RCP8.5
We compare the number of CS occurrences in the historical run  and future climate (2065-2095) (Figure 8).The number of CS occurrences per winter is 10.2 yr −1 in the observational record and 8.2 yr −1 in the ensemble mean for the historical  run.All nine CGCMs underestimate CS numbers compared to the observations, which may be caused by coarse resolution of CGCMs and different characteristics of SAT (see Table S1 for details) [23].The CGCMs also project a decrease in CS occurrences in the future.Four CGCMs in particular (ACCESS1-3, CanESM2, MPI-ESM-LR, and NorESM1-M) predict significant decreases in the number of CS occurrences.In the ensemble mean, the frequency of CS occurrences decreases from 8.2 yr −1 to 7.1 yr −1 under global warming.The results of our sensitivity tests suggest that the weakened meridional SAT gradient in the future may be responsible for the decrease in CSs via weakening of cold advection.S1 for details) [23].The CGCMs also project a decrease in CS occurrences in the future.Four CGCMs in particular (ACCESS1-3, CanESM2, MPI-ESM-LR, and NorESM1-M) predict significant decreases in the number of CS occurrences.In the ensemble mean, the frequency of CS occurrences decreases from 8.2 yr −1 to 7.1 yr −1 under global warming.The results of our sensitivity tests suggest that the weakened meridional SAT gradient in the future may be responsible for the decrease in CSs via weakening of cold advection.
The climatological SAT pattern in the historical run shows a typical winter structure in which high-latitude regions are colder than low latitude regions.While warming is clearly seen over the whole Eurasian continent in the future climate (Figure 9), the warming at high latitudes is more pronounced, particularly over eastern Eurasia continent.In the ensemble mean of the future, the zonal mean meridional SAT gradient along 80 • E to 160 • E is expected to decrease by about 1.5 • C per 10 degrees of latitude (Figure 9a).All the CGCMs project similar decreases in the zonal mean meridional SAT gradient.These changes in the SAT distribution are clearly seen in most regions except for high altitude regions and the Sea of Okhotsk.The reduced meridional SAT gradient leads to relatively warmer advection.The warmer advection is unfavorable to CS development, leading to the projected decreases in CS frequency in the future climate.The climatological SAT pattern in the historical run shows a typical winter structure in which high-latitude regions are colder than low latitude regions.While warming is clearly seen over the whole Eurasian continent in the future climate (Figure 9), the warming at high latitudes is more pronounced, particularly over eastern Eurasia continent.In the ensemble mean of the future, the zonal mean meridional SAT gradient along 80° E to 160° E is expected to decrease by about 1.5 °C per 10 degrees of latitude (Figure 9a).All the CGCMs project similar decreases in the zonal mean meridional SAT gradient.These changes in the SAT distribution are clearly seen in most regions except for high altitude regions and the Sea of Okhotsk.The reduced meridional SAT gradient leads to relatively warmer advection.The warmer advection is unfavorable to CS development, leading to the projected decreases in CS frequency in the future climate.
In addition, circulation changes related to the EAWM and AO can affect SAT patterns in the lower troposphere.The EAWM is related to the circulation structure over East Asia, so we examined the projected changes in circulation patterns in the upper and lower troposphere in the future (Figure 10).The East Asian westerly jet stream is projected to shift poleward in the ensemble mean (Figure 10a).CMCC-CESM, CMCC-CMS, MPI-ESM-LR, and MPI-ESM-MR simulate particularly strong poleward shifts (Figure 10e,f,h,i).The enhancement of the jet stream on the northern flank of jet stream's core is a noticeable feature of weak EAWM years [38].Projected changes in zonal wind at 250 hPa are similar to the pattern in weak EAWM years.In the lower troposphere, zonal wind is projected to intensify over the regions where the East Asian westerly jet intensifies.On the other hand, zonal wind at 850 hPa is projected to decrease over low-latitude regions.This anticyclonic change dominates over the western North Pacific around 40° N as the weak EAWM pattern does.Both the upper and lower level circulation changes imply a climatological weakening of the EAWM and a decrease in the number of CS occurrences in the future due to weakened cold advection.In addition, circulation changes related to the EAWM and AO can affect SAT patterns in the lower troposphere.The EAWM is related to the circulation structure over East Asia, so we examined the projected changes in circulation patterns in the upper and lower troposphere in the future (Figure 10).The East Asian westerly jet stream is projected to shift poleward in the ensemble mean (Figure 10a).CMCC-CESM, CMCC-CMS, MPI-ESM-LR, and MPI-ESM-MR simulate particularly strong poleward shifts (Figure 10e,f,h,i).The enhancement of the jet stream on the northern flank of jet stream's core is a noticeable feature of weak EAWM years [38].Projected changes in zonal wind at 250 hPa are similar to the pattern in weak EAWM years.In the lower troposphere, zonal wind is projected to intensify over the regions where the East Asian westerly jet intensifies.On the other hand, zonal wind at 850 hPa is projected to decrease over low-latitude regions.This anticyclonic change dominates over the western North Pacific around 40 • N as the weak EAWM pattern does.Both the upper and lower level circulation changes imply a climatological weakening of the EAWM and a decrease in the number of CS occurrences in the future due to weakened cold advection.The AO can also affect CS occurrences over East Asia by modulating atmospheric circulation [13,14].All CGCMs successfully simulate the AO pattern, in which negative SLP anomalies are located over the Arctic Ocean and positive SLP anomalies are located over the North Pacific and North Atlantic, as determined by empirical orthogonal function analysis (not shown).Regression analysis of the historical run shows that the SAT gradient is reduced from Mongolia to East Asia in the positive AO phase (Figure 11).The projection value between the AO loading pattern and the mean SLP change in future climate is positive in all the CGCMs, with a value of +0.55 for the ensemble mean.

Atmosphere 2018, 9 ,
x FOR PEER REVIEW 4 of 16 there are expansions of Siberian High (Figure1c), which would satisfy the first criterion, daily SLP exceeding 1030 hPa.Second, the large SAT drops are identified in CS occurrence dates in which SAT anomaly is below daily climatology as third criterion.Termination dates are identified before recovering climatology of SAT, which separates two CSs.In addition, the reference date of all CS composites is the CS occurrence date.In the CMIP5 model results, we calculate the present-day and future SAT anomaly as deviations from the daily climatology for 1975-2005 and 2065-2095, respectively.To allow comparison with the 30-year future run, σ is calculated over only 30 years of the historical run.Thus, we use the same SAT drop criteria to define CSs in the historical and future runs.

Figure 1 .
Figure 1.(a) The SAT domains of East Asia (two dashed boxes; 35° N-40° N, 125° E-130° E and 40° N-45° N, 120° E-125° E) and the domain of the Siberian High (solid box; 35° N-55° N, 90° E-115° E) used to define cold surges.The climatology of observed mean SLP (contour interval of 2 hPa) and SAT (shading) during the winter season (November to March).(b) Times series of daily SAT (solid line) and its climatology (dash line) averaged over 35° N-40° N, 125° E-130° E during January in 2009.Red (blue) dots indicate the occurrence (termination) date of CSs.(c) Time series of daily SLP averaged over the domain of Siberian High during the same period mentioned in Figure 1b.

Figure 1 .
Figure 1.(a) The SAT domains of East Asia (two dashed boxes; 35 • N-40 • N, 125 • E-130 • E and 40 • N-45 • N, 120 • E-125 • E) and the domain of the Siberian High (solid box; 35 • N-55 • N, 90 • E-115 • E) used to define cold surges.The climatology of observed mean SLP (contour interval of 2 hPa) and SAT (shading) during the winter season (November to March).(b) Times series of daily SAT (solid line) and its climatology (dash line) averaged over 35 • N-40 • N, 125 • E-130 • E during January in 2009.Red (blue) dots indicate the occurrence (termination) date of CSs.(c) Time series of daily SLP averaged over the domain of Siberian High during the same period mentioned in Figure 1b.

Figure 2 .
Figure 2. Composite of geopotential height anomalies at 250 hPa (contour interval of 30 m; values significant at the 99% confidence level are represented by thick lines) and at 850 hPa (shading; values significant at the 99% confidence level are represented by black dots) for cold surge occurrence dates identified in the ERA-Interim data, the output of the nine CGCMs for the historical run, and the ensemble mean of the CGCMs for the historical run.The number of CS occurrences for composite are represented on the top-right of each panel.

Figure 3
Figure 3 displays the spatial patterns of the SLP and SAT anomalies during CS occurrences in the reanalysis data and CGCM outputs.The positive SLP anomaly in the inner continent and the

Figure 2 .
Figure 2. Composite of geopotential height anomalies at 250 hPa (contour interval of 30 m; values significant at the 99% confidence level are represented by thick lines) and at 850 hPa (shading; values significant at the 99% confidence level are represented by black dots) for cold surge occurrence dates identified in the ERA-Interim data, the output of the nine CGCMs for the historical run, and the ensemble mean of the CGCMs for the historical run.The number of CS occurrences for composite are represented on the top-right of each panel.

Figure 3
Figure3displays the spatial patterns of the SLP and SAT anomalies during CS occurrences in the reanalysis data and CGCM outputs.The positive SLP anomaly in the inner continent and the negative SLP anomaly in the western North Pacific are well shown (Figure3a).The expansion of the positive SLP anomaly to East Asia is associated with intensification of the Siberian High (not shown).The anomalous northerly located between positive and negative SLP anomalies induces strong cold advection into East Asia across the climatological meridional SAT gradient.The magnitude and location of the SLP and SAT anomaly maxima in the ensemble mean results are in good agreement with those in the observations (Figure3b).In the ACCESS1-3, the positive SLP anomalies are weaker than the negative SLP anomalies, in contrast to the observations and the results from the other CGCMs (Figure3c).BNU-ESM, CanESM2, and GFDL-ESM2G simulate relatively wide and intense cold anomalies (Figure3d,e,h).Overall, the historical simulations capture the westward tilt of the geopotential height patterns and cold SAT anomalies in eastern flank of positive SLP anomalies at 850 hPa, indicating that the CGCMs are able to simulate the general structure of CSs.

Figure 3 .
Figure 3. Composite of SLP anomalies (contour interval of 2 hPa; values significant at the 99% confidence level are represented by thick lines) and SAT anomalies (shading; values significant at the 99% confidence level are represented by black dots) for cold surge occurrence dates identified in the ERA-Interim data, the output of the nine CGCMs for the historical run, and the ensemble mean of the CGCMs for the historical run.

Figure 3 .
Figure 3. Composite of SLP anomalies (contour interval of 2 hPa; values significant at the 99% confidence level are represented by thick lines) and SAT anomalies (shading; values significant at the 99% confidence level are represented by black dots) for cold surge occurrence dates identified in the ERA-Interim data, the output of the nine CGCMs for the historical run, and the ensemble mean of the CGCMs for the historical run.

Figure 4 .
Figure 4. Composites for cold surge occurrence dates of the ERA-Interim (a) temperature anomalies (°C), (b) advection of climatological temperature by anomalous wind (°C day −1 ), (c) advection of anomalous temperature by climatological wind, and (d) advection of anomalous temperature by anomalous wind at 850 hPa.(e-h) are same as (a-d) but for the ensemble mean of the nine CGCMs in the historical run.The spatial correlation coefficients between each advection term and temperature anomalies are shown in the top-right of each panel (values significant at 99% confidence level are represented by black stars).

Figure 4 .
Figure 4. Composites for cold surge occurrence dates of the ERA-Interim (a) temperature anomalies ( • C), (b) advection of climatological temperature by anomalous wind ( • C day −1 ), (c) advection of anomalous temperature by climatological wind, and (d) advection of anomalous temperature by anomalous wind at 850 hPa.(e-h) are same as (a-d) but for the ensemble mean of the nine CGCMs in the historical run.The spatial correlation coefficients between each advection term and temperature anomalies are shown in the top-right of each panel (values significant at 99% confidence level are represented by black stars).

Figure 5 .
Figure 5. Changes in the frequency of cold surge occurrences (yr −1 ) with changes in the criteria for the identification of cold surges (the temperature drop criterion is shown in the white bar and the temperature anomaly criterion is shown in the black bar).The horizontal axis is divided by the standard deviation of the winter daily SAT anomaly.

Figure 5 .
Figure 5. Changes in the frequency of cold surge occurrences (yr −1 ) with changes in the criteria for the identification of cold surges (the temperature drop criterion is shown in the white bar and the temperature anomaly criterion is shown in the black bar).The horizontal axis is divided by the standard deviation of the winter daily SAT anomaly.

Figure 6 .
Figure 6.The interannual correlation coefficient between the occurrence of cold surges and the boreal winter SAT for (a) 1979-2016 in ERA-Interim (values significant at the 95% confidence level are represented by black dots), (b) 1975-2005 in the historical run, and (c) 2065-2095 in ensemble mean of nine CGCMs based on RCP8.5 scenario.Regions with same sign in at least 6 CGCMs are represented by black dots.

Figure 6 .
Figure 6.The interannual correlation coefficient between the occurrence of cold surges and the boreal winter SAT for (a) 1979-2016 in ERA-Interim (values significant at the 95% confidence level are represented by black dots), (b) 1975-2005 in the historical run, and (c) 2065-2095 in ensemble mean of nine CGCMs based on RCP8.5 scenario.Regions with same sign in at least 6 CGCMs are represented by black dots.

Figure 7 .
Figure 7. (a) SAT anomalies (°C) during 1988-1994 from climatology in winter season and composite of daily wind anomalies at 850 hPa for cold surge occurrence dates, (b) time series of the number of winter cold surge occurrences.The dashed lines show the mean number of cold surge occurrences during1979-1987, 1988-1994, and 1995-2016.

Figure 7 .
Figure 7. (a) SAT anomalies ( • C) during 1988-1994 from climatology in winter season and composite of daily wind anomalies at 850 hPa for cold surge occurrence dates, (b) time series of the number of winter cold surge occurrences.The dashed lines show the mean number of cold surge occurrences during1979-1987, 1988-1994, and 1995-2016.

3. 3 .
Projection of Future CS Occurrence Frequency under RCP8.5We compare the number of CS occurrences in the historical run (1975-2005) and future climate (2065-2095) (Figure 8).The number of CS occurrences per winter is 10.2 yr −1 in the observational record and 8.2 yr −1 in the ensemble mean for the historical (1975-2005) run.All nine CGCMs underestimate CS numbers compared to the observations, which may be caused by coarse resolution of CGCMs and different characteristics of SAT (see Table

Figure 9 .
Figure 9.The SAT climatology for the historical run (1975-2005; shown as a contour with intervals of 5 • C).Changes in SAT between 2065-2095 and 1975-2005 (shading; left) from nine CGCMs and their ensemble mean.The climatology of zonal mean SAT along 80 • E to 160 • E (black lines) and its change (red lines) in winter from nine CGCMs and their ensemble mean.

Figure 9 .
Figure 9.The SAT climatology for the historical run (1975-2005; shown as a contour with intervals of 5 °C).Changes in SAT between 2065-2095 and 1975-2005 (shading; left) from nine CGCMs and their ensemble mean.The climatology of zonal mean SAT along 80° E to 160° E (black lines) and its change (red lines) in winter from nine CGCMs and their ensemble mean.

Figure 10 .
Figure 10.The climatology of zonal wind at 250 hPa (contour: intervals of 10 m s −1 ) for 1975-2005 and the changes in zonal wind at 250 hPa (shading) and wind vector at 850 hPa (m s −1 ) in winter from nine CGCMs and their ensemble mean.

Figure 10 .
Figure 10.The climatology of zonal wind at 250 hPa (contour: intervals of 10 m s −1 ) for 1975-2005 and the changes in zonal wind at 250 hPa (shading) and wind vector at 850 hPa (m s −1 ) in winter from nine CGCMs and their ensemble mean.
• .Daily temperature, geopotential height, and wind at 250 and 850 hPa are used to identify CS structures.The monthly mean AO index is obtained from the National Oceanic and Atmospheric Administration Climate Prediction Center.We analyze 388 cold surges detected in cold seasons (November through March) of 1979/80-2016/17.We use the results from nine CGCMs belonging to Coupled Model Intercomparison Project Phase 5 (CMIP5): ACCESS1-3, BNU-ESM, CanESM2, CMCC-CESM, CMCC-CMS, GFDL-ESM2G, MPI-ESM-LR, MPI-ESM-MR, NorESM1-M.Brief model descriptions are given in Table

Table 1 .
Details of CMIP5 models used in this study.