Projecting Spring Consecutive Rainfall Events in the Three Gorges Reservoir based on Triple-Nested Dynamical Downscaling

The Three Gorges Reservoir area (TGR) is a mountainous and highly populated region located in the middle of the Yangtze River basin, central China (Fig. 1a). It often suffers from geological hazards such as collapses and landslides, which block the rivers running to the reservoir, reduce the effective storage of the reservoir and result in huge damages to properties and lives. Consecutive rainfall events (CREs), during which it rains for one week and even longer with a moderate intensity, are powerful triggers for such geological hazards (Corominas and Moya, 1999; Jibson, 2006; Zheng et al., 2020). Therefore, a credible projection of the future change in CREs is of vital importance in evaluating and mitigating the impacts of geological hazards in the TGR.

Figure 1.  (a) Topography of the TGR. The south and north sections are the Wuling Mountains and the Daba Mountains, and the east and west sections are Jianghan Plain and Sichuan Basin. The red dot upstream is Chongqing City, and the purple dot downstream is the Great Dam. (b) The three domains for WRF’s downscaling: the outermost domain has a horizontal resolution of 36 km (Domain1; black box); the middle domain has a horizontal resolution of 12 km (Domain2; red box); and the innermost domain has a horizontal resolution of 4 km (Domain3; yellow box).

Rainfall features will change in their frequency, duration, and intensity in a warmer climate. Of particular importance is the regionalization and intensification of rainfall. An example of this is the increased trends in rainfall amount and intensity observed in East and Northwest China and the decreased trends in the banded southwest–northeast-orientated region (Zhai et al., 2005; Zhou et al., 2016). This can lead to changes in CREs to some extent. Zheng et al. (2020) analysed the future trend of CREs in the TGR using the global climate model (GCM) outputs of phase 5 of the Coupled Model Intercomparison Project (CMIP5) and found a significant increase in the occurrence and intensity of CREs in spring under different greenhouse gas emission scenarios. However, the projection is still fraught with uncertainties: 18 out of the 20 CMIP5 models overestimated the occurrence frequency of CREs with a mean bias of 60.8%; and models differ significantly in terms of the trend and spatial distribution of CREs (Zheng et al., 2020).

One of the major sources of these uncertainties comes from the structural differences among GCMs. Since different GCMs utilize different parameterizations and schemes to account for processes that cannot be resolved by the GCMs’ coarse grids, they often yield significant differences in their rainfall simulations (Chen et al., 2020; Sun and Liang, 2020; Rastogi et al., 2022). Compared to its predecessor, the GCMs of CMIP6 have not exhibited substantial improvement in the simulation of rainfall at the regional scale. For example, the uncertainty range is not narrowed, and systematic biases like the overestimated rainfall frequency are not reduced (Jiang et al., 2020; Xin et al., 2020). As a result, utilizing CMIP6 to project CREs may not provide a meaningful improvement, especially in topographically complex locations like the TGR. Adopting a dynamical downscaling approach based on regional climate models (RCMs) may overcome some of these limitations to improve the projection of CREs.

The dynamical downscaling by RCMs can show more realistic climatological rainfall in specified regions compared with GCMs, because it can compensate for the errors related to the coarse grids in GCMs. RCM downscaling has been widely applied to simulate past and future regional climate. The higher resolution leads to a better representation of finer-scale microphysical and resolved convective processes, as well as the effects of topography and land surface. For example, Yu et al. (2015) and Bao et al. (2015) found that RCM downscaling can eliminate artificial rainfall maxima and correct overestimations of rainfall amount simulated by GCMs. The deficiencies in the annual cycle and frequency of rainfall in GCMs are also reduced well by RCMs (Huang and Gao, 2018; Liang et al., 2019; Gao, 2020; Wang et al., 2021b). Additionally, the multi-model ensemble mean based on the Coordinated Regional Climate Downscaling Experiment has been shown to provide more specific features and less uncertainties in future rainfall in China (Gu et al., 2018; Li et al., 2018).

RCM downscaling can correct for many errors in GCMs, but it is not a panacea. The downscaling of RCMs does not imply a more realistic simulation if there are large uncertainties in the boundary conditions of the driving GCMs (Sato and Xue, 2013; Wang and Kotamarthi, 2015; Rastogi et al., 2022). On the other hand, a realistic simulation of the present climate is a necessary but insufficient test of the accuracy of future climate changes shown by RCM downscaling, since the relation between the accuracy of a present-day simulation and future projection is not guaranteed (Racherla et al., 2012; Singh et al., 2017; Potter et al., 2020). To address these two issues and to provide credibility in downscaling results requires careful selection of realistic driving models and identification of robust physical mechanisms underlying the projected changes (Hall, 2014; Rastogi et al., 2022). Additionally, simulation discrepancies between the driving model and the downscaling model are prevalent in the projected climate change. For example, Bao et al. (2015) found that WRF downscaling exhibited an opposite change in the heavy rainfall amount and intensity in the Yangtze River basin to that of the model used to drive it (GFDL-ESM2G); Ji and Kang (2013) showed that the rainfall rates were different between BCC_CSM1.1 and RegCM4 downscaling over the Tibetan Plateau; Bartók et al. (2017) showed a significant decline in cloudiness for Europe in future based on GCMs, while corresponding regional downscaling experiments showed no significant changes. These discrepancies are issues that need to be focused on in downscaling simulation.

Generally, the horizontal grid spacing of single-nest RCMs is 30–60 km, which is largely determined by the GCM’s resolution (Gao et al., 2011; Im et al., 2015). However, such a resolution does not perform well in the TGR where there is complex terrain and a prominent drop in elevation (Strong et al., 2017; Wang et al., 2017; Huang et al., 2020a; Fu et al., 2023). Downscaling simulation based on a single nest exhibits visible biases in the probability distribution of daily rainfall and fails to reproduce the spatial pattern of extreme rainfall in mountainous areas of central China (Bao et al., 2015; Yu et al., 2015). Thus, a much finer resolution is required, which can be obtained by the multi-nested technique. A double-nested RCM downscaling in the TGR revealed that the inner domain produced rainfall patterns that the parent domain failed to capture (Wu et al., 2012). Using WRF (Weather Research and Forecasting) model triple-nested downscaling, Li et al. (2019b) found that the inner domain exhibited a visible improvement in rainfall distribution and related atmospheric circulation compared to the parent domain.

In the present study, we utilized a WRF triple-nested dynamical downscaling strategy to project future CRE changes over the TGR. To cope with the subsequent flooding of the monsoon season, reservoir discharge in spring is necessary. As the water level drops, the force exerted by the water on the mountain slope changes and seepage at the slope’s base increases. If a prolonged rainfall event occurs at this time, the slope is prone to deformation and can develop into a landslide. Thus, spring is the peak season for CRE-induced geological hazards in the TGR, and the focus of this study.

The remainder of this paper is organized as follows: section 2 describes the data, models, methods and the definition of CREs; section 3 comprehensively validates WRF’s performance in reproducing rainfall and CREs in historical experiments, gives future projections under Shared Socioeconomic Pathway (SSP) scenarios, and assesses the reasonability of the projections; and finally, a summary and discussion are given in section 4.

To assess the model performance in simulating the observed rainfall, a daily gridded gauge precipitation dataset, CN05.1, with a resolution of 0.25° × 0.25°, is employed (Wu and Gao, 2013; Wu et al., 2017). CN05.1 was produced based on 2400 observation stations across China and has been widely used (Sui et al., 2015). Considering the lack of hourly rainfall measurements in CN05.1, we use the China Meteorological Forcing Dataset (CMFD) to evaluate the diurnal cycle. CMFD was derived from a combination of station data, reanalysis datasets and remote sensing products, and has a 3-h temporal resolution and a 0.1° × 0.1° spatial resolution (He et al., 2020). Due to its continuous temporal coverage and consistent quality, CMFD is one of the most widely used climate datasets for China (Peng et al., 2022; Wang et al., 2022). In addition, the fifth major global reanalysis produced by ECMWF (ERA5) is utilized to assess the simulation of atmospheric circulation (Hersbach et al., 2020).

In this study, the RCM used is WRF version 4.0, which was developed at the National Center for Atmospheric Research. It is a mesoscale numerical model consisting of a numerical solver of the fully compressible Eulerian and non-hydrostatic equations (Skamarock et al., 2019). WRF has been widely used in regional climate studies over China, and the results prove it can reproduce the main features of climate and extremes well. Moreover, the performance of the interannual variation provided by WRF is obviously better compared to its driving data (Sato and Xue, 2013; Yu et al., 2015).

Although the biases in GCMs can be mitigated to some extent through the downscaling process, they still potentially impact the downscaling results. Thus, regional downscaling is best driven by GCMs that reasonably simulate the region’s climate. Considering WRF is driven by the GCM’s prognostic variables that are related to rainfall but not by rainfall itself, a realistic simulation of rainfall by WRF is dependent on the accuracy of these prognostic boundary conditions (Sato and Xue, 2013; Rastogi et al., 2022). Hence, we first employ 10 GCMs (at the time, the only 10 models providing 6-h temporal-resolution output) and evaluate the simulated spring low-level winds (850 and 925 hPa), which are highly correlated with spring rainfall in China (Wu et al., 2012; Li et al., 2019b). The overall performance of each model can be judged using a comprehensive rating index (MR), which provides a statistical summary between simulations and ERA5, in terms of their spatial correlation coefficient, root-mean-square difference, and standard deviation difference. The details of MR can be found at Li et al. (2016a). From Table 1, MIROC6 (0.59) exhibits the best MR performance among all 10 models, followed by NorESM2-MM (0.57) and then GFDL-ESM4 (0.56). In pervious works (Kataoka et al., 2020; Tian et al., 2021), MIROC6 performed better in reproducing the seasonal progression of the East Asian monsoon, which is consistent with our results. Therefore, MIROC6 is used for the boundary and initial conditions of WRF. MIROC6 is composed of three sub-models: atmosphere, sea ice–ocean and land. The atmosphere model is based on the Center for Climate System Research, University of Tokyo/National Institute for Environmental Studies (Numaguti et al., 1997) atmospheric general circulation model; the sea ice–ocean model is based on the Center for Climate System Research Ocean Component model (Hasumi, 2006); and the land surface model is based on the Minimal Advanced Treatments of Surface Interaction and Runoff model (Takata et al., 2003). The atmospheric component of MIROC6 has a horizontal resolution of 1.4° × 1.4° (Tatebe et al., 2019).

To further illustrate the accuracy of the driving MIROC6, we provide a basic comparison between the observations and MIROC6 (Fig. 2). In general, the spring rainfall climatology of MIROC6 matches that of CN05.1 with a correlation of 0.64. Additionally, MIROC6 can capture the strength and pattern of the winds well, with a correlation coefficient of 0.74 at 925 hPa and 0.84 at 850 hPa between ERA5 and MIROC6. The strength and pattern of water vapor flux in MIROC6 are consistent with those in ERA5, with a pattern correlation coefficient of 0.88. Generally, the spring rainfall and associated atmospheric conditions in MIROC6 are consistent with those in observations, which lay the basis for achieving a reliable rainfall simulation by using MIROC6 to drive WRF.

Figure 2.  The spatial distribution of spring (a–c) total rainfall amount (units: mm), (d–f) 925-hPa wind speed (units: m s⁻¹), (g–i) 850-hPa wind speed (units: m s⁻¹), and (j–l) water vapor flux in the whole layer below 200 hPa [arrows; units: kg (m⁻¹ s⁻¹)] and its divergence [shaded; units: 10⁻⁵ kg (m⁻² s⁻¹)] in 1995–2014 from observations, MIROC6 and WRF.

A triple-nested downscaling simulation based on WRF is applied in the present study. Figure 1b shows the boundary and model topography of each domain. The outermost-domain (Domain1) grid is centered at (29°N, 110°E) and has 36 km horizontal grid spacing (with 91 × 91 horizontal grid points). To ensure smooth solutions, Domain1 grid cells of WRF closer than five cells from an outer boundary are relaxed towards MIROC6. The innermost-domain (Domain3) grid has 4 km grid spacing (with 127 × 166 horizontal grid points). Regarding the vertical coordinates, we configured 33 terrain-following eta levels from the surface to 50 hPa in each domain. The initial and lateral boundary conditions for Domain1 are derived from the historical simulation and future projection of MIROC6 at 6-h intervals for meridional–zonal wind, specific humidity, air temperature, calculated geopotential height, 2-m air temperature, skin temperature, surface pressure and sea level pressure. Besides, daily sea surface temperature, monthly soil moisture and temperature are also provided by MIROC6.

For WRF downscaling, we use the following configuration: the New Thompson microphysics scheme (Thompson et al., 2008); the Kain–Fritsch convection scheme (Kain, 2004); the Noah Land Surface Model (Niu et al., 2011); the Rapid Radiative Transfer Model longwave radiation scheme; the Dudhia shortwave radiation scheme (Dudhia, 1989); and the Yonsei University boundary layer scheme (Hong et al., 2006). Among them, the microphysics scheme constitutes a key configuration in rainfall modelling. The New Thompson microphysics scheme includes ice, snow, graupel, and their associated processes, which performs well in high-resolution rainfall simulations over the complex terrain of East Asia (Strong et al., 2017; Tiwari et al., 2018; Huang et al., 2020b). Considering the spatial resolutions are sufficiently high for the non-hydrostatic dynamical core of the model to partially resolve sub-mesoscale convective motions (Li et al., 2019a), we therefore negate the convection parameterization of Domain3.

To save computational resource, two time-slice integrations are used to represent the typical present and future climate. Here, we use the end of the 21st century as a typical future climate to emphasize the signals of change in precipitation due to warming (Feng et al., 2011). The two integrations cover the period from 24/25 February to 1 June during 1995–2014 and 2075–2094 (last 20 years of available SSP output of MIROC6), respectively. The first five days (24/25 February to 28/29 February) is considered as a spin-up and is therefore not included in the subsequent analysis. The projection outputs of MIROC6 forced with the SSP2-4.5 and SSP5-8.5 scenarios are employed. SSP2-4.5 is a combined scenario of a medium energy-intensive, socioeconomic developmental path with a rising radiative forcing peaking at 4.5 W m⁻² in 2100. It is comparable to several planned emission pathways required by the mitigation policies of the Paris climate agreement, and is therefore considered a more preferable scenario compared with other SSP scenarios (O'Neill et al., 2016). SSP5-8.5 represents a combined scenario of a high energy-intensive, socioeconomic developmental pathway with strong radiative forcing, which peaks at 8.5 W m⁻² by 2100. Although SSP5-8.5 is the highest emission scenario available in CMIP6, it still assumes emissions well below what the current energy mix would produce in the future (Birkinshaw et al., 2017). Based on the SSP5-8.5 scenario, we can discuss the rainfall changes under extreme warming conditions (Wang and Kotamarthi, 2015).

A non-parametric bootstrapping approach is used to estimate confidence intervals based on the following procedure. First, 500 non-parametric bootstrap samples are generated based on time series obtained by resampling of time steps from the original 20-year dataset of regional averages. Then, quantiles are applied to the 500 estimates of change to obtain 90% confidence intervals (Shao et al., 2012; Hailegeorgis et al., 2013).

To understand the differences in changes in CREs between WRF and MIROC6, we employ a moisture budget analysis (Chou et al., 2013; Ma et al., 2015b; Li et al., 2020). The moisture budget equation is

where P is precipitation; $ \mathbf{V} $ is wind vector; q is specific humidity; E is evaporation, which is calculated from the latent heat fluxes; δ is a residual term including transient eddies, contributions from model interpolation bias, and surface processes; $\left\langle{\;}\right\rangle $ represents a vertical mass integration through the troposphere (below 200 hPa); $ {-}{\partial}_{\text{t}}\left\langle{{q}}\right\rangle $ is the time derivative of q, which implies the change in local water vapor storage; and $ {-}\nabla\cdot \left\langle{{q}\mathbf{V}}\right\rangle $ is the convergence of integrated moisture flux. This can further be divided into two terms according to the mass conservation equation: a horizontal moisture advection term and a vertical moisture advection term. Therefore, Eq. (1) can be reformulated as

To understand the mechanisms responsible for the trend of change, according to Eq. (2), the precipitation changes can be decomposed into

Here, the prime symbol $ {}^{'} $ represents the departure from the historical simulation. In Eq. (3), the changes of vertical moisture advection $ -{\left\langle{\mathrm{\omega }{\partial }_{\mathrm{p}}\mathrm{q}}\right\rangle}^{{{'}}} $ can be further divided as follows:

where $ \stackrel{-}{\omega } $ denotes the climatology in the historical simulation. The first term on the right-hand side of Eq. (4) denotes the changes in q with $ \omega $ unchanged, commonly called the thermodynamic component of the vertical moisture advection term, contributed by the changes in water vapor; the second term denotes the changes in $ \omega $ with q remaining constant, associated with changes in vertical velocity, which is induced by atmospheric circulation changes, and often called the dynamic component; the third term is a nonlinear term that is induced by changes in q and $ \omega $, and is much smaller than the other terms and can be ignored (Seager et al., 2010; Li et al., 2020). Therefore, Eq. (3) can be written as

To analyse the effect of synoptic processes on precipitation, each term of the moisture budget is calculated using the daily mean data.

When the accumulated rainfall amount of one day is equal to or more than 1.0 mm within 24 h, the day is defined as one rainy day. A CRE is set to begin if any of the consecutive conditions occur, as described in Zheng et al. (2020). To quantify the CREs, we use four indices related to geological disasters: the occurrence frequency (OCF), rainy days (TRD), rainfall amount (ACR), and rainfall intensity (INT) (Corominas and Moya, 1999). The OCF describes the frequency of CREs, while the TRD and ACR are the sum of CRE rainy days and accumulated rainfall amounts, respectively. INT describes the mean daily rainfall intensity during CREs, which is not independent of ACR and TRD but equal to ACR divided by TRD (Zheng et al., 2020).

To illustrate the skill of WRF in downscaling rainfall, we first compare WRF’s results for Domain1 with observations during 1995–2014. CN05.1 shows a spring rainfall climatology that decreases from southeast coastal areas to the northwest (Fig. 2a). WRF reproduces the spatial pattern of rainfall well (Fig. 2c). The location of the rainfall maximum in WRF is consistent with observation, indicating that WRF tends to correct the northward bias of the rainfall maximum in MIROC6. For example, the spatial correlation coefficient of rainfall in Domain1 between observation and WRF downscaling (MIROC6) is 0.81 (0.64). Additionally, compared to MIROC6, the root-mean-square error in WRF downscaling is reduced by 19%.

To understand the better performance of WRF, we further evaluate the simulated spring low-level wind and water vapor (Wu et al., 2012; Li et al., 2019b). In general, WRF reproduces the spatial distribution of level-wind throughout East China, especially the wind vectors (Figs. 2f and i). This can be attributed to the realistic boundary conditions provided by MIROC6. Besides, the regional wind features in WRF around the TGR are improved substantially. For example, the enhancement and the deflection of winds on the northwest side of the TGR are reproduced well by WRF, while MIROC6 cannot successfully simulate the regional low-level winds surrounding the TGR (Figs. 2e vs 2d; 2h vs 2g). Additionally, the strength and pattern of water vapor flux in WRF downscaling are consistent with those in ERA5, with a pattern correlation coefficient of 0.90. WRF exhibits a more realistic distribution of divergence than MIROC6 (Fig. 2l). Specifically, in ERA5, the northern boundary of the water vapor convergence zone in Southeast China is located at 32°N. In MIROC6, it is located at 36°N, consistent with the model’s northward bias of the rainfall centre. WRF corrects the northern boundary of water vapor convergence to about 32°N.

Figure 3a shows the monthly evolution of rainfall from observations and simulations in the TGR (Domain3). The two observations indicate that the monthly rainfall amount increases during spring. WRF captures the monthly increases of rainfall, but has a similar positive bias in all three months. However, WRF clearly reduces the overestimation of rainfall in MIROC6 in March and April. The bias of the spring seasonal mean rainfall in WRF is 20.9% relative to CN05.1, while it is 39.8% in MIROC6.

Figure 3.  (a) The monthly evolution (units: mm), (b) the daily probability distribution (%), (c) the diurnal cycle (units: mm h⁻¹), and (d) the frequency of different durations of spring rainfall. The results are for the period 1995–2014 from observations and model simulations over the TGR.

Figure 3b compares the probability distribution of the daily rainfall amount from observations and simulations. For this analysis, data from all grid cells are considered. The two observations show that the probability of the daily rainfall amount decreases as the daily rainfall amount increases. WRF exhibits a better agreement with the observations, generally with less difference than MIROC6 from observations. MIROC6 underestimates the frequency of very light rainfall but overestimates the frequency of the other rainfall bins, which is a common bias in GCMs (Watanabe et al., 2010; Tian et al., 2021; Zheng et al., 2022). For rainfall less than 1.0 mm d⁻¹, the probability is 59.6% and 57.6% in CN05.1 and CMFD respectively, but only 45.8% in MIROC6. WRF corrects the MIROC6 bias and simulates a probability of 61.8%, which is much more consistent with observations. For rainfall between 1.0 and 30.0 mm d⁻¹, the bias ratio of WRF is about −9.6% relative to CN05.1, whereas for MIROC6 it is about 33.3%. MIROC6 and WRF both overestimate the probability of daily rainfall being greater than 30.0 mm d⁻¹. MIROC6 overestimates the probability by 59.4% relative to CN05.1, while WRF overestimates it by 92.6%. There are several reasons to believe that the higher probability of daily rainfall exceeding 30.0 mm d⁻¹ simulated by WRF is realistic. First, being a high-resolution model, WRF can better resolve the mesoscale topography and convective rainfall process, which coarse-resolution models tend to smooth out (Heikkilä et al., 2011; Yu et al., 2015). Second, the interpolation of station data to form the observational gridded dataset causes underestimation in the frequency of extremes compared to the original station gauges. And third, the locations of gauges in the TGR are mostly in valleys and lowland areas, implying a lack rainfall measurements on hilltops or mountain slopes (Argüeso et al., 2013). These reasons can explain why WRF has a higher frequency of heavy rainfall than observations or MIROC6.

Figure 3c gives the diurnal cycle of rainfall in the TGR. CMFD exhibits a peak during overnight hours (0500 LST; LST=UTC+8 hours) and a trough during afternoon/early evening hours (1700–2000 LST). WRF is in good agreement with observation in terms of the peak and trough hours, as well as giving similar hour-by-hour rainfall amounts. In MIROC6, meanwhile, the simulated trough has a 1–2-hour lag compared to CMFD, and the hourly rainfall is clearly overestimated. The bias of the rainfall magnitude in MIROC6 is about two times that in WRF.

Reproducing CREs requires models that can capture the duration of continuous rainfall spells. Figure 3d gives the frequency of rainfall events of different duration. For rainfall events lasting 2–4 days, the frequencies in WRF and MIROC6 are similar to that observed, occurring about eight times per spring. For rainfall events lasting 5–7 days, the frequency in CN05.1 and CMFD is 1.06 and 1.34 per spring, respectively. In WRF, it is 1.30 per spring, which is consistent with CMFD. In contrast, MIROC6 overestimates the frequency as 2.01 per spring. An observed occurrence of a duration more than 8 days is rare, with 0.25 and 0.26 per spring in CN05.1 and CMFD, respectively, and 0.31 per spring in WRF. Meanwhile, MIROC6 shows a visible overestimation of 0.85 per spring.

Figure 4 shows the spatial distribution of rainfall features averaged in 1995–2014 from observations, WRF and MIROC6. For the total rainfall amount (Figs. 4a–d), observations show a peak over the southeast TGR, decreasing northwestward. WRF downscaling successfully reproduces the rainfall pattern, decreasing from southeast to northwest, but it overestimates rainfall over the Wuling Mountains. Compared with WRF, MIROC6 simulates an artificial high rainfall area over the northeastern edge. The biases of rainfall here may arise from the coarse resolution and the systematic precipitation biases of MIROC6 (Zheng et al., 2020). For total rainy days (Figs. 4e–h), observations exhibit a maximum in the south, decreasing northward. It is worth noting that the high-value centres between the two sets of observations show differences, which may lead to uncertainties in the observed spatial pattern of CREs. WRF downscaling captures the observed spatial features, but it underestimates the frequency in the southeast. MIROC6 simulates too many rainy days, especially in the southeast. In observations, the rainfall intensity has a similar spatial pattern as the total rainfall amount. The intensity distribution simulated by WRF is similar to observations, with a spatial correlation coefficient between WRF and CN05.1 (CMFD) of 0.80 (0.87). Meanwhile, MIROC6 simulates an artificial intensity maximum in the northeast and has a low correlation coefficient of 0.68 (0.65) with CN05.1 (CMFD).

Figure 4.  Spatial distributions of spring (a–d) total rainfall amount (units: mm), (e–h) total rainy days (units: d), and (i–l) rainfall intensity (units: mm d⁻¹) over the TGR from observations, WRF, and MIROC6 averaged in 1995–2014. The value in the upper-right gives the regional average.

The spatial distribution of rainfall events lasting for different durations is shown in Fig. 5. There is a high frequency of rainfall events with a duration of 2–4 days throughout the region in WRF, consistent with CN05.1, with a correlation coefficient of 0.53. Meanwhile, MIROC6 fails to reproduce the pattern and shows no correlation with the observations. For rainfall events with a duration of 5–7 days, observations show a pattern of more in the south and less in the north. WRF underestimates the frequency of rainfall in this category in the southeast because of the low number of rainy days (Fig. 4g). The spatial correlation coefficient between CN05.1 and WRF downscaling (MIROC6) is 0.42 (0.31). The two sets of observations show obvious discrepancies when it comes to rainfall lasting longer than 8 days. CN05.1 has two high-value centres along the Wuling Mountains and one high-value centre over the southeast plain, whereas CMFD only has large values in the southwest. In WRF, the maximum centre appears along the Wuling Mountains, while in MIROC6 the frequency is overestimated in the whole region.

Figure 5.  Frequency distribution of rainfall events with a duration of (a–d) 2–4 days, (e–h) 5–7 days, and (i–l) ≥ 8 days over the TGR from observations, WRF and MIROC6 in 1995–2014.

Figure 6 displays the spatial distribution of CREs in observations, WRF, and MIROC6. In CN05.1, the OCF of CREs is 1.67 per spring averaged over the TGR, with a higher value in the south than in the north. The TRD and ACR in CN05.1 is 17.13 days and 136.15 mm, respectively, and the spatial distribution is similar to the OCF. The observed INT exhibits a high value in the southeast and a low value in the northwest. Generally, the distributions of CREs in CN05.1 and CMFD are similar, although there are also some magnitude differences, which are due to the different features in rainfall duration (Fig. 5).

Figure 6.  Spatial distributions of (a, e, i, m) occurrence frequency (OCF), (b, f, j, n) total rainy days (TRD) (units: d), (c, g, k, o) accumulated rainfall amount (ACR) (units: mm), and (d, h, l, p) mean daily rainfall intensity (INT) (units: mm d⁻¹) of spring CREs over the TGR from observations, WRF, and MIROC6 in 1995–2014. The value in the upper-right gives the regional average.

In WRF, the spatial pattern of OCF is generally consistent with the observations, despite an underestimation in the southeast corner. The spatial correlation coefficient of OCF between WRF and CN05.1 is 0.54, which is higher than that between MIROC6 and CN05.1 (0.45). Besides, the average OCF in WRF is 1.90 per spring, which is closer to the observations than MIROC6. MIROC6 obviously overestimates the magnitude with a value of 3.33 per spring. The spatial correlation coefficient between WRF and CN05.1 for TRD/ACR is 0.58/0.62, which is also better than that between MIROC6 and CN05.1 (0.51/0.49). For INT, WRF exhibits a southeast-to-northwest decreasing pattern, which is similar to CN05.1 and CMFD, but the magnitude (9.31 mm d⁻¹) is higher than observed (6.46 mm d⁻¹ and 7.61 mm d⁻¹, respectively).

The above analyses demonstrate that WRF downscaling exhibits a significant improvement in terms of reproducing spring rainfall, as well as CREs, in the TGR. The more accurate simulation of present climate provides confidence in the credibility of future projections.

Figure 7a gives the regionally averaged percentage changes of future monthly rainfall derived from WRF and MIROC6 in 2075–2094 relative to 1995–2014. The spring rainfall in WRF is projected to decrease by −2.6% and −3.9% under SSP2-4.5 and SSP5-8.5, respectively, and the decrease is mainly contributed by the decline in May. In MIROC6, meanwhile, the rainfall change is opposite to that in WRF, increasing by 11.4% and 16.2% under SSP2-4.5 and SSP5-8.5, respectively.

Figure 7.  The percentage changes of (a) monthly rainfall, (b) the daily rainfall probability distribution, (c) the frequency of continuous rainfall, and (d) CRE indices over the TGR in SSP projected runs (2075–2094 averaged) relative to historical runs (1995–2014 averaged) from WRF and MIROC6. Vertical lines and ranges depict the 90% bootstrap confidence intervals.

Figure 7b displays the change in the probability distribution of the daily rainfall amount. In WRF, the probability of rainfall less than 1 mm d⁻¹ is expected to increase by 7.5% and 8.6% under SSP2-4.5 and SSP5-8.5, respectively. However, rainfall between 1 and 30 mm d⁻¹ exhibits an obvious decline (−16.4%/−13.0% under SSP2-4.5/SSP5-8.5). As for rainfall above 30 mm d⁻¹, the probability will increase, according to WRF. The ratio increase of the 30–50 mm bin and >50 mm bin is 3.7%/6.9% and 7.5%/27.4% under SSP2-4.5/SSP5-8.5, respectively, while the increase under SSP2-4.5 is insignificant. In contrast, the change in MIROC6 is different from that in WRF. For rainfall less than 10 mm d⁻¹, MIROC6 exhibits no clear signal, while for rainfall greater than 20 mm d⁻¹ the rainfall frequency is projected to increase profoundly. For example, the ratio increases of the 20–30 mm bin, 30–50 mm bin, and >50 mm bin are 8.1%/10.8%, 48.8%/28.5% and 70.1%/80.3% under SSP2-4.5/SSP5-8.5, respectively. This is consistent with previous studies that found GCMs simulate intensified daily rainfall in a warmer climate (Feng et al., 2011; Gu et al., 2018). It is worth noting that the growth of heavy rainfall in WRF is lower than that in MIROC6. Previous studies have also found such a reduced rainfall change ratio when using WRF downscaling, and they suggest that the change in WRF is more reasonable because of its more proper representation of convective rainfall (Bao et al., 2015; Rastogi et al., 2022).

Figure 8a–d shows the spatial changes of projected rainfall amount over the TGR under two emission scenarios. WRF simulates clear regional features with an increase in the mid-west and a decrease in the east. MIROC6 projects an increase in the entire region with a maximum over the southeast. Next, we check the change in rainfall frequency and intensity (Figs. 8e–l), which is associated with the change in rainfall amount. The rainfall intensity in WRF shows a pronounced increase in most parts of the TGR. In MIROC6, the intensity is projected to increase, especially in the east. As for the changes in rainfall frequency, obvious differences are found between WRF and MIROC6. WRF shows a decreased frequency in the whole region (−4.6/−5.2 days averaged in the TGR under SSP2-4.5/SSP5-8.5), and the decrease is more prominent in the southeast, while the frequency in MIROC6 exhibits slight changes by 0.01/0.84 days under SSP2-4.5/SSP5-8.5 averaged over the TGR. It is worth noting that the change in frequency in MIROC6 is inconsistent between the SSPs. For example, the frequency in MIROC6 decreases in the northeast and increases in the southwest under SSP2-4.5, while it increases over almost the entire region under SSP5-8.5 (Figs. 8g and h); and for the change in the probability distribution of rainfall less than 20 mm, the two SSPs show opposite changes in some intervals (Fig. 7b). This indicates uncertainties in the MIROC6 projections.

Figure 8.  Spatial distributions of future changes (2075–2094 relative to 1995–2014) of (a–d) total rainfall amount (units: mm), (e–h) rainy days (units: d), and (i–l) rainfall intensity (units: mm d⁻¹) over the TGR under the SSP emissions scenarios as projected by WRF and MIROC6. Shading denotes the ratio 0 (no change) does not fall into the 90% bootstrap confidence intervals range.

As WRF’s rainfall frequency decreases, the duration of continuous rainfall also decreases. Figure 7c gives the change in frequency of rainfall events with different durations averaged over the TGR. For rainfall events lasting 2–4 days, WRF suggests that the frequency will decrease significantly, while MIROC6 shows no obvious change. WRF and MIROC6 disagree on the changes in rainfall durations of more than 5 days: a significant decrease is projected by WRF, and the longer the duration, the more pronounced the decrease; while MIROC6 generally exhibits an increase in the frequency of rainfall events lasting longer than 5 days. As for the spatial distribution of change by WRF (not shown), rainfall events lasting 2–4 days show an increase in the Wuling Mountains and a decrease in other parts of the TGR under SSP2-4.5, while a decrease is projected in the entire region under SSP5-8.5. For rainfall events lasting longer than 5 days, the frequency is projected to decrease in most of the TGR under the two SSPs, especially in the Wuling Mountains.

Figure 7d gives the projected changes of CRE indices averaged in the TGR. Overall, the CREs of WRF show a reduction in OCF, TRD and ACR, and a weakening in INT, under SSP2-4.5. Under SSP5-8.5, CREs also decrease, but the magnitude is weaker than under SSP2-4.5. In contrast, MIROC6 exhibits no change in OCF, and an opposite signal in TRD under the two SSPs. However, ACR and INT are projected by MIROC6 to increase, despite the change in INT being insignificant under SSP5-8.5.

Figure 9 gives the spatial distribution of projected CREs. WRF simulates a decrease in OCF over the whole region under the two SSPs (Figs. 9a and e), due to the decreased rainfall frequency. For TRD, WRF projects a prominent decrease in the east and a slight increase in the west under SSP2-4.5. Under SSP5-8.5, WRF exhibits a decrease in the northwest and southeast, and a regional increase in the central TGR. The spatial distributions of ACR and INT change according to both SSPs, generally following those of TRD, which are projected to decrease in most parts of the TGR, with regional increases seen in the central-west area.

Figure 9.  Spatial change distributions (2075–2094 relative to 1995–2014) of CREs’ (a, e, i, m) occurrence frequency (OCF), (b, f, j, n) total rainy days (TRD; units: d), (c, g, k, o) accumulated rainfall amount (ACR; units: mm), and (d, h, l, p) mean daily rainfall intensity (INT; units: mm d⁻¹) over the TGR under the SSP scenarios as projected by WRF and MIROC6. Shading denotes the ratio 0 (no change) does not fall into the 90% bootstrap confidence intervals range.

Compared to WRF, the spatial changes of CREs projected by MIROC6 are quite different. Under SSP2-4.5, OCF and TRD in MIROC6 increase in the southwest and decrease in the east of the TGR, and ACR and INT increase in the whole region. Meanwhile, under SSP5-8.5, MIROC6 projects an increase of OCF in the east and a decrease in west. The spatial distribution of TRD is opposite to that of OCF, which indicates the duration of a single CRE in the west will prolong despite the decreased occurrence. Besides, an increase in ACR is projected by MIROC6. For INT, there is a decrease in the central area and an increase in the southeast.

Previous studies have presented many discrepancies between the driving model and the downscaling model in the projected climate change (Ji and Kang, 2013; Bao et al., 2015; Bartók et al., 2017). In the present study, the future changes of CREs projected by WRF and MIROC6 also show differences. These differences are attributable to the differences in the simulated rainfall frequency. In WRF, the rainfall frequency shows an overall decline in the entire TGR, while MIROC6 shows no obvious change. Thus, the question that remains to be addressed is why the projected changes in rainfall frequency are different between MIROC6 and WRF. Since the changes in atmospheric variables based on SSP2-4.5 and SSP5-8.5 are similar, but the signal is stronger under SSP5-8.5, unless otherwise stated, the results presented hereafter are for SSP5-8.5.

To understand the different changes in the rainfall frequency between WRF and MIROC6, we diagnose the moisture budget over the TGR. Figure 10a shows the climatological rainfall and other moisture terms in ERA5 and simulations in 1995–2014 [Eq. (2)]. ERA5 shows that local rainfall is generally contributed by the vertical moisture advection and evaporation term. The other terms contribute little to the spring rainfall. WRF and MIROC6 show similar contributions in reproduced evaporation. However, they exhibit differences in the vertical moisture advection term. In WRF, the contribution of vertical moisture advection is 3.20 mm d⁻¹, similar to that in observation (2.96 mm d⁻¹), while MIROC6 (5.75 mm d⁻¹) shows obvious overestimation that is about twice as much as observed. This bias indicates that the rainfall simulation and projection in MIROC6 could bear large uncertainties.

Figure 10.  (a) Moisture budget for the rainfall climatology in 1995–2014 and (b) moisture budget changes (2075–2094 relative to 1995–2014) over the TGR. The bars indicate precipitation (P), evaporation (E), local moisture storage (Dq), horizontal moisture advection (Vq), vertical moisture advection (wq), the residual term (R), thermodynamic component (TH), and dynamic component (CD) of the vertical moisture advection term. Different colours represent the results of observations and simulations.

In ERA5, the spatial pattern of the vertical moisture advection term is inhomogeneous (not shown). The high-value centres are located in the Wuling Mountains as well as the southern area, which is similar to the rainfall distribution and implies that the mountainous terrain in the south TGR can enhance the convergence of water vapor flux and then favour rainfall formation. The spatial pattern of the vertical moisture advection term in WRF is similar to the observation, while there is an artificial maximum located in the central TGR in MIROC6.

To further illustrate the relationship between moisture budget terms and synoptic rainfall in the TGR, we calculate the regression fitting of the individual terms to daily rainfall. Specifically, in ERA5, the positive linear correlation between daily rainfall and the vertical moisture advection term is significant (R² = 0.61, p < 0.05), while the relationship in WRF is weaker than observed but still significant (R² = 0.38, p <0.05). This indicates that the daily rainfall frequency is closely tied to the upward motion. The contribution of the other terms is not obvious. The correlation is relatively low (R² = 0.10, p < 0.05) in MIROC6.

To understand the different trends in the frequency between WRF and MIROC6, we analyse the change of each moisture budget term in 2075–2094 relative to 1995–2014, as shown in Fig. 10b [Eq. (5)]. In WRF, the vertical moisture advection term shows a decrease (−0.73 mm d⁻¹), while evaporation and horizontal moisture advection show an increase. Clearly, the decrease in rainfall can be mainly attributed to the decline in the vertical moisture advection term. This indicates that fewer rainfall events can be triggered as upward motion decreases, which leads to a decrease in the rainfall frequency in WRF. We further separate the change in the vertical moisture advection term into two components, and find that the decrease in the dynamic component (−0.98 mm d⁻¹) overwhelms the increase in the thermodynamic component (0.24 mm d⁻¹), resulting in a decline in the vertical moisture advection term. Compared to WRF, the changes in the moisture budget terms in MIROC6 are somewhat different. The vertical moisture advection term shows an obvious increase (1.25 mm d⁻¹), which is responsible for the increased rainfall in MIROC6. To balance the increase in the vertical moisture advection, the horizontal moisture advection exhibits a decrease (−0.73 mm d⁻¹). When we decompose the vertical moisture advection term into dynamic and thermodynamic components, we find that both exhibit an increase (0.84 mm d⁻¹ and 0.36 mm d⁻¹, respectively).

But how can we understand the change in the vertical moisture advection term? It may be related to the difference in the variation of the monsoon. Studies have suggested that the intensified latent heating released by the warming sea surface temperature over the South China Sea could build and reinforce the heat-induced atmospheric circulation, suppress the ascending motion over south-central China, and thus lead to a decreasing trend of rainfall frequency (Li et al., 2016b). This could be responsible for the decreased upward motion in WRF over the TGR. However, it cannot explain the increased vertical moisture advection in MIROC6.

Denson et al. (2021) argued that relative humidity will play an important role in changes to the climatological rainfall frequency at the local scale in a warmer climate. Because atmospheric humidity does not increase at the same rate as temperature, there must be a decrease in light-to-moderate rainfall and the frequency of rainfall events (Mao et al., 2015; Mishra, 2019). Thus, we further check the possible relationship between projected relative humidity and rainfall frequency over the TGR. Generally, WRF and MIROC6 both show a consistent increase in specific humidity (by 3.5 g kg⁻¹ and 3.7 g kg⁻¹, respectively) in the TGR; and as the increase in specific humidity with temperature is slower than the increase in saturation vapour pressure (Wang et al., 2020), a decrease in relative humidity is caused (Fig. 11a and b; decreased by −7.4 % and −5.8 % in WRF and MIROC6, respectively, in the TGR). Cloudiness and associated rainfall are only possible when relative humidity reaches 100% and condenses, which is often induced by air rising and cooling to the dewpoint temperature adiabatically. When the average relative humidity is decreased, the air will be less likely to reach saturation, despite the increased water content, resulting in less cloud formation for an equivalent uplift in energy (Mao et al., 2015; Wang et al., 2020; Denson et al., 2021).

Figure 11.  Spatial change distributions (2075–2094 relative to 1995–2014) of (a, b) relative humidity (%) in the whole layer below 200 hPa, and (c, d, e, f, g, h) cloud cover fraction (%) in different levels under the SSP5-8.5 scenario as projected by MIROC6 and WRF.

To test the above inference, we check the simulated changes in cloud cover fraction (Figs. 11c–h). In general, good agreements are found between the two models: warming leads to a reduction in total cloud cover across eastern China in the future (decrease by −4.8% and −10.1% in MIROC6 and WRF, respectively, over the TGR). When it comes to the clouds at different levels, the decrease in low cloud (below 2 km) is the most obvious, with an average of −3.5% (MIROC6) and −6.8% (WRF) within the TGR. The mid-level cloud (2–6 km) is projected to decrease by −1.6% and −5.5% in MIROC6 and WRF, respectively. In the TGR, the spatial correlation between relative humidity and total cloud cover is 0.87 (0.93), and the temporal correlation is 0.94 (0.96), under SSP5-8.5 (SSP2-4.5). This implies that changes in cloud cover are closely associated with the change in relative humidity.

The spring rainfall of the TGR is primarily triggered by low–mid-level cloud associated with large-scale fronts (Xia et al., 2016). The study of Mishra (2019) showed the occurrence of moderate rainfall can be mainly attributed to mid-level cloud, while light rainfall is associated with low cloud. This indicates that the reduction in cloudiness, especially in low–mid-level cloud, could be responsible for the decrease in light–moderate rainfall in WRF (Fig. 7b). The statistical relationship between low–mid-level cloudiness and rainfall frequency anomalies is significant, with a spatial correlation of 0.54 (0.58), and a temporal correlation of 0.89 (0.90), in the TGR, under SSP5-8.5 (SSP2-4.5). Thus, the projected decrease in rainfall frequency in WRF can be explained by the physical connection to cloudiness and relative humidity. This result is also consistent with the dynamical mechanism described above in section 3.3.1, i.e. the less upward motion would yield lower relative humidity and lower cloud cover.

Despite the cloud cover changes being generally similar in WRF and MIROC6, the change in rainfall frequency in MIROC6, however, cannot be explained by the cloud cover decline. This may be because of the differences in the rainfall schemes. MIROC6 has a limited ability to resolve convective-scale rainfall, because of its coarse resolution. Thus, MIROC6 implements a microphysics scheme and a convective parameterization scheme to characterize the rainfall. Meanwhile, the WRF downscaling approach used in the present study has a convection-resolving resolution (4 km in Domain3), which omits the need for a convection parameterization scheme. We check the frequency change in large-scale rainfall (derived from the microphysics scheme) and convective rainfall (calculated by the convective parameterization scheme) in MIROC6. The projected change of MIROC6 based on the microphysics scheme shows a slight decrease in rainfall frequency under the two SSPs (by −0.27/−0.34 days under SSP2-4.5/SSP5-8.5), which seems consistent with the decrease in WRF. The slight increase in the total rainfall frequency in MIROC6 seems dominated by the increase in the convective rainfall frequency (by 0.28/1.19 days under SSP2-4.5/SSP5-8.5). Thus, the differences in rainfall between WRF and MIROC6 are caused by the convective parameterization scheme of MIROC6. Previous studies have concluded that the convective scheme of MIROC6 (Chikira–Sugiyama) tends to overestimate the change in rainfall frequency in East Asia (Miyakawa et al., 2018; Tatebe et al., 2019), which is similar to our finding. Therefore, the projected increase in rainfall frequency in MIROC6 may not be realistic.

In this study, a triple-nested configuration of WRF is used to downscale the spring rainfall over the TGR for the present (1995–2014) under historical forcing, and future (2075–2094) under SSP scenarios. To obtain a more reasonable driving model of WRF, we first compare the simulated spring low-level winds of 10 GCMs with observations based on an MR index. Generally, MIROC6 exhibits the best performance among the 10 GCMs, and thus it is used for the boundary and initial conditions of WRF. Then, the performance of the downscaled rainfall is validated by comparing observations with the model outputs. The results demonstrate that, compared to MIROC6, WRF improves the simulation of regional rainfall and associated atmospheric circulation. Because of its finer resolution and better represented physical processes, WRF can simulate more realistic rainfall patterns and more regional features than MIROC6. In terms of the monthly evolution, the daily probability distribution, the diurnal cycle, and the duration of continuous rainfall, WRF also shows significant improvement against MIROC6. It generally reproduces more reasonable spatial distribution patterns by successfully eliminating the erroneous rainfall maximum simulated by MIROC6. For the reproduced CREs, WRF displays visible advantages against MIROC6. It presents a more reasonable regional average of CREs. WRF also better captures the observed regional features of CREs, and exhibits higher spatial correlations with observations.

Under the SSP scenarios, WRF projects an increase in rainfall amount in the central-west TGR and a decrease in the east, while the rainfall intensity is projected to become intensified over the entire region. As for the rainfall frequency, WRF suggests a clear decrease. In terms of the projected CREs, WRF suggests an overall decreased frequency over the entire TGR, while it exhibits a regional increase in the duration, rainfall amount, and intensity of CREs in the central-west TGR. These results may imply a higher risk of geological hazards in the central-west TGR in the future as far as one individual CRE event is concerned. Additionally, the TGR is an important water source for the South-to-North Water Diversion Project, so the decreased rainfall in the east may have decision-making and management implications for the project. In MIROC6, the change in rainfall frequency is in contrast to that in WRF and has obvious uncertainties. MIROC6 projects a slight change in frequency and an increase in the strength of CREs.

The inconsistency in the CRE trend between WRF and MIROC6 highlights the importance of detecting discrepancies in rainfall changes between models, and using downscaling to improve regional climate studies. In this study, the inconsistency can be attributed to the difference in rainfall frequency. The moisture budget diagnosis using ERA5 data shows that the rainfall frequency is dominated by the intensity of vertical moisture advection. Such a synoptic contribution is reproduced well by WRF, which is not the case in MIROC6, indicating large uncertainties in the CREs projected by MIROC6. In WRF, the decreased dynamic component of the vertical moisture advection term dominates over the increased thermodynamic component, leading to the weakening of upward motion and the decline in rainfall frequency. In MIROC6, meanwhile, the enhancement of vertical moisture advection causes the increase in rainfall frequency. In addition, WRF shows that warming will cause reduced relative humidity across the TGR, making it harder for the atmosphere to become saturated, thus reducing cloud formation (especially low and mid-level cloud). The decrease in cloud cover can physically explain WRF’s decreased rainfall frequency (Mishra, 2019). However, a similar reduction in cloud cover cannot explain the frequency changes in MIROC6, which is likely due to the convective parameterization scheme of MIROC6 overestimating the increase in rainfall frequency. Moreover, MIROC6 simulates too many rainy days (using a 1 mm d⁻¹ threshold), which is a common systemic bias in state-of-the-art GCMs (Feng et al., 2011; Li et al., 2019a; Zheng et al., 2022). These findings add uncertainties to the simulated rainfall frequency of MIROC6.

Changes in CREs can be useful to qualitatively evaluate the changes of geological hazards in the TGR (Tang et al., 2017). Furthermore, quantitative assessment of landslides based on a combination of CRE changes and field geotechnical parameters can greatly support disaster mitigation and prevention. Using a physical deterministic landslide model is an effective way to assess landslide changes, because it adequately considers the practical rainfall process and quantitatively evaluates the impact of rainfall intensity and duration on landslide stability. Thus, to conduct a quantitative landslide projection, it is necessary to couple RCMs and landslide models. We aim to do this in future work.

It is important to acknowledge some limitations in this study. First, a clear relation between the accuracy of a present-day simulation and a future projection is not guaranteed. Racherla et al. (2012) suggested that the skill of an RCM in simulating the mean climate conditions is not highly related with the skill in capturing future climate change. Second, to save computational resources, only one GCM is used to drive WRF. The bias of MIROC6 is another source of uncertainty because the RCM’s skill is strongly limited by the GCM used to drive it (Bao et al., 2015). Before using a GCM to drive a regional model, correcting the systematic biases of the GCM or employing a multi-model ensemble output as boundary conditions could be a necessary and affordable way to account for the above projection uncertainties (Ma et al., 2015a). Alternatively, directly correcting the biases or adopting machine learning to reduce the biases of downscaling results can also reduce the projection uncertainties (Zhu et al., 2008; Wang et al., 2021a).