Joint analysis of changes in temperature and precipitation on the Loess Plateau during the period 1961–2011

The Loess Plateau is particularly sensitive to climate change owing to its fragile ecological environment and geographic features. Here, we present a comprehensive analysis of the joint probabilistic characteristics and tendencies for bivariate and trivariate precipitation and temperature indices across the plateau, based on copula theory. The results show that the southeast region of the plateau had a higher potential for flooding: the 10-year return levels for the number of days with heavy and very heavy precipitation (R10mm, R20mm) and for the maximum 5-day precipitation value (RX5day) were higher in this region. The northwest region of the plateau, however, had a higher potential for drought, as reflected in the high and increasing 10-year return levels for the number of consecutive dry days (CDD) and the number of days with low precipitation (R1mm). In a joint analysis of precipitation indices, large areas of the Loess Plateau showed a relatively high risk of concurrent extreme precipitation events. However, the risk of concurrent extreme wet and dry events did not increase over the past half century, as demonstrated by nonsignificant changes in the probability of concurrently long CDD and long consecutive wet days (CWD). A trivariate copula analysis showed that some grid locations in the southeast of the plateau had an increasing risk of extreme precipitation events occurring at a high frequency and a high intensity, and forming a large percentage of the annual precipitation. Joint analysis of precipitation and temperature indices showed that the risk of higher temperatures and longer spells of consecutive dry days had increased over the past 50 years in grid locations scattered in the northern and southern regions: there were negative trends in the bivariate return periods for warm days (TX90p) and CDD. In addition, there was a decreased probability of concurrent long spells of consecutive wet days and colder temperatures, as demonstrated by the positive trends in the bivariate return periods for cold days (TX10p) and CWD across large areas of the plateau. Overall, the risk of severe floods and droughts over the Loess Plateau increased over the period 1961–2011.


Introduction
Global warming has resulted in significant changes to climates and environments. Climate change includes not only changes in mean climate but also changes in weather extremes (Fischer and Knutti 2015;Miao et al. 2015). Increases in the number of extreme events have been identified at both regional and global scales (Perkins et al. 2012;Rahmstorf and Coumou 2011). For instance, there have been increases in heavy precipitation events in western and northeastern China and in the low reaches of the Yangtze River basin region, coupled with decreases in light precipitation across most areas (Wang and Fu 2013).
Abstract The Loess Plateau is particularly sensitive to climate change owing to its fragile ecological environment and geographic features. Here, we present a comprehensive analysis of the joint probabilistic characteristics and tendencies for bivariate and trivariate precipitation and temperature indices across the plateau, based on copula theory. The results show that the southeast region of the plateau had a higher potential for flooding: the 10-year return levels for the number of days with heavy and very heavy precipitation (R10mm, R20mm) and for the maximum 5-day precipitation value (RX5day) were higher in this region. The northwest region of the plateau, however, had a higher potential for drought, as reflected in the high and increasing 10-year return levels for the number of consecutive dry days (CDD) and the number of days with low precipitation (R1mm). In a joint analysis of precipitation indices, large areas of the Loess Plateau showed a relatively high risk of concurrent extreme precipitation events. However, the risk of concurrent extreme wet and dry events did not increase over the past half century, as demonstrated by nonsignificant changes in the probability of concurrently long CDD and long consecutive wet days (CWD). A trivariate copula 1 3 The risk of extreme summer heat in eastern China has rapidly increased, with the five hottest summers on record all occurring in the twenty-first century ). Madsen et al. (2014) reviewed a large number of studies on trend analysis and climate change projections for extreme precipitation events and floods in Europe and demonstrated a general increase in extreme precipitation in Europe. Powell and Keim (2015) found that, although there was regionwide warming of extreme minimum temperatures and cooling of extreme maximum temperatures in the southeastern United States for the period 1948-2012, the intensity and magnitude of extreme precipitation events increased overall. These changes inevitably pose risks for natural and human systems (IPCC 2014), particularly with regard to water, agriculture, and food security (Schmidhuber and Tubiello 2007;Wheeler and von Braun 2013).
However, most previous studies have focused on changes in the pattern of a single climate variable, and there has been limited analysis of concurrent climatic extremes. Actually, most of the extreme climate events is the combined results of multiple climate variables. For example, high temperatures combined with reduced precipitation triggers the occurrence of heat waves and droughts. The 2014 California drought is archetypical of an event characterized not only by low precipitation but also by extreme high temperatures (AghaKouchak et al. 2014). A simultaneous increase in the intensity and frequency of extreme precipitation events may increase the risk of flooding. Changes in the interaction of a number of climate variables may have a more severe influence on plants and the environment than changes in any single variable (Estrella and Menzel 2013). There have been some studies that assess the relationship between climate variables and extremes (Estrella and Menzel 2013;Liu et al. 2014;Zhang et al. 2013). Zhang et al. (2013) analyzed the joint probabilistic behavior of precipitation extremes using copulas and suggested that the intensification of precipitation extremes has the potential to increase the probability of occurrence of natural hazards, particularly floods and droughts. Liu et al. (2014) presented the joint probabilistic characteristics and tendencies of bivariate precipitation extremes in southwest China, and showed that Yunnan is characterized by a high risk of drought and an increasing probability of concurrent extreme wet and dry events, with decreasing severity of droughts but increasing severity of floods. Meanwhile, although there is no significant trend in meteorological drought, the concurrence of meteorological droughts and heatwaves shows statistically significant substantial increases across the United States from 1960 to 2010 (Mazdiyasni and AghaKouchak 2015). Joint projections of US East Coast sea level and storm surge also revealed future coastal flood risk will also be further increased forced by sea-level rise and changes in the frequency and intensity of tropical cyclones simultaneously (Little et al. 2015). Wahl et al. (2015) pointed out that changes in the joint distributions of storm surge and heavy precipitation associated with climate variability and change also augment risk of compound flooding in major US cities. Simultaneous occurrences of extreme events may result in a significant and substantial impact on the ecosystem and society. Analyzing changes in concurrent climate extremes is critical to preparing for and mitigating the negative effects of climatic change and variability.
The Loess Plateau (Fig. 1a) in China is the most developed region of loess in the world in terms of extent, thickness and depositional sequence and is also the regions with the most serious soil erosion in the world (Fu 1989). The Loess region covers an area of 0.62 million km 2 ; the typical loess area covering 0.42 million km 2 and the area of severe erosion 0.28 million km 2 (Chen et al. 2007). The Loess is prone to vertical cleavage and its surface soil is soft and loose (Fu 1989), that is susceptible to the forces of wind and water (Chen et al. 2008;Wang et al. 2010b). The loess soil was formed in the Quaternary as a result of wind deposition (Liu 1999). The composition of loess generated by rock and soil forming processes results in substantial amounts of fine sandy loam and silt loam with strong tensile and stress joints so that collapse of the soil may happen after only 1-2 min in water (Shi and Shao 2000). The thick loess is extremely liable to suffer from rain and flow damage during storms (Liu 1999). Based on the 2nd national soil erosion survey, over 41 % area of the Loess Plateau suffers the erosion of moderate and higher degrees (Fig. 1b). The high-intensity storms (Wang et al. 2012), steeply sloping landscapes (Wang et al. 2011), low vegetation cover, and highly erodible soils (Chen et al. 2008;Wang et al. 2010b) make the ecological environment on the Loess Plateau fragile. The terrestrial ecosystems of the Loess Plateau exhibit a comparatively early ecological response to changes in the global climate (Wang et al. 2011). Climate warming accelerates evapotranspiration, aggravates water deficiency, and causes serious soil desiccation and the development of dried soil layer (Wang et al. 2010b). It can negatively affect ecological and hydrological processes and as a result suppresses the vegetation growth (Kong et al. 2015). On the other hand, the potential increase of occurrence of extreme events (flood, drought, etc.) under climate warming would increase the risk of soil and water loss, which directly affects local agricultural and industrial productivity and deplete land resources and degraded the eco-environment in the Loess Plateau (Shi and Shao 2000). Thus, knowledge of climate change over the Loess Plateau would be helpful for ecological restoration of the environment on the plateau. Annual precipitation on the Loess Plateau has decreased and air temperatures have increased over the last 50 years (Bi et al. 2009;Miao et al. 2011;Sun et al. 2015;Wang et al. 2012). These synchronous changes in temperature and precipitation can damage and challenge the Loess Plateau in many ways, including altering water availability, increasing soil erosion, and increasing the occurrence of droughts and floods. Thus, attention should be paid to the changing joint probabilistic characteristics of climate variables over the Loess Plateau.
The objectives of this study were therefore: (1) to use a copula-based analysis to assess the joint probabilistic characteristics and tendencies of bivariate and trivariate precipitation and temperature indices and (2) to evaluate the changes in the probability of droughts, floods, and soil erosion implied by the joint changes in precipitation and temperature.

Study area
The Loess Plateau (35-41ºN, 102-114ºE) is situated at the upper and middle reaches of the Yellow River, and extends from central to northwestern China (Fig. 1). In general, the plateau has a warm or temperate continental monsoon climate with extensive monsoonal influence. The average annual precipitation ranges from 200 mm in the northwest to 750 mm in the southeast, and the rainy season (June-September) accounts for 60-70 % of the total precipitation (Li et al. 2012). Rainfall is mainly from high-intensity storms in the region (Wang et al. 2012). The region spans arid, semiarid, and semihumid zones and is therefore considered to be a semiarid to semihumid transitional zone, which is very sensitive to climate change (Liu and Sang 2013). The Loess Plateau is also one of the most eroded regions in the world because it has highly erodible soils, steep slopes, heavy storms, and low vegetation cover stemming from intensive cultivation and improper land uses (Li et al. 2011). The average annual soil loss rate from the plateau is 5000-10,000 t/km 2 (Zhang et al. 2008).

Data
The observed daily precipitation (PR) datasets were obtained from a gauge-based daily precipitation analysis on a 0.5° × 0.5° grid (Shen et al. 2010). The datasets were developed through the optimal interpolation method involving 2400 meteorological observation stations over mainland China, including 299 stations on the Loess Plateau ( Fig. 1). All report data used in the datasets was subject to quality control including an extremes check and an internal consistency check (Shen and Xiong 2015). The daily maximum temperature (TX) datasets were obtained from a surface temperature 0.5° × 0.5° gridded dataset (V2.0) developed by the National Meteorological Information Center of the China Meteorological Administration (National Meteorological Information Center 2012). Owing to some nonclimatic factors, such as such as changes in observation time, instrumentation, location (altitude and latitude) and nonstandard siting, the number of station over China have changed considerably (Zhou et al. 2004), but remained stable at more than 2000 total stations since 1961. To ensure more observation information involved in, the dataset was constructed since 1961. Although the precipitation and temperature datasets was updated to 2012, there were some missing data in maximum temperature data in 2012. Considering the consistency of involved climate variables over a long timeframe, we finally focused on the period between 1961 and 2011 in this study.  (Sillmann et al. 2013). Nine indices are for extreme precipitation events and two are for extreme temperature events.

Marginal probability distributions
To calculate the return periods and estimate the return levels, the marginal distribution of each variable was first constructed. Several types of distribution, including Gaussian, Poisson, exponential, Weibull, generalized extreme value (GEV), binomial, negative binomial, normal, lognormal, generalized Pareto (GP), and extreme value (EV) distributions were used to fit the indices. The Kolmogorov-Smirnov (K-S) test (Massey 1951) was applied to evaluate the goodness-of-fit of each theoretical distribution to the empirical distributions to assess which distribution type was most appropriate for each index. The higher the p value, the greater is the confidence in the fit of the theoretical distribution. The distributions used to fit each index are listed in Table 2.

Copula methodology
A copula is a joint multivariate distribution in which the marginal distribution is uniform over (0, 1) and which can be used to derive the joint distribution of two or more variables, regardless of their original marginal distributions Leonard et al. 2008). Copulas have the advantage that the marginal behaviors and the dependence structure can be studied separately, and have recently emerged as a practical and efficient analysis method for multivariate events (Wang et al. 2010a). Owing to their flexibility and efficiency, copulas have been successfully and widely applied in the field of hydrology (Salvadori and De Michele 2007). An n-dimensional copula (C)  is the joint probability of univariate marginal distributions F 1 , F 2 , . . . , F n as follows: where u 1 , u 2 , . . . , u n refer to the cumulative distribution functions (CDF) of x 1 , x 2 , . . . , x n . Copulas are categorized into several families, among which the Archimedean (Clayton, Gumbel, Frank, etc.) and elliptical copulas (Gaussian and Student's t, etc.) are the two most commonly applied. In recent years, copulas have become popular for analyzing hydrometeorological extremes. In this study, the Gumbel, Clayton, Frank, Gaussian, and t copulas were considered in an analysis of the joint probability distributions of bivariate and trivariate precipitation extremes. To identify the appropriate copulas, the Akaike information criterion (AIC) (Akaike 1974) and square Euclidean distance (SED) were used. AIC includes two parts: lack-of-fit of the model, which can be obtained either by maximizing the likelihood function of the distribution, or the residual sum of squares of the model, and the unreliability of the model due to the number of model parameters. We estimated AIC based on the residual sum of squares directly as previous studies (Zhang and Singh 2007;Zhang et al. 2013). The AIC is expressed as follows: where k is the sample size, m is the number of parameters, and RSS is the residual sum of squares. The smaller the AIC value, the better is the goodness-of-fit of the copula. The SED is calculated as follows: where C k is the values estimated from empirical copulas, and C k refers to the values from theoretical copulas. Also, the smaller the SED value, the better is the goodness-offit of the copula. In addition, the two-sample Cramér-von Mises test (Anderson 1962;Genest et al. 2009) are used for testing the goodness of fit of theoretical copula. A smaller p value means that we accept the null hypothesis of both samples coming from the same underlying distribution and  the confidence in the fit of the theoretical distribution. The multivariate joint distributions are constructed using the 'copula' package in MATLAB. Therefore, the most appropriate copula functions for different combinations of variables were selected on the basis of the AIC, Cramér-von Mises test and SED values (Table 3). Figure 2 presents the relationship between the empirical copulas and the theoretical copulas. Overall, the theoretical copulas selected for the different index combinations showed good performance because the values from the theoretical copulas were close to those from the empirical copulas.
In general, the different combinations of variables chosen in the study represent different climate conditions. T(CDD, CWD) was applied to analyze the probability of concurrence of long dry periods and long wet periods within a year. T(R95p, R95pTOT) and T(R95p, R95pTOT, PI) were constructed to estimate the joint probability of the frequency, amount, and intensity of extreme precipitation. T(R1mm, R10mm) and T(R1mm, R20mm) were used to verify the concurrent changes in wet, heavy, and very heavy precipitation days within a year. T(CWD, RX5day) was used to reflect the characteristics of the extreme wet events. To analyze the concurrence of extreme temperature and precipitation events, the joint precipitation and temperature indices T(TX90p, CDD) and T(TX10p, CWD) were constructed.

Return periods
The return period can be defined as the average elapsed time or mean interarrival time between occurrences of critical events. Joint return periods can be investigated, depending on the selected copula. For bivariate copulas, the return periods were calculated as follows: and for trivariate copulas, the return periods were calculated as follows:

Trend detection
To analyze the trends in the 10-year return levels and the joint return periods for multiple 10-year return levels, a moving window approach was applied (Ghosh et al. 2012;Liu et al. 2014). Ten-year return levels were generated over 30-year moving windows (30 years is a recognized climate timescale). Then, the nonparametric Mann-Kendall methods (Mann 1945) was applied to estimating the statistical significance of trends of 10-year return levels for the indices at each grid point on the Loess Plateau and the magnitudes of the trends were computed with Sen's slope estimator (Sen 1968). The rankbased nonparametric M-K test is more frequently encountered than parametric statistical tests in analysis of hydrometeorological time series (Yue et al. 2002) and make no assumptions about the probability distribution (Onoz and Bayazit 2003). The methods can examine trends in a time series without requiring normality or linearity. For the joint return periods, the 10 years return levels for each indices were calculated based on the entire time series using the selected distribution at each grid. Next, moving window approach (30 years moving window) was used to generate joint distribution of different indices pairs using copula. Based on the selected copula, joint return periods at 10 years return levels were investigated at each 30 years moving window. Finally, the trends of joint return periods were also estimate using the Mann-Kendall method as that used for the 10 years return levels described above. It should be noted that we assume the climate stationarity when estimating trends of the joint return periods. Figure 3 shows the joint cumulative distribution functions, calculated from the marginal distributions for the (5)

Bivariate analyses
where F x (x), F y (y), and F z (z) are the marginal distributions for the precipitation variables X, Y, and Z, F(x, y) and F(x, y, z) are the corresponding bivariate and trivariate joint probability distributions, and T (X > x, Y > y) and T (X > x, Y > y, Z > z) denote the joint return periods for the specific thresholds being exceeded simultaneously. The change of return period indicate the likelihood of joint occurrence of these two phenomena and the compound risk from extreme climate events. different individual indices and from the copula method for different pairs of indices. The dependence structure and inherent correlations between precipitation indices are clearly reflected in Fig. 3. R95p and R95pTOT, R1mm and R10mm, R1mm and R20mm, and CWD and RX5day had high correlations, with the dots that represent observations concentrated along the diagonals. The corresponding correlation coefficients for these four  Figure 4 shows the contour lines for the multivariate return periods for the different extreme precipitation indices, calculated from the copula equation and Eq. 4. The red dots represent the joint return periods for different combinations of extreme precipitation conditions observed during the period 1961-2011 over the Loess Plateau. The concurrent return periods for the precipitation events mostly ranged from 5 to 200 years. Concurrent events occurring in the upper right corner in Fig. 4 correspond to severe extreme precipitation conditions and have longer return periods. For example, the red dot in the upper right corner in Fig. 4a indicates that the Loess Plateau experienced severe extreme precipitation conditions in that particular year, with a large number of very wet days and a large percentage of the annual precipitation derived from very wet days. Overall, for different concurrent extreme conditions, most observed events showed return periods within 50 years. Figure 5a-h shows the spatial distribution of changing trends in the 10-year return levels for the different precipitation indices. The 10 years return levels indicate the maximum magnitudes of the indices at 90 % occurrence probability. Trends in 10 years return levels can be interpreted as trends in (1/10) % extremes. R95p 10 denotes the 10-year return level for the R95p index, with similar notation used for the other indices. R1mm 10 , R20mm 10 , R10mm 10 , CWD 10 , and RX5day 10 showed similar spatial characteristics, with values that gradually decreased from the southeast to the northwest. The opposite spatial pattern was seen for CDD 10 and R95pTOT 10 , which is consistent with the subarid climate in the northwest region and subhumid climate in the southeast region. R95p 10 did not show obvious regional characteristics across the plateau (Fig. 5a), with most values being between 22 and 30 days.  Fig. 5i-p. For R95p 10 , statistically significant negative trends (p < 0.05) over the past half century were found over almost all regions on the plateau, with positive trends identified in just a handful of grid locations in the southeast region (Fig. 5i). For R95pTOT 10 , the trends were positive in the southeast and middle part of the Loess Plateau, and slightly negative in the northern regions. For R1mm 10 , most grid locations showed negative trends, with rates of change of about −0.4 to −0.2 days/years. Only a few scattered grid locations in the northern, southwest, and southeast of the Loess Plateau showed positive trends, and all but one of these were nonsignificant. Trends for R20mm 10 and R10mm 10 were negative in the central part of the Loess Plateau and positive in the southeast region. CDD 10 exhibited significant positive trends across large areas of the plateau, except in the southwest region. In the southeast region, values for the positive trends were up to 1.0 days/years. For CWD, the 10-year return levels did not show zonal characteristics and the grid locations with positive and negative trends were interspersed. For RX5day 10 , most grid locations had significantly negative trends, especially in the southwest regions, although some locations in the middle of the plateau showed positive trends (Fig. 5p).
The spatial distributions of the joint return periods for different pairs of precipitation indices are presented in Fig. 6a, c, e, g and i. Figure 6b, d, f, h and j show the corresponding trends for the bivariate return periods during the period 1961-2011. Overall, the joint return periods were greater than the return periods for single indices. Across large areas of the plateau, T(R95p 10 , R95pTOT 10 ) was approximately 20 years, with handful of grids scattered in the middle and northwest regions showing larger T(R95p 10 , R95pTOT 10 ) values, indicating a lower probability of concurrent high values for R95p and R95pTOT. Positive trends in T(R95p 10 , R95pTOT 10 ) were found in the northern regions, but trends were negative in the southeast regions. The probability of concurrent high values for R1mm and R10mm, and for R1mm and R20mm, were similar, with T(R1mm 10 , R10mm 10 ) and T(R1mm 10 , R20mm 10 ) both approximately 25 years in most grid locations. Only a few grid locations had larger T(R1mm 10 , R10mm 10 ) and T(R1mm 10 , R20mm 10 ), and hence lower probability of concurrent return. The spatial trend patterns were also similar for T(R1mm 10 , R10mm 10 ) and T(R1mm 10 , R20mm 10 ). The probability of concurrent high values for R1mm and R10mm, and for R1mm and R20mm, decreased during the period 1961-2011, with T(R1mm 10 , R10mm 10 ) and T(R1mm 10 , R20mm 10 ) showing positive trends in most grid Fig. 6 a, c, e, g, i The spatial distribution of bivariate return periods at 10-year return levels for five precipitation index pairs. b, d, f, h, j The corresponding temporal trends over the period 1961-2011. Black dots represent a significance level of p < 0.05. The T(X 10 , Y 10 ) represents the joint return period that the climate indices X and Y are larger than 10-year return levels locations, apart from a few in the southeast and middle of the plateau. The probability of concurrent long CDD and long CWD was low, as indicated by the large T(CDD 10 , CWD 10 ). This suggests that the risk of co-occurrence of droughts and floods within a single year is low on the Loess Plateau. T(CDD 10 , CWD 10 ) did not show any apparent trends areas across regions, but grid locations with positive and negative trends were interspersed across the plateau. Grid locations with a large T(RX5day 10 , CWD 10 ) were found in the southwest of the plateau, but smaller values were observed in most other regions. The probability of concurrent large RX5day and long CWD decreased over the period 1961-2011, with significant positive trends in T(RX5day 10 , CWD 10 ) in many grid locations, especially in the northeast region.

Trivariate copula application
On the Loess Plateau, rainfall is mainly composed of highintensity storms, the loess in this region is generally characterize by a loose texture, large porosity and high erodibility owing to the deposition forms and composition. All contribute to the Loess Plateau having some of the most serious erosion in the world (Shi and Shao 2000;Wang et al. 2012). Analysis of this frequent heavy rainfall is therefore essential. Using trivariate copulas, we assessed the joint behavior of the frequency (R95p), intensity (PI), and annual fraction (R95pTOT) of precipitation in very wet days. Figure 7a shows that the probability of concurrent high values for R95p, R95pTOT, and PI increased from the southeast to the northwest of the Loess Plateau, as indicated by the spatial distribution of T(R95p 10 , R95pTOT 10 , PI 10 ). T(R95p 10 , R95pTOT 10 , PI 10 ) showed positive trends across more than half of the area of the plateau, especially in the eastern and northwest regions, where the rates of change in return period per time shift (time shift ΔT equals to 1 year in this case) were greater than 6 years. However, negative trends were observed in the southeast, indicating an increase over the past 50 years in the probability of extreme precipitation events occurring with high frequency and high intensity, and forming a large percentage of the annual precipitation.

Joint analysis of precipitation and temperature indices
The concurrence of extreme precipitation events and extreme temperature events would strengthen the intensity of the extreme events. For example, high temperatures combined with reduced precipitation would increase the risk of drought. Therefore, we analyzed the bivariate return periods for temperature and precipitation indices. As shown in Fig. 8a, T(TX10p 10 , CWD 10 ) was close to 60 years for most grid locations, with a few scattered locations having even greater bivariate return periods. This indicates that the probability of concurrently high TX10p and long CWD was low. Most grid locations had rates of change in return period per time shift (time shift ΔT equals to 1 year in this case) greater than 1.0 year, indicating that the probability of concurrently high TX10p and long CWD had decreased over the period 1961-2011. A few grid locations located at the edge and middle of the plateau showed negative trends during the period 1961-2011. The probability of concurrently high TX90p and long CDD was similar across most of the plateau, with T(TX90p 10 , CDD 10 ) equaling approximately 56 years in most locations. A few grid locations in the southern regions showed greater values for T(TX90p 10 , CDD 10 ). Some grid locations in the eastern and western regions had increasing trends in T(TX90p 10 , CDD 10 ) over the past 50 years, indicating a decrease in the probability of concurrently high TX90p and long CDD. In the northern and southern regions, negative trends were found, indicating an increase in the probability of concurrently high TX90p and long CDD. In these grid locations, the risk of higher temperatures and longer spells of consecutive dry days increased over the period 1961-2011, which may have resulted in more severe droughts. Fig. 7 a The spatial distribution of trivariate return periods at 10-year return levels for the extreme precipitation indices R95p, R95pTOT, and PI. b The corresponding temporal trends over the period 1961-2011. Black dots represent a significance level of p < 0.05. The T(X 10 , Y 10 ) represents the joint return period that the climate indices X and Y are larger than 10-year return levels

Discussion and conclusions
We presented here a comprehensive analysis of the joint probabilistic characteristics and tendencies of bivariate and trivariate precipitation and temperature indices for the Loess Plateau, calculated using copula theory. For individual indices, the 10-year return levels generally represent the maximum index magnitudes at 90 % occurrence probability (Liu et al. 2014). For eight precipitation indices, the spatial distribution and trends in the 10-year return levels were investigated over the period 1961-2011. The 10-year return levels for the selected indices showed regional features. The southeast region of the Loess Plateau had a greater potential for severe floods, owing to higher R20mm 10 , R10mm 10 , CWD 10 , and RX5day 10 , whereas the northwest region had a greater potential for severe drought because of higher CDD 10 and lower R1mm 10 . The results showed that the value for CDD 10 increased significantly between 1961 and 2011 in large parts of the plateau, especially in the southeast region, indicating a potential trend towards more severe droughts. At same time, significant positive trends in R95p 10 , R20mm 10 , R10mm 10 , and RX5day 10 in some southeastern grid locations demonstrated a potential trend towards more severe floods at those locations. An increase in the severity of floods would to some extent also increase the risk of soil erosion.
For the joint analyses of precipitation indices, the joint return periods were greater than the return periods for single indices. Overall, large areas of the Loess Plateau showed a relatively high risk of concurrent extreme precipitation events, as shown by the relatively low values of 10-year return levels (T(RX5day 10 , CWD 10 ) and T(R1mm 10 , R20mm 10 )) across the plateau. This is mainly due to the inequality of the temporal distribution of precipitation across the year on the Loess Plateau. On the other hand, the probability of co-occurrence of drought and flooding within a single year was low. The risk of concurrent extreme wet and dry events had not increased over the past half century, as shown by the non-significant changes in T(CDD 10 , CWD 10 ) illustrated in Fig. 6h. Over the past 50 years, positive trends in T(R95p 10 , R95pTOT 10 ), T(R1mm 10 , R10mm 10 ), T(R1mm 10 , R20mm 10 ), and T(RX5day 10 , CWD 10 ) were observed in the northern part of the plateau, indicating that the probability of extreme precipitation events with a long duration or high intensity has decreased. Some grid locations in the southeast region showed an increased risk of extreme wet events, as indicated by negative trends in T(R95p 10 , R95pTOT 10 ), T(R1mm, R10mm), T(R1mm, R20mm), and T(RX5day 10 , CWD 10 ). Moreover, trivariate copula analysis for R95p, R95pTOT, and PI showed that these locations also had an increased risk of precipitation events that were even more extreme, with higher frequency and intensity, and forming a larger percentage of the annual precipitation. The Loess Plateau located in the fringe zone for the monsoon and water vapor mainly comes from the East Asian summer monsoon, which result in precipitation change in the area is very sensitive to the change of the East Asian monsoon (Zhou et al. 2010;Wang et al. 2012). In the context of global warming, the monsoon in China has been Fig. 8 a, c The spatial distribution of bivariate return periods at 10-year return levels for joint precipitation and temperature indices. b, d The corresponding temporal trends over the period 1961-2011. Black dots represent a significance level of p < 0.05. The T(X 10 , Y 10 ) represents the joint return period that the climate indices X and Y are larger than 10-year return levels relatively weak since the 1970s, which would explain the positive trends in T(R95p10, R95pTOT10), T(R1mm10, R10mm10), T(R1mm10, R20mm10) in larger part of Loess Plateau to some extent. However, in some girds scattering in the southeast of Loess Plateau, the probability of occurrence of more and stronger extreme precipitation events increased. The results indicate the frequency and intensity of rainstorm were not impacted significant owing to weak monsoon. Some meso-scale convections circulations are also important factors affecting the precipitation formations. In these regions, the small scale, short duration, high intensity precipitation caused by local strong convective conditions is the main form of the rainstorm (Jiao et al. 1999). The convective precipitation responds much more sensitively to temperature increases and increasingly dominates events of extreme precipitation under global warming (Berg et al. 2013). Besides, the Loess plateau is characterized by complex terrain; high gully density and valleys are formed owing to tectonic movement and soil erosion (Shi and Shao 2000). The complex topography would affects the microclimate significantly and result in the spatial heterogeneity of the changing trends in precipitation extremes, such as the existence of hot-day valley (Ghosh et al. 2012).
The Loess Plateau has experienced significant warming over the past 50 years (Wang et al. 2012). Climate warming and the associated rise in extreme temperatures substantially increase the chance of concurrent droughts and heat waves (AghaKouchak et al. 2014;Fischer and Knutti 2015). The most reliable assessment of drought severity and risk is achieved by taking both extreme temperature and extreme precipitation indices into account. We used the copula technique to investigate the combination of extreme temperatures and extreme precipitation events. Our results indicate that the risk of concurrent higher temperatures and longer spells of consecutive dry days was increased in grid locations scattered in the northern and southern regions. This could trigger more severe droughts or heat waves. Concurrent longer spells of consecutive wet days and greater numbers of cold days tended to be less frequent, which may be attributable to the observed warming. Concurrent drought and warmer temperatures would result in an adverse impact on crop yields because the area available to farm would likely reduce owing to the potential reduction in water availability. Moreover, other policies such as afforestation and reforestation efforts would be affected because climate conditions dominate the success of vegetation coverage (Fan et al. 2014a, b). Regional climate changes inevitably influence the growth of plants.
The Loess Plateau is particularly sensitive to climate change owing to its fragile ecological environment and geographic features. Previous studies have reported no significant changes in precipitation extremes in the region (Li et al. 2010;Wang et al. 2012). However, by focusing on probabilistic characteristics, we have demonstrated significant changes in some precipitation extremes for some regions of the plateau. Moreover, from the joint analysis of precipitation and temperature extremes, the risk of severe floods or droughts over the Loess Plateau has increased. The results also showed a non-uniformity of changing trend in spatial scales. In a warming climate, the concurrent decrease of different class precipitation events (R1mm, R10mm, and R20mm) and increase risk of compound extreme warm days and longer dry days would trigger more severe droughts in larger northern region in Loess Plateau. The increasing anabatic evapotranspiration owing to warming and water deficiency would causes serious soil desiccation and the development of dried soil layer (Wang et al. 2010b). In the southeast regions, though vegetation coverage had increased with the afforestation and reforestation efforts in recent years, the higher probability of occurrence of more and stronger extreme precipitation events still posed a great threat to soil and water loss. Thereby, the soil and water conservation should be still a priority for all in Loess Plateau. Management of soil erosion, water availability, and agriculture on the Loess Plateau continues to face many challenges. A comprehensive knowledge of the trends and probabilistic characteristics of precipitation and temperature extremes and how these relate to future changes in climate is essential for the management and mitigation of natural hazards over the Loess Plateau.