Observed Quantile Features of Summertime Temperature Trends over China: Non-Negligible Role of Temperature Variability

Changes in temperature variability can have more serious social and ecological impacts than changes in the mean state of temperature, especially when they are concurrent with global warming. The present study examines the summertime temperatures’ trends over China from the quantile perspective. Through fully investigating the quantile trends (QTs) of the maximum (Tmax) and minimum temperature (Tmin) using the homogenized observation data and quantile regression analysis, we identify evident region-specific quantile features of summertime temperature trends. In most of northern China, the QTs in Tmax and Tmin for all percentiles generally show strong uniform warmings, which are dominated by a warm shift in mean state temperatures. In contrast, the QTs of Tmax in the Yangtze River Basin show distinguishable inter-quantile features, i.e., an increasing tendency of QTs from cooling trends in the lower percentile to warming trends in the higher percentile. Further investigations show that such robust growing QTs of Tmax across quantiles are dominated by the temperature variance. Our results highlight that more attention should be paid to the region-specific dominance of temperature variability in trends and the related causes.


INTRODUCTION
Global warming during the past 50 years is unequivocal, according to the continuous assessments of the Intergovernmental Panel on Climate Change (IPCC) reports (e.g., IPCC, 2012IPCC, , 2013IPCC, , 2018IPCC, , 2021. Such warming is not spread uniformly across the globe; instead, distinctive region-specific features with various magnitudes or even opposite signs are detected. In particular, previous studies have demonstrated that China has also experienced robust warming trends during the past decades from the national mean perspective (e.g., Ding et al., 2007). However, strong regional and seasonal differences of such warming trends have also been reported (e.g., Qian and Qin, 2006;Xu et al., 2017;Li and Zha, 2019). For example,  found that from the long-term (>50a) scale, summertime temperature during the past half-century exhibited a larger warming trend over the north of China, whereas very slight warming or even cooling trends were detected in the south of China. Similarly, this pattern of warming trends is also recognized by several previous studies (Qian and Qin, 2006;Xu et al., 2013;Xu et al., 2017), despite some slight gaps in magnitudes due to different data, methods, and analysis periods. Regarding daily maximum temperature (Tmax) trends, Li et al. (2015) found that the trends of summer Tmax during 1961-2012 over China also showed robust warming in the north and slight warming or cooling in the south, which resembles the above-mentioned warming pattern of mean temperature.
Concurrent with such warming patterns, the trends of hightemperature extreme events (e.g., heatwaves), however, are evidently distinct from that of temperatures, particularly in regions of southern China (e.g., Yangtze River Basin). As revealed by Chen and Zhai (2017), the summertime compound hot events were growing 15-20% (i.e., 2-4 days) per decade, especially in the north China and south China. Hu et al. (2017) investigated the trends of summertime heat days over China and identified a significant increasing trend of heat days in southern China, particularly in the south of the Yangtze River and the southeast coast of China. Meanwhile, evident decreasing trends of heat days were also detected in Henan and Shandong province in their study. Similarly, Qian et al. (2019) also showed that the summertime warm days (TX90P) have significantly increased in western China and East China, whereas slight decreasing trends were coherently detected in some regions of central-eastern China. From the above evidence, it is worth noting that the increasing (decreasing) pattern of heat extremes seems inconsistent with the warming (cooling) pattern of Tmax, especially for South China. This further suggests that the trends of Tmax in different percentiles are probably distinct because most studies used percentile-based indices (e.g., TX90P) to define extreme events. What is the feature and the main cause of the above-mentioned inconsistency (i.e., different trends across percentiles) over China remains unknown.
According to the IPCC SREX report (IPCC, 2012), the variability of extreme temperatures is not only influenced by the mean state shift in temperature, but also related to the changing shape of the temperature distribution (e.g., changes in temperature variance). Furthermore, McKinnon et al. (2016) investigated the relative roles of four statistical moments of temperature distribution (i.e., mean, variance, skewness, and kurtosis) on the summertime temperature trends. Their results showed that the warming trends in summer temperatures over the majority of the Northern Hemisphere were due to shifts in the temperature mean state, but in some specific regions (e.g., mid-to high-latitude Eurasia), changes in temperature variance also make a non-negligible contribution to the trends. As global warming and the concurrent growth of extreme events become more severe, there are growing numbers of studies that highlight the impact of temperature variability on extreme events more widely than the mean state shifts (e.g., Schär et al., 2004;Screen, 2014;Wang et al., 2020;Yang et al., 2021). However, for China, the assessment of relative contribution of temperature mean state shifts and changes in temperature variability to the temperature trends in the context of global warming is still lacking.
In recent years, China has suffered multiple dangerous and deadly heatwaves (e.g., 2013 super summertime heat wave event in southern China; Xia et al., 2016), which caused catastrophic social and economic consequences. Multi-model-based future projections suggest that if global warming persists, deadly heat waves like the 2013 event will become increasingly frequent (e.g., Sun et al., 2014;Sun et al., 2018). Therefore, based on the above gaps, two important scientific questions are raised, namely, what are the quantile features of the summer temperature trends in China? And what is the relative contribution of changes in the mean state and variability of temperatures to the trends? The rest of the article is presented in the following order: Section 2 describes the data and methods used in the study, Section 3 shows the main results of the study, and Section 4 contains the necessary discussions and summaries of the current findings.

Data
The homogenized China daily maximum and minimum surface air temperature (Tmax and Tmin hereafter) from 2,419 meteorological stations are used to investigate the variations of Tmax and Tmin during 1961-2017. In addition, precipitation data from the same homogenized dataset are also used in this paper to investigate the covariability between temperature and precipitation. Those datasets are derived from the China Meteorological Administration (http://data.cma.cn/en). Comparing with the original station records, the potential inhomogeneities (i.e., influences from non-climatic factors) of the data are significantly reduced in those datasets, and the quality of the observed data are further improved. Details of the production of these datasets are referred to Cao et al. (2016) and Yang and Li (2014). To further minimize the inhomogeneities from missing data and stations' irregular geographic distribution, two additional criteria are applied: 1) For meteorological variables, the stations that have any missing data of summer days (i.e., June-July-August) during 1961-2017 are excluded; 2) Following Chai et al. (2018) and , we introduced a re-grid algorithm to redistribute the station data onto a 0.5°× 0.5°longitude-latitude grids; details can be found in .
Finally, the processed meteorological variables are continuous (i.e., no missing data) in time and equal-weighted in space (i.e., 0.5°× 0.5°grid boxes). The illustration of the original stations' distributions and the re-gridded equal-weighted grid boxes can be found in Supplementary Figure S1.

Quantile Regression
To investigate the role of the changing shape of temperature distributions on its trends across different percentiles, we Frontiers in Earth Science | www.frontiersin.org December 2021 | Volume 9 | Article 762563 introduce the quantile regression (QR) method to estimate the quantile trends (QTs). The QR method can be seen as an extension of the method of ordinary least squares (OLS); it can discover more comprehensive relationships between variables in any quantiles (e.g., median). This method has been widely used for assessing the variability of extreme events in recent decades (e.g., McKinnon et al., 2016;Mueller et al., 2016;Yin et al., 2018). A full description of QR method can be found in Koenker and Bassett (1978) and Cade and Noon (2003). The QR source code can be implemented using the python statsmodels (Seabold and Perktold, 2010) package. Detailed information of statsmodels package and QR source code can be found in http:// www.statsmodels.org/stable/_modules/statsmodels/regression/ quantile_regression.html#QuantReg.
In this study, we first remove the seasonal cycle of Tmax and Tmin in each summer between 1961 and 2017 according to the daily climatology of the whole period (i.e., 1961-2017) to obtain the daily Tmax and Tmin anomalies at summer months (i.e., June-July-August). The QR method is then applied to a 57-year (i.e., 1961-2017) time series of summer daily Tmax (Tmin) anomalies for each grid box (i.e., total series length for QR is 57 years × 92 days), and linear trends in Tmax (Tmin) are calculated in 5% steps from the 5th percentile to the 95th percentile. It should be noted that the interpretation of QR results (i.e., QTs) is very similar to the OLS regression, with the difference that instead of predicting the mean of the dependent variable, quantile regression looks at the quantile of the dependent variable. Here, we illustrate with a specific example: suppose we do QR on a Tmax series and obtain its 90th percentile QT of 0.2°C/ 10a. This means that if the trend is translated into a change of extremes, then hot days will plausibly increase in this example. Similarly, the 10th percentile QT of Tmax (Tmin) implies a change in cold days (nights). In contrast to OLS regression, which can only predict the mean of the dependent variable, QR can be applied to any quantile points, so by investigating the QT at multiple percentiles, we can establish a better understanding of the variation in the conditional distribution of the temperature series. Specifically, for QTs at the two tails of the temperature distribution, it is closely related to changes in extreme events.
Furthermore, in order to quantify the contribution of changes in temperature mean state and shape of temperature distribution (STD) to the temperature QTs, we can assume that at a certain percentile, the trend in temperature is only influenced by changes in the temperature mean state and STD. Thus, we can measure the change in temperature QTs due to STD as the difference between the temperature QTs and the mean temperature trend (i.e., Trend STD Trend QT − Trend MEAN ), and the relative contributions of changing STD in temperatures' trends magnitudes can be calculated as: Contrib TrendSTD |Trend STD |/(|Trend STD | + |Trend MEAN |). In addition, some commonly used statistical methods (e.g., linear trend, Empirical Orthogonal Function, Kernel Density Estimation, etc.) are also employed in our study. The significance of the Tmax and Tmin trends in all quantiles are tested using a two-tailed t-test.

Quantile Features of Summertime Tmax and Tmin Trends over China
Previous studies have suggested that changes in STD (e.g., changes in variance, skewness, and kurtosis) can also exert important impacts on temperature extremes (IPCC, 2012;McKinnon et al., 2016;Li et al., 2017). Note that the term "changes in temperature variability" in the following text refers specifically to changes in STD (i.e., variance, kurtosis, and skewness) other than the mean state shift. To investigate whether there are evident differences of Tmax (Tmin) trends among percentiles, Figure 1 shows the selected QTs of summertime Tmax and Tmin at three different percentiles (i.e., 10th, 50th, and 90th percentiles) during 1961-2017 over China. As seen from Figures 1A1-A3, the Tmax QTs clearly show distinguishable features from low to high percentiles, while the Tmin QTs are less variable across different quantiles.
Specifically, the QTs of Tmax at the 10th percentile exhibits consistent significant warming trends (more than 0.2°C/10a) over most regions in the north of China (i.e., to the north of 40°N), especially for North China ( Figure 1A1). In contrast, broad cooling trends are detected in the Yangtze River Basin (YRB) and Huaihe River Basin with significant (p < 0.05) cooling trends of −0.3 to −0.4°C/10a. Compared to the low percentile (10th percentile), the trend of the median (i.e., 50th percentile) of Tmax generally shows a similar warming pattern to the mean Tmax trend over China (e.g., mentioned in Li et al., 2015 and Supplementary Figure S2A). In particular, significant warming trends with a magnitude of 0.2-0.4°C/10a are shown in the north of China (i.e., to the north of 40°N) and coastal areas of South and East China, whereas there is no clear trend in the Central China (slightly warming or cooling with a magnitude less than ±0.1°C/10a; Figure 1A2). For the Tmax at high percentile (i.e., 90th percentile), the pattern of significant warming trends is very similar to that of the median of Tmax, i.e., significant warming trends are detected in the north of China (i.e., to the north of 40°N) and the coastal regions of southeastern China ( Figure 1A2 vs. Figure 1A3). Notably, there are two distinguishable characteristics of Tmax trends at 90th percentile in contrast to that of the median. First, for YRB, the robustness and warming magnitude of Tmax trends at the high 90th percentile is stronger (i.e., more grids show significant warming trends) and higher (i.e., 0.2°C/10a higher) than that at the median. Second, significant cooling trends of Tmax at the high percentile with the magnitude of −0.2 to −0.4°C/10a are also detected in the Huang-Huai-Hai plain, where no clear trends are shown at the median ( Figure 1A2 vs. Figure 1A3). In conclusion, the QTs of Tmax show a largely consistent and significant warming in the north of China, but in the YRB region, this exhibits an asymmetrical feature across quantiles, i.e., an increasing tendency from cooling at the low percentile (i.e., 10th) to warming at the high percentile (i.e., 90th).
To further interpret the current results, we can roughly link the trend of Tmax at the 90th (50th) percentile to the trend of hot extremes (mean Tmax). Therefore, the current results of Tmax Frontiers in Earth Science | www.frontiersin.org December 2021 | Volume 9 | Article 762563 trends at the 50th and 90th percentile generally agree with previous studies regarding the observed trends of mean Tmax and hot extremes during the past 50 years (e.g., Deng et al., 2019;Qian et al., 2019). Moreover, it should be noted that the evident opposite Tmax QTs between the high and low percentiles at the YRB region may be caused by changes in non-mean factors of temperature (e.g., variance and kurtosis), according to McKinnon et al. (2016). This further implies that the changing shape of Tmax distributions (i.e., changes in variance, skewness, and kurtosis) may dominate the trend of hot/cold extremes in the YRB region. Unlike the Tmax, the QTs of Tmin are more consistent among different percentiles ( Figures 1B1-B3), with all QTs showing almost the same trends as the mean Tmin (Supplementary Figure S2B). Taking the trend of the median of Tmin as an example, broad and consistent warming trends are detected nationally, while the warming magnitude in the north of China (0.3-0.5°C/10a to the north of 35°N) is higher than that in the south of China (0.1-0.4°C/10a to the south of 35°N). In addition, the trends of Tmin at 10th and 90th percentiles generally resemble that of the median (Figures 1B1-B3). Comparing the above QT results of Tmin and Tmax, two distinguishable features can be concluded: 1) The QTs of Tmax show an asymmetric feature across percentiles (i.e., an increasing tendency of Tmax trends from low to high percentiles), whereas there are symmetric responses of Tmin trends in all percentiles ( Figure 1B), and 2) the warming trends magnitudes of Tmin in all percentiles are higher than that of Tmax, which agrees with the previous findings regarding the more rapid warming of Tmin than Tmax in recent decades (e.g., Davy et al., 2017). Moreover, this asymmetric trends between Tmax and Tmin, irrespective of the quantile, also implies that the diurnal temperature range (DTR) has shown a clear decreasing trend over the past decades, which has been well recognized in recent studies (Xia, 2013;Xue et al., 2019). In addition, physical reasons for the declining DTR are probably related to the sunshine duration, cloud cover, and large-scale vertical motion anomalies (e.g., Liu et al., 2018).

Stability of Quantile Trends of Tmax and Tmin over the YRB Region
In the previous part, we found that there were distinct responses of QTs between Tmax and Tmin, particularly for the YRB region. However, it is obvious that a trend of temperature can be potentially affected by the start and end year (i.e., the position and length of trend window), which in turn introduce uncertainties and reduce the confidence of the current conclusions. Therefore, in this part, we will fully discuss the stability (i.e., uncertainty) of the observed features of the QTs of Tmax and Tmin over the YRB region. First, we take one particular point (i.e., 120°E, 31°N) in the YRB region as an example to conduct a running window QR analysis (5th to 95th in step of 5%) that covers the window length (WL) from 49 to 57 years and the start year (SY) from 1961 to 1969. Through the running window analysis, the QTs of Tmax and Tmin for 45 different periods will be calculated. It is worth noting that we focus on long-term trends in temperatures in the following analysis, and therefore, only WLs of about 50a or more are shown here. Figures 2A,B exhibit the selecting results (i.e., 10th, 50th, and 90th percentiles) of running window QR analysis for Tmax and Tmin, respectively. A first inspection shows that the QTs of Tmax and Tmin from low to high percentiles within different periods on the whole exhibit similar features to the spatial results of Section 3.1. In detail, consistent cooling trends are detected at the low percentile (i.e., 10th percentile) of Tmax in most of the selected periods (44 periods out of 45; Figure 2A1). For the robustness, significant cooling trends of Tmax in the low percentile with −0.2 to −0.4°C/ 10a are shown in relatively long periods (i.e., 53-57 years), whereas weak and insignificant cooling trends are shown in relatively short periods (49-51 years). The trends of the median of Tmax generally show consistent and significant moderate (0.1-0.3°C/10a) warming in all periods (except for the trends start at 1961; Figure 2A2). Similarly, consistent but stronger warming trends (0.3-0.5°C/10a) are detected in the Tmax at the high percentile (i.e., 90th) for all periods ( Figure 2A3). All the Tmax trends at 90th percentile within different periods are statistically significant (i.e., p < 0.01). For Tmin, running window QT results show uniformly strong and significant warming at every percentile ( Figures 2B1-B3), which indicates that the uniform warmings of Tmin at all quantiles for this specific point have a high degree of confidence. In addition to the above results, it appears that the significance of the 10th percentile cooling QTs of Tmax are not consistent across WLs. In particular for the results after the SY of 1965, significance of Tmax QTs across windows are very weak. For interpreting this finding, we extend the WL (SY) of running window QT analysis to 20-57a  in order to provide a panoramic view (Supplementary Figure S3). Interestingly, in addition to the characteristics we show for the long-term QTs, it is also clear that Tmax QTs in all percentiles exhibit significant interdecadal variations, i.e., the cooling trends of the 60-80s and the significant warming trends of the 80s and onwards (Supplementary Figures  S3A1-A3). These interdecadal variations reflect the effect of internal climatic variability (e.g., Atlantic Multidecadal Oscillation) or other external forcing factors. Thus, this complex multi-scale interdecadal variations superimposed on the long-term QTs we identified will plausibly lead to different degrees of significance for different WLs.
Furthermore, we calculate the statistics of running window QTs for all the grid points at the YRB region to see whether the quantile features of Tmax and Tmin are still pronounced for the whole region. The results of the statistics (i.e., the first and third quartiles, the median, and the kernel density estimation) are shown as a violin plot (Figure 3). In general, the QT characteristics of Tmin and Tmax mentioned above remain  (120°E, 31°N) within the YRB region, unit:°C/10a. The x-axis represents the start year of the analysis period, while the y-axis represents the length of the analysis period. White and black dots represent trends that are statistically significant at 99% and 95% confidence levels, respectively.
Frontiers in Earth Science | www.frontiersin.org December 2021 | Volume 9 | Article 762563 valid in the statistical information ( Figure 3) for the YRB region, i.e., the tendency for the Tmax QTs to increase from low to high percentile and the consistent warming of the Tmin QTs at each percentile. Specifically, the uniform features of Tmin are very stable for the YRB region, with the shape of the QTs distribution (outlines of the violin plot, Figure 3) being almost identical for all percentiles. Tmin QTs exhibit a warming of 0.2°C/10a at the median for all selected percentiles, and the majority (over 80%) of the data in the region show warming trends. For Tmax QTs, the shape of the distribution is also consistently spread across the percentiles, but as the quartiles increase, the overall distribution of Tmax QTs shifts upwards (i.e., shifts towards warming). As a result, the 10th percentile of Tmax shows a clear cooling QT (i.e., ∼−0.1°C/10a at the median) and 75% of the data in the YRB region agree with the cooling, while the 90th percentile of Tmax is closer to Tmin QTs ( Figure 3). In addition, we can also find that the rate of increase in the QTs of Tmax from the low to high percentiles raises rapidly from −0.1°C/10a at the median to about 0.1°C/10a at the median in the low to medium percentiles (i.e., 10th∼50th, Figure 3). However, in the 50th∼90th interval, Tmax QTs increase at less than half the rate of the 10th∼50th interval. Therefore, it can be inferred from the above results that the tendency of the Tmax intra-quantile QTs identified in the YRB region may be largely attributed to changes in the lower percentiles (e.g., 10th).  McKinnon et al., 2016;Li et al., 2017;Tamarin-Brodsky et al., 2020). Hence in this part, we will discuss the relative contributions of changing mean state and STD on the QTs of Tmax and Tmin. Figure 4 shows the spatial features of Tmax and Tmin trends in the 10th and 90th percentiles due to changing STD and their relative contributions. As we expected, the changing STD of Tmax caused a broad cooling (warming) trends in the 10th (90th) percentile over the most parts in the south of China (i.e., to the south of 35°N), especially for the YRB region with a cooling (warming) trend of up to ±0.3°C/10a (Figures 4A1,A2). Furthermore, relative contributions of changing STD on Tmax QTs at 10th and 90th percentiles in all grids are also calculated ( Figures 4B1,B2). Results show that the absolute magnitudes of Tmax QTs due to changing STD contribute over 70% of the whole Tmax trends in both 10th and 90th percentiles over the majority of the YRB region ( Figures 4B1,B2). In addition, further investigations regarding the detailed reasons of the changing STD are also conducted. According to McKinnon et al. (2016), we can roughly express the spatial patterns of temperatures' QTs due to variance, skewness, and kurtosis by the first three modes of Empirical Orthogonal Function (EOF) analysis on the intraquantile variations of temperatures QTs. The EOF results of intra-quantile modes of Tmax trends are shown in Figure 5. Looking at the principal components of the EOF, the first mode shows a clear non-linear trend characteristic and the large variation areas in the spatial distribution are also highly consistent with the variability-dominated YRB regions that we focus on (Figures 5A1,B1). This leading model can represent the effect of changes in temperature variance on the temperature QTs across percentiles based on the findings of McKinnon et al. (2016). Obviously, we can further conclude that the observed quantile features of Tmax trends in YRB due to changing STD are mainly caused by changes in variance of Tmax from the evidence (e.g., the percent variance of the first EOF mode reaches 80.2%) of EOF results. In addition, the variance contributions are smaller FIGURE 3 | Violin plot (showing kernel density estimate of distribution) of running window quantile trends at 10th, 30th, 50th, 70th, and 90th percentiles for all the grid points at the YRB region, with the differences between Tmax (red) and Tmin (blue). Note that the quantile trend analysis for each point is exactly the same as Figure 2. The number of data in each violin is: 45 different running window trends × 318 (265) grid points of Tmax (Tmin) in the YRB region. The meanings of each element within the violin plot can be found in the attached legend box.

The Role of Non-mean Factors in the Tmax and Tmin Trends
Frontiers in Earth Science | www.frontiersin.org December 2021 | Volume 9 | Article 762563 for the second and third EOF modes (i.e., only 10.1% and 3.4% for EOF2 and EOF3), which can also represent the effect of temperature skewness and kurtosis changes on Tmax intraquantile QTs, respectively ( Figures 5A2,B2,A3,B3).
Contrasting with Tmax trends, the magnitudes of Tmin trends and the contributions due to the changing STD are still weak ( Figures 4A3,A4,B3,B4), which is consistent with the abovementioned mean state warming-induced uniformly pattern of Tmin trends. In addition, our current results generally agree with McKinnon et al. (2016), which also shows some hotspots regarding the variance-induced Tmax trends over the north of 30°N.
To sum up, all the above results highlight the non-negligible role of changing STD (especially increasing variance) on Tmax trends during the past decades in the context of global warming, particularly for some variance-sensitive regions (e.g., YRB). These also imply that the summers in the YRB region have become more and more volatile during the past 57 years, i.e., accelerated warming in the high percentiles (i.e., 90th∼99th) of Tmax with rapid increases of extreme hot events (e.g., deadly heatwaves).

DISCUSSION AND CONCLUSION
In the above sections, the quantile features of Tmax and Tmin trends over China during the past half-century are investigated. The evident differences between Tmax and Tmin trends are observed over the YRB region, i.e., "variance-induced region-specific trends patterns of Tmax" versus "mean state shift-dominated uniform warmings of Tmin". In particular, for Tmax, we emphasized the non-negligible role of temperature variance changes in modulating Tmax trends. However, in addition to the above conclusions, there are several additional caveats that should be highlighted and discussed.
First, we discuss whether such quantile features and patterns of Tmax trends are still pronounced in the sub-seasonal scale (i.e., monthly). As a matter of fact, the current results represent the observed robust region-specific variance-induced Tmax trends in the view of seasonal scale. Hence, we further repeated our QR analysis for each month in summer (i.e., June, July, and August) in order to investigate the QTs of temperatures for sub-seasonal scale. The results generally show that the observed differences between Tmax and Tmin as well as the different patterns of Tmax trends between low and high percentiles over the YRB region are still existent and significantly pronounced (Supplementary Figures S4-S6). Despite these similarities, stronger variance-induced distinctive Tmax trends between high and low percentiles are also observed in late summer than in early summer (i.e., August vs. June). In brief, the sub-seasonal features of Tmax trends show not only consistency with seasonal scale but also additional unique different features among each month.
Second, we note that variance-induced robust increasing trends of Tmax only occur in the YRB region, so why is this such a unique regional signature? In summer, the baroclinicity and circulation variability over mid-latitudes of East Asia are much weaker than those in winter half year. Therefore, the temperature variability is less impacted by mid-latitude dynamics, e.g., cold front or cold air intrusion. In contrast, temperature, and particularly Tmax, is more sensitive to radiative processes, involving cloud cover and rainfall. Several studies have highlighted the importance of covariability between mean temperature and precipitation (e.g., Trenberth and Shea, 2005;Miao et al., 2016;Guo et al., 2020). This relationship may also be valid for the coupling between temperature variability and precipitation variability at summertime. During summer, the synoptic-to-subseasonal scale precipitation variation is typically negatively correlated with temperature, especially Tmax, as cloud cover and precipitation apparently reduce incoming shortwave radiation and thus Tmax. Therefore, it comes naturally to speculate that if the precipitation at any location varies strongly, then the Tmax will also exhibit strong variations. We further examined whether the varianceinduced Tmax trends in the YRB region are linked with the strength of intraseasonal variability of precipitation. As we expected, both intra-seasonal Tmax variability ( Figure 6A) and precipitation variability ( Figure 6B) in summer show a very similar pattern of robust increasing trends, and only in the YRB region. This also implies that the variability of summer precipitation in the YRB region is plausibly the physical essence of variance-induced increasing trends of Tmax. The factors that control the combined intra-seasonal variability of temperatureprecipitation over the YRB region, however, require detailed exploration in subsequent studies. Nevertheless, although a quantitative and deterministic analysis for the physical mechanisms is beyond the scope of our current study, some evidence regarding the causes of changing temperature variability has been discussed in several recent studies (e.g., Hsu et al., 2017;Zhu et al., 2020). For instance, Hsu et al. (2017) demonstrated that the boreal summer intraseasonal oscillation could modulate the occurrence of the heatwaves for the monsoon regions. In particular for the YRB region, the 10-30 days boreal summer intraseasonal oscillation is closely linked to the increases in heatwave occurrence. For modeling evidence, a growing number of studies have focused on the dynamic causes (i.e., thermal advections, circulation changes, and land-atmosphere interactions) of the increasing temperature variability over Europe in the context of global warming (e.g., Andrade et al., 2012;Douville et al., 2016;Holmes et al., 2016). These results provide constructive clues and insights for the physical causes of the covariability of Tmaxprecipitation over the YRB region. Frontiers in Earth Science | www.frontiersin.org December 2021 | Volume 9 | Article 762563 8 In addition, it is worth noting that the "trends" of Tmax and Tmin in our study only represent the observed trends during the past 57 years. Whether these quantile features of "trends" are long-term trends or just an increasing phase of some multidecadal scale variations is unknown due to the limitation of observation record. From another perspective, it is more concerned that the "trends" is a signal (from anthnographic forcings) or a noise (internal variability of the climate system). In future studies, we can use the model results (e.g., CMIP5/6 simulations) to investigate such quantile features of Tmax trends or the changing variability as well as the possible mechanisms, particularly for the YRB region. Some modeling evidence has been highlighted regarding the variations of temperature variability in several studies (Gregory and Mitchell, 1995;Rowell, 2005;Fischer and Schär, 2008;Holmes et al., 2016;Chan et al., 2020). For example, Chan et al. (2020) revealed that there are strong couplings between warming and increasing variability in various CMIP5 models in the context of RCP8.5 projections, especially for mid-latitude regions (e.g., Europe, southeastern China). This also implies that higher warming of mean temperature probably accompanies by a higher increase in temperature variability in the high-emission scenario. The compound effects of increasing mean and variability of temperature could significantly increase the risk from potentially deadly heatwaves (Mora et al., 2017;Baldwin et al., 2019).
In summary, the unique quantile features of temperatures' trends in summer over China are shown in the current study. The most important added value of our study is identifying a variancedominated pattern of Tmax trends for the YRB region compared with the previous mean state shift-induced warming pattern. Such compounded increasing variability and the uniform global warming in the mean state could subsequently cause parts of the world exposed to increasing risks of the potential deadly heatwave (Chan et al., 2020). We highlighted that the anthropogenic contributions and the physical mechanisms of the increasing temperature variability should be paid more attention to in future works in order to better plan strategies of adaptation and mitigation solutions for the possible upcoming more intense extremes and disasters.

DATA AVAILABILITY STATEMENT
The data analyzed in this study is subject to the following licenses/ restrictions: We acknowledge China Meteorology Administration (CMA) for providing the homogenized daily maximum and minimum surface air temperature and precipitation data of China. Requests to access these datasets should be directed to China Meteorological Data Service Centre, http://data.cma.cn/en.