Non-stationary dynamics of compound climate extremes: a WRF-CMIP6-GAMLSS framework for southeastern China

Understanding future changes in compound climate extremes (CCEs) is critical for climate risk assessment. Existing research, however, has largely relied on stationary assumptions, overlooking the dynamic evolution of CCEs under non-stationary climate change. To address this gap, this study employs an enhanced Generalized Additive Model for Location, Scale, and Shape (GAMLSS) framework to provide novel insights into the non-stationary characteristics of hot-wet (HW), hot-dry (HD), cold-wet (CW), and cold-dry (CD) extremes under future climate scenarios. We focus on the Minjiang River Basin (MRB) in Southeast China. A high-resolution dataset for CEE detection was generated by dynamical downscaling a bias-corrected CMIP6 dataset, using the Weather Research and Forecasting (WRF) model. Our results indicate that (1) CCEs increase significantly at a rate of 3.55 d per decade under the SSP5-8.5 scenario, with hot extremes (HW and HD) being the dominant contributors. Spatially, the increases exhibit a distinct west to east gradient, peaking in the downstream areas of the MRB. (2) Under the SSP5-8.5 scenario, CCEs exhibit a marked shift from stationary to non-stationary characteristics, with non-stationarity detected in 95.20 % of grid cells. This transition is primarily driven by mean warming, which explains 80.81 % of the change, rather than by variability. (3) The non-stationary results demonstrate that the severity and recurrence frequency of CCEs are systematically underestimated under stationary assumptions. Most CCE types (except for CD) show an increasing recurrence frequency under the SSP5-8.5 scenario. For instance, the frequency of events with a 100-year return period increases at a stronger trend of 3.12 d per decade. This study emphasizes the necessity of updating the frequency changes of CCEs under a non-stationary framework.

1 Introduction

Global warming is leading to more frequent and intense compound climate extremes (CCEs) (Sauter et al., 2023; Liu et al., 2024a; Zhang et al., 2024; You et al., 2025). CCEs have posed severe threats to global social, economic, and ecological systems, with impacts that surpass those of individual extremes in both range and severity (Mukherjee et al., 2023; Zeng et al., 2024; Miao et al., 2024). For example, the Yangtze River Basin in China experienced unprecedented compound hot-dry extremes in August 2022, characterized by record-breaking heatwaves and severe droughts, which directly affected over 50 million people (Jia et al., 2025). The Sixth Intergovernmental Panel on Climate Change (IPCC) report indicated that the probability and intensity of future CCEs are projected to increase (IPCC, 2021). Therefore, a systematic assessment of the future evolution of CCEs is critical for mitigating socio-economic risks and optimizing climate adaptation strategies.

Recent studies have increasingly focused on CCEs, highlighting their growing significance in the context of climate change. Zscheischler et al. (2018) were the first to clearly define the concept of compound events, emphasizing how the interaction of multiple climate and meteorological drivers can amplify extreme impacts. Building on this, Ridder et al. (2022) conducted the first global-scale assessment of the changes in compound events, specifically examining the co-occurrence of heatwaves and drought, extreme winds, and precipitation. Wu et al. (2023) revealed that under warming conditions, the risks associated with global compound pluvial–hot extreme events are projected to be significantly greater in the future than those observed during the historical period. Fang et al. (2025) investigated the future changes of sequential heatwaves and precipitation events as well as concurrent drought and heatwave events in China, with projections indicating an increase in both the frequency and intensity of these events.

While large-scale studies play a crucial role in advancing our understanding of global climate change and extreme events, their practical relevance for disaster risk management and adaptation strategies in medium- and small-scale regions is relatively limited due to their lower spatial and temporal resolution. To overcome this constraint, dynamical downscaling, which utilizes nested high-resolution regional climate models (RCMs), provides a critical technical pathway to investigate climate response mechanisms at fine-scales (Tapiador et al., 2020; Rahimi et al., 2024). Compared to Global Climate Models (GCMs), RCMs offer higher spatial resolution, allowing for more precise simulations of local climate effects induced by topography, such as local convective precipitation, orographic effects, and regional climate heterogeneity (Gilbert et al., 2025). In regions with complex terrain, RCMs are particularly effective at capturing spatial variations of climate variables, such as the differences in wind patterns, precipitation, and their distribution caused by topography in mountainous or basin areas (Imran and Evans, 2025). For example, Byun et al. (2023) assessed the ability of RCMs and GCMs to simulate storm tracks in East Asia, revealing that RCMs are better able to capture high-resolution topography, thereby reducing the biases found in GCMs. Lin et al. (2022) showed that RCMs driven by ERA-Interim reanalysis data are capable of capturing small-scale processes, such as orographic effects, and outperform GCMs in reproducing the large-scale features of the Heat Wave Magnitude Index-daily. Torrez-Rodriguez et al. (2023) also demonstrated that RCMs are better at reproducing the main spatio-temporal characteristics of precipitation in subtropical complex terrain regions. As an advanced convection-permitting RCM, the WRF model significantly enhances the simulation capability for meteorological processes at 1–10 km scales through its fully compressible, non-hydrostatic dynamic core framework (Talbot et al., 2012). This high-resolution simulation capability gives the WRF model a unique advantage in capturing small-scale meteorological phenomena. Zhou et al. (2024) developed a 9 km resolution regional reanalysis dataset covering the Tibetan Plateau based on the WRF model, and demonstrated its superior applicability compared to the fifth generation European Centre for Medium-Range Weather Forecasts Reanalysis (ERA5). Yang et al. (2024) revealed that the WRF model provides better accuracy in simulating snow depth during the cold season in high-elevation regions compared to ERA5-Land.

Additionally, traditional extreme event analyses rely on stationarity assumptions (when analyzing time series data, it is assumed that the statistical properties of the series remain constant over time), presuming that the probability and distributional parameters of climate variables are constant (Sun et al., 2018; Nerantzaki et al., 2023). However, driven by synergistic effects of global warming and anthropogenic forcing, extremes exhibit significant shifts in distributional characteristics (Gao et al., 2018). Therefore, traditional models fail to capture the non-stationary (the statistical properties of a time series change over time and do not remain constant) changes in these extreme events. Many studies have applied the Generalized Additive Models for Location, Scale, and Shape (GAMLSS) (Rigby and Stasinopoulos, 2005) to address non-stationary problems in hydrological and meteorological extremes, enabling updated risk analysis of evolving climate extremes. Lei et al. (2021) investigated the non-stationary changes of extreme precipitation in the Poyang Lake Basin and found that the stationary assumption underestimates the intensity of extreme precipitation in this region. Shao et al. (2022) innovatively proposed the Nonstationary Standardized Runoff Index (NSRI), and the results indicate that, compared to the stationary index, it can more accurately identify drought events. Salarijazi et al. (2023) evaluated the nonstationarity of maximum temperatures in multiple urban areas of Iran and concluded that traditional stationary approaches tend to underestimate the risk of annual maximum temperatures. However, existing non-stationary analyses only focus on individual extremes, and the potential non-stationarity of CCEs has not been established. The comprehensive assessment of future changes in CCEs within a non-stationary framework is also lacking.

Figure 1 Flowchart of CCEs projection in a non-stationary framework.

To address these research gaps, this study adopts a high-resolution approach by integrating the WRF model with GAMLSS. We first perform dynamical downscaling using the WRF model to refine the bias-corrected Coupled Model Intercomparison Project Phase 6 (CMIP6bc) data (1.25°×1.25°) to a 3 km resolution. Based on these high-resolution WRF outputs, CCEs are then identified and used as input for the GAMLSS framework to analyze their non-stationary characteristics. This methodology overcomes the limitations of traditional coarse-resolution models and addresses the shortcomings of stationary assumptions in CCE analysis. By focusing on the Minjiang River Basin (MRB) – a subtropical, monsoon-dominated basin in southeastern China where complex topographic–climatic interactions give rise to high-intensity compound hydroclimatic extremes (Gan et al., 2025; Geng et al., 2024; Wang et al., 2024) – this research aims to examine four types of CCEs: hot-wet (HW), hot-dry (HD), cold-wet (CW), and cold-dry (CD) events. The analysis proceeds as follows (Fig. 1): Supplement Sect. S1 presents the validation of CMIP6bc applicability. Section 3.1 characterizes the spatio-temporal patterns of CCEs under both a middle-of-the-road scenario (SSP2-4.5) and a high-emissions scenario (SSP5-8.5). The non-stationarity detection of CCEs is described in Sect. 3.2. The recurrence frequency changes in CCEs under non-stationary conditions is evaluated in Sect. 3.3. The work establishes a scientific basis for addressing the environmental and climatic challenges posed by CCEs, thereby contributing to effective strategies for regional sustainability and climate resilience.

Figure 2 Study area and model configuration. (a) Topographic features of the MRB (m) and (b) Model configuration with 9 km (D01) and 3 km (D02) nested domains (Zhang et al., 2025). Basemap source: Esri | Powered by Esri, https://services.arcgisonline.com (last access: 3 May 2026).

2 Study region, methods and data

2.1 Study region

The MRB is a complex topographic basin in southeastern coastal China (Fig. 2a). The Minjiang River, the main stream of the basin, drains an area of 60 992 km² – accounting for nearly half Fujian Province's territory. Encompassing three principal tributaries (Jianxi, Futunxi, and Shaxi rivers), the MRB experiences a subtropical monsoon climate characterized by 1700 mm mean annual precipitation and 18 °C mean temperature (Zheng et al., 2023). The basin displays spatio-temporal heterogeneity in precipitation, with flood seasons from April to September that often accompany CCEs. Particularly in the late flood season (July to September), the MRB experiences frequent typhoon-related compound disasters: the upper and middle reaches are commonly affected by typhoon-rainstorm-landslide events, while the lower reaches face high occurrences of typhoon-rainstorm-urban waterlogging and typhoon-rainstorm-flood events (Yang et al., 2025). In 2023, for example, typhoon Doksuri (No. 2305) caused approximately 66 794 people to be affected in Fuzhou, the downstream city of the MRB, with direct economic losses reaching 588 million RMB (Yan et al., 2024). In addition, the region also exhibits a climate characteristic of concurrent rainfall and heat, with CCEs frequently occurring during the warm season, driven by high temperature and heavy rainfall (Sun et al., 2025).

2.2 Data

Obtained from the Science Data Bank, the CMIP6bc dataset serves as the foundation for this investigation (Xu et al., 2021, https://www.scidb.cn, last access: 3 May 2026), which is constructed using the ERA5. This dataset incorporates an 18-model CMIP6 ensemble mean (Supplement Table S1), maintaining both climatological mean and interannual variability statistics while preserving nonlinear temporal trends. Compared with original CMIP6 data, CMIP6bc demonstrates superior performance in extreme event simulation. Despite its widespread application in previous studies (Jamal et al., 2023; Huang et al., 2024; Wu and Zheng, 2023), we conduct a 10-year validation over the MRB. Considering both the reliability of the dataset and the need to optimize computational resources (simulating one year over the MRB requires approximately four days on 80 CPU cores), we select a 10-year historical period (2005–2014) as sufficient to demonstrate the reliability of the CMIP6bc dataset for driving the WRF simulations. ERA5, as a widely used WRF-driven dataset (Arnault et al., 2021; Jiang et al., 2021; Varga and Breuer, 2022; Shang et al., 2022), is utilized here as a reference for the simulation results. Fujian Provincial Meteorological Bureau provides daily precipitation and temperature records from its 30 monitoring stations, which are used as validation data for this study. Future scenarios under SSP2-4.5 and SSP5-8.5 (2025–2065) are employed to project future changes in CCEs.

2.3 WRF model setup

This research utilizes the WRF Version 4.3 with two-domain nested configuration, featuring grid spacings of 9 and 3 km (Fig. 2b). Table 1 summarizes the optimal physics parameterization schemes selected through our comprehensive sensitivity experiments (Sect. S2) (Lin et al., 2023; Zhang et al., 2025). At sufficiently high model resolutions, deep convective processes can be explicitly resolved (Arakawa and Jung, 2011). Therefore, the cumulus parameterization scheme is deactivated in the inner domain (D02) to leverage convection-permitting capability. We first simulate daily precipitation and temperature over the MRB from 1 January 2005 to 31 December 2014, using both CMIP6bc and ERA5 forcing data. Subsequently, future projections from 1 January 2025 to 31 December 2065 are conducted using CMIP6bc under two climate projection scenarios.

Table 1 Settings for WRF model in this study.

2.4 Definition of CCEs

This study considers four types of CCEs: hot-wet events (HW), hot-dry events (HD), cold-wet events (CW) and cold-dry events (CD). We adopt the widely used thresholds (the 90th and 10th percentiles) to identify CCEs (Croitoru et al., 2016; Song et al., 2019; Patel et al., 2024). We first extract daily precipitation (>0.1 mm) and temperature data for each grid during 2025–2065, defining the 90th and 10th percentiles as thresholds to identify hot/cold and wet/dry extremes, respectively. Specifically, we define extreme temperature events as occurring when daily temperatures are higher (hot extremes) or lower (cold extremes) than the threshold. Wet events are characterized by rainfall surpassing the threshold (90th), while dry events are characterized by seven consecutive days without rainfall. For dry conditions defined as seven consecutive days without precipitation, each day within a dry spell is treated as a dry day in the counting of CCEs. We conduct calculations independently for each grid point, applying thresholds specific to each point. Specifically, for interannual variation, we sort the precipitation and temperature data over a 40-year period and determine the thresholds based on the 10th and 90th percentiles to identify CCEs. For seasonal variations, we separately extract the precipitation and temperature data for the summer (JJA) and winter (DJF) seasons, applying the same sorting method to calculate the respective thresholds, thus analyzing the distribution characteristics of CCEs for each season.

2.5 GAMLSS model

GAMLSS is a flexible statistical model used for analyzing distributions with non-stationary characteristics (Rigby and Stasinopoulos, 2005). It extends the traditional generalized linear models (GLMs) and generalized additive models (GAMs) by introducing joint modeling of all distribution parameters (location, scale, and shape). Unlike traditional regression models, GAMLSS effectively characterizes both linear and nonlinear dependencies linking predictors to response variables (Stasinopoulos and Rigby, 2007).

This study employs the semi-parametric GAMLSS, which accommodates parametric terms, nonparametric smooth functions, and random effects within a unified modeling structure (Gao et al., 2018). Consider z independent samples $y_i(i=\mathrm{1},\mathrm{\dots },z)$ following a distribution F_y(y_iθ_i), where the parameter vector ${\mathit{\theta }}_{iT}=({\mathit{\theta }}_{i\mathrm{1}}{\mathit{\theta }}_{i\mathrm{2}},\mathrm{\dots },{\mathit{\theta }}_{ik})$ contains k components representing location (Mn), scale (Var), and shape (skewness and kurtosis), with k normally not exceeding 4. Model selection is performed using Akaike's Information Criterion (AIC) (Akaike, 1974), with the optimal configuration identified through minimum AIC values, and model fitting quality is assessed by the Filliben correlation coefficient (Filliben, 1975). The GAMLSS is formally defined as follows:

$$g_k\left({\mathit{\theta }}_k\right)={\mathrm{\varnothing }}_k{\mathit{\beta }}_k+{\sum }_{j=\mathrm{1}}^{j_k}{h}_{j_k}\left(x_{j_k}\right)$$

where k denotes the indicator of distribution parameters, θ_k is the distribution parameter vector, ∅_k represents n×j_k matrix of covariate variables, β_k is the coefficient vector of length j_k. g_k(⋅) is the link function connecting distribution parameter to linear predictor. $h_{j_k}(\cdot )$ defines how the distribution parameter varies with covariate variable $x_{j_k}$. In the GAMLSS model, time (year) is used as the independent variable (x), and the number of days per year for each type of CCE is treated as the dependent variable (y), thereby enabling the calculation of the non-stationary characteristics of each CCE. We employ two types of GAMLSS models to capture potential changes in the distribution of meteorological variables. The first is the traditional two-parameter location–scale model (mean μ and variance σ), which assumes a fixed distributional shape. The second is a more flexible four-parameter location–scale–shape model (mean μ, variance σ, skewness ν, and kurtosis τ), allowing the distributional shape to evolve over time. These extended models retain the mean and variance parameters while incorporating two additional shape parameters to capture asymmetry and tail behavior. We assess the models by comparing their goodness of fit and select the one that best represents the data distribution for subsequent analyses. Regarding the calculation of stationarity, we consider the CCEs to be stationary if both the mean and variance remain stable. If either the mean or variance shows significant variation, the CCEs are considered non-stationary. Table S2 lists all the distribution functions implemented in our study. The R code for implementing GAMLSS model can be accessed at https://github.com/gamlss-dev/gamlss (last access: 3 May 2026).

Figure 3 Annual Spatio-temporal patterns of CCEs across the MRB from 2025 to 2065.

3 Results

3.1 Spatio-temporal patterns of CCEs under future scenarios

Figure 3a–j illustrate the spatial distribution of annual CCEs over the MRB during 2025–2065. Overall, total CCEs are higher under the SSP5-8.5 scenario (45.48 d) than under the SSP2-4.5 scenario (42.58 d). Dry-related extremes (HD and CD) dominate across the basin, whereas wet-related extremes (HW and CW) are less frequent. Both hot extremes occur more frequently under SSP5-8.5 than SSP2-4.5: HW increases from 0.84 to 1.33 d, and HD rises from 22.69 to 27.45 d. In contrast, cold extremes decrease: CW declines from 3.17 to 2.15 d, and CD from 15.88 to 14.19 d.

Spatially, the two scenarios exhibit similar geographic patterns. The highest frequencies occurring in downstream regions, particularly for HD and HW. HD shows a wider spatial extent, extending into the Futun River and Jian River Basins. Meanwhile, CW and CD display a clear west-to-east increasing gradient, with the highest values concentrated near the Jiufeng Mountains.

Temporally (Fig. 3k–o), CCE trends diverge sharply between scenarios. Under high-emission SSP5-8.5, total CCEs increase significantly at 3.55 d per decade, whereas SSP2-4.5 projects stabilized frequencies. Hot extremes (HW and HD) increase more rapidly under SSP5-8.5, at rates nearly double those under SSP2-4.5. Conversely, the decline in CD is stronger under SSP2-4.5 – approximately 1.7 times that under SSP5-8.5.

Figure 4 Spatial patterns of CCEs in summer (JJA) and winter (DJF) across the MRB from 2025 to 2065.

Building on the analysis of annual changes, we further investigate the seasonal variations of CCEs. As shown in Fig. 4a–j, the spatial distribution of summer CCEs largely resembles the annual pattern, except for HW, which remain infrequent in this season. Nonetheless, a marked increase in HW under the SSP5-8.5 is observed in the downstream MRB. In summary, total CCEs increase during summer under both scenarios, rising more rapidly under SSP5-8.5 (2.26 d per decade) than under SSP2-4.5 (0.79 d per decade). Moreover, differences between scenarios are more pronounced for hot-related extremes (HW and HD) in summer than on the annual scale, whereas cold-related events (CW, CD) show consistent patterns (Fig. 5a–e).

In winter, CCEs exhibit a contrasting spatial distribution (Fig. 4k–t). They occur predominantly in the western MRB, particularly concentrated in the Futun River Basin, with higher frequencies under SSP5-8.5 (11.87 d) than under SSP2-4.5 (11.51 d). This shift is primarily driven by hot-related extremes (HW and HD), which move from the downstream MRB toward western mountainous areas. Meanwhile, cold-related extremes (CW and CD) remain most frequent in the Jiufeng Mountain areas. Winter CCEs show an increasing trend under SSP5-8.5 (0.76 d per decade) but a slight decreasing trend under SSP2-4.5 (Fig. 5f–j). Among individual event types, wet-related extremes (CW and HW) show insignificant trends, whereas HD increases (1.18 and 0.46 d per decade under SSP5-8.5 and SSP2-4.5, respectively) and CD deceases (−0.50 and −0.70 d per decade under the two scenarios).

Figure 5 Temporal patterns of CCEs in summer (JJA) and winter (DJF) across the MRB from 2025 to 2065.

Given the importance of precipitation and temperature as key climate indicators, we calculate their basin-averaged annual values and interannual trends (Supplement Fig. S5). Precipitation varies only slightly, maintaining relatively stable annual fluctuations. In contrast, temperature shows a distinct upward trend, which accelerates under high-emission SSP5-8.5 conditions (0.46 °C per decade). These results suggest that changes in CCEs across the MRB are primarily governed by temperature-driven physical processes – where intensifying hot extremes coincide with declining cold extremes – a conclusion supported by earlier studies (Wu et al., 2020; Zhao et al., 2024; Duan et al., 2024).

Figure 6 Stationary and non-stationary characteristics for CCEs in the MRB (a–e, k–o), percentage of non-stationary and stationary characteristics across five basins (f–j, p–t).

3.2 Non-stationary characteristics of CCEs

To quantify the non-stationary characteristics of CCEs, we analyze variations in both the mean (Mn) and variance (Var) using GAMLSS (Fig. 6). The model demonstrates excellent fitting performance for all indices except HW, as indicated by Filliben coefficients exceeding 0.95 in Fig. S6.

Our results reveal a pronounced shift from stationary to non-stationary behavior in CCEs under SSP5-8.5 compared to SSP2-4.5. This shift is primarily driven by non-stationarity in Mn, which governs 80.81 % of grid cells across the MRB. An additional 11.07 % of grids are influenced by the combined effects of both Mn and Var, primarily distributed in the Shaxi River Basin and the downstream MRB. Dry-related extremes (HD and CD) show a transition from Var to Mn non-stationarity. For HD, the area dominated by Mn non-stationarity expands from 22.65 % to 57.03 % of grids, becoming the prevailing pattern (covering 60 % of grids) in downstream regions where stationarity was previously dominant (54.71 %). Similarly, Mn non-stationarity for CD expands to cover nearly the entire downstream MRB. For wet-related extremes (HW and CW), both Mn and Var non-stationarity increase notably. The area dominated by Mn non-stationarity for HW rises from 34.30 % to 54.02 % of grids, while for CW, it undergoes a more substantial increase from 0.06 % to 23.25 %. Overall, CCEs under SSP5-8.5 display more widespread non-stationarity, with dry-related extremes primarily governed by changes in Mn, whereas wet extremes are influenced by a combination of Mn and Var effects.

Figure 7 Results of Mann-Kendall test for Mn (a–j) and Var (k–t), showing the spatial distribution of Z values. Z values indicate trend significance: $\left|Z\right|>\mathrm{1.65}$ denotes 90 % confidence, while $\left|Z\right|>\mathrm{2.58}$ corresponds to 99 % confidence.

Figure 7 further details the spatial patterns of Mn and Var in CCEs. Compared to Var, Mn exhibits more pronounced and spatially extensive variations. For hot-related extremes (HW and HD), Mn increases significantly across the entire basin under both SSP2-4.5 and SSP5-8.5 (at the 99 % confidence level), indicating that climate warming predominantly amplifies the mean frequency of hot-related extremes rather than their temporal variability. For cold-related extremes, CD exhibits more stronger changes than CW, which are also predominantly driven by reduction in Mn. Under SSP5-8.5, Mn for most CCEs increases significantly across nearly the entire basin, while under SSP2-4.5, it remains relatively stable. In contrast, Var exhibits only minimal changes under both scenarios.

Figure 8 Comparison of non-stationary (NS) and stationary (S) characteristics for CCEs under 20-, 50-, and 100-year return periods (a–j). Spatial distributions of trends in CCEs under 100-year return periods (k–t), 20- and 50-year return period result are provided in Fig. S7. Blank areas indicate grid points that failed to pass the 99 % confidence test.

3.3 Changes in the recurrence frequency of CCEs

To assess changes in the recurrence frequency of CCEs, we fit both stationary and non-stationary models. As shown in Fig. 8a–j, substantial differences emerge between the return periods estimated under stationary versus non-stationary assumptions. Our results reveal that stationary models systematically underestimate the future recurrence frequency of CCEs, whereas non-stationary models reasonably capture their time-evolving characteristics. This underestimation gap widens notably over time.

Under SSP5-8.5, non-stationary projections indicate significant increases in the frequency of CCEs corresponding to 20-, 50-, and 100-year return periods. For instance, total CCEs with a 100-year return period are projected to increase at a rate of 3.12 d per decade. Stationary models consistently underestimate CEE frequency after around 2045. Specifically, HW, HD, and CW all show positive trends in recurrence frequency under SSP5-8.5 with increases of 0.36, 11.65, and 0.30 d per decade, respectively. Under SSP2-4.5, the corresponding increases are lower (0.15, 2.53, and 0.06 d per decade). In contrast, CD exhibits a decreasing trend, with a more rapid decline under the SSP2-4.5 (−3.38 d per decade) than under SSP5-8.5. The 100-year return events display greater sensitivity to climate change, highlighting amplified non-stationary effects on these high-impact extremes.

We further quantify the deterministic trends in recurrence frequency using Empirical Mode Decomposition (EMD, Sect. S1). Figure 8k–t present the variations in the 100-year recurrence frequency (complete results for all return periods are provided in Fig. S7). Overall, CCE frequency shows strong upward trends under SSP5-8.5, with significant spatial heterogeneity. Frequency generally decreases from west to east, with notably high-value areas concentrated in the Shaxi River Basin (exceeding 12 d per decade in some grids). Among all indices, HD exhibits the most severe increase, with basin-wide trends exceeding 8 d per decade. CD is the only index showing a decreasing trend under both scenarios, with a less pronounced decline under SSP5-8.5, suggesting that greenhouse warming has only partially offset the reduction in CD recurrence frequency.

4 Discussion

Although earlier research has emphasized the necessity of analyzing extreme events under non-stationary conditions (Cheng et al., 2014; Byun and Hamlet, 2020; Liu et al., 2024b), the evolution of CCEs within a non-stationary climate remains lacking. Our study develops an innovative non-stationary framework integrating WRF-based dynamical downscaling with GAMLSS to reassess future recurrence frequency of CCEs. The results indicate that traditional stationary models systematically underestimate the frequency of CCEs, highlighting the critical importance of incorporating time-varying risk assessments to avoid misleading projections and to inform robust climate adaptation strategies (Abdelmoaty and Papalexiou, 2023). This innovative framework enables regional-scale reassessment of CCEs which is transferable to other regions.

4.1 Dominance of hot extremes and temperature-driven shifts

The projected increase in CCEs, particularly under SSP5-8.5, aligns with global trends of intensifying hydroclimatic risks under continued warming (Asadieh and Krakauer, 2017; Zhang et al., 2021b; Shu et al., 2024). Our findings indicate that hot-dry extremes dominate both spatially and temporally, increasing at 2.26 d per decade in summer under SSP5-8.5, while cold extremes decline. This pattern is consistent with studies highlighting the rising prevalence of hot-stagnation and hot-dry extremes in East Asia (Yin et al., 2025). In particular, across China, hot–dry extremes are projected to double by 2050 under high-emission scenarios (Yao et al., 2024). Guo et al. (2023) also revealed that under global warming, extreme hot conditions are projected to dominate most regions of China, with some areas experiencing more than 50 extreme hot days. Moreover, there is evidence suggesting that both compound hot–dry and compound hot–wet events are projected to increase significantly under future climate scenarios. Notably, compound hot–dry events exhibit substantially higher frequency, longer duration, and greater intensity than compound hot–wet events (Fang et al., 2025). The reversal between hot and cold extremes has been robustly linked to enhanced radiative forcing from anthropogenic greenhouse gas emissions (Samset et al., 2018; Kramer et al., 2021). Our analysis further reveals that temperature – not precipitation – is the primary driver of CCE changes in the MRB, as evidenced by the strong warming trend (0.46 °C per decade under SSP5-8.5) alongside relatively stable precipitation (Fig. S5). This supports the hypothesis that thermodynamic effects, rather than dynamic ones, dominate mean-state changes in extremes (Horton et al., 2016; Van Der Wiel and Bintanja, 2021). However, it should be emphasized that precipitation responses are often manifested through pronounced changes in intensity and intermittency, rather than in long-term totals alone (Du et al., 2022; Maity and Maity, 2022). In particular, warming can lead to a redistribution of precipitation toward shorter duration but more intense events, even in regions where mean precipitation exhibits statistically insignificant trends (Fowler et al., 2021). For HW events, these changes are closely linked to both thermodynamic and moisture-related processes. As global and regional temperatures rise, the atmosphere's moisture holding capacity increases at an approximate rate of 6–7 % °C⁻¹, consistent with Clausius–Clapeyron scaling (O'Gorman and Muller, 2010; Ali et al., 2021). Previous studies have demonstrated that such moistening effects can substantially intensify CCEs, even in regions where average precipitation remains relatively unchanged (Gimeno et al., 2022; Zhang et al., 2024).

4.2 Non-stationarity: mean shifts outpace variability

A key advance of this study is the explicit detection of non-stationary characteristics in CCEs, which has been largely overlooked in prior compound event analyses. We find that under SSP5-8.5, 95.20 % of grid cells exhibit non-stationarity, predominantly driven by changes in the Mn rather than Var, accounting for 80.81 % of the transitions. This suggests that climate warming amplifies extremes primarily through shifts in baseline intensity – a thermodynamic effect – rather than through increased temporal variability. Similar findings have been reported at global scales, where mean warming dominates changes in extreme temperature distributions (Patel et al., 2024; Nordling et al., 2025). The spatial concentration of Mn-driven non-stationarity in downstream MRB and the Shaxi River Basin may reflect localized warming amplification due to urban heat islands or land-atmosphere feedbacks, a phenomenon noted in other subtropical regions (Gao et al., 2018; Wu et al., 2020).

4.3 Frequency of recurrence systematically underestimated by stationary model

Our comparison between stationary and non-stationary model indicates that the latter captures a significant increase in recurrence risks, particularly for 100-year CCEs, which are projected to rise by 3.12 d per decade under the SSP5-8.5 scenario. Stationary model systematically underestimates recurrence frequency after 2045, and the stronger non-stationary response of 100-year events highlights the heightened vulnerability of high-impact, low-probability extremes. Additionally, the west-to-east gradient in recurrence frequency, with hotspots in the Shaxi River Basin, may be attributed to topographic and land-surface heterogeneity, which modulate local hydroclimatic responses (Zheng et al., 2023). The projected increase in the frequency of 100-year CCEs has important implications for flood control and water resources management in the MRB, particularly in the middle and lower reaches and densely populated urban centers such as Fuzhou. The basin's complex topography and highly heterogeneous land cover amplify the spatial variability of rainfall and runoff, while low-lying urban areas in Fuzhou are particularly vulnerable to both pluvial and fluvial flooding. More frequent HW events are likely to produce intense, concentrated rainfall, rapidly depleting reservoir flood control capacity and increasing the likelihood of spillway activation, while simultaneously exacerbating urban inundation in flood prone districts (He et al., 2024). HD events may reduce river storage, but if followed by subsequent heavy rainfall, they can increase uncertainty in peak flows and place additional stress on flood management systems (Hariharan Sudha et al., 2024). Previous studies have similarly shown that conventional extreme value models fail to capture the escalating severity of extreme events under climate change. Singh et al. (2021) demonstrated that under climate change, precipitation and temperature exhibit a non-stationary dependent structure, and treating them independently can substantially underestimate the occurrence of compound extremes. Feng et al. (2020) compared the recurrence probabilities of floods under stationary and non-stationary conditions and found that, under non-stationary conditions, the annual variability of floods is significantly greater. Xu et al. (2025) evaluated changes in global marine heatwaves and found that stationary models significantly underestimate their frequency, intensity, and duration. In summary, these studies collectively highlight that accounting for non-stationarity is crucial for accurately assessing various extreme events.

4.4 Methodological limitations

Our integrated “bias-corrected CMIP6–WRF dynamical downscaling–GAMLSS” framework represents a significant methodological advancement over the direct use of purely statistical downscaling for projecting CCEs. By resolving mesoscale circulations and explicitly simulating convective processes, our approach more faithfully captures the fine-scale spatiotemporal heterogeneity of precipitation and temperature fields in complex terrain, a capability that statistical methods, reliant on historically derived statistical relationships, fundamentally lack (Gutmann et al., 2012; Rahimi et al., 2024). Nevertheless, certain limitations persist. Even at convection-permitting resolution (3 km), the WRF model exhibits systematic biases in simulating orographic precipitation, a well-documented challenge often stemming from uncertainties in microphysical parameterization schemes and the representation of land-atmosphere energy and moisture exchanges over mountainous regions (Talbot et al., 2012; Zhang et al., 2025). Specifically, the overestimation of precipitation in high-elevation regions by WRF may artificially enhance the frequency of wet-related compound extremes, whereas the underestimation of precipitation in lowland and downstream areas may bias the detection of dry-related compound extremes. Because the identification of CCEs relies on precipitation and temperature thresholds and their co-occurrence, systematic precipitation biases can further influence threshold estimation and event classification, thereby affecting the spatial distribution and frequency of detected CCEs. To address this limitation in future work, we plan to explore machine learning based post-processing approaches (Yin et al., 2021; Xie et al., 2023), in which the residuals between historical WRF simulations and station observations are learned and subsequently transferred to future simulations in a physically constrained method. This approach offers a promising pathway to reduce systematic regional model biases while preserving the internally consistent climate change signal. Furthermore, while statistically robust, our current non-stationary GAMLSS framework employs time merely as a proxy covariate for climate change. While this formulation is widely adopted and effective for detecting long-term trends, time itself serves only as an indirect proxy for the underlying physical processes driving changes in CCEs (Ragno et al., 2019). From a physical perspective, the evolution of compound hot–dry events are governed by a combination of thermodynamic and dynamic mechanisms, including background warming, shifts in large-scale atmospheric circulation, and land–atmosphere feedbacks (Bevacqua et al., 2022; Zhang et al., 2021a; Tian et al., 2024). For example, global mean surface temperature (GMST) can serve as a physically meaningful indicator of anthropogenic thermodynamic forcing, directly linking greenhouse gas increases to enhanced surface heat stress (Gillett et al., 2021). In this subtropical monsoon dominated basin, variations in the intensity, westward extension, and persistence of the Western Pacific Subtropical High (WPSH) can regulate subsidence strength, cloud cover, and surface radiative forcing, thereby favoring persistent hot conditions while suppressing precipitation (Li et al., 2024; An et al., 2025). Changes in the strength and variability of the East Asian summer monsoon (EASM) influence large-scale moisture transport and rainfall timing, affecting both the onset and persistence of dry spells (Park et al., 2020; Dou et al., 2025). In addition, antecedent soil moisture has been shown to influence the persistence and amplification of heat extremes through soil moisture–temperature feedbacks (Jiang and Wang, 2024), further highlighting the value of including land-surface states as covariates. Therefore, incorporating physically based covariates into the GAMLSS framework may therefore improve the interpretation of non-stationarity.

5 Conclusions

Through this intensive case analysis, we establish a transferable framework for assessing the non-stationarity of CCEs. This work advances the understanding of the evolution of CCE recurrence frequency under climate change and offers important perspectives to support adaptive strategies and strengthen disaster risk governance. The main findings are summarized as follows:

1. CCEs increase significantly across the MRB, with a stronger trend under SSP5-8.5 (3.55 d per decade) than under SSP2-4.5. HD extremes dominate both spatially (downstream-focused) and seasonally (summer-peaked), rising at 2.26 d per decade, whereas cold extremes decline. These shifts are primarily temperature-driven, as pronounced warming amplifies hot-related extremes but suppresses cold-related extremes.

2. CCEs shift markedly toward Mn-dominated non-stationarity under SSP5-8.5 scenarios, contrasting with the largely stationarity under SSP2-4.5. Spatially, Mn non-stationarity governs 80.8 % of the MRB under SSP5-8.5, with dry-related extremes (HD and CD) showing the most abrupt transitions. For HD, Mn non-stationarity expands from 22.7 % to 57.0 % of the basin and dominates 60 % of downstream grids – an increase of nearly threefold relative to SSP2-4.5. For CD, Mn-driven shifts cover >90 % of the downstream MRB. Var contributes minimally in both scenarios, confirming that warming amplifies extremes primarily through shifts in mean intensity rather than through increased variability.

3. Non-stationary modeling reveals that stationary approaches systematically underestimate future CCEs recurrence frequency. Under SSP5-8.5, most CCEs types (except CD) exhibit increasing recurrence frequency. Climate change impacts are significantly amplified for 100-year return period events, which rise at 3.12 d per decade – a response highlighting their heightened non-stationary sensitivity. Spatially, recurrence frequency shows a distinct east-to-west gradient, with a significantly increase occurring in the western mountainous areas.