Impacts of anthropogenic aerosols on a snowfall event – a case study in the Guanzhong Basin and its surrounding areas, China

Impacts of anthropogenic aerosols on clouds and snowfall during winter precipitation events remain highly uncertain, particularly under heavy pollution. A winter snowfall event over the Guanzhong Basin (GZB) and its surrounding regions (GZBs), China, has been simulated using a cloud-resolving, fully coupled WRF-Chem model to quantify the respective roles of aerosol–radiation interactions (ARIs) and aerosol–cloud interactions (ACIs). The simulated temporal variation and spatial distribution of air pollutants and precipitation generally agree with the observations in the GZB + GZBs. Sensitivity experiments are performed to evaluate effects of ARIs and ACIs by changing the anthropogenic emissions. The precipitation response to ARIs and ACIs exhibits regional contrast in GZB and GZBs due to different aerosol concentrations. In the GZB, exclusion of ARIs leads to a slight increase in precipitation with increasing emissions, mainly associated with enhanced ice-phase precipitation induced by ACIs. ARIs increase the precipitation in the GZB when emissions increase reaches a threshold, caused by ARI-induced enhancement of relative humidity (RH) which increases ice water path and favors survival of falling ice particles. In contrast, precipitation in the GZBs decreases with increasing emissions, reflecting suppression of liquid-phase precipitation by ACIs and reductions in RH caused by ARIs. In addition, changes in anthropogenic emissions exert limited influence on the spatial distribution of precipitation across the combined GZB–GZBs region. These findings provide process-level insight into how ARIs and ACIs regulate snowfall under polluted conditions, with implications for improving aerosol–precipitation coupling in regional climate and weather models.

1 Introduction

Anthropogenic aerosols can directly and indirectly modulate the dynamics, thermodynamics, and microphysics of the cloud-precipitation process by absorbing or scattering the solar radiation and serving as cloud condensation nuclei (CCN) or ice nuclei (IN) (e.g., Boucher et al., 2013; IPCC, 2013; Huang and Ding, 2021; Zhao et al., 2024), constituting one of major uncertainties in climate prediction and weather forecasting (e.g., Quaas, 2015; Myhre et al., 2013; Makar et al., 2015; IPCC, 2021). Aerosol effects on the solar radiation are referred as to aerosol–radiation interactions (ARIs), which traditionally include the direct and semi-direct effects (e.g., Charlson et al., 1992; Sun and Zhao, 2021; Ackerman et al., 2000). Aerosol effects on CCN and IN and further cloud-precipitation are referred as to aerosol-cloud interactions (ACIs), which is also called aerosol indirect effect (e.g., Mitchell, 1971; Braham, 1974; Twomey, 1977; Albrecht, 1989).

Past studies have proposed many mechanisms about aerosol effects on precipitation, such as the positive precipitation-aerosol relationships attributed to invigoration effect (e.g., Khain et al., 2005; Pinsky et al., 2013; Rosenfeld et al., 2008), the negative aerosol effect on precipitation due to cloud water competition (e.g., Albrecht, 1989; Rosenfeld, 1999), and ARIs that reduce solar radiation down to the surface and modify atmospheric temperature profiles (e.g., Ackerman et al., 2000; Koren et al., 2004). The two mechanisms proposed for the invigoration effect include (1) mixed- or cold-phase invigoration, whereby aerosol induced higher CCN suppress warm-rain formation, resulting in greater lofting of cloud condensate mass and increased fusion heating as the droplets freeze (e.g. Rosenfeld et al., 2008), and (2) warm-phase invigoration, whereby aerosol induced higher CCN increase condensation heating (e.g. Kogan and Martin, 1994; Pinsky et al., 2013). It is worth noting that the invigoration effect is still debatable due to insufficient observational evidence or flawed methodology (e.g., Öktem et al., 2023; Romps et al., 2023; Varble et al., 2023).

However, the net effect of aerosols on clouds and precipitation remains controversial and depends on many factors, including meteorological conditions, specific aerosol loadings, cloud/precipitation types, cloud/precipitation development stages, aerosol composition and size distribution, relative location of aerosol and cloud vertical locations, and orography conditions (e.g. Khain et al., 2008; Tao et al., 2007; Guo et al., 2018, 2014; Zhang et al., 2002; Ackerman et al., 2000; Yang and Li, 2014). First, the aerosol effect on precipitation depends on different types of precipitation. For example, aerosols can enhance convective precipitation but inhibit stratiform precipitation, which lead to intensified flooding in southern China in summer and aggravated drought in East Asia' s middle and low latitudes area especially during winter and spring time (e.g., Menon et al., 2002; Zhang et al., 2017; Huo et al., 2021; Xie et al., 2022; Ryu and Min, 2022; Sun et al., 2022). Second, aerosol effect on precipitation varies with precipitation grade. Studies show that increasing aerosol concentrations can enhance heavy rain but reduce moderate and light rain, which may also reduce the frequency of the moderate and light rain (e.g., Wang et al., 2011; Shao et al., 2022). Third, the specific impact of aerosols on precipitation varies across different phases of precipitation process due to the competition between ARIs and ACIs (e.g., Lee et al., 2016). Furthermore, the influence of aerosols on precipitation also relies on the position of aerosols. As an example, low-level anthropogenic aerosols inhibit the precipitation process while high-level aerosols enhance it (e.g., Bai et al., 2020; Sun et al., 2022). Most recently, Yun et al. (2024) have found that aerosol effects on MCS©\ averaged precipitation is mostly similar between the polluted and clean conditions due to compensations among the subregions with divergent spatiotemporally discretized changes. Wang et al. (2024) have emphasized the aerosol IN effect on cloud microphysics, convection intensity, and regional rainfall distribution. Peng et al. (2025) have pointed out that aerosols exert stronger impacts on convective than stratiform precipitation as mentioned in Sun and Zhao (2021) based on data analysis. Moraglia and Crippa (2026) have also indicated that differences induced by land-use perturbation and aerosols is smaller in stratiform rain than that in convective rain because the absence of strong vertical mixing prevents the activation of more aerosols.

Given the complexity of precipitation phase states, there has been relatively limited exploration into how aerosols influence wintertime precipitation processes, especially snowfall event (e.g., Saleeby et al., 2010; Thompson and Eidhammer, 2014; Guo et al., 2021; Cao et al., 2021; Zhang et al., 2022; Xiao et al., 2022; Pravia-Sarabia et al., 2023). Saleeby et al. (2010) have showed that an increase in aerosol concentrations results in a spillover of snowfall from the windward slope to the leeward slope, but the impact is largely controlled by the synoptic-scale flow and available moisture. Guo et al. (2021) have indicated that increased aerosol number concentrations have less impact on snowfall at the initial period due to the lower supersaturation produced by the urban heat and dry island but increase snowfall in the later period due to higher supersaturation caused by the feedback of the cloud microphysics. Based on sensitivity simulations of a winter extreme mixed-phase storm, Pravia-Sarabia et al. (2023) have indicated that the total precipitation remains basically unchanged, with maximum changes of 5 % while the production of snow is largely altered.

Guanzhong Basin (GZB) is surrounded by the Loess Plateau to the northwest and Qinling Mountain to the south (Fig. 1). Due to fast-growing industries and city expansions, the aerosol loadings inside the basin have been elevated during the past 3 decades, particularly during wintertime (e.g., Bei et al., 2016a, b, 2017b). Aerosol effects on precipitation over the mountain area and inside the GZB have been investigated based on observational data at Mt. Hua's summit and its surrounding area (e.g., Rosenfeld et al., 2007; Yang et al., 2013a, b). Rosenfeld et al. (2007) have found that the precipitation at Mt. Hua in the south of GZB could be decreased by 30 % to 50 % during hazy conditions. Yang et al. (2013a, b) have further pointed out that the decreasing trend of precipitation is correlated well with deterioration of the air pollution at Mt. Hua and inside the GZB, supporting the hypothesis that both aerosol microphysical and radiative effects could reduce precipitation. Most recently, based on sensitivity simulations of a summertime short-term heavy rainfall event under different aerosol scenarios, Bei et al. (2025) have shown that the synergetic effect of ARIs and ACIs consistently decreases the precipitation both in the whole simulation domain and GZB with increasing aerosols, but ARIs play a more important role in the decreasing trend of the precipitation with deterioration of PM pollution as mentioned in Zhou et al. (2024).

Figure 1 (a) WRF-Chem simulation domain with topography and (b) Guanzhong basin (GZB) and its surrounding areas (GZBs). In (b), the black dots denote the meteorological sites with rain gauge and the blue dots denote centers of cities with air pollutants observations and the size of blue circles denotes the number of ambient monitoring sites of cities. The white contour outlines the GZB and the area outside of GZB in (b) is referred to as GZBs.

The present study aims to investigate the synergetic effects of ARIs and ACIs on a snowfall event occurred in the GZB and GZB surrounding areas (GZBs) using a fully coupled cloud-resolving WRF-Chem model. The model and experiment design are described in Sect. 2. The main results and discussion are presented in Sect. 3. The conclusions are given in Sect. 4.

2 Model, data, and experiments design

2.1 WRF-Chem model configuration and experimental design

A specific version of WRF-Chem model, incorporating original contributions by Grell et al. (2005) and modifications by Li et al. (2010, 2011a, b, 2012), is employed to investigate the impact of anthropogenic aerosols on a snowfall event occurred in the GZB and GZBs from 26 to 28 January 2022. We use Goddard shortwave module developed by Chou and Suarez (1999) and Chou et al. (2001) to consider the ARIs. The ACIs are evaluated using a two-moment bulk microphysics scheme developed by Morrison et al. (2009). Detailed model description of the WRF-Chem model, including the aerosol radiative module to calculate aerosol optical properties and the aerosol microphysical module to consider activation of aerosols to CCN and IN, can be found in the Supplement (SI, Sects. S1–S3).

The model is configured with one domain using the grid spacing of 6 km for a total of 500×500 grid points in the east-west and south–north directions, respectively, and the domain is centered at Xi'an (34.25° N, 109° E) (Fig. 1a). Figure 1b shows the study area, which is divided into two parts: GZB and its surrounding areas (GZBs). We use GZB + GZBs to represent the study area. The vertical direction is divided into 35 layers extending from the surface to 50 hPa, with a stretched vertical grid spacing ranging from 30 m near the surface, to 500 m above 2.5 km to achieve a finer vertical resolution within the planetary boundary layer. The meteorological initial and boundary conditions are derived from 6-hourly 1°×1° NCEP-FNL (National Centers for Environmental Prediction final operational global gridded analysis). The chemical initial and boundary conditions are interpolated from the 6 h output of a global chemical transport model for ozone and related chemical tracers (MOZART) (Horowitz et al., 2003).

The model is integrated for a 126 h period from 12:00 UTC 23 January to 18:00 UTC 28 January 2022 and the simulations from 16:00 UTC 25 January to 16:00 UTC 28 January 2022 are analyzed. In order to examine the effect of aerosols with different concentrations on the 3 d snowfall event, we first change the atmospheric aerosol concentrations via scaling the anthropogenic emissions of all species in the model. A set of 19 anthropogenic emission scale factor (AESF) is used in sensitivity simulations, ranging from 2⁻³ to 2³ with an exponential increasing step of 1/3 (Table S2). It is worth noting the use of a wide range of AESFs might not represent realistic emission-control scenarios. In the study, the main purpose of applying a wide range of AESFs is to systematically explore the nonlinear responses of aerosol–cloud–precipitation interactions across aerosol environments ranging from relatively clean to heavily polluted conditions. Wintertime PM_2.5 concentrations in northern China exhibit substantial variability associated with haze development, boundary-layer evolution, and unfavorable meteorological conditions, often varying by more than one order of magnitude between clean and polluted periods (Li et al., 2021). In the GZB, frequent severe winter haze events have been reported under stagnant meteorological conditions and complex basin topography, with PM_2.5 concentrations commonly exceeding 200–300 µg m⁻³ during polluted episodes (Bei et al., 2017b). As shown in Fig. S1 in the Supplement, the average wintertime [PM_2.5] have decreased from about 160 µg m⁻³ in 2013 to 65 µg m⁻³ in 2026 in response to emission-control measures, highlighting the large variability of aerosol loading in this region. Therefore, the selected AESF range is designed to span the observed variability of aerosol conditions and to identify potential threshold and nonlinear behaviors of ARIs and ACIs under different pollution regimes. Then we design two groups of experiments to verify the contribution of ARIs, ACIs and the both to precipitation. In the first group sensitivity simulations, both ARIs and ACIs are considered with the AESF ranging from 2⁻³ to 2³ (hereafter referred as to F_BASE). The second group of sensitivity simulations is the same as the F_BASE, but the ARI effect is not considered (hereafter referred as to F_ARI0). In the F_ARI0 experiment, the aerosol effect on solar radiation is not considered in the shortwave radiation transfer module in the WRF-Chem model, such that aerosols no longer modify atmospheric heating rates and further influence meteorological fields. The model setup is the same for all experiments, except for the anthropogenic emission amplitude and the option of ARIs (on or off). In the F_BASE, the member with the AESF of 1.0 is referred to as the benchmark simulation (CTRL), which is used to validate the model performance.

2.2 Model validation and statistical metrics

Hourly precipitation observational data for the GZB and GZBs stations are provided by the China Meteorological Administration (CMD) service center (http://data.cma.cn/, last access: 3 January 2026), and hourly observations of PM_2.5, O₃, NO₂, and SO₂ are released by the Ministry of Ecology and Environment of China. Three statistical metrics, including the normalized mean bias (NMB), root mean square error (RMSE), and the index of agreement (IOA) are used to evaluate the model performance. The population mean (p-mean hereinafter) of a given variable across all qualified grid points is applied to verify the general impact of aerosols on cloud properties. Detailed description about NMB, RMSE, IOA and p-mean can be found in Sect. S4.

3 Results and discussion

3.1 Snowfall event overview

A heavy snowfall event occurred from 26 to 28 January 2022 is investigated in GZB + GZBs, with the total precipitation of 6.4 mm, and the maximum daily precipitation of 2.3 mm is on 28 January. Meanwhile, despite the occurrence of the snowfall event, the observed near-surface PM_2.5 concentrations ([PM_2.5]) are still high during 36 d, with the average of 108.5 and 74.1 µg m⁻³ in GZB and GZBs, respectively. The snowfall event is basically the result of the interaction between the northerly cold air and southerly and easterly warm air in low-level atmosphere, as a cold front system proceeds from the North China (25 January) to the south-east (26–27 January) and then south-west direction (28 January) mainly along the coast (Fig. S2). GZB + GZBs is located to the southwest (25 January), the west (26 January) and the northwest (27–28 January) side of the cold front, which is not the main area of the precipitation that occurred near the cold front. During 26 to 28 January, there is clear evidence of weak cold air incursion in the GZB + GZBs at 925 hPa (Fig. S2a–c). At 850 and 700 hPa (Fig. S2d–i), the weak cyclonic wind shear and positive vorticity are maintained in the GZB during the event. At 500 hPa (Fig. S2j–l), GZB + GZBs is situated in the front of a very shallow trough, which is favorable for the development of low-level convergence, and eventually leads to the snowfall event in this area. The low-level water vapor is brought to the snowfall area mainly by the southwesterly and easterly airflow from Bay of Bengal and the South China Sea. The snowfall event lasts approximately 72 h from 00:00 LT 26 January to 00:00 LT 29 January 2022. The observed hourly snowfall rate (Fig. 4) shows a multi-peak structure, with the maximum instantaneous rate reaching ∼0.21 mm h⁻¹ around 12:00 LT 28 January. The majority of the precipitation at the surface is in the form of snow, as the near-surface temperature during the event remains between −4 and 0 °C (Fig. S3). Despite the continuous snowfall, the observed PM_2.5 concentrations ([PM_2.5]) remain high, with an average of 108.5 and 74.1 µg m⁻³ in GZB and GZBs (Fig. 2a), respectively, providing a unique polluted environment to investigate aerosol-snowfall interactions.

Figure 2 Pattern comparisons of simulated (color counters) vs. observed (colored dots) average near-surface mass concentrations of (a) PM_2.5, (b) O₃, (c) NO₂ and (d) SO₂ from 26 to 28 January 2022. The white contour outlines the GZB and the area outside of GZB in panels is referred to as GZBs.

Figure 3 Comparison of observed (black dots) and simulated (solid red lines) diurnal profile of near-surface hourly mass concentrations of (a) PM_2.5, (b) O₃, (c) NO₂, and (d) SO₂ averaged at monitoring sites in GZB + GZBs from 26 to 28 January 2022.

3.2 Model evaluations

Figure 2 shows the pattern comparison of simulated PM_2.5, O₃, NO₂, and SO₂ concentrations averaged during the study period with observations at monitoring sites in GZB and GZBs. The model generally reproduces the spatial distribution of air pollutants when comparing to observations. The simulated [PM_2.5] are generally more than 75 µg m⁻³ in the GZB, but the SO₂ level is not high in GZB and GZBs, less than 20 µg m⁻³. The model tends to underestimate O₃ and overestimate NO₂ concentrations in the GZB when comparing to observations. During wintertime, near-surface O₃ concentrations are strongly influenced by NO titration effects, and the overestimation of NO₂ concentrations might partly contribute to the underestimation of O₃ concentrations in the simulations. The NO₂ overestimation might be associated with uncertainties in the anthropogenic emission inventory or PBLH simulations. Figure 3 presents the simulated and observed time series of average hourly mass concentrations of PM_2.5, O₃, NO₂, and SO₂ at monitoring sites in the GZB + GZBs from 26 to 28 January 2022. The model reasonably simulates the temporal variation of PM_2.5, O₃, and NO₂ concentrations, with the IOA exceeding 0.60. However, the model does not show good performance in simulating the SO₂ temporal variation, with an IOA of 0.25. The larger discrepancy for SO₂ mainly results from the sensitivity of point-source SO₂ simulations to meteorological uncertainties. The transport and dispersion of SO₂ from point sources (i.e., power plants and industrial facilities) are highly sensitive to the accuracy of simulated wind fields (Hu et al., 2022). In the GZB, complex terrain and weak synoptic forcing further increase the uncertainty in meteorological simulations (Bei et al., 2017b). Additionally, parameterizations of SO₂ oxidation, including aqueous-phase and heterogeneous pathways, under the persistent high-RH conditions during this snowfall event also introduces further uncertainties (Travis et al., 2025). In general, the model slightly underestimates the air pollutants concentrations against observations, with the NMB ranging from −8.4 % to −2.3 %. These underestimations may partly result from uncertainties in anthropogenic emissions, meteorological fields, and chemical or physical parameterizations in the model, which can influence air pollutant accumulation and secondary formation processes during wintertime (Bei et al., 2017a; Cheng et al., 2021; Li et al., 2025).

Figure 4 provides the temporal variation of precipitation rates averaged at meteorological sites in GZB + GZBs from 26 to 28 January. The model generally reproduces the temporal evolution of the observed hourly precipitation, with NMBs of −3.3 %, −7.4 %, and −1.3 % and IOAs of 0.85, 0.75, and 0.78 over GZB + GZBs, GZB, and GZBs, respectively. The observed five rainrate peaks are successfully reproduced in all regions, but the underestimation of the first and fifth rainrate peak is considerable in GZB + GZBs. Compared with GZB, the model exhibits better overall performance in GZBs, while a larger underestimation is found over GZB during periods of intense precipitation. These underestimations are likely related to the model's sensitivity to initial conditions and its performance during the spin-up and decay stages of the event, respectively, especially given that the precipitation occurs on the periphery of the cold front where large-scale forcing is weak (Liu et al., 2019; Lin et al., 2023). Additionally, the model also produces higher light precipitation during 27–28 January compared to the observations. As shown by Liu et al. (2024), the lack of a shallow convection scheme at gray-zone resolutions (i.e., 6 km) allows excessive lowlevel moisture to accumulate and prematurely convert to light precipitation under stratiform conditions. These localized discrepancies do not compromise the model's overall performance. The pattern comparison of the 3 d accumulative precipitation over the GZB + GZBs is shown in Fig. 5. The model generally replicates the precipitation distribution compared to the observations (Fig. 5a). For all stations in GZB + GZBs, the simulated accumulated precipitation agrees reasonably well with observations, with an NMB of −3.7 % and an IOA of 0.85 (Fig. 5b). When evaluated separately, the model yields NMBs of −10.2 % and −1.3 % over GZB and GZBs, respectively, indicating that the underestimation is mainly concentrated within the basin region. Nevertheless, the relatively high IOAs of 0.89 and 0.85 over GZB and GZBs suggest that the model adequately captures the spatial variability of precipitation in both regions.

Figure 4 Comparison of observed (black dots) and simulated (solid red lines) diurnal profile of hourly rain rate averaged at monitoring sites in (a) GZB + GZBs, (b) GZB, and (c) GZBs from 26 to 28 January 2022.

Figure 5 (a) Pattern comparisons of simulated (color counters) vs. observed (colored dots) accumulative precipitation, and (b) scatter plot of simulated and observed accumulative precipitation in GZB, GZBs, and GZB + GZBs from 26 to 28 January 2022. The white contour outlines the GZB, and the area outside of GZB in (a) is referred to as GZBs.

Figure 6 Average (a) near-surface PM_2.5, (b) black carbon mass concentration, (c) AOD at 550 nm, and (d) absorbing AOD at 550 nm in GZB + GZBs from 26 to 28 January 2022, as a function of the scale factor of anthropogenic emissions.

3.3 Effects of ARIs on meteorological fields in GZB + GZBs

Figure 6a and b provides the variation of [PM_2.5] and black carbon concentrations ([BC]) averaged in the GZB + GZBs during the study period with increasing AESF. [PM_2.5] and [BC] generally increase with increasing AESF linearly. The ARI effect tends to hinder development of planetary boundary layer (PBL) and decrease the PBL height (PBLH), increasing the near-surface concentration of air pollutants. Sensitivity results show that when the AESF increases from 0.125 to 8.0, the PBLH in GZB + GZBs decreases from approximately 450 to 300 m, a reduction of about 34 % based on model output statistics (see Table S2 for detailed values). The quasi-linear relationship between [PM_2.5] and AESF reflects that the variation of [PM_2.5] with the AESF is mainly influenced by the emission of primary aerosols and the precipitation change caused by ARIs and ACIs. In addition, it is worth noting the decreased PBLH is also favorable for the formation of secondary aerosols, including sulfate, nitrate, ammonium, and secondary organic aerosol (SOA). However, secondary aerosols are generally water-soluble, and easily removed by wet deposition or activated to form cloud droplets. The aerosol optical depth (AOD) and absorbing AOD show the similar increasing relationship with the AESF in GZB + GZBs as that of [PM_2.5] and [BC]. The average AOD and AAOD at 550 nm in GZB + GZBs are 0.58 and 0.071 in the CNTL case, respectively, and the single scattering albedo at 550nm is about 0.88, indicating a moderately absorbing component within an otherwise scattering-dominated aerosol population over GZB + GZBs.

Figure 7 Profile variation of average (a) air temperature, (b) vertical velocity, and (c) water vapor over GZB + GZBs from 26 to 28 January 2022 caused by ARIs, as a function of the scale factor of anthropogenic emissions.

ARIs include direct scattering and/or absorbing of incident solar radiation by atmospheric aerosols and the induced adjustments of the surface energy budget, thermodynamic profile and cloudiness (IPCC, 2013). Therefore, aerosols in the atmosphere reduce the solar radiation down to the Earth surface and further sensible heat flux to the atmosphere, lowering the temperature of low-level atmosphere. Furthermore, light-absorbing aerosols also heat the atmosphere. Figure 7a shows the effect of ARIs on the temperature profile averaged in GZB + GZBs during the study period by differentiating the F_BASE and F_ARI0 under various aerosol conditions. It should be noted that the temperature anomalies shown in Fig. 7a represent the total atmospheric response to ARIs, including both the direct aerosol radiative forcing and the subsequent dynamical adjustments, such as changes in vertical motion, heat transport, and boundary-layer mixing. ARIs decrease the air temperature of the low-level atmosphere, and the air temperature decrease is mainly concentrated below 1.8 km and becomes increasingly significant with decreasing height or increasing AESF (Fig. 7a). Previous studies have reported that absorbing aerosols cause a warming effect above the PBL (Ding et al., 2016; Gao et al., 2016; Wilcox et al., 2016; Wu et al., 2025). Consistent with these studies, the ARI-induced temperature response in the present simulation is characterized by cooling below approximately 1.8 km and weak warming aloft (Fig. 7a). The lower-atmospheric cooling is mainly caused by reduced surface solar radiation which decreases the surface sensible heat, whereas the elevated warming is associated with aerosol absorption. The resulting thermal contrast enhances lower-tropospheric stability and suppresses turbulent mixing. Consequently, the average PBLH over GZB + GZBs decreases from about 450 to 300 m with the AESF increasing from 0.125 to 8.0 in the F_BASE, but the warming layer is generally above 1.8 km and can extend to 7–8 km with the AESF exceeding 0.5. It is worth noting that the vertical distribution of aerosol-induced heating does not necessarily coincide with the vertical distribution of absorbing aerosols. As discussed by Wu et al. (2025), radiative heating depends not only on aerosol loading but also on the vertical structure of radiative fluxes, aerosol optical properties, and air density, particularly as well as the decreased surface sensible heat. Therefore, the weak warming above approximately 1.8 km shown in Fig. 7a should be interpreted as the integrated response of aerosol–radiation interactions rather than a direct reflection of aerosol concentration at a given altitude. During the study period, the profile of vertical velocity shows occurrence of downdraft within the PBL and updraft above the PBL (Fig. S4a). The cooling effect of ARIs within the near-surface atmosphere induces a downward motion and meanwhile an upward motion is generated above the PBL, which is resulted from the warming effect caused by absorbing aerosols and enhanced convergence by the ARI-induced downward motion (Fig. 7b). The vertical velocity response shown in Figure 7b does not directly mirror the temperature anomalies at each height. Instead, it reflects a dynamical adjustment to the vertically contrasting radiative forcing. Aerosol scattering cools the lower atmosphere through surface dimming, while aerosol absorption produces weak warming above approximately 1.8 km. The resulting thermal contrast favors ascent above the PBL, while downward motion develops within the boundary layer due to the combined effects of enhanced lower-tropospheric stability and compensating subsidence associated with the ascent aloft. Consequently, ARIs strengthen downward motion in the lower atmosphere and upward motion aloft. These ARI-induced variations in divergence and water vapor flux are shown in Fig. S5. ARIs induce divergence (positive values) below approximately 700 m and convergence (negative values) above this level, with the convergence strengthening as AESF increases. Correspondingly, negative water vapor flux anomalies occur within the divergent layer, indicating moisture export from the lower troposphere, while positive anomalies occur within the convergent layer aloft, indicating moisture accumulation. This circulation pattern contributes to the drying of the PBL. Therefore, the upward motion above the PBL is enhanced by ARIs and more warm and moist air in the lower layer is being transported upward. The ARI-enhanced downward motion within the near-surface atmosphere is subject to transporting more air in the upper layer downward, causing a dry and cold effect, as shown in Fig. 7c.

Figure 8 Average (a) cloud droplet CND, (b) variation of cloud droplet CND due to ARIs, (c) ice crystal CND, and (d) variation of ice crystal CND due to ARIs over GZB + GZBs from 26 to 28 January 2022, as a function of the scale factor of anthropogenic emissions.

3.4 Response of cloud microphysic properties to aerosols

The impact of aerosols on cloud and precipitation is mediated by the combined effects of ACIs and ARIs (Twomey, 1977; Albrecht, 1989). Figure 8a presents the dependence of the column number density (CND) of cloud droplets and ice crystals over GZB + GZBs during the study period on the AESF. Increasing anthropogenic emissions escalate aerosol concentrations in the atmosphere, providing more CCN and IN to activate to form cloud droplets and ice crystals. Without ARIs, the CND of cloud droplets increases monotonically with increasing AESF, and the impact of increasing AESF on the CND of cloud droplets is not very significant, i.e., when the AESF increases from 0.125 to 8.0 or anthropogenic emissions increase by 64 times, the CND increases by 9.2 times. However, the CND of ice crystals is not sensitive to increasing anthropogenic emissions. The different response of cloud droplets and ice crystals to change of anthropogenic emissions is mainly caused by their formation at different height. The cloud water is mainly formed in the atmosphere below 3 km, which is more susceptible to the influence of anthropogenic emissions (Fig. S4b). However, the formation of ice crystal is concentrated in the atmosphere between 4 and 7 km, which is less affected by anthropogenic emissions due to downward motion within the PBL and continuous wet deposition processes (Fig. S4c).

Furthermore, the CND of ice crystals generally decreases with increasing AESF, which is not consistent with the hypothesis that increasing anthropogenic emissions provide more IN to promote the formation of ice crystals. Since the formation of ice crystals is sensitive to relative humidity (RH), the possible reason for the decreasing trend with the AESF is that ACIs alter the RH in the atmosphere due to increasing anthropogenic emissions. In the F_ARI0, the simulation using an AESF of 2⁻³ is referred to as the reference simulation (REF0). For discussion convenience, we define the RH relative to water surface as RHW and the RH relative to ice surface as RHI. Figure S6 presents the impact of ACIs on the RHW and RHI profile averaged in GZB + GZBs during the study period by differentiating all members in F_ARI0 with REF0. ACIs generally decrease RHW and RHI in the atmosphere and the decrease in RHW and RHI becomes significant with increasing AESF. Therefore, decreasing RHI with increasing AESF causes the decreasing trend of the CND with the AESF. In addition, the ACI-induced decrease in RHW provides a physical explanation for the relatively slow increase of cloud droplet CND with AESF (Fig. 8a), as lower RHW makes it more difficult for the additional CCN (from increased emissions) to activate into cloud droplets, thus limiting the rate of CND increase.

Figure 9 Average (a) CWP, (b) variation of CWP due to ARIs, (c) IWP, and (d) variation of IWP due to ARIs over GZB + GZBs from 26 to 28 January 2022, as a function of the scale factor of anthropogenic emissions.

When ARIs are considered, the CND of cloud droplets shows an increasing trend with the AESF when the AESF is less than 2.5, and when the AESF is more than 2.5, the CND decreases with the AESF, which is inconsistent with the Twomey effect (Twomey, 1977). In addition, ARIs decrease the CND at the same AESF, and with the AESF exceeding 2.5, the CND decrease is more than 20 % (Fig. 8b). As shown in Fig. 7, ARIs adjust the profile of air temperature and specific humidity, further modifying the profile of RH. Figure S7a shows impacts of ARIs on the RHW profile averaged in GZB + GZBs during the study period by differentiating the F_BASE and F_ARI0. ARIs decreases RHW in the atmosphere below 3 km where the cloud droplets are formed (Fig. S4b), and the RHW decrease becomes increasingly significant with increasing AESF. Decreased RHW is not favorable for the formation of cloud droplets, so with the same AESF, ARIs decrease the CND of cloud droplets, as shown in Fig. 8b. ARIs progressively decrease RHW with increasing AESF, continuously decreasing activity ability of CCN. Therefore, although increasing anthropogenic emissions provide a large amount of CCN, most of CCN can not be activated due to decreased RH. There exist a threshold AESF (2.5) beyond which the CND commences to decrease (Fig. 8a). For ice crystals, impacts of ARIs on the CND is generally not significant with the AESF less than 1.0 (Fig. 8c), and the CND variation due to ARIs is in the range between -0.5% and 0.5 % (Fig. 8d). When the AESF is more than 1.0, ARIs consistently increase the CND of ice crystals with increasing AESF. The main reason is that with the AESF exceeding 1.0, ARIs increase RHI in the atmosphere between 4 and 7 km, enhancing the formation of ice crystals (Fig. S7b).

In the F_ARI0 without ARIs, cloud water path (CWP) decreases with increasing AESF. Many studies have proposed that increased CCN reduce cloud particle sizes to decrease the efficiency of collision and collection and further hinder auto conversion of cloud water to rainwater, causing more cloud water to exist in the atmosphere (e.g., Li et al., 2008, 2009). The decreasing trend of CWP with the AESF is caused by the decrease in RHW due to ACIs on the one hand (Fig. 9a), as shown in Fig. S7a. On the other hand, the presence of a large amount of small cloud droplets in the atmosphere is also beneficial for the collection of cloud water by ice-phase particles. Therefore, as can be seen from Fig. 9c, although ACIs decrease RHI with increasing AESF, the ice water path (IWP) shows an increasing trend with the AESF. Apparently, ARIs reduce the CWP and the decrease in CWP becomes significant with increasing AESF (Fig. 9a and b), which is mainly caused by the decrease of RHW due to ARIs (Fig. S7a). The impact of ARIs on IWP is insignificant with the AESF less than 1.0, in the range between −1 % and 1 %. When the AESF exceeding 1.0, ARIs increase the IWP due to ARI-induced increase of RHI in the atmosphere between 4 and 7 km (Fig. S7b).

Figure 10 (a) p-mean of R_effc, (b) variation of p-mean of R_effc due to ARIs, (c) p-mean of R_effi, and (d) variation of p-mean of R_effi due to ARIs over GZB + GZBs from 26 to 28 January 2022, as a function of the scale factor of anthropogenic emissions.

Figure 10a shows the variation of p-mean of cloud droplet effective radius (R_effc) with the AESF in F_BASE and F_ARI0. R_effc decreases with increasing AESF in the F_BASE and F_ARI0, which is well consistent with the variation of CND and CWP of cloud droplets. R_effc decreases with increasing AESF in both experiments. This decrease is primarily driven by the simultaneous decrease in CWP (Fig. 9a and b). In the F_ARI0, although CND increases monotonically with AESF (Fig. 8a), the reduction in CWP dominates, leading to a decreasing R_effc. In the F_BASE experiment, R_effc continues to decrease even when CND begins to decline after AESF exceeds 2.5 (Fig. 8a), reflecting the dominant role of CWP reduction. When the AESF exceeds 2.5, ARIs increase R_effc (Fig. 10b), which is caused by the ARI-induced decrease in the CND (Fig. 8b). When ARIs are not considered in the F_ARI0, p-mean of ice crystal effective radius (R_effi) increases with increasing AESF, which is caused by the ACI-induced decrease in CND (Fig. 8c) and increase in IWP of ice crystals (Fig. 9c). However, the increase in R_effi is not significant, i.e., when the AESF increases from 0.125 to 8.0, the R_effi increases by about 2.0 % (Fig. 10c). In addition, the impact of ARIs on the R_effi is insignificant, and the variation in the R_effi is less than 0.2 % (Fig. 10d).

3.5 Aerosol effects on precipitation in GZB and GZBs

Considering the difference in aerosol concentrations between GZB and GZBs, we discuss the impact of aerosols on precipitation in each region separately. Figure 11a provides the variation of the average 3 d accumulative precipitation in the GZB with the AESF. When ARIs are excluded (F_ARI0), the accumulative precipitation exhibits a weak increasing trend with the AESF, which is not consistent with previous studies over the GZB that have reported a suppression effect of ACIs on precipitation (e.g., Yang et al., 2013a, b; Bei et al., 2025). Yang et al. (2013a, b) have correlated the decreasing trend of precipitation inside the GZB with deterioration of the air pollution, revealing the microphysically suppressive effect of increased aerosols on precipitation. Bei et al. (2025) have also shown that ACIs tend to decrease convective precipitation in the GZB based on sensitivity simulations of a summertime short-term heavy rainfall event. In this study, as we have discussed above, increasing anthropogenic emissions decrease CWP and R_effc, lowering the conversion of cloud water to rainwater to decrease the liquid-phase precipitation. However, increasing emissions increase IWP and R_effi, facilitating the formation and growth of snow flakes to enhance the ice-phase precipitation. Figure 12 shows dependence of the liquid-phase and ice-phase precipitation at the ground surface in the GZB with the AESF. Liquid-phase precipitation shows a decreasing trend but ice-phase precipitation shows an increasing trend with the AESF. Overall, liquid-phase precipitation is lower at higher AESF, while ice-phase precipitation is higher. The increase in ice-phase precipitation partly compensates for the decrease in liquid-phase precipitation, leading to a net weak increase in total precipitation over the GZB for higher AESF (Fig. 11a).

Figure 11 Average accumulative precipitation in (a) GZB and (c) GZBs, and variation of precipitation due to ARIs in (b) GZB and (d) GZBs from 26 to 28 January 2022, as a function of the scale factor of anthropogenic emissions.

Figure 12 Average accumulative (a) liquid-phase and (c) ice-phase precipitation, and variation of (b) liquid-phase and (d) ice-phase precipitation due to ARIs in GZB from 26 to 28 January 2022, as a function of the scale factor of anthropogenic emissions.

ARIs decrease the precipitation with the AESF less than 0.6 and increases it when the AESF is more than 0.6 in the GZB (Fig. 11b). ARIs considerably decrease liquid-phase precipitation and the decrease is more than 50 % with the AESF greater than 1.0 (Fig. 12b). Except the decreasing CWP and R_effc with increasing the AESF, the decrease in RHW below 3 km caused by ARIs also contributes to the decreasing trend of liquid-phase precipitation in the GZB (Fig. S8a). However, ARIs generally increases the ice-phase precipitation monotonically with increasing AESF (Fig. 12c). ARI-induced increase in IWP can partially explain the increase in ice-phase precipitation with the AESF exceeding 0.6. Another possible reason is that ARI-induced increase in RHI facilitates survival ice-phase particles when they fall off from clouds (Fig. S8b). The decrease in liquid-phase precipitation exceeds the increase in ice-phase precipitation when the AESF is less than 0.6 and it is opposite with the AESF more than 0.6 (Fig. 12b and d).

Figure 13 Average accumulative (a) liquid-phase and (c) ice-phase precipitation, and variation of (b) liquid-phase and (d) ice-phase precipitation due to ARIs in GZBs from 26 to 28 January 2022, as a function of the scale factor of anthropogenic emissions.

In the GZBs, without ARIs, ACIs decrease precipitation with increasing AESF (Fig. 11c), and ARIs also appreciably decrease the precipitation with the same AESF (Fig. 11c and d), which are generally contrary to those in the GZB (Fig. 11a and b). There are two reasons that can explain the different responses of precipitation to ARIs and ACIs in the GZB and GZBs. In the F_ARI0 without ARIs, the variation of liquid-phase and ice-phase precipitation with increasing AESF in the GZBs is similar to those in the GZB. However, the ratio of liquid-phase precipitation to the total precipitation varies considerably in the GZB and GZBs. The proportion of liquid water path is consistently higher in GZBs than GZB across all emission scenarios (Fig. S9). For example, in the REF0, the ratio in the GZB is 2.5 %, but in the GZBs, it is 9.1 %. To further understand the contrasting precipitation responses between GZB and GZBs, the background thermodynamic and cloud microphysical conditions of the two regions are compared in the F_BASE (Fig. S10). Compared with GZB, GZBs exhibit higher lower-tropospheric relative humidity, substantially larger cloud water content, and smaller ice water content above 6 km. Therefore, aerosol-induced perturbations to liquid cloud microphysics are expected to exert a stronger influence on precipitation in GZBs, providing a physical basis for the contrasting precipitation responses shown in Fig. 11. Meanwhile, when ARIs are excluded, ACIs decrease the liquid-phase precipitation by 50 % when the AESF increases from 0.125 to 8.0 in the GZBs (Fig. 13). Although ACIs increase the ice-phase precipitation with increasing the AESF, the increase is less than the decrease in liquid-phase precipitation, causing the decreasing trend of precipitation with the AESF in the GZBs. In addition, ARIs decrease the liquid-phase and ice-phase precipitation in the GZBs with the same AESF (Fig. 13b and d), This reduction is attributed to the ARI-induced decrease of RHW and RHI in the atmosphere below 3 km (Fig. S11), which would enhance the evaporation of precipitation particles.

Figure 14 presents the 3 d accumulative precipitation distribution under different AESF (0.125, 1.0, and 8.0) in the F_ARI0 and F_BASE. Generally, increasing anthropogenic aerosols do not change the whole precipitation pattern significantly, i.e., the decrease in precipitation from north to south and the occurrence of heavy precipitation in the southeast of GZB + GZBs. However, the variation of different precipitation levels can be observed, for example, the increase in the area with precipitation between 5.0 and 10.0 mm in the GZB and the decrease in the area with precipitation exceeding 10.0 mm in the south. We also classify the simulated daily precipitations into slight (more than 0.0 and less than 0.1 mm d⁻¹), light (0.1–2.5 mm d⁻¹), moderate (2.5–5.0 mm d⁻¹), and heavy (more than 5.0 mm d⁻¹) snowfall according to the Chinese national standard GB/T 28592-2012 Grade of precipitation. In GZB + GZBs, increasing aerosols increase occurrence of slight and light precipitation but decrease moderate and heavy precipitation (Fig. S12).

Figure 14 Distribution of accumulative precipitation from 26 to 28 January 2022 for various scale factor of anthropogenic emissions when ARIs are (a, c, e) excluded and (b, d, f) included. The white contour outlines the GZB. The white contour outlines the GZB and the area outside of GZB in panels is referred to as GZBs.

4 Summary and conclusions

Using a cloud-resolving fully coupled WRF-Chem model, this study investigated a heavy snowfall event occurred in GZB + GZBs from 26 to 28 January 2022, with a focus on disentangling the roles of ACIs and ARIs under different aerosol loadings. The model reasonably reproduces the observed evolution of air pollutants and snowfall, with IOAs exceeding 0.60 for PM_2.5, O₃, and NO₂, and 0.85 for accumulated precipitation. A series of sensitivity simulations covering aerosol loadings from relatively clean to heavily polluted conditions are conducted to quantify the nonlinear responses of cloud microphysics and snowfall to anthropogenic aerosols.

Sensitivity results reveal that ARIs substantially modify the thermodynamic structure of the lower troposphere. Increasing aerosol loading reduces the mean PBLH by approximately 34 %, while producing cooling below 1.8 km and weak warming aloft, thereby enhancing lower-tropospheric stability and inducing secondary downward motion within the PBL and upward motion aloft. These dynamical adjustments redistribute atmospheric moisture, decreasing RHW below 3 km while enhancing RHI in the mixed-phase cloud layer under polluted conditions. Consequently, CND of cloud droplets exhibits a nonlinear response to aerosol loading, with a threshold occurring at an AESF of approximately 2.5 when ARIs are included. Meanwhile, the increase in IWP under polluted conditions is primarily associated with enhanced collection of cloud water by ice crystals.

When ARIs are excluded, increasing anthropogenic emissions raise CND of cloud droplets, but decrease LWP which is caused by ACI-induced decrease of RHW. While ACI-induced decrease in RHI reduces CND of ice crystals with increasing emissions, and the increasing trend of IWP with the AESF is mainly caused by enhanced collection of cloud water by ice crystals.

The response of precipitation to ARIs and ACIs is different in GZB and GZBs. In the GZB, excluding ARIs leads to a weak increase in 3 d accumulated precipitation with increasing emissions. Although increasing aerosols considerably decrease the liquid-phase precipitation, the ACI-induced increase in ice-phase precipitation outweighs the decrease in liquid-phase precipitation. ARIs increase the precipitation in the GZB when the AESF exceeds 0.6, due to ARI-induced increase in RHI which increase IWP and facilitate survival ice-phase particles in the atmosphere. In contrast, precipitation in the GZBs decreases with increasing aerosols, reflecting the dominance of liquid-phase precipitation that is strongly suppressed by ACIs and further reinforced by ARI-induced reductions in RHI. Increasing anthropogenic emissions do not substantially alter the precipitation pattern in GZB + GZBs, but can modify the frequency distribution of precipitation intensity.

Compared with previous studies over the GZB, which have primarily reported aerosol-induced suppression of precipitation through cloud microphysical effects (e.g., Yang et al., 2013a, b; Bei et al., 2025), the present study demonstrates that the response of winter snowfall cannot be explained solely by warm-cloud microphysics. Instead, the competition between liquid-phase suppression and ice-phase enhancement, together with ARI-induced thermodynamic adjustments, determines the net snowfall response. These results highlight that aerosol impacts on winter mixed-phase precipitation depend strongly on cloud phase partitioning and ambient thermodynamic conditions, providing new understanding of aerosol–cloud–radiation coupling under polluted winter environments. The identified regional contrast also suggests that aerosol effects on snowfall may vary considerably within the same synoptic system owing to differences in cloud microphysical characteristics and background humidity, emphasizing the importance of representing coupled aerosol–radiation–cloud processes in numerical weather prediction and regional climate models.

The sensitivity experiment design used to separate ARIs and ACIs is subject to uncertainties arising from the nonlinear coupling between radiation, thermodynamics, and cloud microphysics. The differencing approach does not fully isolate pure ARIs but rather includes contributions from nonlinear interactions between ARIs and ACIs. While this approach is widely used in the previous studies (Bei et al., 2025; Ding et al., 2019; Sun and Zhao, 2021), the quantitative ARIs effects reported here may contain some uncertainty, particularly in regions with strong aerosol-cloud-radiation feedbacks. Nevertheless, the robust regional contrast in precipitation responses that enhancement in the GZB vs. suppression in the GZBs is substantially larger than the estimated uncertainty due to nonlinearities, and the qualitative conclusions of this study remain robust. In addition, the present conclusions are based on a single snowfall event, and their broader applicability should be evaluated using additional winter precipitation cases. Future studies combining long-term observations, multi-case simulations, and more advanced radiative perturbation approaches will further improve the quantitative attribution of ARIs and ACIs and help assess the generality of the threshold behaviors identified in this study.