Special Issue: Ecohydrology

Attribution of trends in meteorological drought during 1960-2016 over the Loess Plateau, China

  • GUO Mengyao , 1, 2 ,
  • SHE Dunxian , 1, 2, * ,
  • ZHANG Liping 1, 2 ,
  • LI Lingcheng 3 ,
  • YANG Zong-Liang 3 ,
  • HONG Si 1, 2
Expand
  • 1. State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan 430072, China
  • 2. Hubei Key Laboratory of Water System Science for Sponge City Construction, Wuhan University, Wuhan 430072, China
  • 3. Jackson School of Geosciences, University of Texas at Austin, Austin, TX 78712-1692, USA
*She Dunxian, E-mail:

Guo Mengyao (1994-), PhD Candidate, specialized in atmosphere-ecosystem interactions. E-mail:

Received date: 2020-09-22

  Accepted date: 2021-05-10

  Online published: 2021-10-25

Supported by

National Natural Science Foundation of China(41877159)

The National Key Research and Development Program of China(2017YFA0603704)

The Scholarship from China Scholarship Council(CSC), No(201906270109)

Abstract

This study uses two forms of the Palmer Drought Severity Index (PDSI), namely the PDSI_TH (potential evapotranspiration estimated-by the Thornthwaite equation) and the PDSI_PM (potential evapotranspiration estimated by the FAO Penman-Monteith equation), to characterize the meteorological drought trends during 1960-2016 in the Loess Plateau (LP) and its four subregions. By designing a series of numerical experiments, we mainly investigated various climatic factors' contributions to the drought trends at annual, summer, and autumn time scales. Overall, the drying trend in the PDSI_TH is much larger than that in the PDSI_PM. The former is more sensitive to air temperature than precipitation, while the latter is the most sensitive to precipitation among all meteorological factors. Increasing temperature results in a decreasing trend (drying) in the PDSI_TH, which is further aggravated by decreasing precipitation, jointly leading to a relatively severe drying trend. For the PDSI_PM that considers more comprehensive climatic factors, the drying trend is partly counteracted by the declining wind speed and solar radiation. Therefore, the PDSI_PM ultimately shows a much smaller drying trend in the past decades.

Cite this article

GUO Mengyao , SHE Dunxian , ZHANG Liping , LI Lingcheng , YANG Zong-Liang , HONG Si . Attribution of trends in meteorological drought during 1960-2016 over the Loess Plateau, China[J]. Journal of Geographical Sciences, 2021 , 31(8) : 1123 -1139 . DOI: 10.1007/s11442-021-1888-y

1 Introduction

Drought, characterized by below-normal available water for a region over an extended period, is a complex and recurring natural hazard that severely impacts the environment and social economy (Wilhite, 2000; Dai, 2011b; Zhang et al., 2019; Li et al., 2020). A substantial increase in terrestrial aridity, mainly driven by anthropogenic global warming, has been revealed by both the observational records and the model simulations since the latter half of the 20th century (Nicholls, 2004; Dai, 2011b; Annamalai et al., 2013; Dai, 2013; Mueller and Zhang, 2016). Moreover, droughts are projected to become more frequent and widespread worldwide in the coming decades (Sheffield and Wood, 2008; Cook et al., 2015; Li et al., 2015), which will challenge the security of water resources, the sustainability of socioeconomic development, ecosystem health, human lives, and so forth. Therefore, characterizing the changing patterns of drought and clarifying their attribution is valuable for understanding the drought propagation mechanism, initiating drought preparedness and mitigation, improving our ability to adapt to climate change, and promoting regional sustainable development.
Drought is usually categorized into four types: meteorological, agricultural, hydrological, and socioeconomic drought (Javadinejad et al., 2020). Meteorological drought, referring to a prolonged period with below-normal precipitation (Palmer, 1965), has attracted much interest because it is easy to be understood and is usually considered as a prerequisite for the other three types of drought (Wilhite, 2000; Dai, 2011b; She and Xia, 2018). Many indices have been developed to assess meteorological drought, such as the Standardized Precipitation Index (SPI) (McKee et al., 1993), the Standardized Precipitation Evapotranspiration Index (SPEI) (Vicente-Serrano et al., 2010), and the Palmer Drought Severity Index (PDSI) (Palmer, 1965). Among them, the PDSI is physically-based and considers precipitation, evaporation, runoff, and soil moisture within the soil water balance modeling, containing a long-term memory of previous moisture conditions for a given month (Bonsal et al., 2013). The PDSI has been widely applied in drought monitoring, assessment, and predictions in many regions worldwide (Dai et al., 2004; Zhang et al., 2019). As a vital input variable in the PDSI algorithm, potential evapotranspiration (PET) was initially estimated by Palmer (1965) using the empirical, temperature-based Thornthwaite equation (denoted as PET_TH) (Thornthwaite, 1948). Alternatively, the Food and Agriculture Organization of the United Nations (FAO) Penman-Monteith equation (denoted as PET_PM) (Allen et al., 1998) is a physically-based PET model, whose parameterization includes the radiative and advective dynamics at the land surface-atmosphere interface (Hobbins et al., 2012; She et al., 2017). Some studies (Hobbins et al., 2008; Sheffield et al., 2012; Zhang et al., 2016) indicate that the drought trends under climate change given by the PDSI based on the PET_TH (PDSI_TH) or the PET_PM (PDSI_PM) could conflict due to the different datasets and PDSI formulation. This discrepancy highlights the importance of comparing different PET estimations in detecting meteorological drought trends in various regions.
The attribution of drought trends is important for better understanding the drought formulation mechanism and the interaction between the atmosphere and ecosystem. Some previous studies have provided a brief view of the causes of drying or wetting trends under climate change based on the SPI, SPEI, and PDSI, which helps drought variability investigation and drought prediction. For example, using the PDSI, Dai et al. (2004) revealed a drying trend over northern and southern Africa, the Middle East, Mongolia, and eastern Australia during 1870-2002 resulting from the combined effect of changes in precipitation and surface temperature; Sun et al. (2016) quantified the contributions of climatic factors to the variations in the SPEI-identified drought for Southwest China using a separating method; Zhang et al. (2016) used both the PDSI_TH and PDSI_PM at annual time scale to assess drought changes over China during 1961-2013 and found that they had quite different responses to warming due to the different PET calculations; Li et al. (2017) found that the PDSI-identified droughts in Central Asia during 1965-2014 are more sensitive to air temperature than to other meteorological variables based on numerical experiments; Zhang et al. (2019) indicated that increases in precipitation variability, temperature, and net radiation along with a decrease in relative humidity resulted in an increasing extreme drought rate in the 21st century in China, but it was partly offset by an increasing trend in precipitation. All these studies provide useful information to support a comprehensive understanding of the attribution of drought variability on annual or even multiannual time scales. However, the attribution of seasonal drought trends, which is vital for agricultural production and the ecological environment, is still limited.
The Loess Plateau (LP), located in north-central China, is one of the most important agricultural regions in China (Liu et al., 2016). Many studies have reported a drying trend over the LP under climate warming (Li et al., 2012; Sun et al., 2015; Liu et al., 2016), which leads to severe water shortages and threatens both the ecological environment and local agriculture. However, research on attributing this drying trend to changing climatic factors is still a gap, leaving it difficult to acquire useful information to help reduce the adverse impacts of drought. The majority of the LP region crops grow under rain-fed conditions, mainly relying on the precipitation falling from June to September (Qin et al., 2013; Liu et al., 2016). Therefore, the attention to drought trends in summer (June-July-August, JJA) and autumn (September-October-November, SON) is also important.
In this study, we assess the meteorological drought trends during the period 1960-2016 over the LP based on the long-term PDSI series at annual, summer, and autumn time scales and then attribute the trends to climatic variables changes. The major objectives of this study are to: (1) respectively investigate the spatial and temporal changes in meteorological drought characterized by the PDSI_TH and the PDSI_PM; (2) provide an integrated comparison of the sensitivity of meteorological drought to contributing climatic factors using the probability density function from kernel density estimation; and (3) assess the contributions of individual climatic variables to the drying/wetting trends by designing a series of numerical experiments.

2 Study area, data, and methodology

2.1 Study area

The LP is located in the north-central part of China and traversed by the upper-middle reaches of the Yellow River, ranging between 33°-42°N and 99°-114°E with an area of roughly 6.4×105 km2 (Figure 1). Approximately, the long-term annual precipitation in this region is 123.3-948.9 mm, and the annual mean temperature is 4.3-14.3°C (Fang et al., 2019). The precipitation in summer and autumn accounts for approximately 70% of annual amount in the form of high-intensity storms, which often caused serious soil erosion, and drought occurs frequently (Zhang et al., 2012; He et al., 2019). Considering the water supply during summer and autumn is crucial for agricultural production and eco-environmental stability, we extend this study to these two seasonal time scales besides the annual time scale.
Figure 1 Location of the Loess Plateau and the 79 meteorological stations and the distribution of the four subregions based on K-means clustering

2.2 Data

The input climatic variables for the period 1960-2016 from the 79 meteorological stations located around the LP region (Figure 1) are provided by the China Meteorological Administration (CMA) (website: http://data.cma.cn) and subject to high-quality control before their release to the public. These data include the daily air temperature (℃) at 2 m height, relative humidity (%) at 2 m height, wind speed (m/s) at 10 m height, sunshine duration (h), air pressure (kPa), and precipitation (mm).

2.3 Methodology

2.3.1 Calculation of the Palmer Drought Severity Index

The PDSI was first proposed by W. C. Palmer in 1965 to evaluate meteorological drought based on a simple two-layer soil water balance model (Palmer, 1965) and has become a very widely used meteorological drought index worldwide (Kogan, 1995; Zou et al., 2005; Zhai et al., 2010; Van Der Schrier et al., 2011; Dash et al., 2012). In the PDSI algorithm, the severity of the meteorological drought is quantified by measuring the degree to which monthly accumulative precipitation departs from its climatically appropriate for existing conditions (CAFEC) value. The simple description of the PDSI calculation is given below; for details, refer to Palmer (1965) and Dai (2011a).
The primary basis for PDSI computation is to obtain monthly hydrologic accounting over the evaluation period, which includes five parameters in a two-layer soil water balance model: precipitation (P), soil moisture loss (L), soil moisture recharge (R), evapotranspiration (ET), and runoff (RO), as well as the potential values of the last four parameters, i.e., potential soil moisture loss (PL), potential soil moisture recharge (PR), potential evapotranspiration (PET), and potential runoff (PRO). Using the two-layer soil water balance model, we can obtain the precipitation under the CAFEC condition, which represents the expected amount of atmospheric water supply during a specific month i, in the following Equation (1):
$\begin{matrix} \widehat{P}=\widehat{ET} \\ \end{matrix}+\widehat{R}+\widehat{RO}-\widehat{L}$
where $\widehat{P}$ is the CAFEC P, which is a key parameter in the PDSI model for measuring the departure of atmospheric moisture from the normal level. Similarly, $\widehat{ET}$, $\widehat{R}$, $\widehat{RO}$, and $\widehat{L}$ are the CAFEC ET, R, RO and L in “normal” weather, respectively, which can be estimated based on the following water balance coefficients:
${{\alpha }_{i}}=\frac{\overline{E{{T}_{i}}}}{\overline{PE{{T}_{i}}}},{{\beta }_{i}}=\frac{\overline{{{R}_{i}}}}{\overline{P{{R}_{i}}}},{{\gamma }_{i}}=\frac{\overline{R{{O}_{i}}}}{\overline{PR{{O}_{i}}}},{{\delta }_{i}}=\frac{\overline{{{L}_{i}}}}{\overline{P{{L}_{i}}}},$
where α, β, γ, and δ are the coefficients of ET, R, RO, and L, respectively; i indicates the month from 1 to 12; each month's mean is averaged over the long-term series for the same month during the evaluation period. With these coefficients that represent the relationship between the actual values and potential values, the CAFEC ET, R, RO, and L can be calculated as
$\begin{matrix} \widehat{ET}=\alpha PET,\widehat{R}=\beta PR,\widehat{RO}=\gamma PRO,\widehat{L}=\delta PL \\ \end{matrix}$
Hence, Equation (1) can be further extended as
$\begin{matrix} \widehat{P}=\alpha PET+\beta PR+\gamma PRO-\delta PL \\ \end{matrix}$
In this study, the PET is calculated respectively by the Thornthwaite and FAO Penman-Monteith equations (Thornthwaite, 1948; Allen et al., 1998):
$\begin{matrix} \begin{matrix} PET\_TH=16{{\left( 10T/H \right)}^{A}}, \\ PET\_PM=\frac{0.408\Delta \left( {{R}_{n}}-G \right)+\gamma \frac{900}{T+273}{{u}_{2}}\left( {{e}_{s}}-{{e}_{a}} \right)}{\Delta +\gamma \left( 1+0.34{{u}_{2}} \right)}, \\ \end{matrix} \\ \end{matrix}$
where T is the monthly averaged air temperature at 2 m height; H is the annual heat index, calculated with $H=\underset{i=1}{\overset{12}{\mathop \sum }}\,{{\left( \frac{{{T}_{i}}}{5} \right)}^{1.514}}$; A is a constant, calculated with $A=6.75\times {{10}^{-7}}{{H}^{3}}-$ $7.71\times {{10}^{-5}}{{H}^{2}}+1.792\times {{10}^{-2}}H+0.49$; Rn is the net radiation; Δ is the slope of the vapor pressure curve; G s the soil heat flux; γ is the psychometric constant; u2 is the wind speed at 2 m height; es is the saturation vapor pressure at 2 m height; ea is the actual vapor pressure at 2 m height.
Then, the P departure (d) from its CAFEC value for each month is calculated as
$\begin{matrix} d=P-\widehat{P} \\ \end{matrix}$
Palmer (1965) converts d to an index of moisture anomaly (Z) by multiplying it with a climatic coefficient K, which represents averages for climate characteristic for that month:
$\begin{matrix} Z=dK \\ \end{matrix}$
The first approximation of K(k) is estimated by the ratio of average moisture demand to average moisture supply:
$\begin{matrix} k=\frac{\overline{PET}+\overline{R}}{\overline{P}+\overline{L}} \\ \end{matrix}$
The final K is modified in accordance with local climatic conditions. The drought severity for a month depends on the moisture anomaly for the current month and on the drought severity for the previous month. As a result, the PDSI value for month t(Xt) is computed as
$\begin{matrix} {{X}_{t}}=p{{X}_{t-1}}+q{{Z}_{t}}, \\ \end{matrix}$
where Xt-1 is the PDSI for the previous month t-1; Zt is the moisture anomaly index for the current month t; and the p and q coefficients are calculated based on local climates.

2.3.2 Quantifying the contributions of climatic variables to the PDSI-identified drought trends

In this study, we employ the numerical experiment method developed by Zhang et al. (2016) to quantify the individual contribution of each climatic variable, i.e., precipitation (P), air temperature (T), wind speed (WS), solar radiation (Rs), and relative humidity (RH), to the PDSI-identified drought trends. This method separates the effect of an individual variable on the drought trend by preserving the original trend in this variable while removing the trends of the other variables (Williams et al., 2015; Li et al., 2017). Table 1 shows the seven experimental cases needed in the study. It should be noted that WS, Rs, and RH cases are only applied in the PDSI_PM as they are excluded in the PDSI_TH calculation. Climatic factors are detrended by the following equation:
$\begin{matrix} {{F}_{det,{{y}_{i}}}}={{F}_{org,{{y}_{i}}}}-\alpha \left( {{y}_{i}}-{{y}_{0}} \right), \\ \end{matrix}$
where yi indicates the year over the evaluation period (1960, 1961, …, 2016 in this study); y0 is the start year; α is the linear trend in the climatic factor series; ${{F}_{org,{{y}_{i}}}}$ is the original climatic factor in yi; ${{F}_{det,{{y}_{i}}}}$ is the detrended climatic factor in the same year. This processing may induce a P, WS, and Rs less than zero or an RH less than zero or greater than one, which is not physically possible. We have modified this by setting any P and Rs less than zero to zero, any WS and RH less than zero to 0.01, and any RH greater than one to 0.99.
Table 1 Framework for the numerical experiments
Case Denotation Computation
Original PDSI_TH_Org PDSI_TH calculated on all original climatic variables
PDSI_PM_Org PDSI_PM calculated on all original climatic variables
Base PDSI_TH_Base PDSI_TH calculated on all detrended climatic variables
PDSI_PM_Base PDSI_PM calculated on all detrended climatic variables
P PDSI_TH_P PDSI_TH calculated on original P and detrended T
PDSI_PM_P PDSI_PM calculated on original P and detrended T, WS, Rs, and RH
T PDSI_TH_T PDSI_TH calculated on original T and detrended P
PDSI_PM_T PDSI_PM calculated on original T and detrended P, WS, Rs, and RH
WS PDSI_PM_WS PDSI_PM calculated on original WS and detrended P, T, Rs, and RH
Rs PDSI_PM_Rs PDSI_PM calculated on original Rs and detrended P, T, WS, and RH
RH PDSI_PM_RH PDSI_PM calculated on original RH and detrended P, T, WS, and Rs
Based on the experimental cases, the individual contribution of one climatic factor's change to drought trends (denoted as $Con_F$) can be quantified by
$\begin{matrix} Con\_F=\partial PDSI\_F/\partial t-\partial PDSI\_Base/\partial t, \\ \end{matrix}$
where $\partial PDSI\_F/\partial t$ is the trend in the PDSI time series in this climatic factor case; $\partial PDSI\_Base/\partial t$ is the trend in the PDSI time series in the base case. According to the dry/wet classification of the PDSI, a positive PDSI trend means a wetting trend, while a negative PDSI trend indicates a drying trend. Therefore, if the $Con\_F$ is positive, this climatic factor contributes to a wetter condition; on the contrary, this climatic factor leads to a drier situation.
We also employ the kernel density estimation method (Węglarczyk, 2018) to describe the probability density functions (PDFs) of the trends in the PDSI time series at all meteorological stations. Comparing the changes of the PDF curves in various cases gives the sensitivity of PDSI-identified droughts to climatic variables (Cook et al., 2015; Zhang et al., 2016; Li et al., 2017).

2.3.3 K-means data clustering algorithm

To better reveal the spatial variability of meteorological drought trends, we divide the whole LP into several subregions based on the location, atmospheric moisture condition, and drought trend at each meteorological station using the K-means (Ghaemi et al., 2009) clustering algorithm. The K-means clustering algorithm, an unsupervised method for data clustering, categorizing items into k groups based on the similarity that is typically measured by squared Euclidean distance. The algorithm is given as follows:
Let $M=\left( {{{\vec{s}}}_{1}},{{{\vec{s}}}_{2}},\ldots,{{{\vec{s}}}_{n}} \right)$ be the set of meteorological stations (n = 79 in this study). Each ${{\vec{s}}_{i}}$ is an attribute vector of the i-th meteorological station. First, k cluster center points are initialized randomly, the so-called mean. Then, each $\vec{s}$ is categorized to its nearest mean, and the mean can then be further updated to the average of the thus-far clustered vectors. This process will be repeated until the mean no longer changes. Finally, the clusters on M are obtained, denoted as $\widehat{M}=\left( \widehat{{{C}_{1}}},\widehat{{{C}_{2}}},\ldots,\widehat{{{C}_{K}}} \right)$, where $\widehat{{{C}_{j}}}$ is the mean of the j-th cluster. The optimal number of subregions is determined by the Akaike information criterion (AIC) (Aho et al., 2014), a standard for evaluating statistical models' goodness of fit. The smaller the AIC value is, the better the model.
Notably, the attribute vector used here consists of latitude, longitude, altitude, annual P, annual PET_TH, and annual PDSI_TH trend. Among them, latitude, longitude, and altitude help cluster the closer stations into a subregion; annual P and PET help cluster the stations with similar atmospheric moisture condition into a subregion; annual PDSI trend helps cluster the stations with similar drought trend into a subregion. Therefore, the clustering algorithm adopted divides the whole LP region into several subregions with similar climate conditions and trends, which is reasonable for studying the drought trends under climate change. Considering some of the input variables in the FAO Penman-Monteith approach are not applicable for the Thornthwaite approach, we use PET_TH (PDSI_TH) instead of PET_PM (PDSI_PM) here to ensure consistency.

3 Results and discussion

3.1 The division of the whole Loess Plateau into subregions

In the subregion number evaluation, AIC obtains its minimum when k equals four. Therefore, the entire LP can be divided into four subregions according to the K-means clustering algorithm, i.e., SR1, SR2, SR3, and SR4, located in the north-central, southwestern, southern, and southeastern parts of the LP, respectively (Figure 1). The station numbers, annual mean P and PET_TH, and the aridity index (AI) in each subregion are shown in Table 2, from which we can see SR2 is relatively humid while SR1 is the driest subregion.
Table 2 Station numbers and climatic condition in each subregion
Stations P (mm) PET (mm) AI (P/PET)
SR1 26 337.8 554.4 0.61
SR2 13 433.5 429.6 1.01
SR3 13 392.2 515.2 0.76
SR4 27 591.7 641.2 0.92

3.2 Temporal changes in climatic factors, PET, and PDSI-identified droughts

Figures 2a-2o shows the temporal changes of the five climatic factors (P, T, WS, Rs, and RH) at annual, summer, and autumn time scales for the whole LP and each subregion from 1960 to 2016. The significance of the linear trends is tested using the Student's t-test (Hox et al., 2017). P (Figures 2a-2c) shows a slightly decreasing trend in the whole LP and most subregions (p > 0.05), which is in good agreement with other studies in the same area (Wang et al., 2012; Sun et al., 2015). Still, a slightly increasing trend is in the annual and summer P in SR2 and in the autumn P in SR1. T has been significantly (p < 0.05) increased in all subregions and at all time scales (Figures 2d-2f), consistent with the widely confirmed climate warming. WS declines significantly (p<0.05) (Figures 2g-2i) in the whole LP. The downward trend in near-surface WS, termed as “stilling”, has been observed in many mid-latitude regions over the past decades (Roderick et al., 2007; McVicar and Roderick, 2010), caused by both natural variation (complicated interactions within the climate system) and anthropogenic factors (e.g., increased terrestrial surface roughness from urbanization). Rs shows negative trends in each subregion and at each time scale (Figures 2j-2l), supporting the results given by Wild (2009) that there is a worldwide “dimming” in the 20th century. Global dimming closely links to the properties of the atmosphere, which have been changed by the increases in aerosols and other air pollutants (Stanhill and Cohen, 2001). For RH, there are slight decreases (p > 0.05, except for annual and summer trends in SR1 and the annual trend in the whole LP) from 1960 to 2016 (Figures 2m-2o). This confirms the result given by Byrne and O'Gorman (2016) that there is a decrease in RH over continents with global warming in climate model simulations, which is because the increasing rate of the moisture transported from ocean to land cannot keep pace with the increasing rate of atmospheric moisture demand over land.
Figure 2 Time series of P, T, WS, Rs, RH, and PET at annual, summer, and autumn time scales for the averages of the whole LP and each subregion (1960-2016). The numbers in the top corner are the changing trends fitted by the linear regression model, with the different colors corresponding to different regions according to the legend. Underlined numbers indicate the trends have passed the Student's t-test of p < 0.05.
Figures 2p-2u display the PET_TH and PET_PM time series during 1960-2016. PET_TH (Figures 2p-2r) shows significantly positive trends (p < 0.05, except for the summer in SR4), with a rate ranging from 0.162 mm/yr2 (autumn in SR2) to 0.848 mm/yr2 (annual in SR3). As T is the sole climatic input for PET_TH, the temporal changes of PET_TH are somewhat similar to those of T (Figures 2d-2f and 2p-2r). PET_PM shows decreasing trends in SR2 and SR4 as well as for the whole region averages, which are quite different from the increasing trends of PET_TH (Figures 2p-2r and 2s-2u). The different climatic inputs could explain the large differences in the changing trends between PET_TH and PET_PM (Sun et al., 2016; Nouri et al., 2017). Generally, the decline in WS and Rs can decrease the aerodynamic and radiative components of PET, ultimately reducing the evaporation rate (Roderick et al., 2007); a higher RH value can often hinder evaporative processes, whereas a lower value favors evaporation (Sun et al., 2016). Take the SR1 in summer as an example. The rising T leads to the increasing PET_TH with the trend of 0.276 mm/yr (Figure 2q). However, although the increase in T (0.021℃/yr) and decrease in RH (-0.066%/yr) can increase the PET_PM, the decline in WS (-0.006 m/s/yr) and Rs (-0.023 MJ/m2d/yr) can decrease the PET_PM; therefore, the PET_PM ultimately shows a decreasing trend of -0.19 mm/yr (Figure 2t) under the joint effect of all four climatic variables (Sun et al., 2016; Nouri et al., 2017). Trenberth et al. (2014) have concluded that using different methods to estimate PET can cause apparently conflicting results in detecting PDSI-identified droughts. Our work further implies that a comparison of trends in the PET_TH and PET_PM is vital for PDSI drought trends analysis.
The time series of the PDSI_TH and PDSI_PM are shown in Figure 3. All series show declining trends, with a rate ranging from -0.091 to -0.02 per year (p < 0.05) for PDSI_TH and from -0.061 to -0.022 per year for PDSI_PM. Since P is the same input for both PDSI forms, different PET estimations should be responsible for the declining trend in the PDSI_TH larger than that in the PDSI_PM. As analyzed previously, incorporating WS, Rs, and WS in the PET_PM calculation leads to the difference between the PET_PM and the PET_TH, and this further results in different trends between the PDSI_PM and the PDSI_TH. Black curves are the time series of the PDSI_TH minus the PDSI_PM, from which the difference in drought severity is displayed. Interestingly, despite the smaller drying trend in the PDSI_PM, the drought severity given by the PDSI_PM is generally higher than that given by the PDSI_TH before the 21st century, which may be the result of a higher PET_PM than PET_TH.
We examine the correlation between the PDSI_PM and PDSI_TH using Pearson's correlation coefficient (r) (Table 3). The high linear correlation indicates the consistency between the PDSI_PM and PDSI_TH variability, which is also reported by the study of Zhang et al. (2016).
Table 3 Pearson's correlation coefficient (r) between the PDSI_TH and PDSI_PM time series during 1960-2016 (all values have passed the Student's t-test of p < 0.001)
Annual Summer Autumn
LP 0.870 0.871 0.901
SR1 0.897 0.895 0.926
SR2 0.769 0.775 0.859
SR3 0.930 0.928 0.929
SR4 0.887 0.888 0.903
Figure 3 Time series of the PDSI_TH, PDSI_PM and PDSI_TH minus the PDSI_PM at annual, summer, and autumn time scales for the averages of the whole LP and each subregion (1960-2016). The numbers in the top right corner are the changing trends fitted by the linear regression model, with the different colors corresponding to different PDSI according to the legend. Underlined numbers indicate the trends have passed the Student's t-test of p < 0.05.

3.3 Attributing the trends in meteorological droughts based on the PDSI

Before attributing the trends in meteorological droughts characterized by the PDSI_PM and PDSI_TH using Equation (11) based on the designed numerical experiments, the efficiency and rationality of this method should be examined through a comparison of Cal_∂PDSI/∂t versus Obs_∂PDSI/∂t at all meteorological stations (Figure 4). Cal_∂PDSI/∂t denotes the calculated PDSI trends, which is the sum of contributions from individual climatic factors; Obs_∂PDSI/∂t denotes the observed PDSI trends, which is the result from subtracting the PDSI_Base trends from the PDSI_Org trends. A better correlation between Cal_∂PDSI/∂t and Obs_∂PDSI/∂t indicates a better fit of the numerical experiments. From Figure 4, Cal_∂PDSI/∂t fits quite well with Obs_∂PDSI/∂t, with the very high Pearson correlation coefficient r at all the time scales for both the PDSI_TH and PDSI_PM. Therefore, the designed numerical experiments are reasonable and reliable, and Equation (11) can be used to quantify the contributions of climatic factors to meteorological droughts.
Figure 4 Fitting results between Cal_∂PDSI/∂t (the calculated PDSI trend, the sum of individual contributions of each climatic factor) and Obs_∂PDSI/∂t (the observed PDSI trend, the result of the PDSI_Org trend minus the PDSI_Base trend) from the 79 meteorological stations. The r values in the top left corner are the Pearson's correlation coefficients between Cal_∂PDSI/∂t and Obs_∂PDSI/∂t.

3.3.1 The sensitivity of the PDSI-identified droughts to climate change

We graphed the PDFs of the PDSI trends from all 79 meteorological stations based on the kernel density estimation (Figure 5). The sensitivity of drought trend to one climatic factor's change is easy to be obtained by comparing the PDF under this factor's case with that under the base case. The PDFs of the PDSI trends in the base cases are roughly normally distributed, with a mean of nearly zero, which is similar to the findings of Zhang et al. (2016) and Zhang et al. (2019). The near-zero mean of the PDSI trends in the base case implies that the PDSI series has no trend when all climatic variables are detrended, proving the detrending method used in the study is effective.
For the PDSI_TH (Figures 5a-5c), the PDFs in the P case move to the left compared to the base case (by -0.018, -0.015, and -0.020 at the annual, summer, and autumn time scales, respectively), which suggests PDSI_TH shows a drier trend due to the decreasing P. In the T case, the PDFs move even further to the left compared to the base case, by -0.026, -0.028, and -0.029 at the annual, summer, and autumn time scales, respectively. Therefore, the drying trend given by PDSI_TH is more sensitive to temperature rising than to precipitation decreasing. However, some studies point out that the PDSI_TH may overrespond to T because of its flawed PET model (Sheffield et al., 2012; Zhang et al., 2016). Notably, the right
Figure 5 The PDFs of the PDSI trends in the base case and the climatic factors' cases from the 79 meteorological stations. The numbers in the top right corner are the mean and standard deviation values, with the different colors indicating different experimental cases corresponding to the legend.
heavy-tailed PDF curves in the P case (Figures 5a-5c) indicate that some stations show a wetter trend due to the increasing P. Besides, the standard deviation (σ) in the P case changes more than that in the T case, which implies that, by comparison, the impact of P on drought trend is more inhomogeneous while the impact of T is more uniform spatially.
Different from the PDSI_TH, for the PDSI_PM (Figures 5d-5f), the PDFs in the P case move more to the left compared with the base case (by -0.015, -0.012, and -0.016 at the annual, summer, and autumn time scales, respectively) than those in the T case (by -0.011 at all time scales). Therefore, the drying trend given by PDSI_PM is more sensitive to precipitation than temperature. Also, the small left movements in the RH case indicate that the drying trend is not very sensitive to RH. In particular, the PDFs in the WS and Rs cases move to the right a little bit compared to the base case, suggesting these two variables cause a slight wetter trend.

3.3.2 The individual contributions of climatic factors to drought trends

The individual contribution of each climatic factor to the PDSI trend is estimated by Equation (11) and shown in Figure 6. From the black bars representing the observed PDSI trends, both the PDSI_PM and PDSI_TH show drying trends over the whole LP, but these trends have high spatial variability. For the PDSI_TH (Figure 6a-c), the drying trend is the largest in SR3, followed by SR4 and SR1, and it is the smallest in SR2. Comparing the individual contributions of P and T, we can find that the contribution of T (red bars) is similar in all subregions, while the contribution of P (blue bars) varies widely in different subregions. In particular, the increasing P in SR2 leads to a wetter trend, but this effect is completed offset by the increasing T, leaving a slight drying trend ultimately. For the PDSI_PM (Figures 6d-6f), the drying trend is the largest in SR3, followed by SR4 and SR2, and it is the smallest in SR1. Both WS and Rs contribute a wetter trend (positive green and purple bars), while such contributions are counteracted by the combined effect of P, T, and RH, ultimately resulting in the drying trends. The significant differences in drying trends between the PDSI_TH and PDSI_PM for the whole LP and all subregions highlight a need for sufficient cautiousness in choosing the PET model when using the PDSI to assess drought.
Unlike most previous studies that focus only on the attribution of drought trends on the annual scale (Zhang et al., 2016; Li et al., 2017), we additionally investigate the drought trends and their contributing factors on the seasonal (summer and autumn) time scale. Among different time scales, climatic factors' contributions to the drought trends are not quite the same, but they have the same negative or positive effects overall. In general, decomposing the annual trends into seasonal time scales gives us more climatic conditions and helps us better understand the drought trend and its attribution.
Figure 6 The individual contributions of climatic factors and the Obs_∂PDSI/∂t (the observed PDSI trend, the result of the PDSI_Org trend minus the PDSI_Base trend) over the LP as a whole and in each subregion

3.4 Uncertainties in the method of the current study

We quantified the contribution of individual climatic factors to the trend in meteorological droughts based on numerical experimental cases; however, there are some uncertainties in this study. They may come from (1) the assumptions of the PDSI model (Palmer, 1965; Alley, 1984), e.g., considering the soil available water capacity (AWC) as the potential P to estimate PRO and simplifying soil evapotranspiration to a process that evapotranspiration does not occur in the soil lower layer until all the moisture evaporates away in the upper layer; (2) the experimental design, e.g., maintaining the trend in one climatic factor while detrending the others may affect the interdependencies of the underlying variables (Zhang et al., 2016); and (3) data processing error, e.g., making an increasing variable slightly decreased in the detrending calculation. Despite the potential uncertainties, the experimental design's rationality has been verified, and the errors in the results are within the acceptable range.
Future study could reduce the uncertainties by the following improvements: (1) using more sophisticated numerical models (e.g., earth system models) and more realistic metrics (e.g., the percentiles of soil moisture, validated) instead of over-simplified drought index to define dry conditions; (2) using more complex detrending method (e.g., additionally removing cyclical components by filter) instead of assuming all variables' time series to follow the simple linear pattern only; (3) excluding the impacts of human activities on the drought trend.

4 Conclusions

The quantitative attribution of the trends in meteorological drought is of great significance for drought forecasting and mitigation, especially for semiarid regions with more droughts expected under future climate change. In this study, we first examined the trends in meteorological drought using the PDSI based on two PET models over the LP in China, then quantified the contribution of climatic factors to the changing trends for drought. Our study provides important insights for attributing the changing aridity under climate change and useful information for regional sustainable development.
The main conclusions are given below:
(1) Based on the different PET estimations using different climatic inputs, the drying trend given by the PDSI_TH is much larger than the PDSI_PM. Nevertheless, variability of the time series of PDSI_TH and PDSI_PM is highly consistent.
(2) The trends in the PDSI_TH are more sensitive to warming than declining precipitation, while the trends in the PDSI_PM are the most sensitive to the decreasing precipitation. The disparity is caused by the integrated effect of different variables in the PDSI_PM calculation.
(3) The drought trends have high spatial variability among different subregions, mainly contributed by the inhomogeneous precipitation changes over the whole LP region.
(4) For the PDSI_TH, the increasing temperature result in a drying trend; this trend is further aggravated by decreasing precipitation, ultimately resulting in a more severe drying trend. For the PDSI_PM, declining wind speed and solar radiation cause a wetting trend, but this trend is counteracted by the combined effect of precipitation, temperature, and relative humidity, still resulting in a drying trend ultimately.
[1]
Aho K, Derryberry D W, Peterson T, 2014. Model selection for ecologists: The worldviews of AIC and BIC. Ecology, 95(3):631-636.

DOI

[2]
Allen R G, Pereira L S, Raes D et al., 1998. Crop Evapotranspiration: Guidelines for Computing Crop Water Requirements. FAO Irrigation and Drainage Paper 56. FAO, Rome, 300(9):D05109.

[3]
Alley W M, 1984. The Palmer Drought Severity Index: Limitations and assumptions. Journal of Climate and Applied Meteorology, 23(7):1100-1109.

DOI

[4]
Annamalai H, Hafner J, Sooraj K et al., 2013. Global warming shifts the monsoon circulation, drying South Asia. Journal of Climate, 26(9):2701-2718.

DOI

[5]
Bonsal B R, Aider R, Gachon P et al., 2013. An assessment of Canadian prairie drought: Past, present, and future. Climate Dynamics, 41(2):501-516.

DOI

[6]
Byrne M P, O'gorman P A, 2016. Understanding decreases in land relative humidity with global warming: Conceptual model and GCM simulations. Journal of Climate, 29(24):9045-9061.

DOI

[7]
Cook B I, Ault T R, Smerdon J E, 2015. Unprecedented 21st century drought risk in the American Southwest and Central Plains. Science Advances, 1(1):e1400082.

DOI

[8]
Dai A G, 2011a. Characteristics and trends in various forms of the Palmer Drought Severity Index during 1900-2008. Journal of Geophysical Research:Atmospheres, 116(D12), doi: 10.1029/2010JD015541.

DOI

[9]
Dai A G, 2011b. Drought under global warming: A review. Wiley Interdisciplinary Reviews:Climate Change, 2(1):45-65.

[10]
Dai A G, 2013. Increasing drought under global warming in observations and models. Nature Climate Change, 3(1):52-58.

DOI

[11]
Dai A G, Trenberth K E, Qian T, 2004. A global dataset of Palmer Drought Severity Index for 1870-2002: Relationship with soil moisture and effects of surface warming. Journal of Hydrometeorology, 5(6):1117-1130.

DOI

[12]
Dash B K, M R, Fahima K et al., 2012. Characteristics of meteorological drought in Bangladesh. Natural Hazards, 64(2):1461-1474.

DOI

[13]
Fang W, Huang S Z, Huang Q et al., 2019. Probabilistic assessment of remote sensing-based terrestrial vegetation vulnerability to drought stress of the Loess Plateau in China. Remote Sensing of Environment, 232:111290.

DOI

[14]
Ghaemi R, Sulaiman M N, Ibrahim H et al., 2009. A survey: Clustering ensembles techniques. World Academy of Science, Engineering and Technology, 50:636-645.

[15]
He G H, Zhao Y, Wang J H et al., 2019. Attribution analysis based on Budyko hypothesis for land evapotranspiration change in the Loess Plateau, China. Journal of Arid Land, 11(6):939-954.

DOI

[16]
Hobbins M, Wood A, Streubel D et al., 2012. What drives the variability of evaporative demand across the conterminous United States? Journal of Hydrometeorology, 13(4):1195-1214.

DOI

[17]
Hobbins M T, Dai A G, Roderick M L et al., 2008. Revisiting the parameterization of potential evaporation as a driver of long-term water balance trends. Geophysical Research Letters, 35(12), doi: 10.1029/2008GL033840.

DOI

[18]
Hox J J, Moerbeek M, Van De Schoot R, 2017. Multilevel Analysis: Techniques and Applications. 3rd ed. New York: Routledge, doi: 10.4324/9781315650982.

DOI

[19]
Javadinejad S, Dara R, Jafary F, 2020. Evaluation of hydro-meteorological drought indices for characterizing historical and future droughts and their impact on groundwater. Resources Environment and Information Engineering, 2(1):71-83.

DOI

[20]
Kogan F N, 1995. Droughts of the late 1980s in the United States as derived from NOAA polar-orbiting satellite data. Bulletin of the American Meteorological Society, 76(5):655-668.

DOI

[21]
Li L C, She D X, Zheng H et al., 2020. Elucidating diverse drought characteristics from two meteorological drought indices (SPI and SPEI) in China. Journal of Hydrometeorology, 21(7):1513-1530.

DOI

[22]
Li L C, Zhang L P, Xia J et al., 2015. Implications of modelled climate and land cover changes on runoff in the middle route of the South to North Water Transfer Project in China. Water Resources Management, 29(8):2563-2579.

DOI

[23]
Li Z, Chen Y N, Fang G H et al., 2017. Multivariate assessment and attribution of droughts in Central Asia. Scientific Reports, 7(1):1-12.

DOI

[24]
Li Z, Zheng F L, Liu W Z, 2012. Spatiotemporal characteristics of reference evapotranspiration during 1961-2009 and its projected changes during 2011-2099 on the Loess Plateau of China. Agricultural and Forest Meteorology, 154:147-155.

[25]
Liu Z P, Wang Y Q, Shao M A et al., 2016. Spatiotemporal analysis of multiscalar drought characteristics across the Loess Plateau of China. Journal of Hydrology, 534:281-299.

DOI

[26]
McKee T B, Doesken N J, Kleist J, 1993. The relationship of drought frequency and duration to time scales. In: Proceedings of the 8th Conference on Applied Climatology, 17(22):2179-2183.

[27]
McVicar T R, Roderick M L, 2010. Atmospheric science: Winds of change. Nature Geoscience, 3(11):747-748.

DOI

[28]
Mueller B, Zhang X B, 2016. Causes of drying trends in northern hemispheric land areas in reconstructed soil moisture data. Climatic Change, 134(1/2):255-267.

DOI

[29]
Nicholls N, 2004. The changing nature of Australian droughts. Climatic Change, 63(3):323-336.

DOI

[30]
Nouri M, Homaee M, Bannayan M, 2017. Quantitative trend, sensitivity and contribution analyses of reference evapotranspiration in some arid environments under climate change. Water Resources Management, 31(7):2207-2224.

DOI

[31]
Palmer W C, 1965. Meteorological drought. Research Paper No.45. Washington DC: US Department of Commerce, Weather Bureau, 59.

[32]
Qin S H, Li L L, Wang D et al., 2013. Effects of limited supplemental irrigation with catchment rainfall on rain-fed potato in semi-arid areas on the western Loess Plateau, China. American Journal of Potato Research, 90(1):33-42.

DOI

[33]
Roderick M L, Rotstayn L D, Farquhar G D et al., 2007. On the attribution of changing pan evaporation. Geophysical Research Letters, 34(17):251-270.

[34]
She D X, Xia J, 2018. Copulas-based drought characteristics analysis and risk assessment across the Loess Plateau of China. Water Resources Management, 32(2):547-564.

DOI

[35]
She D X, Xia J, Zhang Y Y, 2017. Changes in reference evapotranspiration and its driving factors in the middle reaches of Yellow River Basin, China. Science of the Total Environment, 607:1151-1162.

[36]
Sheffield J, Wood E F, 2008. Projected changes in drought occurrence under future global warming from multi-model, multi-scenario, IPCC AR4 simulations. Climate Dynamics, 31(1):79-105.

DOI

[37]
Sheffield J, Wood E F, Roderick M L, 2012. Little change in global drought over the past 60 years. Nature, 491(7424):435-438.

DOI

[38]
Stanhill G, Cohen S, 2001. Global dimming: A review of the evidence for a widespread and significant reduction in global radiation with discussion of its probable causes and possible agricultural consequences. Agricultural and Forest Meteorology, 107(4):255-278.

DOI

[39]
Sun Q H, Miao C Y, Duan Q Y et al., 2015. Temperature and precipitation changes over the Loess Plateau between 1961 and 2011, based on high-density gauge observations. Global and Planetary Change, 132:1-10.

DOI

[40]
Sun S L, Chen H S, Wang G J et al., 2016. Shift in potential evapotranspiration and its implications for dryness/wetness over Southwest China. Journal of Geophysical Research:Atmospheres, 121(16):9342-9355.

[41]
Thornthwaite C W, 1948. An approach toward a rational classification of climate. Geographical Review, 38(1):55-94.

DOI

[42]
Trenberth K E, Dai A G, Van Der Schrier G et al., 2014. Global warming and changes in drought. Nature Climate Change, 4(1):17-22.

DOI

[43]
Van Der Schrier G, Jones P D, Briffa K R, 2011. The sensitivity of the PDSI to the Thornthwaite and Penman-Monteith parameterizations for potential evapotranspiration. Journal of Geophysical Research, 116(D3), doi: 10.1029/2010JD015001.

DOI

[44]
Vicente-Serrano S M, Beguería S, López-Moreno J I, 2010. A multiscalar drought index sensitive to global warming: The standardized precipitation evapotranspiration index. Journal of Climate, 23(7):1696-1718.

DOI

[45]
Wang Q X, Fan X H, Qin Z D et al., 2012. Change trends of temperature and precipitation in the Loess Plateau region of China, 1961-2010. Global and Planetary Change, 92:138-147.

[46]
Węglarczyk S, 2018. Kernel density estimation and its application. In: ITM Web of Conferences, doi: 10.1051/itmconf/2018230003.

DOI

[47]
Wild M, 2009. Global dimming and brightening: A review. Journal of Geophysical Research, 114(D10), doi: 10.1029/2008JD011470.

DOI

[48]
Wilhite D A, 2000. Drought as a natural hazard:Concepts and definitions. In: Wilhite D A (ed). Drought:A Global Assessment, Vol.1. London: Routledge,3-18.

[49]
Williams A P, Seager R, Abatzoglou J T et al., 2015. Contribution of anthropogenic warming to California drought during 2012-2014. Geophysical Research Letters, 42(16):6819-6828.

DOI

[50]
Zhai J Q, Su B D, Krysanova V et al., 2010. Spatial variation and trends in PDSI and SPI indices and their relation to streamflow in 10 large regions of China. Journal of Climate, 23(3):649-663.

DOI

[51]
Zhang B Q, Wu P T, Zhao X N et al., 2012. Drought variation trends in different subregions of the Chinese Loess Plateau over the past four decades. Agricultural Water Management, 115:167-177.

DOI

[52]
Zhang J, Sun F B, Lai W L et al., 2019. Attributing changes in future extreme droughts based on PDSI in China. Journal of Hydrology, 573:607-615.

DOI

[53]
Zhang J, Sun F B, Xu J J et al., 2016. Dependence of trends in and sensitivity of drought over China (1961-2013)on potential evaporation model. Geophysical Research Letters, 43(1):206-213.

DOI

[54]
Zou X K, Zhai P M, Zhang Q, 2005. Variations in droughts over China:1951-2003. Geophysical Research Letters, 32(4), doi: 10.1029/2004GL021853.

DOI

Outlines

/