Global bias-corrected data produced by CESM version 1 that participated in the CMIP5 was obtained from the website of the Research Data Archive (https://rda.ucar.edu/datasets/ds316.1/ #!description). The dataset was provided in the Intermediate File Format specific to WRF and contained all the variables needed for the initial and boundary conditions for WRF simulations. This dataset had already been bias-corrected using the European Centre for Medium-Range Weather Forecasts (ECMWF) Interim Reanalysis (ERA-Interim) fields for 1981–2005, following the method of Bruyère et al. (2014). All the variables had 26 pressure levels and were provided in files at 6-hourly intervals. Files were available for a simulation for the present (1951–2005) and three future scenarios (RCP4.5, RCP6.0, and RCP8.5) spanning 2006–2100.

We ran the WRF model to generate hourly meteorological fields of the “Basecase” (2006–2015) and the three RCP cases (RCP4.5, RCP6.0, RCP8.5) for the winter months (December, January, and February) of 2046–2055. The simulation domain covered the entirety of China with a spatial resolution of 36 km × 36 km (Figure 1). The physical options used to drive the WRF simulations are listed in Table 1 (Hu et al., 2016), including the new Thompson microphysics scheme (Thompson et al., 2008), Rapid Radiative Transfer Model (RRTM) longwave radiation (Mlawer et al., 1997), and Yonsei University PBL parameterization (Hong et al., 2006). Due to corrupted data caused by a modelling bug in CESM, bias-corrected data was unavailable for RCP2.6. Thus, we chose three scenarios (RCP4.5, RCP6.0, and RCP8.5), spanning the winter of 2046–2055, to drive WRF v3.8 for future analyses. Results from the CESM model using the RCP6.0 emission scenario were used as the meteorological boundary and initial fields of the WRF model for the 2006–2015 periods. The WRF model was configured with 44 vertical layers from surface to 50 hPa (Supplementary Figure S1). Because of the vast territory of China, climate and air quality vary in the different regions. To better compare and summarize changes in different regions, we divided China into seven regions based on differences in climate and administration (Figure 1), including East China (E), North China (N), Central China (C), South China (S), Southwest China (SW), Northwest China (NW), and Northeast China (NE).

We adopt the ASI definition proposed by Wang et al. (2018), which is shown in Eq. 1. This parameterization is derived using PBLH and WS10 data to fit the normalized daily PM_2.5 concentrations during the winter. This could better explain haze events and has been used in other studies (Zhang et al., 2019; Gao et al., 2020). For each grid cell (i, j), if the daily cumulative precipitation is <1 mm and the daily average PBLH and WS10 satisfy Eq. 1, the air is considered stagnant and the given day is considered to be an air stagnation day (ASD). We define an ASE as an event in which air stagnation lasts for 3 days or more. P B L H i , j < 0.759 ∗ exp ( − 0.6 ∗ W S 10 i , j ) + 0.264 (1)

Evaluation of model performance is an important step to verify the reliability of modelling results and to establish confidence for further application. The WRF simulations were validated using the fifth-generation European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis data for global climate and weather (ERA5). Monthly averaged data (Hersbach et al., 2019) from ERA5 for 2006 to 2015 were obtained from the website https://cds.climate.copernicus.eu/cdsapp#!/dataset/. In this study, statistical metrics and benchmarks recommended by the United States EPA were used. The normalized mean deviation (NMB), normalized mean error (NME), and the Pearson correlation coefficient (R) of the observed and simulated values were calculated to quantitatively evaluate model performance (Emery et al., 2017). The corresponding equations are shown below: N M B = ∑ ( P j − O j ) O j × 100 (2) N M E = ∑ | P j − O j | O j × 100 (3) R = ∑ [ ( P j − P ¯ ) × ( O j − O ¯ ) ] ∑ ( P j − P ¯ ) 2 × ∑ ( O j − O ¯ ) 2 (4)

Figure 2 shows the spatial distribution of wintertime (DJF) meteorological fields from the Basecase and the corresponding ERA5 reanalysis data. Detailed statistical metrics are listed in Table 2. Note that PBLH, WS10, and T2 were averaged from 2006 to 2015, while PCP was calculated as the average of cumulative winter precipitation. As can be seen from Figure 2, the simulated PBLH in the Tibetan Plateau and its surrounding areas was significantly underestimated, which might be due to the complex terrain in the region (Wang et al., 2021). Combined with Table 2, the mean PBLH within the model domain was 490.4 m in the ERA5 data, while it was 477.92 m according to WRF. The WRF model underestimated PBLH with a mean bias of −12.48 m, while the RMSE was as large as 172.22 m. This indicated that model performances may vary a lot across the different regions. Nevertheless, the simulated PBLH showed a strong positive correlation of 0.84 with ERA5 data, suggesting that the overall model performance was acceptable. The WRF model successfully captured the general pattern of PCP as compared to the ERA5 data, but underestimated PCP in southern China. Previous studies have reported that ERA5 precipitation data are relatively higher than the ground observations over China in winters, with a mean bias of 35% (Jiang et al., 2021; Jiao et al., 2021), which partly explains the discrepancy between WRF simulations and ERA5 data. A comparison of the simulated PCP with the Tropical Rainfall Measuring Mission (TRMM) satellite data for 2006–2015 (Supplementary Figure S2) revealed that the WRF-simulated PCP had similar spatial patterns to that of TRMM data. The WRF-simulated mean PCP within the model domain was 15.0 mm, which was close to that of TRMM data (15.5 mm), giving us confidence in the simulated PCP.

The WS10 and T2 (the temperature at 2 m) predicted by WRF showed relatively consistent spatial patterns with those from ERA5 data. WS10 was slightly overestimated with a mean bias of 0.95 m/s. The ME and RMSE of WS10 were both less than 2 m/s, and thus, met the performance criteria suggested by Emery et al. (2001). The model performance for T2 was the best among the selected meteorological factors with correlation coefficients as high as 0.99. To explore the differences in modelling results before and after downscaling, we compared the CESM output with ERA5 data and WRF simulations (Supplementary Figure S3; Supplementary Table S1). Because bias-corrected CESM output does not provide PBLH and PCP (Bruyère et al., 2014; Bruyère et al., 2015), only projected WS10 and T2 obtained. As Supplementary Figure S3 shows, both WRF and CESM exhibited good spatial agreement with the ERA5 data. As compared to CESM, WRF could reproduce more regional scale details for both WS10 and T2, especially in the Tibetan Plateau and the Sichuan Basin. Previous studies have indicated that CESM cannot capture topography-induced temperature patterns (Liu et al., 2013; Chen et al., 2018), which is consistent with our results. The biases between CESM and ERA5 were 1.07 m/s and 0.95 K for WS10 and T2, respectively, which was slightly larger than those between WRF and ERA5.

The decadal mean values of PBLH, WS10, T2, and PCP under the three RCP scenarios were calculated and compared with the Basecase to obtain their spatial and temporal variations due to climate change. Note that PCP was calculated as the average cumulative winter precipitation. Figure 3 shows the WRF-projected changes in PBLH, WS10, T2, and PCP under RCP4.5, RCP6.0, and RCP8.5 scenarios relative to the present period; relative changes in these parameters in the different regions are depicted in Figure 4. Wintertime PBLH at the national scale was projected to increase for the 2046–2055 period under the three RCP scenarios. The increase in PBLH was more obvious for the RCP6.0 scenario (15.1 m), followed by RCP8.5 (6.4 m) and RCP4.5 (2.1 m). A significant increase mainly occurred in the Tibetan Plateau and some parts of North and Northeast China under the three RCP scenarios. For Central, East, and South China, the WRF-RCP4.5 and WRF-RCP8.5 projected generally shallow PBLH, while WRF-RCP6.0 showed a weak increase.

The projected changes in PCP in China were more varied among different scenarios. The WRF-RCP4.5 projected wetter conditions in South and East China and slightly drier conditions in Central, North, and Northeast China for the 2046–2055 period. Under RCP6.0, PCP showed a decreasing trend in most of China except Northwest and Northeast China. A significant decrease was projected for Central China with a regional mean of 26.65% (−14.13 mm). Under RCP8.5, however, wetter conditions were projected for Central, East, and South China, with a maximum increase of 23.80% (15.1 mm). Significant decreases were found for Southwest and Northeast China, with the regional mean value changes of −2.35 and −1.38 mm, respectively. The PCP in North China showed a consistent decreasing trend under the three RCP scenarios, however, the changes were different for other regions.

For WS10, the national mean value for the 2046–2055 period remained largely unchanged relative to the present 2006–2015 period under the RCP4.5 and RCP6.0 scenarios, while it increased slightly under RCP8.5. The projected WS10 showed an increasing trend under all three scenarios over Southwest China, while the largest increase was found under RCP8.5 (5.93%). For other regions, the projected changes in WS10 were less than 2.02% under the RCP4.5 and RCP6.0 scenarios. As can be seen from Figure 3, the projected wintertime temperatures showed a significantly increasing trend over China under all three scenarios. Due to the highest anthropogenic greenhouse gas emissions in RCP8.5, the projected T2 showed the greatest increase, followed by RCP6.0 and RCP4.5. The projected changes in temperature presented obvious regional differences. The largest warming was found in the Tibetan Plateau, with a maximum of 3.6°C increase in average temperature under RCP8.5. The regional mean relative changes in T2 ranged from 0.18% to 0.56%, and the maximum occurred in Southwest China (0.56%).

The CESM-projected WS10 and T2 distributions over China under the three scenarios for the 2046–2055 period were also compared with distributions for the 2006–2015 period (Supplementary Figures S4, S5). Generally, there are no obvious differences between the CESM and WRF projections. The WRF-projected mean values of changes in WS10 and T2 under the three scenarios were comparable to those obtained from CESM; however, the latter covered a wider range, which is attributable to the finer spatial resolution of WRF that allowed it to better capture regional-scale climate patterns. As compared to the CESM simulations, the WRF projections showed a decrease in WS10 in South and Northeast China under the RCP6.0 and RCP8.5 scenarios.

The spatial distributions of ASD and ASE in the Basecase and their changes under the three RCPs over China are shown in Figure 5. The projected total ASD during the 2006–2015 period was 167 days, while the changes during the 2046–2055 period under the three RCPs aere 6 days (RCP4.5), −9 days (RCP6.0), and −7 days (RCP8.5), respectively. The maximum (minimum) changes projected were 55 (−43), 37 (−83), and 49 (−68) days under RCP4.5, RCP6.0, and RCP8.5, respectively. Under the RCP4.5 scenario, the projected number of ASD increased in most of China, while it decreased significantly in Northwest, Southwest, and Northwest China under the RCP6.0 and RCP8.5 scenarios. Significant reductions in the number of ASD were also found in Sichuan Basin under RCP8.5, with a minimum of about −68 days. The projected ASE exhibited an increasing trend from south to north during the 2006–2015 period, with a national mean value of 56.9. The changes in ASE under the three scenarios showed similar spatial distributions to those of ASD. Positive changes were mainly projected for East and Central China, and the eastern part of Northwest and Southwest China under RCP4.5. Negative changes were found in most of China under the RCP6.0 and RCP8.5 scenarios. At the national scale, the projected ASE decreased in the future for the three scenarios, with the maximum decrease occurring in RCP8.5 (−4.0), followed by RCP6.0 (−2.9) and RCP4.5 (−0.5).

The probability distribution of ASE and the changes in ASE over the 2046–2055 period relative to 2006–2015 was calculated using the total ASE during the winter of the 10 years for each grid cell in China; the distribution is shown in Figure 6 to identify the changes of air stagnation under the different scenarios. The numbers of ASE occurrence in China during the winter of the 10 years were mainly concentrated in 30–75 events (for points with a frequency greater than 1%). The projected frequency of ASE occurrence exhibited a bimodal distribution, with a major peak at around 65–75 times and a minor peak at around 40–55 times. Both the peaks shifted to a less frequent ASE under the three RCP scenarios, which means that the number of wintertime ASE occurrences in China may decrease in the future. As can be seen from Figure 6B, the changes in ASE in the 2046–2055 period relative to the 2006–2015 period follow a normal distribution under the three RCP scenarios. The mean changes of winter ASE are −1, −3, and −4 times for RCP4.5, RCP6.0, and RCP8.5, respectively. This is related to the fact that the mean value of PBLH shows an increasing trend in the future, while the mean values of WS10 and PCP change slightly.

Long-lasting air stagnation events can increase the risk of public exposure to high concentrations of air pollutants. Figure 7 shows the relative frequency distribution of the duration of ASD under the different scenarios for all grid cells over China. For example, assuming that an ASE would last for 3 days, the occurrence under different scenarios was in the order of Basecase (28%) > RCP6.0 (26%) > RCP8.5 (23%) > RCP4.5 (22%). As can be seen from Figure 7, the longest duration of winter ASE in China was projected to be 14 days. The frequency of ASE lasting 2–4 days was the highest (27%–29%) in the Basecase, but with the increase in duration, the proportion of Basecase gradually decreased. The proportion of the RCP scenario gradually increases, especially under the RCP8.5 scenario, which means that the probability of the occurrence of long-lasting ASE will increase in the future.

The changes in ASE lasted for more than 5 days for the 2046–2055 period relative to the present 2006–2015 period in different regions are given in Figure 8. Generally, the occurrence of long-lasting ASE in different regions of China increased by 3–11 times in the mid-century decade. The largest increase was projected unser the RCP4.5 scenario for most regions. Central China was found to have the largest increase among the different regions, followed by East China and Northeast China. The lowest increase was projected for South China and was majorly ascribed to a significant increase in PCP.