Introduction

It is unequivocal that our planet is warming, as authenticated by both observed records and climate models (e.g., Collins et al. 2013). Even the COVID-19 pandemic, which occurred along with a sudden reduction in greenhouse gas emissions, has not seemed to change the trend of global warming (Forster et al. 2020). As a vital component of the climate system, precipitation also changes in its intensity and frequency with rising temperatures, thus bringing challenges to watershed adaptation/elimination management (Miralha et al. 2020; Zhang et al. 2019a). Moreover, recent attention has been given to the temporal unevenness of precipitation because this unevenness can change the potential risks, e.g., both the rainfall duration and antecedent dry period can influence the first-flush effect of nonpoint source pollution (NPS; e.g., Chen 2020; Feng et al. 2013).

NPS refers to accumulated pollutants (e.g., fertilizers) in a watershed entering the downstream receiving water body in a wide, dispersed, and trace form under the flushing and transport of precipitation, thus leading to eutrophication, damaging the aquatic environment, and threatening water security (Miralha et al. 2020; Mosley 2015). Precipitation plays an important role in the processing of NPS, as it is not only the key driving force of NPS (in the way of precipitation runoff and/or soil erosion), but it can also affect the physicochemical properties of soils (e.g., moisture content, pH, C/N ratio, temperature), reshaping the surface water flux, groundwater flux, and ecohydrological flux accordingly (Ouyang et al. 2019; Shortridge 2019). Precipitation is also the dominant input of both physics-based and statistics-based NPS models, indicating the importance of precipitation characteristics to NPS (Miralha et al. 2020; Varekar et al. 2021). Although the precipitation amount is considered a preferred feature of precipitation, precipitation unevenness, which describes the extent to which falling precipitation is asymmetrical over time, can also affect NPS (Nobre et al. 2020; Qiu et al. 2018). For instance, droughts can disrupt nonpoint source pollution processes and increase the internal processes of the nutrient cycle but also provide accumulated nonpoint nutrients for post-drought rainfall (Mosley 2015). Therefore, it is meaningful to study future changes in precipitation unevenness in addition to changes in the precipitation amount before conducting future NPS simulations.

Global climate models (GCMs) are major sources of information for climate projections (Guo et al. 2020a). Nevertheless, compared to observations, GCMs exhibit systematic biases due to the imperfect understanding of physical climatic processes, deficient parameterizations of subgrid processes, misrepresentations of ocean/land-atmosphere feedbacks, and approximation of numerical solutions (Nguyen et al. 2019). Although regional climate models (RCMs) nested in GCMs improve the resolution of GCMs, systematic biases still remain, limiting the direct use of raw GCM or RCM outputs. To reduce these systematic biases and provide more realistic climate simulations, bias correction (BC) techniques have become an imperative postprocessing procedure for almost all watershed climate impact studies (Guo et al. 2019; Nguyen et al. 2020; Zhang et al. 2019a). Statistical BC techniques are widely used in climate change research; in these methods, the statistical relationships between observed climate variables and their corresponding historical simulations are transplanted to the projected periods of the model (Albek et al. 2019; Ouyang et al. 2019; Shalby et al. 2020). Among statistical BC methods, evolving quantile mapping-based methods (QMs) show impressive performances, as they can be easily implemented to adjust not only the mean and variance but also any quantile of the variable of interest (Wang and Chen 2014). However, almost all QMs are univariate BCs that only reproduce the marginal distributions of climate variables. In other words, QMs cannot correct the dependence structure among variables, time scales, and stations (Li and Babovic 2019; Mehrotra and Sharma 2019; Vrac 2018).

To reconstruct this dependence structure, a few different multivariate BCs (MBCs) have been proposed based on QMs over the past two decades. Among these MBCs, rank-order sampling methods have been noticed due to their simple rules and ease of implementation. Vrac and Friederichs (2015) introduced a real-world dependence structure successfully by applying the Schaake shuffle (Clark et al. 2004) method directly to QM outputs. Similarly, by combining the distribution-free shuffle approach, Guo et al. (2019) and Eum et al. (2020) rebuilt the observed precipitation-temperature dependence by reordering the temperature sequence based on the intervariable correlations of the observed time series. Moreover, Vrac (2018) expanded the empirical copula to a multidimensional version (e.g., multisite and/or multivariate) to boost the stochasticity and variability in the corrected results. However, the above methods more or less assume that the projected dependence structure is consistent with that of the observations, as these methods employ either an observed copula function or Cholesky decomposition. Recently, Mehrotra and Sharma (2019) proposed a new rank-order sampling method for lessening the stationary hypothesis of the dependence structure and named it the three-dimensional BC (3DBC). 3DBC can reflect the evolution of the dependence structure that might occur in the future indicated by climate models. Thus, this method can provide more approximate realistic material for future climate impact studies; however, it is rarely used.

Large-scale atmospheric circulation patterns and local physical features are the main factors affecting precipitation unevenness. Atmospheric circulation patterns represent background fields for studying regional precipitation unevenness, while local physical features such as the precipitation grade structure and near-surface temperature play direct roles in regional changes in precipitation unevenness (Dey and Mujumdar 2019; Rajah et al. 2014; Wang et al. 2019). Although different indicators are selected to characterize precipitation unevenness, such as information entropy (Feng et al. 2013; Konapala et al. 2020), the Gini index, or the concentration index (Rajah et al. 2014; Shortridge 2019; Wang et al. 2019), nearly all these indicators are closely related to rainfall intensity. Hence, from a regional perspective, large-scale circulation patterns adjust the regional precipitation unevenness by modifying the precipitation grade structure. Most studies have focused on the impact of extreme precipitation due to its disproportionate consequences, but low-intensity precipitation (LIP), which impacts the carbon and nitrogen cycles in terrestrial ecosystems, such as crop growth (Shortridge 2019) and ecosystem evolution (Romero et al. 2020), can also regulate precipitation unevenness (Rajah et al. 2014).

To obtain prior information on NPS from precipitation, this study investigated future changes in both amount and unevenness of precipitation under a more approximate real-world configuration using 3DBC and then quantified the relationship between LIP/high-intensity precipitation (HIP) and precipitation unevenness conditioned on rising temperatures. The objectives of this study are as follows: (1) establishing a two-stage BC model with a corrected dependence structure, (2) projecting the dynamic future changes in precipitation, and (3) quantifying the influence of LIP/HIP on the unevenness of the rainfall distribution conditioned on temperature.

Study area and data

Study area

The Daqing River Basin (DQRB) is one of the five subbasins in the Hai River Basin (Fig. 1). The drainage area of the DQRB is approximately 43,060 km2 spanning between 113.66–117.72° E longitude and 38.09–40.08° N latitude (Li and Babovic 2019). It is mountainous in the western region, hilly in the central area, and consists of plains in the eastern region. The easternmost part is adjacent to Bohai Bay. Controlled by tropical and midlatitude climate systems, the precipitation in the DQRB varies significantly, making it one of the areas with the most uneven seasonal variations (Guo et al. 2020b; Konapala et al. 2020; Zhang et al. 2019b). The DQRB has a typical temperate continental monsoon climate with annual precipitation of 500–700 mm; approximately 80% of the precipitation occurs during the months of June–September. The plain portion is densely populated, including parts of Beijing and Tianjin, the first- and second-largest cities in North China, and is an important grain production base with high fertilizer application. Furthermore, an above-average warming rate compared to other regions of China and rapid urbanization makes the DQRB attractive for studying future climate change (Ren et al. 2012).

Fig. 1
figure 1

Position of the a Hai River Basin in China, b DQRB in the Hai River Basin, and c weather stations and HadGEM3-RA grids used in the DQRB. The elevation data was sourced from 90-m SRTM data (http://www.gscloud.cn/).

Full size image

Data

Similar to Miralha et al. (2020), this study used both observed and climate model-simulated daily precipitation (Pr), maximum temperature (Tmax), and minimum temperature (Tmin) data considering their availability, robustness, and impact on NPS. The observed data were taken from the China Meteorological Science Data Sharing Service Network (http://data.cma.cn) and covered the period from 1981 to 2015. A total of five standard weather stations within the DQRB, as well as five stations situated around it, were selected, and all the data were quality-controlled (Table S1). The climate model simulations were taken from the Coordinated Downscaling Experiment-East Asia (CORDEX-EA) project (http://cordex-ea.climate.go.kr/cordex). The CORDEX-EA project provides five RCMs for the DQRB that are driven by boundary and initial conditions obtained from historical HadGEM2-AO simulations. Focusing on precipitation, Park et al. (2016) found that the Hadley Centre Global Environmental Model version 3 regional climate model (HadGEM3-RA) has a slight advantage in simulating precipitation extremes in East Asia among CORDEX-EA historical runs. Using a score-based method, Zhang et al. (2019a) also declared that HadGEM3-RA is better than other CORDEX-EA runs in simulating precipitation as well as temperature. Moreover, HadGEM3-RA provides a longer output (1981–2100) than other RCMs, meaning that climate change signals can be captured to a large extent. Therefore, given that the results of BC seriously depend on the RCM outputs, HadGEM3-RA was selected as the raw model in this study due to its good performance simulating precipitation extremes and its long-time-series simulations. HadGEM3-RA was developed by the National Institute of Meteorological Research, Korea, and has a spatial resolution of 0.44°. For future projections, two experiments are run based on the Representative Concentration Pathways (RCPs) 4.5 and 8.5. A comprehensive introduction to HadGEM3-RA can be found in Davies et al. (2005).

The HadGEM3-RA grid cells closest to the weather stations were selected for the simulations (Fig. 1). Additionally, the projected periods of both RCP4.5 and RCP8.5 were divided into near-future (from 2031 to 2065, denoted as a) and far-future (from 2066 to 2100, denoted as b) periods as 3DBC worked on equal-length series. It should be noted that the historical simulation of HadGEM3-RA ended in 2005; accordingly, the HadGEM3-RA RCP4.5 experiment from 2006 to 2015 was treated as an extended historical simulation after two-sample Kolmogorov-Smirnov tests were conducted for each station. The results showed that all extended historical simulations had the same distribution as the historical simulations at the significance level of 5%. Thus, the sequences from 1981 to 2015 were regarded as the current simulated period. Additionally, HadGEM3-RA has a 360-day calendar per year. Although this does not disturb the QM process, five extra days were added from the 36th day every 72 days in a year for the following calculations (Vanderkelen et al. 2018).

Methodology

Description of two-stage BC

As run-free methods, QMs only statistically adjust the distribution of raw simulations to that of observations under the stationary hypothesis and cannot correct unreasonable physical processes in simulations (Chen et al. 2020; Maraun et al. 2017). Despite these limitations, QMs have been applied in many studies to reduce the biases of raw models (Eum et al. 2020; Guo et al. 2020a; Li and Babovic 2019). In this study, a two-stage BC was applied to HadGEM3-RA for the Pr, Tmax, and Tmin data; in the first stage, equidistant cumulative distribution function matching (EDCDFm), a QM variant, was used to adjust the marginal distributions, and in the second stage, 3DBC was employed to reproduce the dependence structure. Brief introductions to these methods are given below, and their specific steps are displayed in the Supplementary Material (Fig. S1).

EDCDFm assumes that the difference between observations and simulations at any percentile during the reference period can also apply to the projected period (Li et al. 2010; Wang and Chen 2014). Thus, the problem of extrapolating beyond historical values can be avoided; moreover, the nonstationarity relative to the reference climate is also adaptive to some extent. The EDCDFm can be written as follows:

where x denotes the daily values for the climate variable; F and F−1 represent the cumulative distribution function (CDF) and the inverse CDF of the climate variable, respectively; subscripts o and m denote the observation and model, respectively; and subscripts f and c denote the future and current periods, respectively.

The accuracy of the EDCDFm results has a vital impact on the second stage, as these results are the inputs of 3DBC. Thus, considering seasonal differences, a 31-day sliding window was adopted, and the nonparametric method called diffusion-based kernel density estimation (Botev et al. 2010) was applied to fit the CDFs of Pr, Tmax, and Tmin after comparison with the gamma, Weibull, and traditional kernel density functions. Specifically, for Pr, daily values greater than 0.02 mm were regarded as rainy days to remove the drizzle effect; the singularity stochastic removal technique (Vrac et al. 2016), which replaces the nonrainy days with tiny values, was used to match the occurrences. A fivefold cross-validation, with each fold containing seven consecutive years (1981–1987, 1988–1994, 1995–2001, 2002–2008, and 2009–2015), was considered for the validation.

3DBC is an improved rank-order sampling method. It holds the advantage of the traditional empirical copula method: by assuming the weather system of each calendar day is similar, 3DBC performs independently in each dimension. This means that 3DBC can be applied to any number of variables, stations, and time scales. Most importantly, replacing the first-order coefficient of autocorrelation (Lag 1) of the corrected QM rank space with that of projected simulations, 3DBC incorporates the dynamic evolution of the dependence structure as indicated by the projected climate model.

To validate the two-stage BC model, the corrected outputs of the current period (1981–2015) were compared with the corresponding observations, including the mean, standard deviation (Std), Lag 1, intervvariable correlations (VC), and intersstation correlations (SC). Although the two-stage BC was performed at the daily scale, the monthly and annual results aggregated by the corrected daily outputs were also analyzed considering low-frequency variability and subsequent studies. The evolution of the dependence structure between the current period and the scenario that likely had the most discrepancy, RCP 8.5b, was examined by comparing their theoretical and practical calculations using the mean absolute error (MAE) (Equations 3 and 4):

$$ MAE=\frac{1}{N}\sum \limits_{n=1}^N\left|{R}_{BC}^s-{R}_n^s\right| $$
(3)
$$ {\displaystyle \begin{array}{l}{R}_n^s={R}_o^s+{R}_f^s-{R}_c^s\\ {}\kern1em =\frac{\frac{1+{R}_f^s}{1+{R}_c^s}\left(1+{R}_o^s\right)-\frac{1-{R}_f^s}{1-{R}_c^s}\left(1-{R}_o^s\right)}{\frac{1+{R}_f^s}{1+{R}_c^s}\left(1+{R}_o^s\right)+\frac{1-{R}_f^s}{1-{R}_c^s}\left(1-{R}_o^s\right)}\end{array}} $$
(4)

where N is the number of stations; Rns denotes the theoretical corrected future correlation following Bárdossy and Pegram (2012); \( {R}_{BC}^s \) represents the corresponding practical corrected correlation; Ro, Rc, and Rf represent the statistics of observations, corrected current simulations, and corrected future simulations, respectively; superscript s represents statistics and can be Lag 1, VC, or SC; and subscript BC denotes the correction method, which can be QM or 3DBC.

Daily dry/wet spells, which refer to continuous days with or without rain, were selected to assess the abilities of the two BC methods to reproduce intermittent characteristics (Li et al. 2016b). The spatial continuity ratio (CR), which reflects the spatial intermittence of the rainfall field, was also taken as an intermittent indicator (Wilks 1998). It reveals the link between precipitation amounts and occurrences by calculating the ratio of the precipitation mean at station i depending on whether station j is wet or dry (Equation 5). The CR is small for highly correlated stations and large for weakly correlated stations.

$$ CR=\frac{E\left({\mathit{\Pr}}_i\left|{\mathit{\Pr}}_i>0\right..02\cap {\mathit{\Pr}}_j\le 0.02\right)}{E\left({\mathit{\Pr}}_i\left|{\mathit{\Pr}}_i>0\right..02\cap {\mathit{\Pr}}_j>0.02\right)} $$
(5)

Relative entropy

Introduced by Feng et al. (2013), different information theory metrics have been used to quantify precipitation unevenness (Konapala et al. 2020; Pascale et al. 2015). Relative entropy (RE), defined as the distance between the daily/monthly rainfall distribution and the uniform distribution in 1 year, is used in this study (Dey and Mujumdar 2019). RE is a nonparametric measurement that can capture higher-order statistics and deal with non-Gaussian distributions such as precipitation distributions. Notably, the autocorrelation of precipitation is not considered for RE. In the following analysis, the RE determined by daily precipitation (RED) is used to articulate the annual precipitation unevenness because it contains differently graded rainfall data, and the RE determined by monthly precipitation (REM) is supplemented to determine the rainfall seasonality. RE can be expressed as follows:

$$ RE=\sum \limits_{k=1}^{1/U}{p}_k{\log}_2\left(\frac{p_k}{U}\right) $$
(6)

where pk is the proportion of the precipitation of the k ‐ th day/month to the annual precipitation and U equals 1/365 or 1/12 for daily or monthly RE, respectively.

According to the left and right limits of RE, its domain is [0, log2(1/U)]; a higher RE value denotes higher unevenness, while a lower RE value denotes lower unevenness. RE only considers changes in proportion, providing a relatively independent dimension for precipitation unevenness. Specifically, if the annual effective rainy days/months, Nrainy, is considerably different from zero, it can be calculated using RE as Nrainy = 1/U · 2−RE (Pascale et al. 2015).

Partial correlation coefficient

It is of great significance to survey the relationship between RE and LIP/HIP to reveal the dynamics of precipitation unevenness in the DQRB. Ignoring the effect of temperature may bias the results because the temperature influences the intensity and frequency of rainfall. In this study, the partial correlation coefficient (Kendall 1977), which measures the dependence of two random variables from several interrelated random variables, is used to assess the correlation between RE and LIP/HIP in the context of changing temperatures. The partial correlation coefficient can be expressed as follows:

$$ {\rho}_{1,2\left|3,\dots, m\right.}=-\frac{P_{12}}{\sqrt{P_{11}{P}_{22}}} $$
(7)

where Pij denotes the algebraic cofactor of ‖ρij‖m × m constituted by the Spearman rank correlation matrix of m random variables.

Results

BC validation

Evaluating the current-day performance of 3DBC is conducive to providing reliable constrained projections that are at least partially explained (Padrón et al. 2019). Fig. 2 shows the observed, raw, and bias-corrected daily, monthly, annual statistics compared with the corresponding observed statistics. In Fig. 2, each point represents one station or one pair of stations for SC. All variables were scaled for displaying purposes following the rule under which the deviation of real data is not visually reduced. Considering that the observations were truncated, the Pearson correlation coefficient was used to calculate the correlation degree of the rainfall data (Zhou et al. 2019). It can be seen that 3DBC and QM can correct the mean values equivalently (Fig. 2a–c), as well as the daily Std values (Fig. 2d). This is not surprising, as these statistics are rank-order independent. Two overall advantages of 3DBC were observed compared with QM: (1) 3DBC can not only correct the daily scale bias (column 1 of Fig. 2) but can also correct the monthly and annual-scale biases (columns 2–3 of Fig. 2); and (2) 3DBC can further improve the dependence structure of the corrected QM outputs (Fig. 2g–o). These results were expected, as 3DBC is designed to correct additional dependence structures. Furthermore, for situation in which the bias of HadGEM3-RA was very small, the results of the two BCs also showed little difference. It was noted that the performance of 3DBC was not better than that of QM or even of HadGEM3-RA for several situations. Similar phenomena were observed in Nguyen et al. (2020) and may have been caused by natural variability and/or model sensitivity (Chen et al. 2020). Overall, the results showed that 3DBC is able to reproduce the first- and second-order observation moments at daily, monthly, and annual scales as well as dependence structures.

Fig. 2
figure 2

Scatter plots of observed statistics (x-axis) and HadGEM3-RA-, QM-, and 3DBC-simulated statistics (y-axis) for the Pr, Tmax, and Tmin data from 1981 to 2015. The statistics are the a–c mean, d–f standard deviation, g–i first-order coefficient of autocorrelation (Lag 1), j–l intervvariable correlations (VCs), and m–o intersstation correlations (SCs). The left, middle, and right columns represent daily, monthly, and annual scales, respectively. Each point denotes a site or a pair of sites for SC.

Full size image

Intermittence is a characteristic of daily precipitation that differs from other continuous climate variables, and it has a significant effect on organisms from lower trophic levels and on plant growth (Romero et al. 2020). Intermittent comparisons of the corrected, raw, and observed results are displayed in Fig. 3. The total numbers of dry/wet spells with different durations at all stations over the period from 1981 to 2015 were calculated, and the results showed that 3DBC (MAE=561) performed slightly better than QM (MAE=500) in modeling the dry spells; both of these methods were able to reproduce the observed dry spells (Fig. 3a). Wet spells of both short and long durations were better modeled by 3DBC than by QM (Fig. 3b). The CR values for all station pairs were used to evaluate the ability of the two BC methods to correct the spatial intermittence (Fig. 3c). It was evident that both 3DBC and QM added considerable value to the raw rainfall field; specifically, 3DBC performed better than QM and only had a slight overestimation compared to the observations. As shown in Fig. 2, this may be partly explained by the fact that 3DBC can better correct all temporally scaled Lag 1 values similarly to the first-order Markov model in many weather generators, which represents precipitation occurrences. The analysis shown in Fig. 3 reveals that 3DBC is capable of depicting the intermittent characteristics of daily precipitation.

Fig. 3
figure 3

Comparison of simulated (by HadGEM3-RA, QM, and 3DBC) and observed (Obs) total a dry and b wet spell counts at all stations and c the continuity ratios (CRs) of all station pairs from 1981 to 2015. The MAE values with different colors represent the mean absolute count errors between the simulated spells (by HadGEM3-RA, QM, and 3DBC) and the observed spells.

Full size image

The main strength of 3DBC lies in its capability to project climate variables with a dynamic, evolving dependence structure. Assuming that the statistical data are additive in nature and the bias and climate change signals are independent of each other without historical biases for corrected sequences, the practical corrected dependence structure changes can be calculated (Mehrotra and Sharma 2019). A comparison of the practical and theoretical dependence structures obtained using the QM and 3DBC outputs for RCP8.5b is exhibited in Fig. 4. The 3DBC method presented a relatively beneficial performance in the daily, monthly, and annual Lag 1 values of Pr, Tmax, and Tmin, except for the annual Pr values (Fig. 4a). Furthermore, 3DBC generally led to improved estimates of VC change signals, especially for annual Pr-Tmin. Additionally, for the daily Pr-Tmax, monthly Pr-Tmax, and monthly Pr-Tmin VCs, QM had an advantage over 3DBC, but the differences were small, and their MAE values were all less than 0.08 (Fig. 4b). The results also showed that 3DBC provided lower MAE values for the SCs than QM with respect to all variables and all temporal scales (Fig. 4c). Thus, it can be concluded from Fig. 4 that 3DBC performs better than QM in capturing changes in the dependence structure.

Fig. 4
figure 4

Comparison of practical and theoretical (a) Lag1; (b) VC; and (c) SC using the QM and 3DBC outputs between the base scenario (1981–2015) and RCP8.5b scenario.

Full size image

Spatiotemporal variation in future precipitation

Providing projected climate information that is close to real conditions is beneficial for watershed management (Li et al. 2019b); thus, the outputs of 3DBC described in the “BC validation” section were further analyzed. Both the spatiotemporal variations in the rainfall amount and the precipitation unevenness variations were investigated; the historical observations (1981–2015) were taken as the base scenario, and the near-/far-future RCP4.5/8.5 simulations were taken as the projected scenarios. Comprehensive descriptions of the annual average features of precipitation and temperature at each station under different projected scenarios are listed in Table S2.

Diagrams describing the interannual variations in annual precipitation and RED in the DQRB under different scenarios are presented in Fig. 5. The annual precipitation output in the four projected scenarios increased compared with the base period. The increase was largest under RCP8.5b, at nearly 30% higher than the base period value. The annual precipitation in the other three scenarios ranged from 480 to 620 mm, showing little difference. Under all projected scenarios, precipitation experienced a slightly decreasing trend (p > 0.05) similar to that observed during the base period using the linear regression model (Fig. 5a). Although the precipitation increased, the RED values of all four projected scenarios ranged from 3.52 to 3.58, showing a slight decrease compared with the base scenario (p > 0.05). Furthermore, the differences in RED among all scenarios were small, which also meant that the differences in effective rainy days were small. In the base scenario, the precipitation unevenness trend was negative, which was also found by Guo et al. (2020b) despite their use of different indicators. A similar finding was obtained by Chen (2020), who observed that the trend was not significant (p> 0.1). Contrary to the trend in the base period, slightly positive trends (p > 0.05) were observed in all projected HadGEM3-RA scenarios (Fig. 5b). These small differences indicated the necessity of analyzing the spatial distribution of precipitation unevenness.

Fig. 5
figure 5

Interannual variations in annual a precipitation and b RED in the DQRB under different scenarios. Each sector represents a scenario, and the red dotted lines denote linear slopes (mm/10a for precipitation/10a for RED). RED daily relative entropy.

Full size image

With global warming, rain belts may be redistributed, which may cause changes in spatial precipitation patterns. Fig. 6 shows the spatial differences in precipitation and RED calculated by inverse distance interpolation under the typical RCP8.5 scenario. The spatial pattern of the base period precipitation is scarce in the northwest and abundant in the southeast, and the precipitation is greatest in the western and eastern marginal regions (Fig. 6a). Compared with the base period, the spatial patterns of precipitation under the projected scenarios exhibited no remarkable change (Fig. 6b–c). Consistent with the results of Yang et al. (2019), the precipitation of the whole basin increased; the increased range of precipitation in the northwestern mountainous area was 0–80 mm, and that in the southeast plain was 100–211 mm. Similar to the spatial distribution of the precipitation amounts, the spatial distribution of RED also presented the characteristics of mountain and plain differentiation; the RED values in the main western mountainous area were between 2.92 and 3.50, while these values were between 3.50 and 3.92 in the eastern plain area (Fig. 6d). Under the projected scenarios, the precipitation unevenness decreased in most areas of the DQRB, which was consistent with Fig. 5, while the precipitation unevenness in the western marginal of the basin increased slightly (Fig. 6e–f). However, annual-scale research may not be sufficient to explain the change characteristics of the DQRB; thus, monthly precipitation variations in the DQRB were further studied.

Fig. 6
figure 6

Spatial distributions of a–c annual precipitation and d–f RED between the base scenario and typical RCP8.5 scenario. The first column represents the means of the base period; the second and third columns represent the typical anomalies under the RCP8.5 scenario for the near- and far-future periods, respectively.

Full size image

The general variations in monthly precipitation under different scenarios are presented in Fig. 7. Since monsoons may arrive at the DQRB within a timespan of 1–2 months, the 12 months within each year were reclassified according to their precipitation amounts. Months with rainfall values > 75% of the 12 months were named extreme wet months (EW); those with rainfall values ≤ 25% were defined as extreme dry months (ED); and the months that fell within the 25–50% and 50–75% percentiles were divided into normal dry months (ND) and normal wet months (NW), respectively. To show the interannual dynamics of precipitation, the months within each year were sorted in descending order from the most precipitation to the least (Fig. 7a–c). Taking RCP4.5b and RCP8.5b as examples, the ED became drier, and the EW became wetter under the projected scenarios compared with the base period. The rainfall in EW accounted for more than 70% of the annual precipitation under all projected scenarios. These results also showed that the maximum monthly precipitation trend in EW and all ED months decreased under the projected scenarios (including under the scenarios that are not shown). As floods and droughts tend to occur in EW and ED months, the frequency of calendar months marked as ED/EW in each year under all scenarios is shown in Fig. 7d–e. The projected ED tended to occur in December and January, while EW tended to occur between June and August. This indicated the likely expansion of precipitation seasonality in the future, as was supported by the REM results (Fig. 7f). Increased precipitation seasonality is often accompanied by high probabilities of flood and drought risks. As a supplement, changes in the daily precipitation grade structure showed that the frequency of extreme precipitation increased, while the number of rainy days and the frequency of light rain (LR, < 10 mm/day) underwent slight changes except under RCP8.5b (Fig. S2). The enhanced precipitation seasonality (Fig. 7) poses challenges to the water resource management of the DQRB.

Fig. 7
figure 7

General variations in monthly precipitation from the a observations (1981–2015) and under b RCP4.5b and c RCP8.5b. Each column represents the monthly series in descending order from greatest to least in a year. EW, NW, ND, and ED are the abbreviations of extreme wet, normal wet, normal dry, and extreme dry months, respectively. The left/right arrow represents negative/positive linear trends, and the single asterisk and triple asterisks represent p < 0.05 and p < 0.001, respectively. The frequencies of calendar months marked as d ED and e EW of each year are compared, as are the f yearly average REM values under all scenarios. The error bar represents the standard deviation of REM, the monthly relative entropy.

Full size image

The dry/wet spell characteristics of the projected scenarios are shown in Fig. 8 and Table S3. All scenarios and almost all stations displayed no significant trends (p < 0.05) with regard to the mean annual total duration of wet spells (Table S3). Based on Fig. 3a–b, the total number of different dry/wet spells at all stations under each projected scenario is shown in Fig. 8. As seen in Fig. 8, 2–7-day dry spells played a dominant role in the DQRB. Compared with the base period, the number of dry spells with durations of 1 day and > 29 days increased, while the other durations decreased. The 1-day and 2–7-day wet spells were dominant in the base period, but their frequencies in the future decreased except under RCP8.5b, while that of ≥ 8-day wet spells increased under all scenarios. The long dry/wet spells indicated that the risk of drought and flooding in the DQRB may increase in the future (Li et al. 2016b).

Fig. 8
figure 8

Comparison between the total numbers of different (a) dry spells and (b) wet spells at all stations under the projected scenarios and during the base period.

Full size image

The relationship between LIP/HIP and precipitation unevenness

It is meaningful to study the relationship between LIP/HIP and RED under a warming climate because their changing frequencies and intensities respond to rising temperatures. In this study, annual accumulated rainfall amounts greater than the 95th quantile of the rainfall sequence and the annual accumulated LR are defined as HIP and LIP, respectively (Chen 2020). We still took winter rainfall into account because winter rainfall events are important for winter wheat, the main crop of the DQRB. As the amount and frequency of snowfall are very small in the study region, the difference between snowfall and rainfall is ignored. A comparison of LIP and HIP between the 3DBC-derived results and observations during 1981–2015 is presented in Fig. S3. It is shown that 3DBC simulated LIP/HIP well. Because the average annual Tmax and Tmin values at all stations were highly correlated (r > 0.88), their means were assumed to be the annual average temperature T. This also means that the T used in this paper was a local value rather than a global value, as our purpose was to consider the influence of the regional climate system more on the precipitation structure than on the trend. In the current study, the future temperature of the DQRB is argued to be higher than that of the base scenario under all future scenarios (Table S2). In the following analysis, we interpolated the station data to a 0.25°×0.25° resolution to enhance the spatial balance; accordingly, a total of 64 grids were obtained (De Caceres et al. 2018).

The partial correlation coefficients between RED and HIP/LIP conditioned on T, i.e., \( {\rho}_{RE^D, HIP\left|T\right.} \) and \( {\rho}_{RE^D, LIP\left|T\right.} \), are shown for all grids under different scenarios in Fig. 9a–b. The \( {\rho}_{RE^D, LIP\left|T\right.} \) values were less than −0.2, while the \( {\rho}_{RE^D, HIP\left|T\right.} \) values tended to be positive (except under RCP4.5b) but also had a certain small probability of being negative. The proportions of significant grids (p < 0.05) were 78–100% for \( {\rho}_{RE^D, LIP\left|T\right.} \); for \( {\rho}_{RE^D, HIP\left|T\right.} \), the proportions under RCP4.5b were 0, and those under the other scenarios were 10–20%. As LIP and HIP can simultaneously influence the RED, the partial correlation coefficients between RED and HIP/LIP conditioned on T and LIP/HIP, i.e., \( {\rho}_{RE^D, HIP\left|T\right., LIP} \) and \( {\rho}_{RE^D, LIP\left|T\right., HIP} \), were further studied (Fig. 9c–d). The absolute values of \( {\rho}_{RE^D, HIP\left|T\right., LIP} \) and \( {\rho}_{RE^D, LIP\left|T\right., HIP} \) were higher than those of \( {\rho}_{RE^D, HIP\left|T\right.} \) and \( {\rho}_{RE^D, LIP\left|T\right.} \), indicating the overall discrepant effects of the LIP and HIP. As additional evidence, the proportions of significant grids (p < 0.05) also increased compared with those of \( {\rho}_{RE^D, HIP\left|T\right., LIP} \) and \( {\rho}_{RE^D, LIP\left|T\right., HIP} \). It was also discovered that the significant grids of \( {\rho}_{RE^D, HIP\left|T\right., LIP} \) and \( {\rho}_{RE^D, HIP\left|T\right.} \) were all positive. However, no regularity was found in the correlation coefficients between the projected and base scenarios, which stressed the complexity of their mechanisms. Very similar results can be seen in station data-based Tables S4 and S5.

Fig. 9
figure 9

Partial correlation coefficient intervals of a RED, LIP | T; b RED, HIP | T; c RED, LIP | T, HIP; and d RED, HIP | T, LIP under all scenarios. The colors of the median circles represent the proportions of significant grids (p < 0.05).

Full size image

Discussion

Providing comprehensive projected climate information that matches the corresponding basin characteristics is helpful for nonpoint source management under rising temperatures. According to the HadGEM3-RA model corrected by 3DBC, the future spatiotemporal variations in precipitation in the DQRB, including the precipitation magnitude and unevenness, were evaluated in this study, as was the relationship between HIP/LIP and precipitation unevenness.

Predictive ability of 3DBC

The validation results for the historical period show that, based on the QM outputs in which only the intensity of HadGEM3-RA was corrected, 3DBC can further correct the overall dependence structure of HadGEM3-RA. Similar to many MBCs mentioned in the “Introduction” section, 3DBC is able to correct the dependence structure of HadGEM3-RA, supporting its applicability in the DQRB (Figs. 2–4). In particular, consistent with the results of Mehrotra and Sharma (2019), the current study showed that 3DBC was able to reflect the evolution of the dependence structure, which was believed to be the result of adapting the Lag 1 values of the climate variables.

However, the preset assumptions and lack of physical processes likely to cause statistical BC, including QMs and MBCs, are disputed in many situations (Chen et al. 2020; Maraun et al. 2017). Nguyen et al. (2019) reviewed these concerns and concluded that statistical BCs were expected to remain attractive to end users to narrow the biases until better methods became available. Another dispute is whether the dependence structure should be stable over time (François et al. 2020). Many researchers have suggested that VCs and SCs change obviously with climate warming (e.g., Guo et al. 2020a). Accordingly, it seems reasonable to consider that the dependence structure is unstable. Nevertheless, from another perspective, SCs are largely constrained by the physical features of special regions, so the dependence structure is likely to be stable over time (Vrac 2018). Even under biased behavior, HadGEM3-RA, which follows simulated physical processes, honestly reflects the response of the studied region to different forcings (Bárdossy and Pegram 2012). Thus, the dependence structure information drew support from HadGEM3-RA in the current study.

In general, 3DBC shows advantages over QM in terms of the dependence structure, which may inspire simulations involving multiple climate variables. Recent studies have suggested that MBCs perform better than QMs in hydrological models (Eum et al. 2020; Guo et al. 2020a). This makes 3DBC attractive for projected simulations of NPS, as 3DBC performs in a dimensionally independent manner and can specifically reflect climate change signals. For example, the mobility of heavy metals is closely related to both temperature and precipitation (Liu et al. 2020; Punia 2021); thus, 3DBC may allow the response of this process to climate change to be modeled. Moreover, 3DBC is likely to function for latitude/elevation-based freeze-thaw processes in watersheds in which freezing influenced by temperature and precipitation changes alter the water conductivity and erosion resistance of the soil (Ouyang et al. 2018). This study also emphasized the necessity of comparing the verification results of QM and MBCs. Our results suggest that using MBCs such as 3DBC sometimes worsens the prediction results. This suggestion may be helpful for analyzing the results of Miralha et al. (2020), wherein QM performed better than MBCs in the NPS model.

Precipitation predictions under studied scenarios

Unlike the general “dry-get-drier, wet-get-wetter” paradigm under a warming climate, precipitation was found to decline during past decades and increase in future decades throughout the DQRB (Figs. 5 and 6), similar to the results obtained in previous studies (Song et al. 2015; Yuan et al. 2019; Li et al. 2016a; Li et al. 2018b). This result is considered to be related to the Western Pacific subtropical high, which determines the intensity of the East Asian monsoon (Wang et al. 2019); i.e., moisture transport and convergence with weak large-scale circulation are insufficient to deliver abundant moisture, while the strengthened and unstable East Asian summer monsoon north of 30° N enhances the availability of water vapor. Future strengthened large-scale circulation patterns are also detected in HadGEM3-RA (not shown). Moreover, along with higher temperatures, an improved atmospheric water vapor capacity can also intensify precipitation (Chou et al. 2012; Li et al. 2016a). The precipitation amount trends under all projected scenarios were found to be similar to that of the baseline scenario (Fig. 5a). Bucchignani et al. (2014) showed that the precipitation trend under the RCP4.5 scenario was negative, but it tended to be positive under the RCP8.5 scenario using a climate model. The precipitation amount trends at the watershed scale were influenced by the considered scale, period, and method and still reflect open questions; thus, we cautiously retain these related results for comparisons with future studies. This does not affect further analysis, as the scenarios are taken as analytical units.

The results suggest that in the future DQRB, precipitation unevenness will decrease slightly in almost all areas, while rainfall seasonality will be strengthened (Figs. 5–7). The RED range in the current area was very similar to that in the central part of the Indian continent (Dey and Mujumdar 2019), which may be because of their similar monsoon characteristics. Generally, accelerating temperature rises will lead to faster development in extreme rainfall than the average, meaning that the precipitation unevenness will intensify (Chen 2020; Collins et al. 2013; Pendergrass and Knutti 2018). However, the current results were inconsistent, and inconsistencies were also found in some long-term observation records in other regions (Guo et al. 2020b; Wang et al. 2019). One possible reason for this result may be that the temporal correlation was not considered in previously used precipitation unevenness indexes. Contrary to the RED results, the REM results predicted that precipitation will fall more unevenly in the future (Fig. 7). This interesting result was also obtained by Guo et al. (2020b), who found that daily precipitation concentrations and monthly precipitation concentrations showed different sensitivities to droughts. This highlights the need to develop new precipitation unevenness indexes based on event risks, which will inevitably result in another layer of complications (Pendergrass and Knutti 2018).

Studying the future characteristics of precipitation can provide some practical implications for NPS simulations as well as some restoration strategies. There is usually a positive correlation between the pollution load and annual precipitation (Qiu et al. 2018); thus, the results of this study imply that the pollution load is likely to increase, especially over the eastern plain region, as shown in Fig. 5a and Fig. 6a–c. However, this linear relationship can be disturbed in part by the temporal distribution of precipitation (Qiu et al. 2018; Varekar et al. 2021). Thus, it is essential to additionally notice that the impact of precipitation unevenness on NPS is event-based. With an increase in rainfall seasonality (Fig. 7), on the one hand, paying more attention to the related sensitive parameters may improve the simulations as the corresponding physical processes change (Xie et al. 2017); on the other hand, the frequencies of long-duration dry/wet spells in the DQRB are expected to increase in the future (Fig. 8). Given that prolonged dry spells increase the hydraulic residence time, which may worsen the water quality, and that prolonged wet spells may make the pollution load concentrate in the receiving water body, the likelihood of eutrophication is expected to increase. Strengthening the monitoring of precipitation and optimizing the best management practices (e.g., rational fertilization) in the watershed, especially in the critical source area during prolonged wet spells, are expected to reduce the pollution load.

Influence of LIP/HIP on RE D

The results of this study indicate that both LIP and HIP can affect the precipitation unevenness conditioned on temperature; in particular, HIP tends to increase the precipitation unevenness degree, while LIP tends to decrease the degree of precipitation unevenness. The RED index is related to the degree of asymmetry between the probability distribution of annual daily precipitation and a uniform distribution, thus focusing on the accumulated contribution of variously graded daily precipitation values to the unevenness. However, the significance between RED and HIP is weak. Several possible factors may explain this outcome. (1) The spatiotemporal distribution of extreme precipitation, which greatly influences RED, is related to elevation (Zheng et al. 2017). Similarly, stations with significant correlation coefficients were mainly concentrated in the southeast region of the study area at low elevations, including station numbers 54606, 54527, and 54623. (2) The Clausius-Clapeyron relation, which states that extreme precipitation is expected to increase by 6–7% per degree increase in centigrade, is restricted by the exponential moisture demand at high temperatures. In the Hai River Basin, Wang et al. (2018) found that the Clausius-Clapeyron relation tended to be negative when the temperature was higher than 23 °C. (3) Soil moisture, which adjusts land-atmosphere interactions and cycles by redistributing near-surface energy and water, may also affect extreme precipitation (Li et al. 2019a). Recent NPS studies regard precipitation unevenness as an influencing factor of NPS (Nobre et al. 2020; Varekar et al. 2021). Therefore, distinguishing the contribution of HIP/LIP to precipitation unevenness may reduce the uncertainty of the export load, considering that precipitation with sufficient energy is the main driving force of pollutant migration.

There are several limitations to this research. HadGEM3-RA was selected in the current study because it had a relatively good historical performance in terms of climate change predictions. Nevertheless, future climate change is so complex that many factors, including raw models, internal climate variabilities, BC methods, interpolations, and dataset resolutions, can influence the prediction results (Ouyang et al. 2019). It is more convincing to apply finer climate models with different realizations to our future work, although only a few RCMs are publicly available. Some studies show that subdaily scale analyses are more consistent with physical processes than studies conducted on daily scale (Chen 2020; Pendergrass and Knutti 2018). Thus, combining disaggregation methods such as segment resampling (Li et al. 2018a) may provide more robust results. Furthermore, teleconnections and sunspot activity were not considered in this study. Taking these factors into consideration would allow a more in-depth study of precipitation unevenness.

Conclusions

Providing rainfall prediction information according to the characteristics of special regions is critical to ensuring the environmental quality of entire watersheds. In this study, first, 3DBC was used to correct the HadGEM3-RA model. Then, using the corrected results, the amount and unevenness of precipitation in the DQRB were analyzed. Finally, the partial correlation coefficient was used to quantify the correlations between LIP/HIP and precipitation unevenness. The main conclusions are summarized as follows.

  1. 1.

    Compared with univariate QM, two-stage 3DBC can not only correct the means and standard deviations output by HadGEM3-RA but can also generally correct the dependence structure as a whole and retain the evolution information of the dependence structure.

  2. 2.

    The precipitation of the DQRB is expected to increase from 495.3 mm in the base period to 555.9–642.6 mm under the future scenarios. The spatial heterogeneity of precipitation in the DQRB is expected to increase in the future, and the precipitation in the central and eastern regions will evidently increase.

  3. 3.

    The future precipitation unevenness of the DQRB is projected to decrease slightly from 3.60 during the base period to 3.52–3.58 under the projected scenarios, and precipitation is expected to become more concentrated between June and August. Spatially, the future heterogeneity of the precipitation unevenness distribution is weakened, with almost all areas experiencing declines except the western margin of the study area.

  4. 4.

    Both HIP and LIP can adjust precipitation unevenness conditioned on temperature; HIP tends to increase unevenness, while LIP tends to reduce unevenness.