A Detectable Anthropogenic Shift Toward Intensified Summer Hot Drought Events Over Northeastern China

This study investigates the influence of external forcings on the various summer hot drought events (SHDEs) over northeastern China (NEC). SHDEs are represented by the probability‐based index (PI), which considers precipitation and temperature anomalies. The results show that SHDEs over NEC increased from 1961 to 2005, and the experiments for historical forcing (ALL), increased greenhouse gases emission forcing (GHG), and anthropogenic forcing (ANT) can largely reproduce the spatial and temporal features of the trends of SHDEs over NEC. Based on the optimal fingerprinting method, the impact of increased anthropogenic activities can be detected at the 90% confidence level. In addition, the attributable changes of PI in response to GHG and ANT forcings resemble the observation, implying that the increasing trends of SHDEs over NEC are primarily attributed to the increased anthropogenic activity. Furthermore, the occurrence probability of SHDEs over NEC will be further increased under different Representative Concentration Pathways in the future. Additional strict control regulations on GHG emissions are thus suggested to mitigate its impact on regional climate changes.


Introduction
Extreme climate and weather (extreme precipitation, flood, heat waves, drought, etc.) have great impacts on the economy, agriculture, and human lives and are of great concern to the public and the government (Chen et al., 2012;Qiu et al., 2017;Wang et al., 2016;. Because northeastern China (NEC) is the granary of the country, droughts often lead to a serious deficiency in water supply and a severe reduction in crop production, which have further influences on the quality of human life and result in numerous economic losses. For example, the extreme summer hot drought event that occurred in 2016 over NEC led to severe yield reductions and economic losses reaching up to CNY15.61 billion (Li, Chen, Wang, Sun, & Ma, 2018). Thus, it is highly important to understand the changing characteristics of summer hot drought events (SHDEs) over NEC and its associated possible mechanisms.
In general, the SHDEs over NEC are characterized by higher temperature and less precipitation. To date, numerous studies have focused on the changing features of the summer precipitation and surface air temperature (SAT) over NEC on different timescales. In terms of the interannual variations, changes of the sea ice concentration over the Arctic Ocean and Greenland , the north Atlantic SST Wu et al., 2011;Zhao et al., 2019), the Indian subcontinent (Zhang, Chen, & Zhao 2019) and the Indian Ocean (Gao et al., 2014;Sun et al., 2019), the Tibetan heating , the Pacific Ocean (Zhu et al., 2019), and the weakened East Asian summer monsoon (Sun et al., 2017) may be linked to the variations of summer precipitation and SAT over NEC. On the interdecadal timescale, summer precipitation is observed to significantly decrease since the end of the 1970s (Ding et al., 2008;Liu et al., 2011;Sun & Wang, 2015), which might be partially linked to the abrupt decrease of preceding winter and spring snows over the Tibetan Plateau (Ding et al., 2009). In addition, the summer precipitation over NEC also experienced an interdecadal decrease after the end of the 1990s (Zhao et al., 2018). Similar to changes in precipitation, the summer SAT over NEC experienced a shift toward a warm phase in the early 1990s (Chen & Lu, 2014). Additionally,  suggested that the Eurasian warming may cause the extreme drought in NEC in the recent two decades. For the long-term trend, many studies have suggested that the summer precipitation over NEC has a decreasing trend (Qian & Lin, 2005;Zhai et al., 2005) while SAT presents an increasing trend (Feng et al., 2015) over the past half-century. Hence, previous studies have deeply investigated the features and causal factors of changes in summer precipitation and SAT over NEC on different timescales from the perspective of the internal climate system (Wang & Zeng, 2018).
In terms of the changes in drought events over NEC, Yang et al. (2016) indicated that drought (represented by the Standardized Precipitation Index) over NEC has become increasingly severe, mainly resulting from the decreased precipitation. When considering the impact of precipitation and temperature simultaneously, Zou et al. (2005) suggested that the drying trend over NEC was significant during the period from 1951 to 2003 based on the Palmer Drought Severity Index. Yu et al. (2014) also indicated that the persistent drought event over NEC became increasingly severe during the period from 1951 to 2010 based on the Standardized Precipitation Evapotranspiration Index (SPEI). In addition, Chen and Sun (2015a) suggested that the drought events have increased over NEC on the centennial scale, and the severe drought events were closely related to higher air temperature based on the results of SPEI. In general, all of these indices coincidently presented an increasing trend in summer drought events over NEC. In addition to these drought indices, another index based the concept of the multivariate copula is introduced (AghaKouchak et al., Cheng et al., 2016;Michele et al., 2005;Salvadori & De Michele, 2010). This drought index uses a probability-based index that considers the joint effect of precipitation deficiency and high SAT, which is easier to identify the severity of hot drought event and has been used to investigate the influence of the internal natural forcing (Barents Sea ice decline in March) on the variations of SHDEs over NEC on the interannual scale (Li, Chen, Wang, Sun, & Ma, 2018). Therefore, this probability-based index could be further applied to investigate the long-term changes in SHDEs over NEC as well as the potential contributions of the strengthened human activities.
Human influence has been detected in the changes in drought events over different regions around the world, such as central and eastern China (Ma, Zhou, Stone, et al., 2017), the Tibetan Plateau , southern China (Chen & Sun, 2017b), East Africa (Lott et al., 2013), and California (Diffenbaugh et al., 2015;Seager et al., 2017;Wang et al., 2014;Williams et al., 2015). However, the role of anthropogenic activity in the increased summer hot droughts over NEC remains unclear. Therefore, the focuses of this study are listed as follows: What is the changing characteristic of summer hot droughts over NEC? Whether human activities have imposed detectable impacts on the increased SHDEs over NEC? How severe the SHDEs over NEC will be in the future?
The outline of this study is listed as follows. Section 2.1 shows the data sets used in this study. Section 2.2 describes the methods, including the multivariate copula method, the ensemble empirical mode decomposition (EEMD) method, and the optimal fingerprinting analyses. Section 3 first presents the analyses regarding the performance of different external forcing simulations from the Coupled Model Intercomparison Project Phase 5 (CMIP5) models to reproduce SHDEs over NEC. After that, the detection and attribution analyses of changes in NEC in response to different external forcings are conducted. Besides, the future changes of SHDEs under different Representative Concentration Pathways (RCP) are further discussed. Section 4 provides detailed discussions concerning anthropogenic activity to influence the increased SHDEs over NEC, and Section 5 finally presents a brief conclusion.

Data
An observational gridded monthly precipitation and SAT data set (referred to as CN05.1) developed by Wu and Gao (2013) is applied in this study. CN05.1 was produced by interpolation of data from 2416 stations across China during the period from 1961 to 2016 using the "anomaly approach" method (Wu & Gao, 2013;Xu et al., 2009), which has a horizontal resolution of 1°× 1°.
To conduct optimal fingerprinting analyses, the monthly precipitation and SAT data sets from CMIP5 are also used and are obtained from the website of https://www.ipcc-data.org/sim/gcm_monthly/AR5/ Reference-Archive.html. To investigate the influences of external forcings on the variations of SHDEs over NEC, multi-model simulations from different external forcings during the period from 1961 to 2005 are used, including historical anthropogenic plus natural forcing (ALL), natural forcing (NAT), and greenhouse gases forcing (GHG). In this study, the anthropogenic impact (ANT) is calculated as ALL minus NAT (Zhang et al., 2013) and that the other anthropogenic forcing (OA) such as aerosol emissions and land use can be obtained as ANT minus GHG . For the purpose of detection and attribution analyses, the preindustrial unforced control (CTL) simulations are also selected. Here, five chunks of nonoverlapping 45-year data sets (225 years in total) of precipitation and SAT are selected for each model. In addition, the simulations under the RCP4.5 and the RCP8.5 scenarios are used to investigate future changes in SHDEs over NEC. In summary, 14 climate models that include the simulations of ALL, GHG, NAT, CTL, and future scenarios are selected. The detailed information on these models is listed in Table 1. To give equivalent weight to each model, the ensemble means for each model is obtained and then the multi-model ensemble median (MME) of the chosen models are calculated when needed . For convenience, all of the data sets derived from CMIP5 are bilinearly interpolated onto 1°× 1°grids as CN05.1, and then the seasonal mean (June, July, and August) of those data sets is calculated.

Methods 2.2.1. Multivariate Copula Method
In our study, the multivariate copula method is used to identify the SHDEs over NEC. The traditional copula method is proposed by Salvadori and De Michele (2010), in which the joint cumulative distribution function of precipitation (X 1 ) and temperature (X 2 ) can be defined as where F 1 (x 1 ) = P(X 1 ≤ x 1 ) and F 2 (x 2 ) = P(X 2 ≤ x 2 ) are the relevant cumulative distribution functions. Here, x 1 and x 2 represent the thresholds of precipitation and temperature and the value of F(x 1 , x 2 ) can be understood as the case that precipitation and temperature are higher bounded (corresponding to low precipitation and low temperature).
Further, Salvadori et al. (2013) proposed the joint survival method to better calculate the return period. Based on this method, the condition of precipitation deficiency and high temperature can be considered simultaneously. This joint survival function is defined as where In the case of low precipitation and high temperature, X 1 represents the original precipitation time series multiplied by −1, and x 1 represents the threshold of precipitation  multiplied by −1 correspondingly. Here, the cumulative distribution functions of summer precipitation (which has been multiplied by −1) and SAT over NEC are calculated based on the kernel smoothing density estimator. Then, the corresponding PI can be calculated using the t copula method, which is used to fit the cumulative probability of temperature and precipitation (https://ww2.mathworks.cn/help/stats/copularnd.html). This t copula method considers the joint survival function of precipitation ( F1 X 1 ð ÞÞand temperature ( F2 X 2 ð Þ) and hence can identify the hot drought events. Accordingly, the values of PI vary between zero and unity, representing the severity of hot drought events. Relatively high (low) thresholds for precipitation (x 1 ) and temperature (x 2 ) would result in a relatively small (large) value of PI, which represents a more (less or no) severe hot drought event. Thus, a period (e.g., 1995-2005) with increased small PI events denotes a period with increased hot drought events (supporting information, Figure S1).
Based on the PI, the corresponding survival Kendall's return period (κx ) can be obtained according to Salvadori et al. (2013), where a large return period suggests a severe small-probability hot drought event. Given a certain pair of X 1 and X 2 , there exists a unique isoline that shares the same probability t (t ∈ PI, t ∈ [0,1]; see the isolines in Figure 1). The corresponding survival Kendall's return period (κx) is defined as where t is a certain probability level of PI, and K̄t ð Þ is obtained by Therefore, the probability level of PI is negatively correlated with κx, and a small value of PI corresponds to a severe hot drought event (longer return period, see the pentagram in Figure 1). Figure 1 shows an illustration of the concurrent return periods of precipitation deficiency and high temperature over NEC in JJA during 1961JJA during -2005 Here, the blue dots give the anomalous precipitation and temperature for each year, and the isolines represent the corresponding survival Kendall's return period. For example, the hot drought event occurred in 2000 is the severest year (with the highest temperature level and the second lowest precipitation level) during the period of 1961-2005, which has the lowest PI value (0.015; Figure S1) with a corresponding survival Kendall's return period of 60 years ( Figure 1). The severity of hot drought events is easily identified based on PI. Thus, the PI is used to analyze the hot drought event in the following sections.

Ensemble Empirical Mode Decomposition Method
In this study, the EEMD method (Wu & Huang, 2009) is employed to analyze the changes in PI at different timescales. Based on this method, the components with truly physical meaning can be extracted from the signal. The steps to obtain different components of c j (t) is summarized as follows: (1) Add a random white noise signal to the original series; (2) decompose the series into several intrinsic mode functions (IMFs); (3) repeat steps (1) and (2) with different white noise series; and (4) obtain the IMF components (Wu and Wang, 2009). According to equation (5), the MME of the PI series can be divided into six IMF components from a high frequency to a low frequency of c j (t) in this study where r n (t) is the residual term and n denotes the total number of IMF components (n = 6). Here, the IMFs are oscillatory functions with varying amplitude and frequency, and the number of extrema and zerocrossings must either be equal or differ at most by one for a single IMF. Based on the EEMD method, the variability part of PI is obtained by summing the first five components (e.g., high-frequency signals, annual signals, and interdecadal signals), while the last component is recognized as the trend part.

Optimal Fingerprinting Method
To detect and attribute changes in hot drought events over NEC (represented by PI), the optimal fingerprinting method is introduced (Zhang et al., 2013). The optimal fingerprinting method is based on the concept of the total least square method (Min et al., 2011;Ribes et al., 2013). In equation (6), y represents the observational historical PI values, x is the climate response to the external forcing considered (the PI values simulated by the MME of a certain external forcing including ALL, GHG, NAT, ANT, and OA), β is the scaling factor that adjusts the magnitudes of the fingerprints to best match the observation, and ε is the Gaussian random residual term relevant to the internal variability. To conduct the fingerprinting analysis, the CTL simulations with 30 segments of 45-year nonoverlapping PI data sets are divided into two samples (six models with five chunks, respectively, which are chosen based on the criteria in Section 3.1). Half of these segments are used to prewhitening the data, and the other half segments are used to calculate the 5-95% uncertainty of regression coefficient relating to β by regressing the observational PI onto model simulated PI series (x) based on the MME of ALL, GHG, NAT, ANT, and OA. More detailed information can be found in Allen and Stott (2003).
To investigate the attributable changes in PI, trends for the MME of PI are first multiplied by the corresponding scaling factors (5-95% marginal of β), and then they are further multiplied by 45 (corresponding to the 45-year period) for each individual external forcings. The observed changes in PI are calculated by multiplying the trend of the observed PI (the 90% uncertainty range of trend is estimated by the total least square method) by 45 (corresponding to the 45-year period).

Performance Evaluations in the Simulation of SHDEs
Before detecting and attributing historical changes in SHDEs over NEC, the model capability in reproducing trends of summer precipitation and SAT over NEC is evaluated. Figure 2 shows the observational trends of  It seems that there exists large diversity in the long-term trend of the SAT and precipitation over NEC among different models, and these biases might be influenced by the internal climate variability (Santiago et al., 2016), the internal model variability (Giorgi & Bi, 2000), the choices of the model schemes , and so on. Thus, six models (bcc-csm1-1, CNRM-CM5, CSIRO-Mk3-6-0, GFDL-CM3, GFDL-ESM 2M, and HadGEM2-ES) that simultaneously well reproduce the trend of SAT and precipitation over NEC (see the bold lettering models in Table 1) are chosen for further analyses. Figure 4 shows the observation and the historical simulations for the trends of PI from 1961 to 2005 over NEC. SHDEs over NEC tended to increase significantly (corresponding to decreasing PI). The increasing trend of SHDEs over NEC is significant at the 90% confidence level based on the Mann-Kendall nonparametric test, largely affected by the increasing trend of surface temperature (Figure 2a) rather than precipitation (the trend of summer precipitation over most parts of NEC is statistically insignificant, which cannot pass the 90% confidence level). Given that SHDEs over NEC often cause great economic and agricultural losses, the features of SHDEs over NEC and the possible influences from different external forcings are further examined in the following analysis. Overall, the decreasing trend of PI over NEC (corresponding to an increase in SHDEs) can be captured by most of the models (significant at the 90% confidence level based on the Mann-Kendall nonparametric test). The spatial pattern correlations between the observation and the individual ALL experiment are larger than 0.83. In addition, the spatial correlation for the trend of PI between MME and the observation reaches 0.95, which is higher than most of the models. The numerical values in each panel show the spatial correlation between the observation and models. The dotted regions for observation and individual models suggest that the trend is significant at 90% confidence level based on the Mann-Kendall nonparametric test, and the dotted regions in MME indicate that more than 70% of the six individual models have the same sign as the MME.

10.1029/2019EA000836
Earth and Space Science LI ET AL. that the trend is significant at 90% confidence level based on the Mann-Kendall nonparametric test, and the dotted regions in MME indicate that more than 70% of the six individual models share the same sign as the MME.

Earth and Space Science
Similar to the ALL forcing, the results based on the ANT forcing simulate well the decreasing trend of PI over NEC (Figures 5a-5g). In particular, the spatial correlations between observation and individual models are all greater than 0.65. Additionally, the spatial correlation between MME and observation reaches 0.93, largely resembling the observation. Consistent with historical simulation and the ANT forcing, the GHG forcing also presents a good resemblance to observation in terms of the trends of PI over NEC (Figures 5h-5n). The spatial correlation between individual models and observation are greater than 0.73, with a spatial correlation between MME and observation reaches 0.95. Unlike the ALL, GHG, and ANT forcings, the results based on the NAT ( Figure S2) and OA forcings ( Figure S3) show less similarity to the observation. For the NAT forcing, most of the individual models and the corresponding MME present an increasing trend of PI, which is opposite to the observation. For the OA forcing, almost none of the individual models and MME show a spatial pattern similar to the observation, except for CSIRO-Mk3-6-0, which has a decreasing trend of PI with a spatial correlation of 0.89 with observation.
In addition to the spatial patterns of PI trends during the period from 1961 to 2005, we also investigate the spatially averaged temporal variations of PI over NEC using the cosine of latitudes as weight. Based on the EEMD method, the original temporal series (Figure 6a) can be divided into two components, including the trend part ( Figure 6b) and the variability part (Figure 6c). The results suggest that the PI over NEC has interannual and interdecadal variations (Figures 6a and 6c) and a decreasing trend (Figure 6b) during the period from 1961 to 2005. The decreasing trend of PI can be satisfactorily captured by the MMEs from simulations of the ALL, GHG, and ANT forcings. However, increasing trends of PI occur for the MME of the NAT and OA forcings, which are opposite to the observation (Figure 6b). These results are consistent with the above findings concerning the spatial distribution of trends in PI. In terms of the interannual and interdecadal variations of SHDEs over NEC, all of the simulations fail to reproduce these observed changing characteristics (Figure 6c). Thus, the possible impacts of these external forcings on the long-term changes of SHDEs over NEC are further explored in the following analyses.
The above evidences indicate that the region of NEC has experienced increased SHDEs during the period from 1961 to 2005. The historical simulation and external forcings such as GHG and ANT can reproduce the change features of SHDEs over NEC, whereas the other external forcings, including NAT and OA, fail to reproduce these features. In the following section, the optimal fingerprinting method is applied to explore the extent to which these external forcings can exert influences on the changes of SHDEs over NEC.

Detection and Attribution of SHDEs Over NEC
In the following, the optimal fingerprinting method is used to investigate the influence of the given external forcings on changes in SHDEs over NEC. According to Section 2.2.3, the observational PI series averaged over NEC are regressed onto the MME of the six chosen models with respect to individual external forcings (ALL, GHG, ANT, NAT, and OA). Figure 7a shows the best estimations and the 90% uncertainty ranges of scaling factors for the trend component of PI series in response to different external forcings. For detection, the 90% (5-95%) uncertainty range of the scaling factors should exclude zero, and its best estimation should be close to unity (Allen & Stott, 2003). The 90% uncertainty ranges of the ALL, GHG, and ANT forcings are larger than zero, suggesting that these signals are detectable according to the above principles. The p values for them are larger than 0.1 based on the residual consistency test (significant at 90% confidence level), suggesting a good fit of these regression models compared with the observation (Ribes et al., 2013). Best estimates for the ALL and ANT forcings are close to unity, indicating that they have good resemblances to the observed PI trends. However, the signals of the NAT and OA forcings fail to be detected, and their negative scaling factors suggest the opposite effect on the changes in SHDEs over NEC, which are consistent with the above analyses in Section 3.1. Both the NAT forcing and the OA forcing lack a primary component to be a predictor; thus, these two regression models are invalid Zhang et al., 2013).
Figure 7b further shows the fingerprints of PI variations in response to external forcings. Results suggest that none of these external forcings can be detected according to the above principles. Similar results can also be found for the original time series ( Figure S4). Consequently, the influence of human activities (GHG and ANT) and historical simulation are detectable in terms of the increasing trend of SHDEs over NEC. However, the interannual and interdecadal variations of SHDEs over NEC are not detected by external forcings, which indeed might be linked to internal variability (e.g., Li, Chen, Wang, Sun, & Ma, 2018).
To estimate the increased SHDEs attributed to different forcings, we further calculate the attributable PI changes based on the method suggested in Section 2.2.3. Figure 7c shows the results of the attributable changes in the trend part of PI. In this study, PI is observed to decrease by −0.21 (−0.33 to −0.09), indicating that the SHDEs over NEC increased over the past half-century. Similar changes can also be found for the GHG, ANT, and ALL simulations. PI is estimated to decrease by −0.16 (−0.31 to −0.01) for the GHG forcing, by −0.15 (−0.30 to −0.01) for the ANT forcing, and by −0.17 (−0.34 to −0.01) for the ALL forcing. These attributable changes in PI are consistent with the observation. However, the median attributable changes in PI for the NAT and OA forcings are quite small, with values of −0.04 and −0.08, respectively. In addition, the 90% confidence intervals for the NAT and OA forcings are not consistent with the observation. Therefore, the effects of the NAT and OA forcings with respect to influence on changes in PI trend are rather small, but the attributable changes in PI trend with the inclusion of human activity (ALL, GHG, and ANT) well resemble the observation, suggesting that the human influence could noticeably increase the occurrences of SHDEs over NEC.

Future Changes in Summer SHDEs Over NEC
The above analyses suggest that anthropogenic activity may increase the occurrences of SHDEs over NEC. Besides, previous studies also indicated that there might be more extreme events  and more drought events across China in the future (Chen & Sun, 2017a;Dai & Zhao, 2017;Li, Chen, Wang, & Yu, 2018), which have great impacts on the security of socioeconomic systems . The issue of how the severity and the occurrence probability of SHDEs over NEC change due to the continuous GHG emissions in the future is thus investigated under different emission scenarios of RCP4.5 and RCP8.5. Figure 8 shows the temporal series of the spatially averaged PI series over NEC during the period from 1961 to 2100. Here, we combine the historical data sets (from 1961 to 2005) and the future data sets (from 2006 to 2100) under RCP4.5 and RCP8.5 scenarios together. The results based on the MME of six chosen models show an evident decreasing trend of PI in the future, suggesting an increasing occurrence of SHDEs over NEC (corresponding to the lower value of PI) in the future. The occurrence probability for SHDEs over NEC mainly presents a higher value in the RCP8.5 scenario than that in the RCP4.5 scenario. Specifically, the severe SHDEs (with PI less than 0.1, see the dashed line in Figure 8) are predicted to occur once every 6-7 years by the end of the 21st century (2080-2099) based on the MME of the six chosen models in the RCP4.5 scenario, while it will constantly occur (all of the PI values will be less than 0.1) in the RCP8.5 scenario by the end of the 21st century (2080)(2081)(2082)(2083)(2084)(2085)(2086)(2087)(2088)(2089)(2090)(2091)(2092)(2093)(2094)(2095)(2096)(2097)(2098)(2099). Further analyses suggest that the historical 20-year return period events during 1961-2005 (with 45-year period) will occur once every seven years in the RCP4.5 scenario during 2055-2099 (with 45-year period), and it will occur once every 3 years in the RCP8.5 scenario during 2055-2099 (with 45-year period) calculated by equations (3) and (4).

Discussion
To verify our results further, another popular drought index known as SPEI (Vicente-Serrano et al., 2010) is also used in this study. SPEI shows a good capability to monitor drought events over China and is also a type of probability-based index (Chen & Sun, 2015b). Based on the observational temporal series of SPEI, a decreasing trend of summer drought events over NEC ( Figure S5b) is evident, accompanied by interdecadal and interannual variability (Figures S5a and S5c). The decreasing trend of SPEI (corresponding to an increasing trend of drought events over NEC) can be captured by the ALL, GHG, and ANT forcings, while no significant trends are shown for NAT and OA ( Figure S5b). These results are similar to those based on PI ( Figure 6). Based on the optimal fingerprinting method, we observed that the ALL, GHG, and ANT forcings can be detected in terms of the trend of SPEI ( Figure S6a), with all of the p values greater than 0.1 based on the residual consistency test. However, the NAT and OA forcings still are not detected in terms of the series of SPEI ( Figure S6a). Besides, the variability part fails to be detected ( Figure S6b) in response to external forcings, for that all of the corresponding p values are smaller than 0.1 based on the residual consistency test. These findings are also consistent with the results of the PI (Figure 7). In terms of the attributable changes in SPEI ( Figure S7), the results also suggest that the ALL, GHG, and ANT forcings satisfactorily reproduce the observational changes in SPEI, whereas the other two forcings fail. Thus, the increased SHDEs over NEC are a robust response to the strengthened anthropogenic activities.

Conclusion
In summary, the SHDEs over NEC increased significantly during the past century, and the historical simulation and the external forcings, including GHG, ANT, and ALL, can reproduce this feature from both spatial distribution and temporal variations. Based on the optimal fingerprinting method, the ALL, ANT, and GHG forcings can be detected. Additionally, the attributable changes in PI trend for the ALL, ANT, and GHG forcings well capture the observational changes in PI. On the other hand, the results suggest that the NAT and OA forcings have almost no contribution to the observational changes in PI trend. Therefore, the anthropogenic activity may be an important factor that influences the increased SHDEs over NEC where greenhouse Figure 8. Temporal series of PI averaged over NEC from 1961 to 2100 during the historical period (black) and during the future period based on RCP4.5 (blue) and RCP8.5 (red). The colored solid lines indicate the MME of the six models, and the corresponding shadings indicate the model spreads. The black dashed line represents the value of 0.1. gas emissions may be the dominant factor for the detected response. Furthermore, the region of NEC is expected to face increased and more severe SHDEs in the future.