Global Sensitivity Analysis of the Standardized Precipitation Evapotranspiration Index at Different Time Scales in Jilin Province, China

: The Standardized Precipitation Evapotranspiration Index (SPEI) has been widely applied, due to its multi-scalar features and the ability to identify different drought types. However, its sensitivity to climatic variables still remains unclear, especially at different time scales. Therefore, this study investigates the sensitivity of SPEI to average temperature ( T mean ), diurnal temperature ranges ( T delta ), relative humidity ( RH ), solar radiation ( Rs ), wind speed ( U 2 ), geothermal flux ( G ), and precipitation ( P ) from 1957 to 2017 using the extended Fourier Amplitude Sensitivity Test at different time scales in Jilin Province, China. Jilin Province experienced a significant rise in T mean , and a sharp decrease in T delta , Rs , and U 2 . P is undoubtedly the most influential factor to the SPEI among the meteorological variables, which explained 59.9% − 97.9% of the total variability, especially during the main crop growing season (from May to September). While T mean , RH , or U 2 observably affect the SPEI and cannot be neglected during the nongrowing season. In terms of spatial distribution, the SPEI was mainly affected by P in the eastern region, while it was also influenced by T mean , RH , and U 2 as well in the western region. The sensitivity of the SPEI differs in time scales: P > T mean > RH > U 2 > Rs > G > T delta (1 to 6 month), P > U 2 > RH ≈ T mean > G > Rs > T delta (7 to 18 month), and P > U 2 > G > T mean > RH > Rs > T delta (more than 24 month time scale), respectively. The results have the potential to provide a reference for agricultural production and management in Jilin Province, China.


Introduction
Increasing global temperatures are associated with a rise in the frequency and intensity of extreme weather events, including severe drought and high temperatures [1], particularly in Eurasia [2][3][4] and Africa [5,6]. This has resulted in regional drought research becoming more intricate. In order to accurately monitor drought, scholars and government departments around the world have proposed a series of drought indicators, including the Chinese Z Index (CZI), the Palmer Drought Severity Index (PDSI), the Water Stress Index (WSI) and the Standardized Precipitation Index (SPI). Among them, the widely used SPI is simple to calculate, flexible, and effective. However, indices such as the SPI and CZI fail to account for evapotranspiration. Moreover, recent studies have demonstrated that the increased evapotranspiration caused by rising temperature cannot be ignored [7,8]. Thus, the use of more comprehensive indicators is becoming more popular. In particular, the SPEI, put forward by Vicente-Serrano et al. (2010), is not only highly sensitive to changes in evaporation, but is also characterized by a simple calculation and its ability to consider multiple time scales. Due to a broad range of meteorological variables in the context of global warming, the SPEI is able to comprehensively identify different drought types [7,[9][10][11][12]. Thus, the SPEI has proved to be highly effective in the detection, monitoring, and exploration of drought [7].
Numerous studies have applied the SPEI in order to evaluate drought characteristics and atmospheric drought mechanisms, and to perform drought assessments and monitoring [9,[13][14][15][16][17][18]. For example, Chen et al. [8] used the SPI and SPEI to analyze the characteristics of drought occurrence in Liaoning Province from 1960 to 2015. They determined that the SPEI was able to better explain yield variation in major cereal crops (39%-78%) due to its accounting of P and reference evapotranspiration (ETo). Peña-Gallardo et al. [19] quantified drought using the SPEI and observed that inter-annual changes in crop yields were strongly correlated to drought time scales. Over the past 50 years, empirical models have been developed to estimate ETo based on different combinations of climatic variables, an important parameter also used for the SPEI calculation, such as the Hargreaves [20], the Priestley-Taylor [21], and the FAO-56 Penman-Montieth (PM) [22] methods. Among these, the FAO-56 PM method was identified as superior for the SPEI calculation [23][24][25][26], due to its ability to accurately estimate drought trends [8].
Sensitivity analysis is effective in the identification of the dominant mechanisms in model behavior [27] and any relationship between variables [28]. It is able to analyze the main role of meteorological variables in regional ETo assessments [29]. Thus, sensitivity analysis of drought indices to meteorological variables can be applied to analyze the relationship between meteorological variables and drought events [30]. In recent years, the PM equation has been widely used in the sensitivity analysis of ETo [31][32][33]. Yao et al. [34] observed distinct sensitivities of P and reference evapotranspiration (ETo) to the SPEI in Xinjiang, and drought condition was mainly determined by the ETo in the arid zone. Vicente-Serrano et al. [35] also pointed out that precipitation and ETo had the same sensitivity to SPEI. Despite the large amount of research on the SPEI [36,37], the contribution of each meteorological variable to SPEI fluctuations, as well as the sensitivity of the SPEI to meteorological variables, still remain unclear.
In order to provide a more in-depth investigation of meteorological variable sensitivities to the SPEI, this study applied the Extended Fourier Amplitude Sensitivity Test (EFAST) method [38] to analyze the SPEI global sensitivity. More specifically, Tmean, Tdelta, RH, Rs, U2, G, and P at different time scales were used as the independent variables, and the SPEI was used as the dependent variables. The objectives of this study are to a) identify the existence of a possible tendency of each of the climatic variables; b) explore the sensitivity of SPEI to Tmean, Tdelta, RH, Rs, U2, G, and P; and c) identify the response of the sensitivity of the SPEI to time scales.

Description of Study Site
Jilin Province (121°38'-131°19' E, 40°50'-46°19' N) is located in the central part of northeast China, with a total area of 187,4000 km 2 ( Figure 1). It exhibits a high altitude in the southeast and low altitude in the northwest. Jilin Province can be divided into two regions: the mountainous terrain in the east and plains in the mid-west [39]. The province has an annual Tmean within the range of 2-6 °C, and a frost-free period for approximately 100-160 d. Mean annual P are higher in the mountains and lower in the western plains, with the average annual P of the whole province being 400-600 mm, more than 80% of which is received in summer.

Meteorological Data
Meteorological data were collected from the China surface climate data daily dataset (V3.0) released by the National Meteorological Information Center of China (http://data.cma.cn/). The surface climate data daily dataset covering 28 meteorological stations ( Figure 1) includes daily precipitation, air temperature (maximum, minimum, and average), atmospheric pressure, relative humidity, wind speed at 10 m above ground, and sunshine hours in Jilin Province from 1957 to 2017.
The Inverse Distance Weighted (IDW) method was used to interpolate the daily meteorological data to account for any gaps in the dataset. In addition, the homogeneity and reliability of the meteorological data were tested and controlled by adopting the method described by Liu et al. [40].

Reference Evapotranspiration
The Penman-Monteith reference evapotranspiration method was selected to estimate seasonal ETo on the basis of monthly interval in this study: where ETo is monthly average daily reference evapotranspiration, mm d -1 ; T is monthly average daily temperature, o C; G is monthly soil heat flux, MJ m -2 d -1 ; Rn is the monthly average daily net radiation, MJ m -2 d -1 , computed from actual duration of daily sunshine hours; γ is the psychrometric constant, kPa o C -1 ; es is monthly average of the daily mean saturation vapor pressure (kPa), ea is monthly average of the daily actual vapor pressure (kPa), Δ is the slope of saturation vapor pressure curve, kPa o C -1 ; and ws is monthly average daily wind speed at height 2 m, m s -1 .
Since daily maximum and minimum air temperature are often correlated, the average daily temperature and diurnal temperature ranges are calculated by averaging the mean of the daily maximum and minimum temperature, the difference of the daily maximum and minimum temperature, respectively [41]. The monthly or other time-scale accumulative climatic variables, P and Rs, were estimated by summing the daily data, while the monthly or other timescale mean climatic variables, Tmean, Tdelta, RH, and U2, were estimated by averaging the daily data. When assuming a constant soil heat capacity of 2.1 MJ m -3 o C -1 and an appropriate soil depth of 2 m, the soil heat flux at different time periods G can be approximated as [22]: where Ti and Ti-1 is the monthly or other timescale mean temperature ( o C) at time i, i-1; m is the number of months within a specific time scale.

Trend Analysis
The Mann-Kendall method (MK) [42] was adopted for the detection of trends for the climatic variables. Prior to this, the trend-free prewhitening method was applied in order to eliminate the autocorrelation of data sequences [43]. Following the MK test, the Theil-Sen method [40,44] was used to calculate the magnitude of trends. The MK trend test and Theil-Sen estimator were both implemented by the R package RKT [45].

Standardized Precipitation Evapotranspiration Index (SPEI)
The SPEI was proposed by Vicente-Serrano et al. [7] for identifying drought episodes. The SPEI was computed using the difference between P and ETo as the water surplus/deficit indicator. The calculation procedures for the SPEI are as follows: (1) The water surplus/deficit equation where i P is monthly precipitation, EToi is potential evapotranspiration and i D is water surplus/deficit at month i.
(2) log-logistic distribution Assuming a three-parameter log-logistic distribution, the probability density function of Di was calculated by the following equation: while the cumulative probability density function was determined as follows: where α, β, and γ are the shape, scale, and position parameters, respectively.
(3) The SPEI index series was subsequently derived from the fitted values through the transform of the standardized normal distribution. The equation used for the calculation is as follows: and P is the probability of exceeding a determination value of D. In particular, for p ≤ 0.5, and p = 1-F(x), while for p > 0.5, P =1-P and the sign of the SPEI is reversed. The variables C0, C1, C2, d1, d2, and d3 are constants. A detailed description of the calculation method can be found in [46] and [7], whereby the specific values of C0 = 2.515517, C1 = 0.802853, C2 = 0.010328, d1 = 1.432788, d2 = 0.189269, and d3 = 0.001308 can be derived.

Extended Fourier Amplitude Sensitivity Test (EFAST)
The theory of global sensitivity analysis (GSA) was initially proposed by Cukier et al. [47][48][49]. All GSA variables change simultaneously within the value ranges. Hence, the GSA always identify not only the impact of each variable on the variance of the response variable, but the interaction effects of all variables. It can also overcome the shortcomings of local sensitivity analysis, such as calculating one variable at a time, analyzing the change of one variable while maintaining other variables constant, and the lack of interaction effect analysis [50,51]. Numerous GSA methods have been applied in the literature, such as the Extended Fourier Amplitude Sensitivity Analysis (EFAST) [50], the Sobol method [52], and the Morris method [53]. In particular, EFAST is one of the most powerful techniques, and has been recommended by many scholars [38,50]. The underlying method of EFAST global sensitivity analysis is to allow for the input variables of the model to oscillate at different frequencies. In addition, EFAST is able to assess the importance of each variable by evaluating the Fourier coefficients of the model outputs. The basic theory can be described as follows: Define = ( , , … , , and variable X as the k-dimensional space of input variable, Here, s is the scale variable, ωi is the frequency variable associated with variable xi and Gi is the transform function. Assume that φi is a phase-shift chosen uniformly within the range [0, 2π], which allows for a change in the starting point of the scheme. Then we have Based on Lu and Mohanty [54], the variable function takes the form, where 1 i F − is the inverse cumulative distribution function of the probability distribution function of Xi. Thus the values of xi are computed as follows: where ∀ is the partial variance generated by the interaction of all orders between X and other , and 0 < P < 1 . Variable oscillates periodically at a frequency of i ω , function f is the periodic function of s, with the period of s taken as less than 2π. Following this, f can be expanded into a Fourier series in the integer frequency domain as follows: where l A and l B are Fourier coefficients defined as: The spectrum of the Fourier series expansion is defined as The partial variance of Y attributed to the main effect of Xi is thus obtained by: In practice the infinite summation is not convenient. The main influence of Xi can be approximated, as described by Mesgouez et al. [50], as follows: where M is denoted as the interference factor and is generally chosen within the range of [4,6]. Xu and Gertner [55] denoted the total effect of Xi as the maximum value of i ω . In addition, ( ( |~ , the sum of the partial variances, is calculated by taking the sum of the frequency spectrum elements for all frequencies between 1 and Note that we assume that the interaction effects i ω are negligible at these frequencies. The total and partial variances of Y attributed to the total effect of are thus obtained as follows: However, improved results can be obtained by using 2 1 38]. In general, the effective simulation of partial variance is a function of the interference factor M and the frequency set ~i ω . In this study, the R package "sensitivity" available on the comprehensive r archive network (CRAN) was used [56]. Figure 2 presents the spatial distribution of the trend of meteorological variables in Jilin Province, China, over the past 60 years, including Tmean, Tdelta, RH, Rs, U2, G, and P. The Mann-Kendall test demonstrated that Tmean, Tdelta, RH, Rs, U2, and ETo each exhibited a significant change (p < 0.05) over the past 60 years. Fluctuations were observed for annual P and G, although there was no significant trend detected (p < 0.05).

Meteorological Variable Trend Analysis
The Tmean in Jilin Province showed a significant increase (p < 0.05) over the study period ( Figure  2A), which is consistent with the global warming general trend. In particular, the maximum increase in Tmean was observed at Dunhua station (0.35 °C/10 a), and the minimum at Nong'an station (0.11 °C (10a) -1 ). The average increasing rate of Tmean was 0.26 °C (10a) -1 , indicating that Tmean increased by 1.56 °C during the past 60 years in Jilin Province, China. Despite the rise in mean temperature, Tdelta declined ( Figure 2B). The Tdelta in 89% of the stations demonstrated a decreasing trend, of which 82% was significant (p < 0.05), as shown in Figure 2B. Averaging the selected 28 stations, Tdelta decreases at the rate of 0.19 °C (10a) -1 , and decreased by 1.15 °C over the past 60 years. The only three exhibiting an upward trend were at Changbai (significant increase), Liaoyuan, and Nong'an. These results indicate that the Tmean significantly increases over the study period in Jilin Province, while the increases in maximum and minimum temperature are not synchronized, showing a significant negative trend in Tdelta.
The RH in Jilin Province generally showed a downward trend, with a mean tendency rate of -0.47% (10a) -1 . Among the 28 meteorological stations, 26 stations exhibited a decreasing trend in the RH, and 69.2% of the stations significantly decreased (p < 0.05). Spatially, the maximum positive trend in RH trend (0.08% (10a) -1 ) was witnessed in Wangqing station in eastern Jilin, and the lowest negative trend was detected in Qiangguo (-0.92% (10a) -1 ) ( Figure 2C). The annual Rs and average U2 generally demonstrated a significant downward trend, with the mean tendency rate being -52.1 MJ m -2 d -1 (10a) -1 and -0.12 m s -1 (10a) -1 , respectively ( Figure 2D,E). The Rs and U2 decreased by 313 MJ m -2 d -1 (10a) -1 and 0.72 m s -1 (10a) -1 , respectively, over the past 60 years. Among the selected 28 meteorological sites, both Rs in 100% of the stations and U2 in 89.3% of the stations showed a province-wide decreasing trend, which was consistent with the findings of Wang et al. [57]. This trend was significant at 71.4% of the stations for Rs and 96.0% for U2, respectively (p < 0.05). The U2 exhibited an evident spatial variation tend, with a decreasing trend observed from the southeast to northwest. The maximum positive trend reached 0.03 m s -1 (10a) -1 at Erdao station in southeast Jilin, while the minimum negative trend (-0.30 m s -1 (10a) -1 ) was observed at Fuyu station. Changes in G and annual precipitation over the past 60 years were relatively small compared to those of annual Rs and average U2 ( Figure 2F,G). Though the G of most stations (67.9%) exhibited a downward trend and the annual precipitation of 53.6% of the stations showed an upward trend, but the changes were not significant at the p < 0.05 significant level. There was a clear spatial distributing in the trend in the ETo ( Figure 2H). In particular, the ETo observed in southeast Jilin exhibited an upward trend (e.g., 12.0 mm (10a) -1 in Erdao, 109 mm (10a) -1 in Changbai, and mm (10a) -1 in Jingyu, p < 0.05). However, northwest Jilin was dominantly characterized by a negative trend (-19.1 mm (10a) -1 in Fuyu, -8.6 mm (10a) -1 in Qianan and -5.4 mm (10a) -1 in Nong'an (p < 0.05)). This may be attributed to higher temperatures, lower P, and greater sunshine levels in northwest Jilin compared to the southeast. On average, over the past 60 years, Jilin Province experienced a decreasing trend in ETo at a tendency rate of 1.96 mm (10a) -1 , with 71.4% of stations identified as having a decreasing trend, among which 20% decline at the p < 0.05 significant level.

Monthly Scale SPEI Sensitivity Analysis
In order to further clarify the sensitivity of the SPEI, Tmean, Tdelta, RH, Rs, U2, G, and P at 1 month time scale were selected as the independent variable, and the SPEI at monthly time scales were used as the dependent variable for the EFAST method. The monthly range values of the meteorological variables selected for the model at each meteorological site are reported in Table 1 according to the historical meteorological data from 1957 to 2017 in Jilin Province, China. Moreover, the first-order sensitivity and total sensitivity indices for the meteorological variables at different months are shown in Figure 3.
The first-order sensitivity index represents the sensitivity of the SPEI to each of the selected meteorological variables (main effect), while the total sensitivity index denotes the overall contribution, including the main effect and their interaction with the other inputs [50]. When averaged over the entire Jilin Province, the first-order sensitivity indices of Tmean, Tdelta, RH, Rs, U2, G, and P are consistent with the corresponding total sensitivity indices (with a negligible difference) in different months over the past 60 years. This suggests that the main effects of the selected meteorological variables explained most of the SPEI total variability, with a relatively low contribution explained by their interactions. To the seven main effects (Tmean, Tdelta, RH, Rs, U2, G, and P), the sensitivity of the SPEI differs in different months. The main sensitivity meteorological variables (the first-order sensitivity coefficients ≥ 5%) for January were P, Tmean, RH, and U2, with the first-order sensitivity coefficients being 65.0%, 12.9%, 9.8%, and 5.6%, respectively. The impact of Tdelta, Rs, and G on the SPEI are negligible due to a relatively low amount of sensitivity coefficient (< 5%). Similarly, the main sensitivity meteorological variables for the other months were P (59.9%), Tmean (16.0%), and RH (10.8%) in February, P (60.2%), Tmean (20.1%), and RH (8.6%) in March, P (71.8%), Tmean (9.1%), RH (8.8%), and U2 (5.6%) in April, P (87.1%) in May, P (92.3%) in June, P (97.9%) in July, P (95.3%) in August, P (94.9%) in September, P (88.8%) and RH (5.0%) in October, P (80.3%), RH (8.3%), and Tmean (5.5%) in November, P (76.3%), Tmean (7.8%), and RH (6.0%) in December. P is undoubtedly the most influential one to the SPEI among the seven selected meteorological variables, which explained 59.9%-97.9% of the SPEI total variability. In particular, P was the only sensitive variable identified during the main crop growing season (from May to September) in Jilin Province, China. However, in other months, Tmean, RH, or U2 observably affect the SPEI and cannot be neglected (the first-order sensitivity coefficients ≥ 5%). Averaged over all the 12 months, the SPEI was the most affected by P (80.6%), followed by Tmean (6.6%), RH (5.3%), U2 (2.8%), Rs (1.0%), and Tdelta (0.1%).   Figure 4 and Figure 5 provide the spatial distributions of the first-order and total sensitivity indices of the SPEI to Tmean, Tdelta, RH, Rs, U2, G, and P, respectively, at the 1 month time scale over the past 60 years in Jilin Province, China. The spatial distribution of the first-order sensitivity indices is generally consistent with that of the total sensitivity indices. Therefore, the first-order sensitivity indices are only discussed.
The spatial distribution of the first-order sensitivity index to Tmean shows a significant upward trend from the southeast to northwest ( Figure 4A). The sensitivity index ranges from 1.6% to 5.2%, with the minimum 1.6% located at Linjiang station, while the SPEI in northwest Jilin was much more sensitive to Tmean (10%), with the maximum of 14.8% at Baicheng station). Similarly, the first-order sensitivity to Tdelta also exhibited an upward trend from the southeast to the northwest ( Figure 4B). However, the first-order sensitivities to Tdelta were all observed to be less than 0.5%, thus suggesting that this variable has a limited impact on the SPEI.  The spatial distribution of the first-order sensitivity to RH exhibited a similar trend to that of Tmean, increasing from the southeast to northwest (0.32%-13.73%) ( Figure 4C). More specifically, the sensitivities within the southeast region were less than 5%, while in the northwest region, values were greater than 5%, where RH is one of the main sensitivity meteorological variable affecting the SPEI. There was no evident spatial distribution in the sensitivity for Rs and G ( Figures 4D,F), with more than 96% of meteorological stations showing a relatively low amount of sensitivity coefficient (<5%), therefore having no significant effect on the SPEI. The sensitivities to U2 exhibited an upward spatial distribution trend from the southeast to northwest (0.55%-5.84%) ( Figure 4E). However, with the exception of Baicheng, Tongyu, and Daan stations in northwest Jilin, the sensitivity indices of the SPEI in the other stations do not show a relatively low response to U2. The SPEI in Jilin Province was observed to be most sensitive to P. The sensitivity of P at all stations were more than 61.2% ( Figure  4G), while those in the southeast were more than 80%. The spatial distribution of this variable presented a gradient downward trend from the southeast to the northwest (94.9%-61.2%), contrasting with the results of Tmean, RH, and U2.

SPEI Sensitivity Analysis at Different Time Scales
Dynamic changes in the first-order and total sensitivity indices of the SPEI to Tmean, Tdelta, RH, Rs, U2, G, and P, respectively, as the time scale increases, are shown in Figure 6 and Figure 7. At the 1 month time scale, the first-order sensitivity values were from high to low as follows: P (80.6%) > Tmean (6.6%) > RH (5.3%) > U2 (2.8%) > Rs (1.0%) > G (0.5%) > Tdelta (0.1%). As the time scale increased, the sensitivity of Tmean experienced a rapid decline, and was followed by a slight increase. In particular, at the 1 to 12 month scale, the sensitivity to Tmean decreased rapidly from 6.6% to 3.0%, while at the 12 to 36 month scale, the Tmean increased slowly from 3.0% to 3.8% and subsequently stabilized. The sensitivity of Tdelta was generally observed to be stable over the time scales, and within the range of 0-0.1%. The sensitivity of RH exhibited a similar trend to that of Tmean, with a decrease of 5.3% to 2.8% corresponding to the time scale increase from 1 to 12 month scale, and followed by a slight increase with further time scale increases up to 36 months. The sensitivity of Rs exhibited an initial increase with time scale (1 to 4 months), followed by a decrease (5 to 12 month scale), and finally stabilized (after 13 month scale). However, observed changes were relatively slight, ranging from 0.5% to 1.5% ( Figure 6A). Both the first-order average sensitivities to U2 and G exhibited a steady trend in variation with the time scale increases. At the time scale of 1 month, the sensitivity indices were 2.8% and 0.5%, and increased to 6.6% and 5.7%, respectively, at time scale of 36 months. The sensitivity index to P fluctuated between 78.4% and the maximum value of 87.2% (11 to 12 months) over all time scales, initially rising then declining. These values are much higher than those of other meteorological variables. The results suggest that the time scale has a significant impact on the sensitivity of the SPEI to meteorological variables in Jilin Province.
In summary, at the short time scale of 1 to 6 months, the order of the SPEI sensitivity was observed as follows: P > Tmean > RH > U2 > Rs > G > Tdelta. At the time scale of 7 to 18 months, the sequence was observed as: P > U2 > RH ≈Tmean > G > Rs > Tdelta. For more than 24 months, the sequence was as follows: P > U2 > G > Tmean > RH > Rs > Tdelta.

Discussion
Numerous in-depth studies have recently been conducted in order to symmetrically evaluate SPEI characteristics over varying temporal and spatial scales [58,59]. However, investigations regarding the sensitivity of the SPEI are still relatively unclear. Previous findings regarding the sensitivity of the SPEI have been reported by Yao et al. [34] and Serrano et al. [35], but the variables are all limited to P and ETo. This may be caused by the calculation of the SPEI giving a high frequency of 0 values [18] when the water is relatively balanced, which makes the nondimensional relative sensitivity coefficients ( SPEI x x SPEI ∂ ∂  ) incalculable due to the "1/0 problem" when using the conventional-partial-derivative-based sensitivity analysis (a local sensitivity analysis method) [60,61]. In this study, a global sensitivity analysis, EFAST model, was therefore applied. The sensitivity of SPEI to meteorological variables at different temporal and spatial scales, as well as at different time scales, was discussed. P is undoubtedly the most influential variable to the SPEI over time and space ( Figure 3, Figure  4, and Figure 5). Our results suggested that P explains 59.9%-97.9% of the SPEI total variability in time from January to December, and 61.2%-94.9% in space with a gradient downward trend from the southeast to the northwest. In particular, during the main crop growing season (from May to September), P was the only sensitive variable identified (the relative sensitivity coefficient ranges from 87.1% to 97.9%), with all the other meteorological variables contributing a negligible amount of sensitivity coefficient (<5%) (Figure 3). In Tibet, China, Li et al. [62], using a multiple regression model, also reported that precipitation has the most significant influence on the variation of the SPEI. In addition, Zhang et al. [33], using the self-calibrating Palmer Drought Severity Index, confirmed our findings, and concluded that the sensitivity of precipitation was higher than that of ETo. Our finding, the high contribution rate of P, would also be a potential evidence and explanation as to why the P-based indices like the SPI show a high level of consistency with the SPEI, and a considerable adaptability in most times and places, such as the findings of Yao et al. [34], Wang et al. [57], Chen et al. [8], and Przybylak et al. [30].
P is undoubtedly the most influential variable to the SPEI, while Tmean, RH, or U2 observably affect the SPEI and cannot be neglected (Figure 3, Figure 4, and Figure 5). Sensitivity analysis models only involved P or P and ETo may fail to consider the impact of Tmean, Tdelta, RH, Rs, U2, and G on drought in the context of global warming [12,35]. Our results using the EFAST model suggest that the SPEI was the most affected by P (80.6%), followed by Tmean (6.6%), RH (5.3%), U2 (2.8%), Rs (1.0%), and Tdelta (0.1%) for the Jilin Province as a whole. During the main crop growing season (from May to September), P was the only sensitive variable identified, but Tmean, RH, or U2 observably affect the SPEI and cannot be neglected (the sensitivity coefficients ≥ 5%) during the nongrowing season ( Figure  3). This result suggests that the SPEI may have an advantage over P-based indices such as the SPI and CZI, especially in the context of global warming, because other meteorological variables are not stationary, which made the assumptions of the SPI invalid [8]. The sensitivity of the SPEI to climatic variables differs not only in months, but also in spatial distribution [57]. Our findings further suggest that, in Jilin Province, the SPEI was mainly affected by P in the eastern region, while it was also influenced by Tmean, RH, and U2 as well in the western region (Figure 4 and Figure 5). This may be due to relatively lower precipitation in the western plains where the ET0 makes more contribution to water balance (P-ETo), and thus explains a higher amount of the SPEI total variability. Regarding the order of sensitivities, Wang et al. [32] analyzed the impact of drought evolution on climate variables using the standardized regression coefficient and concluded that the most sensitive variable in south China, north China and the Qinghai-Tibet plateau was P, followed by U2, Tmean, RH, and Rs. Moreover, in northwest China, the most sensitive variable was identified as U2, followed by P, Tmean, RH, and Rs. These results are similar to those of our study, yet the order of sensitivities differs. This may be attributed to differences in the study areas and sensitivity models. Numerous studies have demonstrated that the SPEI is both multi-scale and comprehensive. In particular, in the context of global warming, considering Tmean, Tdelta, RH, Rs, U2, G, and P, the SPEI is able to effectively identify different drought types [7,[9][10][11]. However, the effects of Tmean, Tdelta, RH, Rs, U2, G, and P on the SPEI still remain unclear, especially at different time scales. Thus, we further investigated the sensitivity of SPEI with respect to meteorological variables at different time scales in Jilin Province. The specific time scales investigated were 1 to 36 months. The model was reconstructed in conjunction with the SPEI, and EFAST was applied for global sensitivity analysis on the SPEI. The EFAST model show that the sensitivity of the SPEI differs in time scales. At the time scale of 1 to 6 months, SPEI are the most sensitive to P, followed by Tmean, RH, U2, Rs, G, and Tdelta. At the time scale of 7 to 18 months, P is the most sensitive climatic variable to the SPEI, followed by U2, RH, Tmean, G, Rs, Tdelta. The order changes to P > U2 > G > Tmean > RH > Rs > Tdelta when the time scale increases up to more than 24 months. The significant response of SPEI to time scales is potentially caused by the multi-scalar phenomenon of the SPEI or drought. The SPEI at the 1 to 6 month time scale is most appropriate to monitor drought impacts with regard to meteorological and agroecological systems, while the SPEI at the 7 to 18 month time scale is appropriate for the drought identification for hydrological systems [8,9,63]. Some other authors also report that 1, 3, and 24 month time scales are used as a good proxy for meteorological, agricultural, and hydrological [30].
With more and more notable effects of global warming, the climate change impact on the SPEI is very complex. Our work focuses on a preliminary analysis of a possible tendency of climatic variables and the sensitivity of the SPEI in response to them, as well as the multi-scalar characteristics using the FAO-56 PM methods. The impact of the topographic feature on the SPEI and variations in SPEI resulting from different ETo calculation methods, not considered in this study, should be investigated in future work.

Conclusion
Jilin Province experienced a significant rise in Tmean, and a sharp decrease in Tdelta, Rs, and U2 in the past 60 years. P is undoubtedly the most influential one to the SPEI among the seven selected meteorological variables, which explained 59.9%-97.9% of the total variability, especially during the main crop growing season (from May to September), while Tmean, RH, or U2 observably affect the SPEI and cannot be neglected during the nongrowing season. In terms of spatial distribution, the SPEI was mainly affected by P in the eastern region, while it was also influenced by Tmean, RH, and U2 as well in the western region. The sensitivity of the SPEI differs in time scales: P > Tmean > RH > U2 > Rs > G > Tdelta (1 to 6 month), P > U2 > RH ≈Tmean > G > Rs > Tdelta (7 to 18 month), and P > U2 > G > Tmean > RH > Rs > Tdelta (more than 24 month time scale), respectively. The results have the potential to provide a reference for agricultural production and management in Jilin Province, China.