Multiscale Spatio-Temporal Changes of Precipitation Extremes in Beijing-Tianjin-Hebei Region, China during 1958–2017

: In this study, based on daily precipitation records during 1958–2017 from 28 meteorological stations in the Beijing-Tianjin-Hebei (BTH) region, the spatio-temporal variations in precipitation extremes deﬁned by twelve indices are analyzed by the methods of linear regression, Mann-Kendall test and continuous wavelet transform. The results showed that the spatial patterns of all the indices except for consecutive dry days (CDD) and consecutive wet days (CWD) were similar to that of annual total precipitation with the high values in the east and the low value in the west. Regionally averaged precipitation extremes were characterized by decreasing trends, of which ﬁve indices (i.e., very heavy precipitation days (R50), very wet precipitation (R95p), extreme wet precipitation (R99p), max one-day precipitation (R × 1day), and max ﬁve-day precipitation (R × 5day)) exhibited signiﬁcantly decreasing trends at 5% level. From monthly and seasonal scale, almost all of the highest values in R × 1day and R × 5day occurred in summer, especially in July and August due to the impacts of East Asian monsoon climate on inter-annual uneven distribution of precipitation. The signiﬁcant decreasing trends in annual R × 1day and R × 5day were mainly caused by the signiﬁcant descend in summer. Besides, the possible associations between precipitation extremes and large-scale climate anomalies (e.g., ENSO (El Niño Southern Oscillation), NAO (North Atlantic Oscillation), IOD (Indian Ocean Dipole), and PDO (Paciﬁc Decadal Oscillation)) were also investigated using the correlation analysis. The results showed that the precipitation extremes were signiﬁcantly inﬂuenced by ENSO with one-year ahead, and the converse correlations between the precipitation extremes and climate indices with one-year ahead and 0-year ahead were observed. Moreover, all the indices show signiﬁcant two- to four-year periodic oscillation during the entire period of 1958–2017, and most of indices show signiﬁcant four- to eight-year periodic oscillation during certain periods. The inﬂuences of climate anomalies on precipitation extremes were composed by di ﬀ erent periodic components, with most of higher correlations occurring in low-frequency components.


Introduction
With the global hydrological cycle accelerating subsequently in the background of climate warming [1], extreme climate events have become more frequent [2], which has a significant impact on human society and even causes serious losses to people's lives and property [3][4][5]. Correspondingly, extreme events and related societal impacts are becoming an increasingly interesting area and has also attracted great attention, especially in the field of meteorology [6]. The Expert Team on Climate and found most extreme precipitation indices showed an insignificant decreasing trend with exception in R × 1day and R × 5day. Shen et al. [38] summarized the spatial and temporal variations of extreme precipitation over the BTH region using the percentile method and 174 station data during 1981-2016, and they found the spatial patterns of extreme precipitation indices except for extreme precipitation days are similar and extreme precipitation events frequently occurred in July, August, and October. Additionally, for the neighboring regions of BTH, Du et al. [39] reported that alterations in extreme precipitation mainly occurred north of 38 • N, and decreasing trends in extreme precipitation were detected at most stations in the Haihe River basin. Zhou et al. [40] reported that all cities in the BTH region showed decline trends in extreme precipitation during 1960-2014, with a significant decreasing trend in megacities based on linear regression method.
Additionally, variations of precipitation extremes are caused by many forcing factors, such as large-scale atmospheric circulation, topography, and local climate [14,32,[41][42][43][44][45][46][47][48][49]. For example, Xiao et al. [41] analyzed the teleconnections between precipitation extremes and El Niño Southern Oscillation (ENSO), North Atlantic Oscillation (NAO), Indian Ocean Dipole (IOD), and Pacific Decadal Oscillation (PDO) using the regression method, and found that regional responses of precipitation extremes to climate indices varied at different stages. Chen et al. [42] investigated the possible association between temporal patterns of precipitation extreme over China with large-scale ocean-atmospheric indices (e.g., NAO, PDO, Arctic Oscillation (AO), NINO3 (Niño3 sea surface temperature (SST) index), NINO3.4 (Niño3.4 SST index), and other SST indices) based on grid daily precipitation data during 1963-2013, and found that precipitation extremes in four different seasons were influenced by the different combinations of large-scale climate indices with different time lags, and their regional relations are complex when the climate indices were at different phases. Xue et al. [10] discussed the long-term trends of precipitation extremes and underlying mechanisms in U.S. Great Basin during 1951-2013, and found that the effects of ENSO and PDO on precipitation and precipitation extremes were low, while the intensifying high-pressure system over the western North Pacific in winter and increasing water vapor and air temperature of the atmosphere over the basin in winter and spring were mainly responsible for the upward trends of precipitation extremes. With advances in the field of synoptic climatology, it has been found that ocean-atmosphere interactions are not chaotic or random, it is therefore possible to identify the relationship between ocean-atmospheric modes and hydrological variables. Related studies revealed that the ENSO, PDO, and other large-scale atmospheric circulation are likely to be important factors in the variability of precipitation extremes on both spatial and temporal scales. And these climate anomalies will be helpful for prediction/forecasting of precipitation extremes. However, few efforts have been made on the possible linkage between precipitation extremes over BTH and the large-scale atmospheric circulation. Thus, to better understanding the changes of precipitation extremes, the associations between regional precipitation extremes and large-scale climate anomalies should be further discussed.
Although many researchers have conducted analyses regarding the changes of precipitation extremes, the long-term trends with different scales over the BTH region and the possible teleconnection with large-scale climate anomalies still need to be strengthened. Additionally, assessing the temporal trends and spatial patterns of precipitation extremes in BTH will provide information for managing water resources effectively and better responding the extreme events. Therefore, the main objectives of our study are (1) to investigate the spatial-temporal variability of extreme precipitation at different scales over the BTH region using several trend detecting and periodicity analysis methods based on the 28 meteorological stations during 1958-2017, and (2) to identify the dominant climate factors in explaining extreme precipitation changes by investigating the statistical association between temporal patterns of extreme precipitation signals with large-scale ocean-atmosphere indices using correlation analysis.

Study Area
The study area (Figure 1a), which is located in 113 • 04 -119 • 53 E, 36 • 01 -42 • 37 N, has a total area of about 217,185 km 2 . It includes Beijing municipality, Tianjin municipality, and other 11 cities in Hebei Province (the BTH region), also known as Jing-jin-ji Metropolitan Region or Jing-Jin-Ji (JJJ). As shown in Figure 1c, the spatial pattern of elevation exhibits the obvious feature of "high in the northwest and low in the southeast". The eastern and southern part of BTH are adjacent to the Yellow River and Bohai Sea, with elevation less than 500 m. Taihang and Yan mountains are in western and northern part of the BTH region, which constitute 53.6% of the total area. The monsoon climate has obvious seasonal variations characterized by hot, rainy summers, and cold, dry winters. Meteorological disasters in this region are frequent and have caused serious damages to human life especially in recent years. Drought often strikes in the spring and flooding occurs in the summer. From the annual precipitation totals, BTH as a whole is commonly regarded as a semi-humid region, with a mean annual total precipitation of 540 mm. Figure 1b showed the spatial distribution of annual mean precipitation over BTH during 1958-2017. The annual mean precipitation increased from the west to the east in the BTH region, ranging from 393 mm to 723.4 mm, with the highest value (Zunhua station (No. 13 in Figure 1)) located in the belt between the mountain and plain areas. especially in recent years. Drought often strikes in the spring and flooding occurs in the summer. From the annual precipitation totals, BTH as a whole is commonly regarded as a semi-humid region, with a mean annual total precipitation of 540 mm. Figure 1b showed the spatial distribution of annual mean precipitation over BTH during 1958-2017. The annual mean precipitation increased from the west to the east in the BTH region, ranging from 393 mm to 723.4 mm, with the highest value (Zunhua station (No. 13 in Figure 1)) located in the belt between the mountain and plain areas.

Precipitation Data
The daily precipitation dataset [51] of 28 meteorological stations (Table 1 and Figure 1) provided by the China Meteorological Administration (CMA) was used to extract time series of extreme precipitation. The weather stations have a reasonable spatial coverage, consequently being able to emphasize the regional behavior and variability of the precipitation extremes in the BTH region. The geographical coordinates of the stations are shown in Table 1. Before constructing the precipitation indices series, data quality control and homogeneity assessment were performed on the initial data set of daily precipitation amounts from each station. Based on the daily data, we calculated the monthly, seasonal and annual precipitation to do the trend analysis.  Table 1. (b) Annual mean precipitation spatial distribution during the period of 1958-2017. The trends and slopes were calculated by the linear regression method, with the significant changes filled red five stars at 0.05 level. (c) Digital elevation model (DEM) data in the study area. DEM data set [50] is provided by Geospatial Data Cloud site, Computer Network Information Center, and the Chinese Academy of Sciences.

Precipitation Data
The daily precipitation dataset [51] of 28 meteorological stations (Table 1 and Figure 1) provided by the China Meteorological Administration (CMA) was used to extract time series of extreme precipitation. The weather stations have a reasonable spatial coverage, consequently being able to emphasize the regional behavior and variability of the precipitation extremes in the BTH region. The geographical coordinates of the stations are shown in Table 1. Before constructing the precipitation indices series, data quality control and homogeneity assessment were performed on the initial data set of daily precipitation amounts from each station. Based on the daily data, we calculated the monthly, seasonal and annual precipitation to do the trend analysis.

Climate Indices
In this study, ENSO, NAO, IOD, and PDO were chosen to analyze their effects on precipitation extremes in the BTH region. ENSO is indicated by the NINO3.4 index, derived from sea surface temperature anomaly estimated in the Niño 3.4 region (5 • N-5 • S, 120 • -170 • W) [52]. NAO is defined as a meridional dipole in the atmospheric pressure with centers of actions near the Azores and Iceland [53]. A higher NAO index is associated with westerlies stronger than average in the mid-latitude North Atlantic, while the opposite is true for negative values of the index [54]. IOD is defined as the SST anomaly difference between the western equatorial Indian Ocean (50 • -70 • E, 10 • S-10 • N) and the south eastern equatorial Indian Ocean (90 • -110 • E, 10 • S-0 • N), referred to as Dipole Mode Index [55]. PDO is most frequently referred to as a long-lived El Niño like pattern of the Pacific climate variability which is considered as the leading principle component of monthly sea surface temperature anomalies in the North Pacific Ocean, poleward of 20 • N [56]. In this study, annual NINO3.4 index is calculated based on the average value from July to December, and the NAO is estimated from December to March [41]. The annual PDO and IOD indices are averaged from January to December based on the monthly series. Four seasons such as spring (March-May), summer (June-August), autumn (September-November), and winter (December-February) were also analyzed in the study. The seasonal climate indices are averaged from the monthly series.

Definition of Extreme Precipitation Indices
The ETCCDI defined 27 core extreme indices based on daily temperature and precipitation. The RClimDex [57] developed by Xuebin Zhang and Yang Feng at Climate Research Division, was used for data quality control and extreme indices calculation. In this study, we employed a set of 12 indices (listed in Table 2) related to extreme precipitation, including two user-defined thresholds (R25 and R50) according to the Standard for Grade Classification of Precipitation Intensity of China released by the China Meteorological Administration (Light rain: Daily precipitation <10 mm; Moderate rain: 10 mm≤ daily precipitation <25 mm; Heavy rain: 25 mm≤ daily precipitation <50 mm; Torrential rain: 50 mm ≤ daily precipitation <100 mm; Extraordinary rainstorm: Daily precipitation ≥100 mm). In this study, the baseline period of 1961-1990 was used for the estimation of threshold values for the percentile-based indices. These indices can reflect the changes of extreme precipitation in different aspects. Some of these indicators alone or together with other indices were previously widely used to assess climate changes of the extreme precipitations in different regions of the world [1,11,33,41].

R10
Annual number of days with more than 10 mm/day days R20 Annual number of days with more than 20 mm/day days R25 * Annual number of days with more than 25 mm/day days R50 * Annual number of days with more than 50 mm/day days CDD Maximum number of consecutive dry days 1 days CWD Maximum number of consecutive wet days 2 days R95p Annual total precipitation when daily precipitation >95th percentile mm R99p Annual total precipitation when daily precipitation >99th percentile mm R × 1day Annual, seasonal and monthly maximum one-day precipitation mm R × 5day Annual, seasonal and monthly maximum five-days precipitation mm SDII Annual total precipitation divided by the number of wet days in the year mm/d PRCPTOT Annual total amount of precipitation cumulated in wet days mm * 25 mm and 50 mm were the thresholds defined by the authors. 1 Dry days are those days when daily precipitation was < 1 mm; 2 Wet days are those days when daily precipitation was ≥ 1 mm.

Statistical Analysis
In this study, temporal trends of indices for each station were determined using the least-squares linear regression method and Mann-Kendall test, which were widely used in the meteorological and hydrological fields [10,30,37,58]. Furthermore, the Spearman's rank correlation and two-tailed t-test were used to test the significance of the relation between precipitation extremes and annual total precipitation and climate indices in Sections 3.3 and 3.4.
Linear regression is a parametric method used to obtain the slope of hydro-meteorological variables over time, which can be represented as The slope b can be used as an indicator of trend and is calculated as where y i is a climatic factor (e.g., precipitation extremes in this study), x i is time, and n is the length of the time series. The linear trends are also considered to be statistically significant on the 0.05 level using Pearson correlation analysis. For a given time series Y = (y 1 , y 2 , . . . , y n ), the Mann-Kendall test statistic S is defined by where y i and y j are the sequential data values, and the function sgn(y) is defined as The statistic S is approximately normally distributed with the mean E(S) = 0 and variance as where t i is considered as the number of ties up to sample i. The standardized normal test statistic Z is given by A positive (negative) value of Z signifies an upward (downward) trend [44]. In this study, we chose a significance level of 0.05 for the Mann-Kendall test.

Wavelet Analysis
Wavelet analysis was widely applied to examine multi-temporal scale features of precipitation extremes or other climate events [59][60][61][62][63][64][65]. Wavelet transform (WT) decomposes time series into time-frequency space and identifies the dominant modes of variability and the characteristics how these modes vary in time [66]. Continuous wavelet transform (CWT) is appropriate for extracting a wide range of possible dominant frequencies from geophysical and hydroclimate time series [67,68]. In this study, CWT with Morlet wavelet was used to obtain information about the periodic structure of precipitation extremes, which can be defined as [68] ψ(η) = π −1/4 e iωη−0.5η 2 (7) where ω is the dimensionless frequency, and η is the dimensionless time parameter. The wavelet is stretched in time (t) by varying its scale (s), so that η = s/t. It is possible to obtain a picture depicting the change of amplitude with scale and its variation with time by varying the wavelet scale s and translating along the localized time index η.
For a given wavelet ψ 0 (η), it was assumed that Yj is a time series of length N (Yj, j = 1, 2, . . . , N) with equal time spacing ∆t. The continuous wavelet transform of a discrete sequence Yj is defined as a convolution of Yj with the scaled and translated wavelet ψ 0 (η): where asterisk (*) indicates the complex conjugate.

Spatial Variation of Extreme Precipitation Events
The spatial distribution of precipitation extremes over the BTH region from 1958 to 2017 was shown in Figure 2. In general, the spatial distributions of precipitation extremes apart from consecutive dry days (CDD) and consecutive wet days (CWD) were similar to that of annual precipitation in BTH, with the highest values in the northeast and the lowest in the northwest (Figure 1b). In terms of PRCPTOT, the highest value (709.8 mm in Zunhua station (No. 13)) was observed in the northeast while the lowest value (375.1 mm in Zhangbei station (No. 1)) was found in the northwest of BTH. Similarly, the SDII varied from 7.2 mm/day (Zhangbei station (No.1)) to 13.8 mm/day (Zunhua station (No. 13)), with the high values (12-14 mm/day) occurring on the east coast around the Bohai Bay. Spatial patterns in R95p and R99p were similar, with the largest being observed in the northeast and the lowest in the northwest. The indices in R × 1day and R × 5day increased from the northwest to the southeast, with the highest value of 105.8 mm and 153 mm in Zhangbei station (No. 1). However, the low value (43.5 mm) in R × 1day also occurred in Qinhuangdao station, which was only about half value of the nearest stations (80-90 mm). On the other hand, the value of R × 5day in Qinhuangdao Hebei province) were larger than 80 days. The spatial pattern of CWD was most unique, with the high value in the north and south and the low value in the middle of the region. Specially, the lowest value of CWD occurred in the eastern coast around the Bohai Bay with the value ranging from 3.4 to 3.7 days. In terms of R10, R20, R25, and R50, the spatial patterns were similar, with a decrease from the east to the west. Overall, the spatial patterns of precipitation extremes found in this study agree well with the findings from Mei et al. [33]. The landforms in the northern and western parts of BTH are of a mountainous nature, which act as orographic barriers for the advancing air masses and limit precipitation distribution [69].

Spatial Variation of Extreme Precipitation Events
The spatial distribution of precipitation extremes over the BTH region from 1958 to 2017 was shown in Figure 2. In general, the spatial distributions of precipitation extremes apart from consecutive dry days (CDD) and consecutive wet days (CWD) were similar to that of annual precipitation in BTH, with the highest values in the northeast and the lowest in the northwest (Figure 1b). In terms of PRCPTOT, the highest value (709.8 mm in Zunhua station (No. 13)) was observed in the northeast while the lowest value (375.1 mm in Zhangbei station (No. 1)) was found in the northwest of BTH. Similarly, the SDII varied from 7.2 mm/day (Zhangbei station (No.1)) to 13.8 mm/day (Zunhua station (No. 13)), with the high values (12-14 mm/day) occurring on the east coast around the Bohai Bay. Spatial patterns in R95p and R99p were similar, with the largest being observed in the northeast and the lowest in the northwest. The indices in R × 1day and R × 5day increased from the northwest to the southeast, with the highest value of 105.8mm and 153mm in Zhangbei station (No. 1). However, the low value (43.5 mm) in R × 1day also occurred in Qinhuangdao station, which was only about half value of the nearest stations (80-90 mm). On the other hand, the value of R × 5day in Qinhuangdao (143 mm) match with its surroundings (130-150 mm). The mean values of CDD were greater in the northwest of BTH than those in the southeast, with the highest value of 94. Specially, the lowest value of CWD occurred in the eastern coast around the Bohai Bay with the value ranging from 3.4 to 3.7 days. In terms of R10, R20, R25, and R50, the spatial patterns were similar, with a decrease from the east to the west. Overall, the spatial patterns of precipitation extremes found in this study agree well with the findings from Mei et al. [33]. The landforms in the northern and western parts of BTH are of a mountainous nature, which act as orographic barriers for the advancing air masses and limit precipitation distribution [69].

Annual Scale Trends in Extreme Precipitation Indices
To analyze overall changes in extreme precipitation of BTH in 1958-2017, Figure 3 shows the average annual time series of BTH for all indices, while Table 3 shows the number of individual stations detected with positive or negative trends for the indices analyzed in the study. The spatial patterns of these results are shown in Figures 4 and 5 based on the Mann-Kendall and linear regression method, respectively. Overall, in these indices, negative trends dominate over positive trends, which mean that extreme precipitation over BTH shows more decreasing than increasing trends in the study period.   Table 3. Results of the statistical tests for regional extreme precipitation indices and the number of stations with positive or negative trends over the BTH region from 1958 to 2017.    and NSI represent non-significant decreasing and increasing trend, respectively. NT means no any change trend.  Table 3. Results of the statistical tests for regional extreme precipitation indices and the number of stations with positive or negative trends over the BTH region from 1958 to 2017. Over 1958-2017, all the indices showed a decreasing trend based on the linear regression method, while only five indices (R50, R95p, R99p, R × 1day, and R × 5day) showed a statistically significant decreasing trend (p ≤ 0.05) during this period, at a rate of 0.1 day/decade, 8.5 mm/decade, 4.78 mm/decade, 2.65 mm/decade, and 5.05 mm/decade, respectively. As for the average time series of PRCPTOT, it exhibited a non-significant (p = 0.183) decreasing trend at a rate of 10.27 mm/decade, whereas the average SDII time series of BTH had decreased by 0.09 mm/d per decade. The CDD, ranging from 45 days to 119 days with a mean value of 77.8 days, had a decreasing trend at a rate of 1.43 days/decade from 1958 to 2017. The decreasing trend (p = 0.06) for CWD was more significant than that for CDD (p = 0.337) although the decreasing rate of CWD is relatively lower (0.07 day/decade). The regionally averaged occurrence of heavy precipitation days (R10) and heavier precipitation days (R20 and R25) have decreasing trends, nonetheless, they were not statistically significant. The decreasing rates in R10, R20, and R25 were 0.09, 0.12, and 0.13 day/decade, respectively. Similar trends can be also concluded from the Mann-Kendall trend test, as shown in Table 3. All the indices except R10 had a decreasing trend, but also only five indices (i.e., R50, R95p, R99p, R × 1day, and R × 5day) showed a statistically significant decreasing at the 0.05 level.

Mann-Kendall Trend Test Linear Regression Trend
Overall, the preponderance of evidence indicated that the decline trends in indices were more frequent than the increasing ones although considerably few decreasing trends exceeding 95% confidence level during 1958-2017. However, there are slight discrepancies in trends between the Mann-Kendall test and the linear regression. Compared with Mann-Kendall test, the greater number of stations (34 versus 23) with significantly decreasing trends are found based on linear regression. Additionally, based on the Mann-Kendall test, about 53.6% (46.4%) of the stations showed upward (downward) trends for R10, while none of these trends are statistically significant. On the contrary, most of stations (71.4%) for R10 showed negative trends based on the result of linear regression method. For the other indices, the similar trends can be found based on the two methods. Generally, although most of locations exhibited downward trends for the precipitation extremes, less than 25% of them showed significant decreasing trends in the BTH region during 1958-2017.

Monthly and Seasonal Changes in Precipitation Extremes
In this study, R × 1day and R × 5day were used to analyze the monthly and seasonal changes in precipitation extremes, as shown in Figure 6. It was found that the monthly mean values in R × 1day and R × 5day during the period of 1958-2017 varied from 1.6 to 56.1 mm, and 2.4 to 88.7 mm, respectively. The wettest season in BTH is summer (June-July-August) due to effects of eastern Asian Monsoon climate, of which the precipitation amount accounts for about 70% of annual totals [67]. Obviously, almost all of the highest values in R × 1day and R × 5day occurred in summer (both about 98.3% of totals), especially in July (37 and 36, respectively) and August (19 and 22, respectively). And the lowest values frequently occurred in winter (December-January-February). On average, the highest values of R × 1day and R × 5day occurred in July (28.6~123.8 mm and 41~182.3 mm, respectively) and the lowest in January (0.02~9.7 mm and 0.1~14.1 mm). The trends of monthly and seasonal time series were also detected by the Mann-Kendall test, as shown in Figure 6e. In terms of R × 1day, a significant increasing trend was found in May while a significant decreasing trend occurred in July and August, which may result the trend of significant upward in spring and significant downward in summer. Other monthly and seasonal series also had increasing or decreasing trends but are not significant at the level of 0.05. Similarly, a significant increasing (decreasing) trend in May (July and summer) was also found in R × 5day. However, a significant upward trend for R × 5day was also detected in June. Overall, the decreasing trend in summer (especially in July and August) may be the cause the decreasing trend in annual series. Furthermore, the statistical results of trends for all the stations are shown in Figure 6f. An increasing trend was dominant in spring (100% of locations), autumn (75% and 71.4%) and winter (96.4% and 89.3%) for R × 1day and R × 5day, respectively, while about 85.7% of locations showed a downward trend in summer for both indices. Certainly, the number of significantly increasing or decreasing trends (SI or SD) was less than that of non-significantly increasing or decreasing trends (NSI or NSD) for all the seasons. Analogously, the monthly series can also be classified by two groups except February and March. One was dominated by increasing trends, including April, May, June, September, October, and December, and the other was leaded by decreasing trends, including January, July, August, and November. In February, the stations with increasing or decreasing trends both accounted for about half of the totals in R × 1day (50%) and R × 5day (53.6% and 46.4%). Interestingly, in March, more stations showed decline trends in R × 1day, whereas the opposite results were found in R × 5day. In conclusion, the decline trends in annual series may be caused by significant decreasing trends in summer, especially in July and August. Although significant increasing trends were also found in several months or seasons, their effects on the trends in annual series were less than that in summer, because most of values in the annual series occurred in summer.
Atmosphere 2019, 10, x FOR PEER R EVIEW 13 of 30 Figure 6. Time series of monthly (a,b) and seasonal (c,d) precipitation extremes in R×1day (a,c) and R × 5day (b,d) in BTH during the period of 1958-2017. The regional trends (e) and statistical results for all the stations (f) in monthly and seasonal precipitation extremes using the Mann-Kendall test. NT means no trends, SD and SI mean significantly decreasing or increasing trends, and NSD and NSI mean non-significantly decreasing or increasing trends, respectively.

Decadal Variability in Precipitation Extremes
We  Figure 7. The detailed information of precipitation extremes in each subperiod can be found in Table S1. It was found that all the maximum values of decadal mean indices for sub-period occurred in T1, except for CDD (T2) and R10 (T6), while all the minimum values except CDD (T3) occurred in T5. The decadal changes in CDD suggested that the BTH appeared more  monthly (a,b) and seasonal (c,d) precipitation extremes in R × 1day (a,c) and R × 5day (b,d) in BTH during the period of 1958-2017. The regional trends (e) and statistical results for all the stations (f) in monthly and seasonal precipitation extremes using the Mann-Kendall test. NT means no trends, SD and SI mean significantly decreasing or increasing trends, and NSD and NSI mean non-significantly decreasing or increasing trends, respectively.

Decadal Variability in Precipitation Extremes
We  Figure 7. The detailed information of precipitation extremes in each sub-period can be found in Table S1. It was found that all the maximum values of decadal mean indices for sub-period occurred in T1, except for CDD (T2) and R10 (T6), while all the minimum values except CDD (T3) occurred in T5. The decadal changes in CDD suggested that the BTH appeared more dryness in T2 because of the higher consequence dry days, and the discrepancies ranged from 9.5 to 19.3 days (about 12.2-24.8% of mean value in CDD during 1958-2017) compared with the other sub-periods. The results from the other indices suggested that the extreme precipitation events or flood events may be more acute in T1 over BTH, whereas the relatively less probability or weights of extreme precipitation events occurred in T5. Figure 7 also showed that there were greater variations in R50, R95p, and R99p compared with other indices, which mean the higher spatial variations. Additionally, it was worth noting that all the indices except CDD experienced clear fluctuations among these sub-periods with a complex trend of decreasing (from T1 to T3)-increasing (from T3 to T4)-decreasing (from T4 to T5)-increasing (from T5 to T6).   1958-1967, 1968-1977, 1978-1987, 1988-1997, 1998-2007 In order to further analyze inter-decadal changes for precipitation extremes, Figure 8 demonstrated decadal changes for all the extreme precipitation indices. In general, the spatial patterns for all the indices in the six sub-periods were similar with slight variations. That is, the spatial distributions in each sub-period had no any significant change, which agreed well with that in the entire period of 1958-2017. As mentioned previously, the decadal changes in CDD were different from the other indices. About 42.9% of stations in the CDD in T2 was larger than 90 days, covering about 54.5% of BTH, while only one station in T3 was larger than 90 days. In T3 sub-period, the lowest value of the CDD located in the east of BTH, less than 60 days. In terms of CWD, in T2, T5 and T6 sub-periods, the lower values occurred in the middle of BTH, while the higher values located in the southern and northern areas of Hebei province. However, in T1 and T4, the lowest value of CWD located in the eastern part of the BTH region. In T1 sub-period, about 35.7% of locations in PRCPTOT was larger than 620 mm, mainly located in the eastern BTH and the south of Hebei province. However, none of them was higher than 600 mm in T5, and PRCPTOT for all the stations except Zunhua (No. 13) and Qinglong (No. 14) stations (located in the eastern BTH) was less than 550 mm. The similar decadal variations were found in R × 1day and R × 5day. The decadal averages of R × 1day and R × 5day had a same pattern of "high value in the southeast and low value in the northwest" for different sub-periods. The similar decadal variations and spatial patterns can be obtained for the indices of precipitation days (e.g., R10, R20, R25, and R50). In terms of R95p and R99p, the spatial patterns are different, with the lowest value occurring in the middle of the BTH region in T5 and in the western mountain areas during the other sub-periods. For the index of SDII, there was no difference among the six sub-periods, with the same spatial patterns of "high value in the southeast,  1958-1967, 1968-1977, 1978-1987, 1988-1997, 1998-2007, and 2008-2017). The upper and lower limits of the box indicate the 75th and 25th percentile value among these stations. The horizontal line in the box indicates the median and the whiskers show the range of these stations. Black solid squares are the mean values of all the stations, connected by the red lines. Open circles represent the outliers.
In order to further analyze inter-decadal changes for precipitation extremes, Figure 8 demonstrated decadal changes for all the extreme precipitation indices. In general, the spatial patterns for all the indices in the six sub-periods were similar with slight variations. That is, the spatial distributions in each sub-period had no any significant change, which agreed well with that in the entire period of 1958-2017. As mentioned previously, the decadal changes in CDD were different from the other indices. About 42.9% of stations in the CDD in T2 was larger than 90 days, covering about 54.5% of BTH, while only one station in T3 was larger than 90 days. In T3 sub-period, the lowest value of the CDD located in the east of BTH, less than 60 days. In terms of CWD, in T2, T5 and T6 sub-periods, the lower values occurred in the middle of BTH, while the higher values located in the southern and northern areas of Hebei province. However, in T1 and T4, the lowest value of CWD located in the eastern part of the BTH region. In T1 sub-period, about 35.7% of locations in PRCPTOT was larger than 620 mm, mainly located in the eastern BTH and the south of Hebei province. However, none of them was higher than 600 mm in T5, and PRCPTOT for all the stations except Zunhua (No. 13) and Qinglong (No. 14) stations (located in the eastern BTH) was less than 550 mm. The similar decadal variations were found in R × 1day and R × 5day. The decadal averages of R × 1day and R × 5day had a same pattern of "high value in the southeast and low value in the northwest" for different sub-periods. The similar decadal variations and spatial patterns can be obtained for the indices of precipitation days (e.g., R10, R20, R25, and R50). In terms of R95p and R99p, the spatial patterns are different, with the lowest value occurring in the middle of the BTH region in T5 and in the western mountain areas during the other sub-periods. For the index of SDII, there was no difference among the six sub-periods, with the same spatial patterns of "high value in the southeast, and low value in the northwest".  1958-1967, 1968-1977, 1978-1987, 1988-1997, 1998-2007, and 2008-2017). The decadal precipitation extremes were computed as the average value over each 10-year period.

Periodicity Analysis of Extreme Precipitation Indices
The wavelet power spectra of the regional average series of precipitation extremes over the BTH region were shown in Figure 9. We can observe significant periodic oscillations in all indices. For all the indices, the significant two-to five-year periods at the 5% level were predominant during 1958-2017. Furthermore, the CDD wavelet power spectrum showed 6-10-year modulations of variation, which occurred in the 1970s to the mid-1980s. Similarly, the significant inter-decadal (approximately 8-11 years) oscillations were active from the 1970s to the 1980s. For SDII, the power was broadly distributed, with a period of 2-11 years from the 1990s to 2000, and with a period of 8-11 years in the 1980s. Additionally, continuous wavelet spectrum for R95p also showed the significant wavelet power at 8-10 years during 1982-1993. The global wavelet power spectrum (GWPS) is an equalweighted average of all the local wavelet power spectra for each scale, which shows dominant scales with no temporal transformation ( Figure 9). Regional averaged GWPS shows a statistically significant low-frequency oscillation with a period of 2-9 years in CDD, a period of 2-12 years in CWD, a period of 2-8 years in R10, R20, R25, and RPCPTOT, a period of 2-10 years in SDII, a period of 2-14 years in R95p, a period of 2-15 years in R50 and R99p, and a period of 2-16 years in R × 1day and R × 5day. According to the local wavelet spectra (Figure 9), inter-annual (2-4 years and 4-8 years) and interdecadal (8-16 years) time scales are significantly above the red noise background spectra over certain time periods. Therefore, these three scale bands (2-4 years, 4-8 years, and 8-16 years) were selected to compute the scale-averaged wavelet power (SAWP), as shown in Figure 10. The results showed that significant inter-annual (2-4 years) oscillations for all the indices almost occurred during the entire period of 1958-2017. In terms of CDD, significant inter-annual (4-8 years) oscillations occurred in the 1970s and 1980s, while inter-decadal oscillations were active from 1966 to 1982. For the index of CWD, significant inter-annual (4-8 years) fluctuations occurred in the periods of 1958-1984 and inter-decadal oscillations were active from 1966 to 1996. PRCPTOT showed an obvious inter-annual (4-8 years) oscillations during the entire period and inter-decadal oscillations during 1974-1986. For R10, R20, R25 and R50, the significant inter-annual oscillations were almost active in the entire period, but inter-decadal oscillations occurred in 1980s for R20 and R25, and in 1990s for R50. Both R95p and R99p showed significant inter-annual fluctuations in 1960s-1990s, and obvious inter-decadal oscillations in 1980s-1990s. Inter-annual oscillations for R × 1day and R × 5day were active from 2009-2017 and inter-decadal oscillations occurred in 1970s. Moreover, the extra inter-annual oscillations for R × 5day were also significant during 1958-1986. SDII showed relatively longer inter-decadal oscillations during 1973-2007.  (1958-1967, 1968-1977, 1978-1987, 1988-1997, 1998-2007, and 2008-2017). The decadal precipitation extremes were computed as the average value over each 10-year period.

Periodicity Analysis of Extreme Precipitation Indices
The wavelet power spectra of the regional average series of precipitation extremes over the BTH region were shown in Figure 9. We can observe significant periodic oscillations in all indices. For all the indices, the significant two-to five-year periods at the 5% level were predominant during 1958-2017. Furthermore, the CDD wavelet power spectrum showed 6-10-year modulations of variation, which occurred in the 1970s to the mid-1980s. Similarly, the significant inter-decadal (approximately 8-11 years) oscillations were active from the 1970s to the 1980s. For SDII, the power was broadly distributed, with a period of 2-11 years from the 1990s to 2000, and with a period of 8-11 years in the 1980s. Additionally, continuous wavelet spectrum for R95p also showed the significant wavelet power at 8-10 years during 1982-1993. The global wavelet power spectrum (GWPS) is an equal-weighted average of all the local wavelet power spectra for each scale, which shows dominant scales with no temporal transformation ( Figure 9). Regional averaged GWPS shows a statistically significant low-frequency oscillation with a period of 2-9 years in CDD, a period of 2-12 years in CWD, a period of 2-8 years in R10, R20, R25, and RPCPTOT, a period of 2-10 years in SDII, a period of 2-14 years in R95p, a period of 2-15 years in R50 and R99p, and a period of 2-16 years in R × 1day and R × 5day. According to the local wavelet spectra (Figure 9), inter-annual (2-4 years and 4-8 years) and inter-decadal (8-16 years) time scales are significantly above the red noise background spectra over certain time periods. Therefore, these three scale bands (2-4 years, 4-8 years, and 8-16 years) were selected to compute the scale-averaged wavelet power (SAWP), as shown in Figure 10. The results showed that of time series of regional average precipitation extremes over the BTH region. The thick black contours designate the 95% confidence level of local power against red noise, and the cone of influence (COI) where edge effects might distort is shown as the dashed U-shape line. The dashed line in the global wavelet power spectrum is the 5% significance level, using a red-noise background spectrum test. of time series of regional average precipitation extremes over the BTH region. The thick black contours designate the 95% confidence level of local power against red noise, and the cone of influence (COI) where edge effects might distort is shown as the dashed U-shape line. The dashed line in the global wavelet power spectrum is the 5% significance level, using a red-noise background spectrum test.

Relationship Between Precipitation Extremes and Annual Total Precipitation
Previous studies have indicated that the annual total precipitation correlates well with extreme precipitation [33,[70][71][72]. In this study, we estimate the Spearman's rank correlations between the extreme precipitation indices and annual total precipitation, as shown in Figure 11. All the indices except CDD showed significantly positive correlations at a level of α = 0.05. The correlation coefficients between extreme indices and total precipitation exceeded 0.6, except for CDD, which implied that these precipitation extremes had significant correlations with annual precipitation, especially the PRCPTOT, R10, R20, and R25 (correlation coefficient >0.90). Therefore, the decrease in precipitation extremes over BTH reflect the decrease in annual total precipitation to some extent, which match with the finding of previous studies [69]. In reverse, the decline in precipitation extremes mainly caused by the descend in annual total precipitation. Figure 11 also illustrates the correlation coefficients among the extreme precipitation indices. Except for CDD, all the indices have a significant correlation with most of the coefficients higher than 0.5 (p < 0.01). The weak negative correlations between CDD and the other extreme precipitation indices were also found in Figure 11. Atmosphere 2019, 10, x FOR PEER R EVIEW 18 of 30

Relationship Between Precipitation Extremes and Annual Total Precipitation
Previous studies have indicated that the annual total precipitation correlates well with extreme precipitation [33,[70][71][72]. In this study, we estimate the Spearman's rank correlations between the extreme precipitation indices and annual total precipitation, as shown in Figure 11. All the indices except CDD showed significantly positive correlations at a level of α = 0.05. The correlation coefficients between extreme indices and total precipitation exceeded 0.6, except for CDD, which implied that these precipitation extremes had significant correlations with annual precipitation, especially the PRCPTOT, R10, R20, and R25 (correlation coefficient >0.90). Therefore, the decrease in precipitation extremes over BTH reflect the decrease in annual total precipitation to some extent, which match with the finding of previous studies [69]. In reverse, the decline in precipitation extremes mainly caused by the descend in annual total precipitation. Figure 11 also illustrates the correlation coefficients among the extreme precipitation indices. Except for CDD, all the indices have a significant correlation with most of the coefficients higher than 0.5 (p < 0.01). The weak negative correlations between CDD and the other extreme precipitation indices were also found in Figure 11.

Possible Linkage to Large-Scale Atmospheric Circulation Patterns
As discussed in previous parts, the large-scale circulation patterns have affected the climate in China. As large-scale ocean-atmospheric phenomena may influence precipitation extremes in one year and the next year [41], the effects of one-year and 0-year ahead ENSO, NAO, PDO, and IOD on the annual precipitation extremes were investigated in this study to fully identify their relationships. The indices of climate variability in previous year (NINO3.4-1, NAO-1, IOD-1, and PDO-1) (i.e., time Figure 11. Scatter matrix and correlation coefficients of regional annual precipitation and the extreme precipitation indices in BTH from 1958 to 2017. Bold font means the value is statistical significance at a level of α = 0.05.

Possible Linkage to Large-Scale Atmospheric Circulation Patterns
As discussed in previous parts, the large-scale circulation patterns have affected the climate in China. As large-scale ocean-atmospheric phenomena may influence precipitation extremes in one year and the next year [41], the effects of one-year and 0-year ahead ENSO, NAO, PDO, and IOD on the annual precipitation extremes were investigated in this study to fully identify their relationships.  Table 4. We found that most of extreme precipitation indices in the BTH region were significantly influenced by ENSO with one-year ahead (NINO3.4-1). The p-value for the correlations between NINO3.4-1 and CDD, PRCPTOT, R10, R20, and R25 is less than 0.01. However, almost of all the regional averaged indices were not significantly influenced by the other climate indices, especially for NAO-1, POD-1, and POD-0 with none of extreme indices being significant correlation at 5% and 10% significant levels. Interestingly, almost all the extreme indices have a converse correlation with the climate indices in previous and present year. For example, CDD is positively correlated with NAO-0 but it has a negative correlation with NAO-1, as shown in Table 4. Table 4. Spearman's rank correlation coefficients between regional extreme precipitation indices and climate variables and the number of stations for positive or negative correlation coefficient between precipitation extremes and climate variables. Additionally, for each station, the correlation coefficients are illustrated in Figure 12. It can be seen that the areas with positive correlations between precipitation extremes and atmospheric circulation indices with one-year ahead were larger than those with 0-year ahead. For instance, about 54.  Note: Correlation coefficients significant at 5% level are bolded. The mark * represents the significant correlation at 10% level. In this study, monthly and seasonal R × 1day and R × 5day were also used to analyze their relationships associated with the climate indices, as shown in Figure 13. In terms of monthly R × 1day, Figure 13a showed that significant correlation occurred in April (NAO-1), May (PDO-1, PDO-0), July (NAO-1), August (NINO3.4-1, NINO3.4-0, PDO-0), and November (NINO3.4-0 and IOD-0) at a significance level of α = 0.05. A similar conclusion can be seen for monthly R × 5day in above five months (Figure 13b). Additionally, R × 5day in September showed a significantly correlated with NINO3.4 in previous year. In terms of seasonal series, extreme precipitation in summer and autumn are significantly influenced by large-scale atmospheric circulation. For example, R × 1day and R × 5day in summer are significantly correlated with the NINO3.4-1, and R×5day in summer is also significantly influenced by NINO3.4-0. R × 1day and R × 5day in autumn are significantly correlated with the IOD-1, and R×5day in summer is also significantly influenced by NINO3.4-1. In this study, monthly and seasonal R × 1day and R × 5day were also used to analyze their relationships associated with the climate indices, as shown in Figure 13. In terms of monthly R × 1day, Figure 13a showed that significant correlation occurred in April (NAO-1), May (PDO-1, PDO-0), July (NAO-1), August (NINO3.4-1, NINO3.4-0, PDO-0), and November (NINO3.4-0 and IOD-0) at a significance level of α = 0.05. A similar conclusion can be seen for monthly R × 5day in above five months (Figure 13b). Additionally, R × 5day in September showed a significantly correlated with NINO3.4 in previous year. In terms of seasonal series, extreme precipitation in summer and autumn are significantly influenced by large-scale atmospheric circulation. For example, R × 1day and R × 5day in summer are significantly correlated with the NINO3.4-1, and R × 5day in summer is also significantly influenced by NINO3.4-0. R × 1day and R × 5day in autumn are significantly correlated with the IOD-1, and R × 5day in summer is also significantly influenced by NINO3.4-1.
To fully understand their teleconnections, the periodic structure of precipitation extremes and large-scale atmospheric circulation indices are considered and estimated in this study. Based on the length of the original series, we have chosen a suitable scale to construct the periodic components of the original series in the wavelet transforms. As Figure 11 showed, precipitation extremes have the significant oscillation with a scale of 2-16 years during a certain period. Therefore, we selected 2, 4, 8, and 16 year scales to develop the periodic components of the precipitation extremes and climate indices. In this study, we applied the CWT approach to decompose the original series into four identical resolutions, i.e., two-year periodicity (C1), four-year periodicity (C2), eight-year periodicity (C3), and 16-year periodicity (C4). The Spearman's correlations between the decomposed series (C1, C2, C3, and C4) of precipitation extremes and that of the large-scale atmospheric circulation indices were calculated, as shown in Table 5. We found that almost half of precipitation extremes were significantly correlated with NINO3.4-0 at different periodic components, but the significantly correlations between precipitation extremes and NINO3.4-1 mainly occurred in C1 and C4. The visibly opposite correlations between precipitation extremes and NINO3.4 with one-year ahead and 0-year ahead were shown for C1, whereas the same correlations were found for C4. All the extreme indices except for R99p and R × 1day showed a significantly positive correlation with the low-frequency component (C4) of NAO, and nine of extreme indices also showed a significantly positive correlation with NAO-1 for C2. Most of precipitation indices were significantly correlated with PDO for the low-frequency components (C3 and C4), accounting for 87.5% (PDO-0) and 66.7% (PDO-1). Moreover, almost all the indices for C1 showed significantly positive correlation with PDO-1 except for CDD and R50. In terms of IOD, eight of extreme precipitation indices at C4 were significantly correlated with IOD-0 and IOD-1, and almost all of them showed a significantly positive correlation with IOD-1 for C2. Overall, the influences of climate indices on regional precipitation extremes were composed by different periodic components, with the higher correlations mainly occurring at a low-frequency component (e.g., C4). However, the discrepancies of correlations between precipitation extremes and climate indices with 0-year and one-year ahead mainly occurred at a high-frequency component (e.g., C1). To fully understand their teleconnections, the periodic structure of precipitation extremes and large-scale atmospheric circulation indices are considered and estimated in this study. Based on the length of the original series, we have chosen a suitable scale to construct the periodic components of the original series in the wavelet transforms. As Figure 11 showed, precipitation extremes have the significant oscillation with a scale of 2-16 years during a certain period. Therefore, we selected 2, 4, 8, and 16 year scales to develop the periodic components of the precipitation extremes and climate indices. In this study, we applied the CWT approach to decompose the original series into four identical resolutions, i.e., two-year periodicity (C1), four-year periodicity (C2), eight-year periodicity (C3), and 16-year periodicity (C4). The Spearman's correlations between the decomposed series (C1, C2, C3, and C4) of precipitation extremes and that of the large-scale atmospheric circulation indices Figure 13. Spearman's rank correlation coefficients between (a) monthly R × 1day, (b) monthly R × 5day, and (c) seasonal R × 1day and R × 5day, and atmospheric circulation indices during 1958-2017. Black dashed lines mean the positive or negative correlation coefficients with statistical significance at 5% level.

Discussion
Previous studies have reported that extreme precipitation has changed over the past few decades, but there are apparent discrepancies in regarding these changes over different regions in China, including BTH or North China [73]. Overall, the temporal trends in our work are similar to previous studies in reporting decreasing trends in precipitation extremes over the past few decades in BTH [33,36,40] and the Circum-Bohai-Sea region in the north of China [63,74]. Certainly, there are some discrepancies between our work and previous studies. For instance, Mei et al. [33] found that the significantly decreasing trends were only shown in R × 1day and R × 5day during 1960-2013. But we found that more indices (e.g., R50, R95p, and R99p) exhibited significant decreasing trends in BTH during 1958-2017 ( Figure 3 and Table 3). In terms of spatial patterns of precipitation extremes, our results are also consistent with the previous studies in BTH or neighboring regions. Overall, for most of precipitation extremes except for CDD and CWD, the highest values located in the eastern BTH region, while the lowest values occurred in the western part (Figure 2), which was mainly influenced by the terrain (Figure 1a) with the high values in the northwest and the low in the southeast of BTH. Spatial patterns of precipitation extremes are similar to that of annual total precipitation (Figure 1b and [69]) with high correlations (Figure 11), which was consistent with previous studies [6,33]. The correlations between the trend magnitudes of annual total precipitation and precipitation extremes for all the stations were also estimated, as shown in Figure 14. All the extreme precipitation indices were positively correlated with annual total precipitation, which is also different from that found in Mei et al. [33]. In Mei's work, they found the annual trend magnitudes of CDD and CWD were negatively correlated with that of annual total precipitation (ATP), with the correlations being significant for CDD. However, in our study, all the precipitation extremes are statistically significant at the 0.05 level except for CDD and CWD. The discrepancies may be caused by the data with different time periods used in the two studies. Furthermore, we also analyzed variations of the monthly, seasonal and decadal precipitation extremes, which has been reported in the other regions [64,75,76] but not in the study area [33,36,40]. As all known, North China has a monsoon climate with salient four seasons and about 60-80% of total precipitation amount in summer. Thus, it is helpful to understand the annual variations by the analysis of monthly or seasonal changes in precipitation extremes. We found that the decreasing trends in annual precipitation extremes were mainly caused by the descend in summer, especially in July and August ( Figure 6). Although the significantly increasing trend can be found in many months or seasons ( Figure 6), their increases have not remedied the descends of annual R × 1day and R × 5day over BTH.
As previous mentioned, changes in precipitation extremes are related to the large-scale circulation change [41,42,[77][78][79]. The periodicity characteristics of precipitation extremes will be helpful to understand their associations or linkages with climate indices, which has been investigated in the other regions in previous studies [40,[60][61][62][63][64][65]. Our results showed that the significant two-to five-year periods at the 5% level were predominant for all the indices during the entire period, which are roughly consistent with the previous work of Wang et al. [63]. As pointed by Moron et al. [77], there are appropriate evidences that both the quasi-biennial (2-3 years) and quasi-quadrennial (4-6 years) peaks of climate variables are related to ENSO. That is, the inter-annual (2-5 years) oscillation in extreme precipitation over BTH may be influenced by ENSO, which was also verified based by the correlation analysis (Tables 4 and 5). Moreover, quasi-decadal oscillations of precipitation extremes were also detected. Furthermore, to better understand their relationships between precipitation extremes and climate anomalies, the quantitative analyses should be conducted. However, the correlations between precipitation extremes and climate indices over BTH or north China has not been examined or fully explained in the previous studies. Thus, we computed the correlations between the precipitation extremes and large-scale climate anomalies to reveal their relationships. Our results showed that most of extreme precipitation indices in the BTH region were significantly influenced by ENSO with one-year ahead (Table 4 and Figure 10), which are consistent with the finding in Xiao et al. [41] and Gao et al. [79]. Besides, we also found that almost all the extreme indices exist a converse correlation with the climate indices in one-year ahead and 0-year ahead. That is, the relationships of annual extreme precipitation indices with the climate indices are different when the climate indices are at different stages [41]. Moreover, the precipitation extremes in different months were influenced by different climate indices or more than one index ( Figure 11). To some extent, our results also showed that precipitation extremes at any scale are influenced by many climate indices together. As previous studies mentioned, the climate of eastern China is dominated by the East Asian Monsoon, which is significantly influenced by ENSO [80,81]. Meanwhile, the East Asian Monsoon is also influenced by NAO, PDO, and IOD, and then affect the regional patterns of precipitation or precipitation extremes. Our findings also provide an illustration for this conclusion. As previous mentioned, changes in precipitation extremes are related to the large-scale circulation change [41,42,[77][78][79]. The periodicity characteristics of precipitation extremes will be helpful to understand their associations or linkages with climate indices, which has been investigated in the other regions in previous studies [40,[60][61][62][63][64][65]. Our results showed that the significant two-to five-year periods at the 5% level were predominant for all the indices during the entire period, which are roughly consistent with the previous work of Wang et al. [63]. As pointed by Moron et al. [77], there are appropriate evidences that both the quasi-biennial (2-3 years) and quasi-quadrennial (4-6 years) peaks of climate variables are related to ENSO. That is, the inter-annual (2-5 years) oscillation in extreme precipitation over BTH may be influenced by ENSO, which was also verified based by the correlation analysis (Table 4 and Table 5). Moreover, quasi-decadal oscillations of precipitation extremes were also detected. Furthermore, to better understand their relationships between precipitation extremes and climate anomalies, the quantitative analyses should be conducted. However, the correlations between precipitation extremes and climate indices over BTH or north China has not been examined or fully explained in the previous studies. Thus, we computed the correlations between the precipitation extremes and large-scale climate anomalies to reveal their relationships. Our results showed that most of extreme precipitation indices in the BTH region were significantly influenced by ENSO with one-year ahead (Table 4 and Figure 10), which are consistent with the finding in Xiao et al. [41] and Gao et al. [79]. Besides, we also found that almost all the extreme indices exist a converse correlation with the climate indices in one-year ahead and 0-year ahead. That is, the relationships of annual extreme precipitation indices with the climate indices are different when the climate indices are at different stages [41]. Moreover, the precipitation extremes Overall, in this study, the spatio-temporal variation and characteristics of precipitation extremes in the BTH region were studied using various methods, and useful information was obtained, which is potentially is useful for natural disaster prevention and mitigation. In this study, we use the simple linear regression and Mann-Kendall test to detect the trends of precipitation extremes. Linear regression method, one of parametric tests, can be capable of quantifying the change in the data, which assumes the time series and the errors following a normal distribution. A normality test of regional precipitation extremes is detected using one-sample Kolmogorov-Smirnov test method, which clearly indicates that all the data can be regarded as normally distributed. Mann-Kendall trend test, one of the widely used non-parametric tests, being a function of the ranks of the observations rather than their actual values, is not affected by the actual distribution of the data and is less sensitive to outliers. Thus, Mann-Kendall test, as well as other non-parametric trend tests, is more suitable for detecting trends in hydrological time series, which are usually skewed and may be contaminated with outliers [82]. A common requirement of both parametric and non-parametric trend tests is that the data be independent [82,83]. In this study, the possible effects of serial correlation have been not considered. However, in the future work, the serial correlation should be detected and eliminated by pre-whitening of the data when these methods are used to detect the trends of hydrologic series. Additionally, the observed data with the sparseness of meteorological gauges (only 28 national-level meteorological stations) may have led to uncertainties and limitation in our results and conclusions, especially in spatial variability [84], which need to be improved in the future work. The higher temporal (e.g., hourly time series) and spatial (e.g., 316 meteorological stations used in Zhao et al. [69]) resolutions would lead to more reliable and specific results. Moreover, the selection of indices of climate extremes is also an important step to assess their variations. A number of indices were suggested, each aimed at depicting a certain aspect of climate extremes. In this study, the 12 precipitation indices were chosen from the indicators recommended by the ETCCDMI. From our findings, several indices showed similar results, especially for these indices belongings to same properties. For instance, the results of spatial patterns and temporal trends in R10, R20, and R25 are almost similar (Table 3, Figures 3-5). Thus, in order to quickly and efficiently analyze the variations of precipitation extremes, these similar indicators can be eliminated based on the different demands. One representative index can be chosen for each type of precipitation extremes, such as the percentile-based threshold, and fix-based threshold.
In this study, the possible relationship between precipitation extremes and climate indices is also discussed, including the in-phase and out-of-phase linkages between them. Past studies have revealed that the relationships between precipitation extremes and different teleconnection phases tend to be asymmetrical and generally do not portray linear behavior [85]. Thus, the nonlinear and different influences associated with different phases of these teleconnections have important implications on this issue. Aside from the above two aspects, the multiple-variables jointly relationships between precipitation extremes and climate indices should be also further investigated in the future work. Moreover, the mechanistic understanding of how these climate drivers impact extreme precipitation in BTH through forcing large-scale atmospheric circulation changes need to be further investigated, and the analyses of the physical causes of these variability will form the basis of our future research.

Summary and Conclusions
In this work, attempts are made to investigate the long-term changes of precipitation extremes at different scales. The analyses of spatial-temporal variations in precipitation extremes across BTH are conducted with 12 indices based on daily precipitation records from 28 meteorological stations. Moreover, the possible relationships between extreme precipitation events and large-scale atmospheric circulation patterns are also analyzed. Based on the above results and discussion, our mainly findings are summarized as follows: (1) The spatial distributions of precipitation extremes apart from CDD and CWD were similar to the spatial patterns of annual precipitation in BTH, with the highest in the northeast and the lowest in the northwest. The mean values of CDD were greater in the northwest of BTH than those in the southeast. The spatial pattern of CWD was most unique, with the high value in the north and south and the low value in the middle of the region.
(2) Regionally averaged precipitation extremes were characterized by decreasing trends, of which five indices (R50, R95p, R99p, R × 1day, and R × 5day) exhibited statistically significant decreasing trends at a level of 0.05. From monthly and seasonal scale, almost all of the high values in R × 1day and R × 5day occurred in summer, especially in July and August due to the impacts of East Asian monsoon climate on inter-annual uneven distribution of precipitation. The significant decreasing trends in annual R × 1day and R × 5day were mainly caused by the significant decreasing trend in summer.
(3) The extreme precipitation events may be more acute in the first decade (T1, 1958-1967) over those areas, whereas the relatively less probability or weights of extreme precipitation events occurred in T5 (1998-2007). All the indices except CDD experienced clear fluctuations among these sub-periods with a complex trend of decreasing (from T1 to T3)-increasing (from T3 to T4)-decreasing (from T4 to T5)-increasing (from T5 to T6).
(4) For all the indices, the significant two-to five-year periods at the 5% level were predominant during the entire period of 1958-2017 based on the wavelet analysis. Specially, there were a statistically significant low-frequency oscillation with a period of 2-9 years in CDD, a period of 2-12 years in CWD,