Controlling Factors of Historical Variation of Winter Tibetan Plateau Snow Cover Revealed by Large‐Ensemble Experiments

This study investigated the factors controlling the interannual variability of Tibetan Plateau snow cover (TPSC) in winter using large‐ensemble simulations. Composite analysis reveals that years with high TPSC are associated with a positive Arctic Oscillation (AO)‐like pattern in which a strengthened subtropical westerly enhances zonal water vapor flux, resulting in greater precipitation over the Tibetan Plateau. Observations and simulations commonly show that TPSC tends to be high during the positive phase of AO, and vice versa. Meanwhile, the magnitude of the contribution of tropical sea surface temperature (SST) to TPSC is different based on observations and simulations. Observations indicate that TPSC is insensitive to El Niño/La Niña events, whereas the simulated ensemble mean TPSC is correlated strongly with Niño 3.4 SST index. In El Niño years, higher specific humidity and more numerous synoptic disturbances are found around the Tibetan Plateau. As several ensemble members show features close to those observed, the observed insensitivity of TPSC to El Niño–Southern Oscillation (ENSO) events is considered attributable to the limited sample size of the observations. This is corroborated in the analysis for a relatively short period, in which the observed TPSC is sensitive to ENSO. Comparison of the results of the historical simulations with those of the simulations without warming suggests that global warming has reduced the snow‐to‐rain ratio over the Tibetan Plateau, which is compensated by the overall increase in precipitation. Consequently, the impact of global warming on TPSC was negligibly weak until the 2000s.

The controlling factors of winter TPSC have been discussed in several previous studies. Based on satellite-derived measurements, Z. Wang et al. (2019) reported that the interannual variation of central-eastern TPSC is explained by the Arctic Oscillation (AO), that is, AO can induce an anomalous circulation pattern that controls winter TPSC. A similar linkage between the Arctic and in situ snow depth and snow cover days on the TP was suggested by both Lü et al. (2008) and Zhang et al. (2019). Meanwhile, Jiang et al. (2019) analyzed in situ snow depth and argued that the effect of oceanic forcing such as the El Niño-Southern Oscillation (ENSO) on TPSC is not negligible. This feature was supported by several other observation-based studies (e.g., Bao & You, 2019;Shaman & Tziperman, 2005;Wang & Xu, 2018). However, these studies were based on in situ or satellite-derived measurements acquired over periods of up to 48 years. This could limit the availability of strong AO and ENSO cases, which would hamper statistical analysis of the response of TPSC to such large-scale forcing. Furthermore, the consequences of the combined impact of AO and ENSO will remain unresolved unless adequate cases of AO and ENSO become available. In this study, we analyzed the output of a large-ensemble general circulation model experiment to show that AO and ENSO are the dominant large-scale controlling factors leading to anomalous winter TPSC. During 1955During -1996, the TP experienced significant warming at a rate substantially higher than that of the midlatitudes in the Northern Hemisphere . In contrast to the consensus regarding this strong warming, validation of an observed trend in TPSC remains controversial. Based on statistical analyses of observational data, some studies have suggested that the winter snow cover over the TP and the Himalayas has declined significantly in recent years (e.g., Li, Su, et al., 2018;Ma et al., 2019;T. Wang et al., 2013), whereas other studies have reported the opposite trend (e.g., Tan et al., 2019;X. Wang et al., 2017). Thus, through analysis of a set of large-ensemble experiments with and without consideration of anthropogenic effects, we also investigated how global warming might alter TPSC and related hydrometeorological regimes.

Observation Data
This study used data of daily snow cover over the TP provided by Chen et al. (2018). This data set is called the "composite long-term gap-filled Tibetan Plateau daily snow cover extent record," which provides satellite-based 5-km-mesh snow cover data during 1982-2016. It is recognized that the satellite data contain uncertainty regarding estimation of the snowpack (Brown et al., 2010;Pulliainen et al., 2020;Yang et al., 2015). Thus, verification was conducted using in situ observation data from 57 stations. The correlation coefficient of the interannual variation during 1982-2013 is 0.77 (p < 0.01; see Figures S1 and S2), which means that the satellite-based data capture the main features of the interannual variation of TPSC.
Gridded (0.5° × 0.5°) monthly 2-m air temperature data were obtained from the air temperature database of the University of Delaware (USA), and monthly precipitation data (2.5° × 2.5°) were obtained from the Global Precipitation Climatology Project (GPCP V2.3; Adler et al., 2003). Circulation and snow-hydrological parameters were derived from the Japanese 55-year Reanalysis data sets (JRA-55; Kobayashi et al., 2015). The analyzed variables comprise geopotential height, wind speed, specific humidity, surface snowfall water equivalent, surface total precipitation, and total column water vapor flux with a spatial resolution of 1.25° × 1.25°. The ERA5 reanalysis data set (Hersbach et al., 2020) with a 0.25° × 0.25° spatial resolution and NCEP/NCAR reanalysis data set (Kalnay et al., 1996) with a 2.5° × 2.5° spatial resolution is also used to validate the results. We confirmed that the results remain robust irrespective with the selection of the reanalysis data set. Sea surface temperature (SST) data were obtained from the Centennial Observation-Based Estimates of SST Version 2 (COBE-SST2; Hirahara et al., 2014). The data mentioned here were averaged over our target period of midwinter to early spring (i.e., from January through March; JFM), and it was these three-monthly averaged fields that were mainly analyzed in this study.

Simulation Data
As mentioned in Section 1, the observation and reanalysis data cover a period of only a few tens of years. This limitation might induce a distorted interpretation of the relationship between TPSC and its controlling factors. Hence, in addition to the observation data, large-ensemble simulations were also used in this study. We used the monthly output from two 100-member ensemble simulations obtained from a database called the "Database for Policy Decision-Making for Future Climate Change" (d4PDF; Mizuta et al., 2017). The simulations used considered the historical climate and a counterfactual historical climate for the period 1951-2010. The historical climate experiment (HIST) was forced by observed SST and sea ice from the COBE-SST2 at the lower boundary conditions and by observed external forcing agents such as greenhouse gases (GHGs), aerosols, and ozone. The counterfactual historical experiment (NAT) was forced by SST excluding the long-term trend, while GHGs, aerosols, and ozone were fixed to values corresponding to the preindustrial period. The ensemble member spread represents the effects of the internal variability of the atmosphere and land coupled system, while the ensemble means can be considered to represent the forced response to SST and other external forcings. Here, the HIST experiment is analyzed comprehensively, and the impact of global warming is investigated through comparison of the ensemble means of the HIST and NAT experiments. The readers may anticipate that the AGCM has limitations in simulating snow cover over complex terrains. We therefore carefully discussed the validity of simulated TPSC in d4PDF product in Section 4.2.

Indices
To assess the variability of TPSC, a TPSC fraction (hereafter, TPSCF) is defined as the percentage of the snow-covered area relative to the entire area of the TP (i.e., the area with an elevation of >3,000 m above mean sea level). For each of the observations and the model, years with high and low values of TPSCF were selected by applying threshold values corresponding to ±0.5 times the standard deviation of the interannual TPSCF variation. To highlight the atmospheric and oceanic conditions of high and low TPSCF years, composite analysis is conducted for the selected years in Section 3.2.
The AO index was calculated as the first principal component score of the empirical orthogonal function analysis of monthly mean sea level pressure anomalies poleward of 20°N (Wallace & Gutzler, 1981). To diagnose ENSO status, the Niño index was defined as the first principal component score of the empirical orthogonal function analysis of monthly SST anomalies over 30°S-30°N, 120°E−50°W (Ashok et al., 2007). It should be noted that in computing these indices, year-round monthly mean anomaly data were used with latitudinal weighting, and the anomaly was defined as the departure from the climatology in each member during 1951-2010. The calculated Niño indices in HIST show reasonable values in comparison with the observations, and the AO indices in HIST were also reasonable in comparison with those obtained from the JRA-55, ERA5, and NCEP/NCAR reanalysis data set (see Figure S3).

Model Performance in HIST
Observed TPSCF shows large interannual variation (Figure 1; black line) in the range of 28.47%-45.76%. The interannual variation of observed TPSCF looks markedly different to that of the ensemble mean of HIST during 1982-2010, as evidenced by their low correlation (r = 0.11). This low correlation suggests that historical TPSCF is strongly controlled by atmospheric internal variability rather than by anthropogenic and/or natural forcing. This is because the HIST was driven by using only SST, GHGs, and other external forcings, and thereby the ensemble mean is likely to represent the forced response of TPSCF to these forcings. In reality, TPSCF will also be strongly affected by atmospheric internal variability, such as AO, as evidenced by the large ensemble spread in Figure 1. Meanwhile, the ensemble mean of HIST for low (e.g., the year of 1999) and high (e.g., the year of 1998) cases deviates significantly from the long-term average TPSCF, indicating that SST also influences TPSCF (Figure 1). The similarity between the ensemble mean TPSCF of HIST and NAT implies that the impact of global warming on midwinter TPSCF was negligibly weak until 2010. This interpretation is consistent with the observations for which the linear trend is insignificant at 95% confidence level (p = 0.085) for the period 1982-2010, although significant decline is detected for the period extended to recent years (1982-2016, p < 0.01). The influence of global warming on TPSCF will be revisited in Section 4.1.
To examine the spatial variation of atmospheric patterns related to TPSCF, high and low TPSCF years were selected by applying thresholds of ±0.5 times the standard deviation of the JFM-mean TPSCF. For the observations, 12 and 9 years were selected as representative of high and low TPSCF years, respectively, while 2033 and 1683 years were selected correspondingly for HIST. The spatial distribution of snow cover over the TP is displayed in Figure 2. Generally, the snow cover fraction is high in western and central-eastern parts of the TP; the difference in snow cover fraction between high and low years is evident in these areas (Figures 2a-2c). The averaged observed TPSCF value for high and low years is 43.28% and 32.37%, respectively. It can be seen that TPSCF is overestimated in HIST (Figures 2d and 2e). The bias is larger over the western TP where the terrain is more complex, which probably reflects the inadequacy of the model resolution (60 km) to simulate snow processes. Nevertheless, the difference in TPSCF between high (82.12%) and low  (c and f) difference between high and low years. TPSCF value is shown in the upper-right corner of each panel. Stippled regions indicate that snow cover fraction difference is statistically significant exceeding the 95% confidence level using the two-sided Student's t-test.
(70.83%) years, corresponding to the range of one standard deviation of interannual TPSCF variation, is 11.29% (Figure 2f), which is very close to that of the observations (10.91%; Figure 2c). Hence, it is considered that TPSCF was simulated reasonably well in the HIST experiment in terms of the range of interannual variation.
Quantitative comparison between the observations and HIST was conducted to further validate the simulated local hydroclimatic processes. The hydroclimatic parameters were averaged over the TP region where ground elevation is >3,000 m. The difference between high and low TP-SCF years was calculated for the following hydroclimatic parameters: 2-m air temperature, precipitation, and convergence of water vapor flux. The results show that the ranges of interannual variation for these parameters are generally simulated well in HIST, except HIST underestimates the difference in water vapor convergence (Table 1). Meanwhile, lower temperature in simulated high TPSCF years causes a higher snow to rain ratio (0.34%) than that in low TPSCF years.

TPSC and Its Linkage to Large-Scale Factors
We investigated the large-scale atmospheric circulation anomalies that characterize high and low TPSCF years. In high TPSCF years, a large-scale wavy pattern is found at 500 hPa in the high-latitudes with significant positive anomalies over the northern U.S., Atlantic-Europe, and Eastern Siberia and significant negative anomalies over Arctic regions and areas extending from the Middle East to the TP (Figures 3a and 3d). Water vapor flux exhibits an enhanced cyclonic circulation around the TP in high TPSCF years (Figures 3b and 3e), carrying more water vapor to the south of the TP. It is speculated that the enhanced water vapor transport causes more snowfall over the TP. This is similar to the claim of Bao and You (2019) Note. Values significant at the 95% and 99% levels are denoted by ** and ***, respectively. Figure 3. Composite maps of (a) geopotential height at 500 hPa, (b) precipitation (color) and total column water vapor flux (arrow), and (c) SST with respect to difference between high and low years of observed JFM TPSCF. (d-f) Same as (a-c) but for HIST. Note that high and low years for HIST are based on TPSCF of individual members for (d and e) but on the ensemble mean TPSCF for (f). Stippled regions and thick black vectors denote differences that are statistically significant exceeding the 95% confidence level using the two-sided Student's t-test.

Table 1 Comparison of JFM Hydroclimatic Parameters Averaged Over the TP Region (>3,000 m) Between High and Low TPSCF Years for the Observations and HIST
that the westerly jet is one of the crucial causes of snow depth variation over the TP. The circulation difference over the midlatitudes of the Eurasian continent resembles a typical AO-induced pattern (see Figure S4). This indicates that TPSCF is likely modulated by AO through change in the moisture transport, for which both the observations and HIST show consistent patterns. Meanwhile, there is disagreement between the observations and HIST regarding SST. The observations suggest that SST in the eastern tropical Pacific does not vary in accordance with TPSCF ( Figure 3c). However, in HIST, the SST difference between high and low TPSCF years is significant, reflecting the El Niño-like SST pattern for high TPSCF years. The SST in the Indian Ocean is also significantly higher for high TPSCF years.
Latitude-height cross sections of specific humidity averaged over 66°-104°E, based on observations and HIST, are presented in Figures 4a and 4b, respectively. In both cases, higher specific humidity is found over and around the TP in El Niño years than in La Niña years. The higher humidity during El Niño years is evident from the surface to the midtroposphere over South Asia (Figures 4a and 4b). The reanalysis data show that storm activity during El Niño years tends to be high around western and southern parts of the TP, as demonstrated by the greater anomaly in the kinematic energy of the transient eddy at 200 hPa ( Figure 4c). These results provide physical support for our earlier suggestion that TPSC tends to be higher in El Niño years. Shaman and Tziperman (2005) also highlighted that upper-tropospheric potential vorticity reaching the TP is greater during El Niño years. Through this process, snowfall over the TP is likely to be increased in El Niño years.
Based on the above analysis, the potential impacts of AO and ENSO on TPSCF can be proposed (i.e., Figures 3, 4 and S4). The role of AO is found to be dominant in terms of both observations and HIST. The positive AO-like atmospheric circulation, which strengthens the subtropical jet, leads to the enhancement of the water vapor flux toward the TP and increased regional precipitation. Meanwhile, higher specific humidity and greater storm activities are found in El Niño years. Another important result is the disagreement in the SST patterns between the observations and HIST (Figures 3c and 3f). As this difference could be attributable to the limited number of observations, the analysis is extended to a probabilistic approach. In the following subsection, assessment of the impacts of AO and ENSO is presented through analysis of individual ensemble members.

Impact of AO and ENSO on TPSC
Modulation of the TPSCF anomaly for AO and ENSO phases was investigated to elucidate their respective impact on TPSCF. The probability distribution of TPSCF anomalies under different AO and ENSO conditions is presented in Figure 5. The TPSCF anomalies in HIST are defined as deviations from the climatological TPSCF in each member. The probability distribution of TPSCF anomalies for the observations shifts to the positive TPSCF side with increase in the AO index ( Figure 5a). Positive TPSCF anomalies appear preferentially when the AO index is >0.5, and vice versa for negative TPSCF anomalies. Similarly, the TP-SCF anomalies in HIST present the same feature, although the shift of the probability distribution is smaller than that for the observations. The weaker dependency of TPSCF in HIST to AO is presumably attributable to the larger sample size used in HIST, which enabled the collection of many AO cases. Conversely, the observed relationship between TPSCF and AO is thought to emphasize the limited availability of cases in the study period .
The probability distributions of observed TPSCF for El Niño and La Niña events appear to overlap. Hence, the observations indicate that TPSCF is insensitive to ENSO phase, although ENSO is known to bring anomalous snowfall over the TP through generation of abnormal vorticity (Shaman & Tziperman, 2005;Wang & Xu, 2018). In contrast with our observational results, HIST shows that the probability distribution of simulated TPSCF differs markedly between El Niño and La Niña years. This gap between the observations and HIST is also recognized in the composite analysis (Figures 3c and 3f). To examine whether the limited number of observation samples could distort the impact of ENSO on TPSC in the observations, the analysis was extended to each member. There are members that have a probability distribution similar to that of the observations, that is, TPSCF anomalies are insensitive to ENSO phase. We confirmed that a majority of the ensemble members exhibit a shift of TPSCF to the higher side in El Niño years (see Figure 6 for selected members and Figure S5 for all members). Therefore, the seemingly different response of TPSCF anomalies to ENSO between the observations and HIST is probably due to the limited number of cases in the observations rather than model bias. It is concluded that the observations tend to reflect stronger linkage to AO and weaker linkage to ENSO owing to the limited number of cases. Analysis of the large-ensemble members reveals that TPSCF is sensitive to both AO and ENSO, but with stronger impact from ENSO than AO.
The large ensemble in d4PDF enables analysis of the combined influence of AO and ENSO on TPSCF.
Owing to the insufficient length of the observations, it is difficult to attribute the TPSCF anomaly to the combined impact of AO and ENSO. The low TPSCF years appear slightly concentrated in the third quadrant (i.e., negative AO and La Niña years) (Figure 7a). For HIST, many samples of the TPSCF value can be obtained for an arbitrary combination of AO and Niño indices. The TPSCF anomalies accounting for the combination of AO and Niño indices were computed using all members of the HIST experiment for the period 1951-2010 (i.e., in total: 60 years × 100 members = 6,000 samples; Figures 7b and 7c). It can be seen from Figure 7 that strong positive TPSCF appears when positive AO coincides with El Niño. Similarly, strong negative TPSCF anomalies can be found when negative AO coincides with La Niña. Meanwhile, in the El Niño condition, TPSCF anomalies tend to be larger than zero, even when coupled with negative AO, and vice versa for the La Niña condition. Hence, oceanic forcing (i.e., El Niño/La Niña) is considered to have greater influence than atmospheric internal variability (i.e., AO) on TPSCF. Regarding the analytical results of large-ensemble experiments, the hydroclimate properties behind the collaborative effect of AO and ENSO on TPSC are further investigated. In both observation and HIST, El Niño can bring greater precipitation to TP, especially when it coincides with positive AO, and vice versa for La Niña condition (Figure 8a). Moreover, neither ENSO nor AO seems to modulate snow cover by regulating air temperature over TP, since the detrended 2-m air temperature anomalies distributed independently to Niño index and AO index in both observations and HIST (Figure 8b). Therefore, the collaborative effect is mainly determined through the precipitation changes rather than the air temperature changes. The results would support the collaborative impacts of AO and ENSO we proposed, considering that the enhancement of moisture and water vapor fluxes has a strong influence on TP in their positive and warm phases.

Impact of Global Warming
It can be seen from Figure 1 that the interannual variations of the ensemble mean TPSCF for HIST and NAT are very similar. As the difference between these experiments is the external forcing agents, this result suggests that the impact of global warming on winter TPSCF is weaker than that of natural variability.  Indeed, the difference in the JFM averaged warming rate (HIST minus NAT) over the TP is 0.14°C/decade (p < 0.01), which is slightly lower than the trend of observed air temperature (0.20°C/decade, p < 0.01) during 1951-2010. The reason behind the weak sensitivity of winter TPSC to global warming in the model might be complicated. The form of precipitation, rain or snow, depends largely on ambient air temperature (Screen & Simmonds, 2012). Higher atmospheric temperatures reduce the snow-to-rain ratio (Bintanja & Andry, 2017;Räisänen, 2008), which tends to reduce snowfall. However, it has been reported that increasing air temperatures globally tend to increase the atmospheric water vapor content and moisture transport, allowing more precipitation to occur (e.g., Salzmann, 2016). Consequently, increased air temperature also has the effect of increasing precipitation, which could lead to greater amounts of snowfall. Owing to these compensating mechanisms, it is unclear whether a warming climate suppresses or enhances the total amount of snowfall over the TP. To investigate the response of snowfall over the TP to air temperature, the snowfall fraction to total precipitation was calculated. The snowfall fraction is defined as the ratio of the area-mean JFM snowfall to total precipitation over the TP. During 1958-2010, JRA55 and HIST show a similar feature in the response of the snowfall fraction to increasing 2-m air temperature ( Figure 9). We confirmed this relationship holds true for ERA5 reanalysis. For a 1°C increase in 2-m air temperature, the snowfall fraction  decreases by 0.46% and 0.40% for JRA55 and HIST, respectively. Hence, it is considered that HIST captures this relationship well.
The ensemble means of precipitable water, total precipitation, snowfall fraction, and snowfall in HIST and NAT were compared to examine the response of hydrological components to a warming climate over the TP. It was found that precipitable water and total precipitation both increase with increasing air temperature (Figures 10b and 10c, respectively). As discussed earlier, the temperature rise significantly reduces the snow-to-rain ratio ( Figure 10a); therefore, more precipitation falls as rain. The reduced snowfall is compensated by the increase in total precipitation, which in turn leads to greater snowfall in HIST than NAT (Figure 10d). This mechanism helps ensure that a certain amount of precipitation still falls as snow, even under the effects of global warming. This finding confirms the results of in situ analysis by Deng et al. (2017).
A strong decline in observed TPSCF is found in recent years (Figure 1). It is presumably not caused by AO or ENSO effects because recent low TPSCF years fall into the four different quadrants in Figure 7a, meaning that recent AO and ENSO phases are not fixed. After removing the linear trend of observed TPSCF, the composited geopotential height at 500 hPa still has a negative anomaly centered over the Arctic region and another negative anomaly branch that extends from the Middle East to the TP (figure not shown), that is, the circulation pattern leading to high TPSCF is robust despite the decline of TPSCF in recent years. As the contribution of global warming to TPSCF during 1951-2010 is negligible (Figure 1), it is unlikely that the recent decline of TPSC is driven by global warming. Recent studies highlighted that light-absorbing aerosols deposited on the snow cover over the TP are found to affect the regional climate through snow albedo feedback (e.g., Yao et al., 2019;Ma et al., 2019). The snowmelt efficiency of the snow cover fraction induced by black carbon deposition on the snowpack is equally as high as that induced by CO 2 increase in January and March (Qian et al., 2011). Moreover, the dust-induced snowmelt efficacy is also comparably high with that of black carbon in the Himalayas (Sarangi et al., 2020). These aerosol-linked impacts might be captured by the observation data, resulting in reduction of the observed TPSCF after 2010 (Figure 1; e.g., Li, Guo, et al., 2018). As sophisticated aerosol processes were not incorporated in the global climate model used for the HIST experiment (Yukimoto et al., 2012), HIST might underestimate the effect of anthropogenic activities. Our understanding of the role of aerosol deposition on snow remains limited; therefore, effort will be needed to improve model performance toward precise assessment of the impact of aerosols.

Assessment of the Simulated TPSC
We noticed that the simulated TPSC in HIST contains some biases, as mentioned in Section 3.1 (Figures 2d  and 2e). To further assess the model biases, three additional observed snow cover data sets, as well as the data set provided by Chen et al. (2018;hereafter C18), are analyzed. The data used are Northern Hemisphere daily 5-km snow cover extent product (named as JASMES; Hori et al., 2017), daily snow cover data at a 24-km resolution obtained from the Interactive Multi-Sensor Snow and Ice Mapping System (IMS; Helfrich et al., 2007), and Northern Hemisphere EASE-Grid 2.0 Weekly Snow Cover and Sea Ice Extent, Version 4.0 (NSIDC-0046; Brodzik & Armstrong, 2013). Observational snow cover data sets were converted to the 60km spatial resolution to match with the model grids. A comparison of the overlapped period of 1982-2010 JFM climatology TPSC was performed between the observational snow cover data sets and 100 ensemble members. For IMS, the climatology was defined for the period 1999-2010. The analysis was performed for 1054 pixels/grids over the Tibetan Plateau (i.e., the area with an elevation of >3,000 m above mean sea level).
The Taylor diagram is presented in Figure S6 to summarize the comparison of snow cover climatology. The spatial correlation of JASMES, IMS, and NSIDC-0046 against C18 is 0.56, 0.65, and 0.69 with normalized root mean square error (NRMSE) of 0.33, 0.28, and 0.24, respectively. This result implies the observational snow cover data sets contain uncertainties. The correlation coefficients of spatial TPSC distribution between observation and 100 ensemble members are all above 0.26, with a median value of 0.29 ( Figure S6). The normalized standard deviations of climatological TPSC in ensemble members fall within the range of observational products indicating the magnitude of the simulated and observed TPSC variations are close to each other. The NRMSE in simulations ranging from 0.50 to 0.53 is greater than those in observations, which denotes the spatial contrast of simulated TPSC is higher in simulations than observation (Figures 2a, 2b, 2d and 2e).
To prove that the model has the potential to produce a reasonable magnitude of the interannual variations, the variance of interannual TPSC variation was drawn as Figure 11. The figure shows although there is a gap among data sets and models, the HIST experiment can produce an essential part of the TPSC variations, especially for central-eastern TP, where the magnitude of interannual variance is close to the observational snow cover data set. As for western TP, the simulated interannual variance of snow cover is similar to JASMES and NSIDC-0046, but it is different from C18 and IMS.
In conclusion, the simulated TPSC shows reasonable interannual variability, which appears to be comparable to the observed values, especially in the central-eastern part of the TP. There is a noticeable limitation of simulated snow cover in terms of its climatological distribution. Meanwhile, the observational snow cover data sets also contain uncertainties. This means our understanding of the real snow cover distribution is yet limited, and hence, the intercomparison of observational product will be needed. Figure 11. Spatial distribution of JFM snow cover fraction (climatology; contours) and its interannual standard deviation (color shading). The observation snow cover data sets were calculated from 1982 to 2016, except that 1999-2016 was used for IMS due to the data availability. The simulated snow cover was calculated from 6,000 samples in HIST (60 years with 100 members).

TPSC and Its Linkage to SST Forcing
As presented in Section 3.3, the simulated TPSC is found to be sensitive to ENSO. In this section, the maximum covariance analysis (MCA) is conducted using interannual variation of the ensemble mean JFM SST and TPSC to study how much proportion of total variance can be explained by ENSO. The MCA (also referred to as singular value decomposition) looks for patterns in two spatially and temporally varying data sets. The first MCA mode for TPSC exhibits a pronounced positive anomaly pattern over the whole TP, especially over the eastern and central TP region, coupled with an El Niño-like SST pattern (Figure 12b). The squared covariance fraction of the first mode is 84.88%, suggesting that ENSO is among the strongest mode that co-varies with TPSC and accounts remarkably for the total variance. The correlation coefficient of principle components between the two fields is 0.72 (p < 0.01) (Figure 12c). The snow cover fraction over TP corresponding to the first mode shows significant positive values (Figure 12a). The principal components of TPSC and SST in the first MCA mode are highly correlated with Niño 3.4 index, with the slopes of linear regression against the Niño 3.4 index being 1.05 for SST and 0.70 for TPSC (p < 0.01; Figure 13), respectively. This result provides evidence that in HIST, the TPSC has a strong linkage with ENSO.
As for the observed ENSO-TPSC relationship, a composite analysis of JFM TPSC was conducted for highlighting the regional snow cover distribution inside of the TP in strong El Niño years and strong La Niña years during 1982-2016. The selected years were denoted in Figure S7. The three snow cover products (C18, JASMES, and NSIDC-0046) were analyzed, and IMS was not used because of its short time coverage. The snow cover anomalies in response to ENSO show a dipole pattern in TP with less snow cover fraction in the northern part of TP and more snow cover fraction in the southern part of TP in El Niño years. The anomaly of snow cover fraction averaged over whole TP is nearly zero, which is a common feature for all three observational products ( Figure S8). The result suggests that TPSCF might be insensitive to ENSO, and it seems to contradict what HIST shows. The following analysis will demonstrate that this gap is related to the period of analysis.
The same analysis as presented above was conducted for a slightly shorter period of 1982-2010 for which the HIST data are available. The result shows the snow cover fraction over whole TP is turned to be higher in El Niño years than in La Niña years, with a stronger positive signal in the southern part of TP than that in the northern part of TP ( Figure S9). These results indicate that the insensitivity of observed TPSCF (or north-south dipole pattern) to ENSO during 1982-2016 is largely affected by the three ENSO events after 2010 ( Figure S7). This possibly implies that the TPSC-ENSO connection in observation has changed in recent years. The 15-year running correlation coefficients between Niño 3.4 index and snow cover fractions Finally, probability distributions comparing observed and simulated TPSC's responses to ENSO are evaluated like Figure 5 but separately for whole TP and southern TP ( Figure 15). The observed snow cover fraction over the whole TP is sensitive to ENSO during 1982-2010 ( Figure 15a). Through the comparison of observation-based probability distributions between Figures 15a and 5b, the TPSC-ENSO linkage is dependent on the period of analysis. The observed ENSO signal is more pronounced over southern TP (Figure 15b). The  simulations also show greater ENSO's impact of snow cover over southern TP than the whole TP, indicating a strong regional dependence of EN-SO's impact. The detailed comparison of the ENSO's impact on regional snow cover anomalies in TP again confirmed the benefit of large-ensemble experiments because the uncertainty increases with narrowing the target area, which makes the detection of climatic forcing impacts more difficult in observation. Finally, the influence of ENSO and AO on TPSC is illustrated in Figure 16, which also includes the magnitudes of their impacts on TPSCF based on the large-ensemble experiments.

Conclusions
This study investigated the controlling factors that lead to TPSC variation. First, the atmospheric circulation patterns leading to years with high/low TPSCF were examined. The observed high TPSCF years are characterized by a positive AO-like circulation pattern that intensifies the subtropical jet and enhances water vapor transport toward the TP, forming conditions favorable for winter snowfall on the TP. El Niño supplies more humidity to the TP and induces more storm activity to the south of the TP, which also brings more snowfall. To ensure the robustness of the mechanism suggested by the observations, the results of large-ensemble experiments were analyzed. In the high TPSC years in HIST, a similar circulation pattern was found over the Eurasian continent, with negative geopotential height anomalies over the Arctic and regions from the Middle East to the TP.
Using a large-ensemble data set, the probability distribution of TPSCF was created for different AO and ENSO phases. For both observations and HIST, positive TPSCF anomalies appear preferentially when the AO index is positive and higher, and vice versa. Despite the insensitivity of the observed TPSCF to ENSO phase, the simulated TPSCF in HIST showed a higher probability of enhanced TPSCF for El Niño years. This disagreement between the observations and HIST is likely attributable to the small sample size of the observations. The probability distribution of  TPSCF against ENSO phase was examined for individual ensemble members in HIST. Interestingly, several members showed that TPSCF anomalies are insensitive to El Niño/La Niña events, similar to the observations. Thus, this discrepancy regarding the impact of ENSO on TPSCF between the observations and the model is likely attributable to the limited sample size of the observations. Moreover, the combined impact of AO and ENSO on TPSCF was quantified. It was found that TPSCF anomalies are greater and positive under strong positive AO and El Niño conditions, and vice versa for strong negative TPSCF anomalies. This is strongly linked to greater precipitation in the El Niño condition than La Niña, in particular when El Niño coincides with positive AO events. The MCA analysis for HIST shows that the ENSO is the most dominant oceanic forcing for the interannual variation of JFM TPSC. Their relationship, however, appears to be weakened and even reversed in the recent observation, the reason behind which is an issue for future studies.
The thermodynamic relationship between precipitation and temperature was compared between HIST and NAT. We found that historical global warming has reduced the snow-to-rain ratio; however, precipitable water and total precipitation have tended to increase. Owing to this compensating effect, the impact of global warming on TPSCF is negligible in the simulations. Comprehensive understanding of the key controlling factors that regulate TPSC during midwinter to early spring is beneficial for improving the predictability of TPSC, even the East Asian summer monsoon in following seasons. This topic merits further examination of associated environmental issues in and around the TP that might reflect TPSC. Variation of TPSC has considerable influence on high-elevation human activities through water availability and flood management at the regional scale. Therefore, further investigation is needed on local-scale changes in snow cover on the TP or other high-elevation regions, which might be strongly affected by elevation or topography. The large-ensemble simulations we used may have limitations in simulating snowpack over high-altitude and complex terrains. The assessment of the AGCM by comparing it with regional models (such as those available in CORDEX; Gutowski et al., 2016) will be crucial in the future. Given the high impacts of ENSO and AO on TPSC, such assessment can desirably cover long duration. Additional high-resolution observations and numerical simulations of snow cover will be required to further explore the linkage between TPSC and large-scale climate variability modes.