Anthropogenic warming of Tibetan Plateau and constrained future projection

Serving as ‘the water tower of Asia’, the Tibetan Plateau (TP) supplies water resources to more than 1.4 billion people. It is warming more rapidly than the global average over the past decades, affecting regional hydrological cycle and ecosystem services. However, the anthropogenic (ANT) influence remains unknown. Here we assessed the human contribution to the observed TP warming based on coupled climate simulations and an optimal fingerprinting detection and attribution analysis. We show that the observed rapid warming on the TP (1.23 °C over 1961–2005) is attributable to human influence, and particularly, to the greenhouse gases with a contribution of 1.37 °C by the best estimate, which was slightly offset by anthropogenic aerosols. As the multi-model ensemble tends to underestimate the ANT warming trend, the constraint from the attribution results suggests an even warmer future on the TP than previously expected, implying further increased geohazard risks in the Asian water tower.


Introduction
The Tibetan Plateau (TP), known as 'the roof of the world' or 'the 3rd pole of the earth' , contains the largest volumes of ice outside the Arctic and Antarctic. Feeding water to many Asian rivers, including the Yangtze River, the Yellow River, and the Ganges River, it has supplied water resources to more than 1.4 billion people, serving as 'the water tower of Asia' (Xu et al 2008, Immerzeel et al 2010, Immerzeel and Bierkens 2012, Biemans et al 2019. The TP also plays an important role in the regional to global climate and hydrological cycle, via interactions with atmospheric circulations such as the Asian monsoon system (Zhang et al 2015(Zhang et al , 2017 and midlatitude westerlies (Schiemann et al 2008, Wang et al 2017.
The TP is undergoing rapid climate change due to its high sensitivity to global warming. Observations have revealed a significant warming trend on the TP since the 1950s with a range from 0.16 • C to 0.67 • C per decade during different periods, which is faster than that over the northern hemisphere and the northern midlatitudes (e.g. Liu and Chen 2000, Rangwala et al 2009, You et al 2010, Duan and Xiao 2015, Kuang and Jiao 2016, Yao et al 2019. The warming trend on the TP exhibited strong seasonal and spatial variations. The greatest warming is observed in winter (e.g. Liu and Chen 2000, You et al 2010, Cai et al 2017. Spatially, larger warming rates are found in the northern part of the TP (e.g. Yang et al 2011, You et al 2015. In addition, a distinct feature is the elevation-dependent warming over and around the TP, as manifested by the amplified warming rates with elevation mainly since the 1990s (Pepin et al 2015, Cai et al 2017, Yao et al 2019. The observed warming of the TP is directly associated with the altered surface energy balance mainly related to the snow-albedo feedback and cloud-radiation interactions. On one hand, the absorbed solar radiation has increased due to the decreased snow cover through the snow-albedo feedback (Rangwala et al 2010), and to the decreased daytime cloud over the southern TP (Duan and Wu 2006, Duan andXiao 2015). Meanwhile, the increased downward longwave radiation at the surface has also contributed to the warming as a result of the increased surface water vapor especially during cold months (Rangwala et al 2009), and the increased nocturnal cloud over the northern TP which enhances atmospheric counter-radiation (Duan and Wu 2006, Duan andXiao 2015).
The multidecadal temperature variations on the TP are affected by both internal climate variability and external forcings. In terms of internal variability, the Atlantic Multidecadal Oscillation (AMO) has been demonstrated to contribute to temperature variations on the TP on the multidecadal timescale, as suggested by proxy records dating back to the last millennium, instrumental records, and model simulations (Wang et al 2009, Shi et al 2017. Specifically, the warm (cold) phases of the AMO are associated with warm (cold) periods of the TP, as well as over East Asia (Wang et al 2009, Shi et al 2017. This teleconnection between the AMO and temperature variations over East Asia is possibly associated with the mid-latitude westerly anomalies and the propagation of Rossby waves related to the AMO (Wang et al 2009(Wang et al , 2013. Other internal variability modes such as the Pacific Decadal Oscillation (PDO) exert weaker influences on the temperature over the TP (Wang et al 2013(Wang et al , 2014. From the perspective of external forcings, anthropogenic (ANT) influence, particularly the greenhouse gas (GHG) forcing, has been demonstrated as the dominant driver of the observed warming trends on the global to regional scales based on formal detection and attribution analyses with high confidence (e.g. Bindoff et al 2013). Nevertheless, this effect has not been established rigorously for the warming over the TP. A recent study has reported the detectable contribution of ANT forcing in extreme temperature changes, including the increasing warm extremes and decreasing cold extremes, on the eastern TP since 1958 (Yin et al 2019). Nevertheless, the influence of different components of ANT forcing (e.g. GHG and aerosols;Duan et al 2006), as well as the natural (NAT) forcing, in the TP warming remains ambiguous. Thus, in this study we aim to assess the quantitative contributions of different external forcings in the observed warming trend on the TP using a formal optimal fingerprinting detection and attribution analysis. Specifically, we aim to address the following scientific questions: (a) Whether and to what extent can the observed TP warming be attributable to human activities? (b) How does the attribution of the past warming imply for the future of the TP?

Observations
We used a gridded monthly surface air temperature (SAT) dataset from the National Meteorological Information Center of China developed from over 2400 observing stations, which covers 1961 to the present (0.5 • × 0.5 • ; Xu et al 2009, Wu andGao 2013). There are ∼100 observing stations on the TP, most of which distributed in the eastern TP (figure 1 in Wu and Gao 2013).
The TP is outlined with an elevation of 2500 m. Detection and attribution analyses were performed for the central-eastern TP (28 • -40 • N, 85 • -105 • E; blue box in figure 2(a)) where the observations are relatively sufficient and reliable.

Model data
We used monthly near-surface SAT data from 14 CMIP5 models (Taylor et al 2012). These models were selected as they provide historical simulations with ANT forcing alone, GHG forcing alone, or anthropogenic aerosol (AA) forcing alone. A total number of 217 historical simulations were used, including 71 simulations with all external (ALL) forcings, 42 simulations with natural (NAT)-only (solar and volcanic activities) forcing, 29 simulations with ANT forcing, 42 simulations with GHG forcing, and 33 simulations with AA forcing. In addition, a total of 8438 years of pre-industrial control simulations from these models were used to estimate internal variability (supplementary table 1 (available online at stacks.iop.org/ERL/16/044039/ mmedia)). As the majority of the separate forcing simulations cover 1850-2005, detection and attribution analyses were performed for 1961-2005. Furthermore, SAT projections in the 21st century under the RCP4.5 and RCP8.5 emission scenarios from 32 CMIP5 models were used to derive observationally constrained projections (supplementary table 2).
The CMIP5 ensemble reasonably reproduces the seasonal cycle and spatial distribution of near-surface SAT in climatology on the TP, with relatively cold temperatures on the western TP where the elevation is the highest (supplementary figure 1). There exists a long-standing cold bias on the TP in model simulations, which is related to biases in surface energy processes (Chen et al 2017).

Optimal fingerprinting detection and attribution
We employed the optimal fingerprinting method based on a generalized multivariate linear regression model Tett 1999, Allen andStott 2003). It is expressed as y = β (X − υ) + µ, where the observed variations (y) is the sum of scaled forced responses or fingerprints (X) estimated from multimodel ensembles and the noise in observations (µ). The sampling noise (υ) arises from a finite ensemble simulations. The total least squares regression was used here. The scaling factor β adjusts the signal magnitude so that the modeled response best matches the observations. A scaling factor with its confidence interval significantly larger than zero indicates the detectable effect of a certain external forcing; additionally, a scaling factor consistent with unity indicates that the modeled response is consistent with the observed changes.
To improve the signal-to-noise ratio, the regression was conducted for the regional mean SAT and for 5 year non-overlapping means over 1961-2005. To enhance robustness of the results, the detection and attribution analysis was also performed with regional information included, by dividing the TP into three subregions based on different warming rates (see division in figure 2(a)).
The estimates of natural climate variability (noise covariance) were based on the control simulations as well as intra-ensemble variability in separate forcing runs. Model-simulated variability is checked against observed residual variance using a standard residual consistency test (Allen and Tett 1999). The detection results are reliable only when internal variability is sufficiently represented in models (Zhang et al 2006). To account for a potential underestimation of variability in models, we also conducted the analyses with noise estimates doubled (as in Polson et al 2013, Polson and Hegerl 2017). Please see supplementary text 1 for details.

Observationally constrained projections
To account for the model overestimation or underestimation of the responses to ANT forcings, we used the attribution result as an observational constraint to the future projections of anthropogenic warming. The observationally constrained projections were derived by multiplying the ensemble mean projections of SAT with the scaling factor of the ANT forcing from the two-signal analyses, as in previous studies (Sun et al 2014).

Results
In the observation, the TP has exhibited a warming trend of 0.27 • C/10 year over 1961-2005, which is statistically significant at the 1% level (black curve in figure 1). The observed overall warming trend is reproduced well under ALL forcings with both NAT and ANT influences, with a significant warming trend of 0.22 • C/10 year (red curve in figure 1(a)). However, the simulations with NAT forcings alone fail in reproducing the observed warming trend (grey curve in figure 1(a)). This suggests that the warming trend in ALL forcing runs is dominated by ANT forcing (orange curve in figure 1(b), with a significant trend of 0.19 • C/10 year), mostly from the GHG forcing (green curve in figure 1(b), with a significant trend of 0.30 • C/10 year), and is partly offset by the AA forcing (blue curve in figure 1(b), with a significant cooling trend of 0.11 • C/10 year). Here we assume that the effect of land use change on the TP is negligible, and that the effect of ANT forcing is dominated by GHG and AA forcings, given that human induced land use changes in the TP are relatively small compared to the other parts of the globe (Foley et al 2005, Venter et al 2016, Song et al 2018. The warming of the TP exhibited distinct spatial differences (figure 2(a)). In particular, the northern TP has warmed faster than the southern TP. A warming center is located in the Qaidam Basin year. The warming pattern in the ALL forcing runs (figure 2(b)) resembles the observations (figure 2(a)) and is mostly due to ANT forcing (figure 2(d)), with a negligible contribution from NAT forcing (figure 2(c)). In particular, the ANT warming is dominated by GHG forcing (figure 2(e)), and is partly offset by AA forcing (figure 2(f)).
To objectively test for the presence of ANT or NAT forcing included responses in the observed warming, we used an optimal fingerprinting detection and attribution analysis to quantify the relative contributions of individual external forcings (Hasselmann 1993, Allen and Tett 1999, Allen and Stott 2003, Zhang et al 2006, Stone et al 2009, Hegerl and Zwiers 2011. In the one-signal analyses, we regressed the observed SAT anomalies onto the simulated anomalies under individual external forcings. The results confirm that external forcings have played substantial roles in driving the observed warming trend over the TP (figure 3(a)). The fingerprint of ALL forcings is detected in the observations, as the scaling factor β for ALL forcing including its 90% confidence interval is significantly larger than zero. The best estimate of β is consistent with 1, and thus the model-simulated responses of the SAT to ALL forcings are consistent with the observed warming magnitude.
To understand which forcing components dominate the detected response in ALL forcings, the observed SAT anomalies over the TP were then regressed onto the individual forcing responses ( figure 3(a)). The scaling factor β for ANT is significantly larger than 0, and the corresponding residual is consistent with internal variability, indicating that the responses to ANT forcing alone are detected in the observed TP warming, although the response is slightly underestimated by the models (with β slightly larger than 1). The GHG forcing is also detected in the observed changes, with β significantly larger than 0 and close to 1 and passing the residual consistency test.
On the other hand, the AA forcing response is not detected in the observed changes, with the 90% confidence interval of β being negative (i.e. the changes are opposite that of the observations). For the NAT forcing, the 90% confidence interval of β includes zero when the modeled internal variability is doubled, and the corresponding residual is inconsistent with internal variability, suggesting that the modeled SAT response to NAT forcing cannot account for the observed changes ( figure 3(a)). Consequently, it is concluded that the detected SAT responses in ALL and ANT forcings are dominated by GHG forcing.
There may exist correlations between signals. To distinguish the influence of one forcing from the others, we further conducted two-signal detection analyses for ANT-NAT, and GHG-AA forcings (figures 3(b)-(c)). The origin (0, 0) is outside the 90% joint confidence regions for both ANT-NAT and GHG-AA, suggesting that the combined effect of ANT-NAT, and GHG-AA can be detected. Specifically, the detectable effect of ANT forcing can be separated from NAT forcing, with β larger than zero for ANT forcing ( figure 3(b)). In addition, the detectable effect of GHG forcing can be separated from AA forcing as evidenced by the β larger than zero for GHG forcing (figure 3(c)). The uncertainty ellipsoid for GHG and AA are strongly tilted, indicating that the scaling factors for these two signals are correlated, and that the warming induced by GHG may be offset by the cooling induced by AA ( figure 3(c)). The results of two-signal analyses are robust with different numbers of EOFs retained (supplementary figure 2). A three-signal analysis further confirmed the separately detectable effect of GHG from those of NAT and AA, indicating a robust attribution to GHG forcing (supplementary figure 3).
To take into account the regional features, we further divided the TP into three subregions, viz. the north-eastern TP, the central-southern TP, and the south-eastern TP according to different warming rates (see blue boxes in figure 2(a)). The observed warming over the subregions are represented well by ALL forcing simulations (supplementary figure  4). We then re-conducted the optimal fingerprinting detection and attribution analyses by taking into account the spatial patterns. The one-signal analyses indicate that the ALL, ANT, and GHG fingerprints are detected over a wide range of EOF truncations, while the signals of NAT and AA are not detected (supplementary figure 5). The two-signal and threesignal analyses further indicate that the ANT forcing is separately detectable from NAT forcing, and that the GHG forcing is detectable from the AA and NAT forcings (supplementary figures 6 and 7). Based on the rigorous detection and attribution analyses, we conclude that ANT forcing, dominated by the GHG forcing, has had an attributable influence in the observed warming over the TP during 1961TP during -2005 The influence of ANT forcing, particularly the GHG forcing, on the observed warming can be consistently detected in all the seasons, whereas the AA forcing has led to a slight cooling trend all year round (figure 4). Quantitatively, among the seasons, the warming trend in winter (0.40 • C/10 year) is approximately twice as large as that in summer (0.20 • C/10 year) in the observation (figure 4). Such a seasonality is reproduced by the ANT and GHG forcing simulations, although the seasonal difference is underestimated. For example, the warming trend in response to the GHG forcing is 0.32 • C/10 year in winter and 0.28 • C/10 year in summer. In particular, this is due to the underestimated warming trend in winter in response to the ANT and GHG forcings in models, as indicated by the corresponding scaling factors which are generally larger than 1 ( figure 4(l)).
To what extent can the warming be attributed to the forcing components? We further estimated the warming attributable to different individual external forcings, by multiplying the warming trends in the ALL, NAT, ANT, GHG, and AA signals by their respective estimates of scaling factors ( figure 3(d)). confidence intervals in one-signal analysis at a truncation of seven EOFs, for simulations with all external forcings (ALL), natural (NAT) forcings alone, anthropogenic (ANT) forcings alone, greenhouse gases (GHGs) alone, and anthropogenic aerosol (AA) alone. The thin error bars indicate 90% confidence intervals when model-simulated internal variability is doubled. The residual consistency test is passed in all cases except for NAT (indicated with the cross symbol). (b), (c) Scaling factors, their 90% confidence intervals (pink bars for x-axis variables and blue bars for y-axis variables), and 90% joint confidence regions (dashed ellipsoids) obtained from two-signal analyses at a truncation of seven EOFs. The points (1, 1) denoted by stars indicate consistency between model-simulated responses and the observed changes. (d) Attributable warming (unit: K) due to ALL forcings from the one-signal analysis, due to NAT and ANT forcings from the two-signal analysis, and due to GHG and AA forcings from the three-signal analysis, along with their 90% confidence intervals. The solid line indicates the observed warming magnitude over 1961-2005. Of the observed annual-mean warming of 1.23 • C during 1961-2005, 0.93 • C (90% confidence interval 0.51 • C-1.35 • C) can be attributed to ALL forcing. Of the three individual components of ALL forcings, the GHG would have warmed the TP by 1.37 • C (0.12 • C-2.79 • C). This attributable warming is offset by a slight cooling due to AA, although its effect is not detected. The contribution of NAT forcing is negligible.
We note that the models tend to underestimate the annual-mean SAT responses to the ANT forcing, as the best estimates of scaling factor for ANT are larger than 1 in both the one-signal and two-signal analyses (figures 3(a) and (b)). This mainly results from the underestimed warming trend in autumn and winter (figures 4(g)-(l)). This implies that the multi-model ensemble of CMIP5 may underestimate the future anthropogenic warming on the TP. Thus, we applied the attribution results as a constraint to future projections, by multiplying the ensemble mean responses of SAT under the RCP4.5 and RCP8.5 emission scenarios with the scaling factor of the ANT forcing from the two-signal analyses. This constraint assumes that the scaling factor based on historical observations and historical ANT forcings remains reasonable for the future (Sun et al 2014), although  aerosols are expected to reduce in the future. The observationally constrained projections suggest an even warmer future on the TP than previously expected (figure 5). Under the RCP8.5 (RCP4.5) scenario, in terms of the annual mean, the TP would warm by 5.95 • C (2.99 • C) by the end of the 21st century (2081-2100) under the best estimate of constraint, which is 0.63 • C (0.32 • C) warmer than the raw output. Even for the mid-term (2041-2060) projection, the TP would warm by 3.00 • C (2.25 • C) under the best estimate of constraint, which is 0.32 • C (0.24 • C) warmer than the raw estimate. Hence, the observational constraint indicates a substantially warmer TP in the mid-to long-term futures than previously expected.

Concluding remarks
Based on the rigorous optimal fingerprinting analyses, the observed warming of the TP (1.23 • C during 1961-2005) is detected and attributed to human influence, dominated by the GHG forcing. Specifically, the GHG forcing has warmed the TP by ∼1.37 • C by the best estimate, which was slightly offset by the AAs. The attributable influence of ANT forcing, particularly the GHG forcing, on the observed warming holds across all the seasons.
The attribution results suggest that the current generation of the state-of-the-art GCMs tend to underestimate the SAT responses over the TP to ANT forcing. Using the attribution results as an observational constraint to future projections, an exacerbated warming is expected over the TP. In particular, by the end of the 21st century (2081-2100), the constrained projections suggest a 0.63 • C (0.32 • C) greater warming of the TP than previously expected under RCP8.5 (RCP4.5) scenario. This implies a greater loss of glacier mass (Yao et al 2012, Lutz et al 2014, Huss and Hock 2018, Best 2019) and further increased geohazard risks like landslides, debris flows, glacial lake outbursts, and snow and ice disasters in this region (Zhang et al 2017, Gao et al 2019, Immerzeel et al 2020. Thus, timely and effective mitigation and adaptation activities are in urgent need for the Asian water tower, on which billions of peoplerely. The detectable and attributable influence of ANT forcing, particularly the GHG forcing, on the global scale warming over the past century has been widely recognized (e. In comparison, the GHGs induced warming is much higher on the TP, by approximately 0.30 • C/10 year, as shown in our results (over a similar period of 1961-2005, with the observed warming of 0.27 • C/10 year). Thus, the greater observed warming on the TP than the global average is largely induced by GHGs.
We note that, in addition to external forcings, regional climate change can also be affected by internal climate variability. Our results show that external forcings have largely accounted for the observed TP warming in recent decades, suggesting a weak influence from internal variability. We acknowledge that multidecadal internal variability (such as the AMO) may have played a role in the observed TP warming. Moreover, the interactions between external forcings and internal climate variability may also affect regional climate (Lu et al 2014). Further dedicated research is needed to disentangle and quantify the contributions of external forcings and internal variability, as well as their interactions, in the observed TP warming, for example, using pacemaker experiments in which observed internal variability modes (e.g. AMO and PDO) are prescribed, such as the GMMIP experiments in the Coupled Model Intercomparison Project Phase 6 (Zhou et al 2016).

Data availability statement
The CMIP5 model data that support the findings of this study are openly available at the following URL: https://esgf-node.llnl.gov/search/cmip5/. The observational SAT data is available from http:// data.cma.cn/en.