Introduction

Aerosols exert critical influences on the climate system through direct and indirect radiative forcing. Although great progress has been made during the past decades in understanding the important effects of aerosols on climate, they are still among the largest uncertainties in understanding and projecting climate changes1,2,3,4,5. It is challenging to accurately quantify the aerosol impacts on the climate system because of insufficient understanding of the integrated aerosol radiative effects and drastic variations of aerosols in both temporal and spatial scales6,7,8,9,10.

Aerosols interact with the planetary boundary layer (PBL) through aerosol radiative forcing4,11,12. Aerosol-PBL interactions (API) are reported as one of the key mechanisms in the formation of severe air pollution5,13,14. In particular, aerosols are known to inhibit the development of the PBL by regulating surface sensible heat (SH)15,16,17,18. A shallow PBL would, in turn, favor the accumulation of aerosols, leading to the well-recognized positive feedback in the API11,13,19,20. However, the PBL entrainment, a process representing the turbulent exchange of air masses and heat fluxes between the PBL and the free atmosphere21,22, has not been well understood in the context of API. Because turbulent fluxes within the PBL cannot be solved mathematically22, the entrainment problem and the possible feedback with aerosol, as noted in this study, reflect the fundamental complexity of a multi-scale chaotic system.

As one of the cornerstones in parameterizing the PBL23,24, a linear scheme has been used to characterize the vertical heat fluxes and the relationship between surface SH fluxes (\(F_{H_S}\)) and entrainment fluxes (\(F_{H_{zi}}\)) in popular models such as the Weather Research and Forecasting25 and the European Centre for Medium-Range Weather Forecasts (ECMWF24,26). The following widely used parameterization was initially proposed by Ball27, followed by Betts28 and Deardorff29:

$$A_i = - F_{H_{zi}}/F_{H_S},$$
(1)

where Ai is called the entrainment parameter. In this way, the entrainment fluxes can be derived from the surface fluxes21,25. While the entrainment parameter is set as an empirical constant in the previously mentioned studies, we argue that Ai can vary considerably with aerosol properties and vertical distributions. Aerosols may come into play in two ways. First, aerosol can reduce the surface sensible heat by diminishing sunlight, thus weakening entrainment and shallowing the PBL. Second, aerosols can modify the thermodynamic structure in the PBL and cause a deviation from the linear structure of heat fluxes, modifying the entrainment. The former effect has been widely recognized13,16,20, whereas the latter remains largely uncertain, motivating this study.

Using comprehensive field observations made in Beijing, China, we examined the response of the entrainment process to aerosols. In particular, we obtained the entrainment rate from the lidar-derived PBL growth. In harmony with variations in aerosol radiative forcing, drastic changes in entrainment rates are noted for different aerosol loadings. Besides its loading, aerosol absorptions and vertical distributions can exert strong and nonlinear impacts on the vertical distribution of buoyancy, affecting the entrainment process. We propose a mechanism of aerosol-entrainment coupling to demonstrate the notable interactions between aerosols and the entrainment process in the vertical. As a critical part of API, such an aerosol-entrainment coupling leads to a much greater sensitivity of observed entrainment rates to aerosols than suggested by traditional calculations.

Results

Entrainment rates for variable pollution levels

Severe air pollution is generally associated with a shallow PBL19,30. API has been considered as a key factor in enhancing the surface pollution level19,20. We present the temporal evolution of the PBL height (PBLH) during the growth period for a clean case (21 August 2019) and a polluted case (1 April 2018) (Fig. 1). The two cases have dramatic differences in PBLHs. The PBL deepens more rapidly when the atmosphere is clean in terms of the mass of fine particles with an aerodynamic diameter less than 2.5 μm (PM2.5 = 12.1 μg m−3 on average) than when the atmosphere is polluted (PM2.5 = 261.8 μg m−3 on average). Even though SH is considerably higher in the polluted case, the PBLH growth rate is much lower, suggesting that other factors play critical roles in suppressing the PBL growth. By altering heating fluxes in the vertical, aerosols may be a key factor dictating the seemingly contradictory relationship between SH and PBL development.

Fig. 1: PBL evolution during the growth period.
figure 1

Lidar backscatter profiles for a the clean case (date: 21 August 2019) and b the polluted case (date: 1 April 2018) during the growth period over Beijing, China. Backscatter is presented as a normalized signal on a log-scale, in arbitrary units. Black dots mark the positions of the PBL top. Red dots and red lines indicate the corresponding measurements of sensible heat and the mean value, respectively. Pink dash lines indicate variations in the PBLH obtained from ERA-5 reanalysis data. The mean values of PM2.5 during these periods are given in each panel. The mean values of observed entrainment rates (\(w_e\)) are presented at the top of each panel. Hereafter, we use millimeter per second as the unit for the entrainment rate (1 mm s−1 = 10−3 m s−1).

Full size image

ECMWF 5 Re-Analysis (ERA-5) data show large overestimations of the PBLH under polluted conditions (Fig. 1b and Supplementary Fig. 1). The systematic differences between PBLH observations and reanalysis data can indicate the impact of unresolved processes in the ECMWF model. Note that the ERA-5 PBLHs are consistent with observed values under clean conditions (Fig. 1a). Under heavily polluted conditions, ERA-5 PBLHs are higher than observed values by more than 60%. Since ERA-5 data only take account of the climatology mean of aerosols loading31,32, notable aerosol radiative effects associated with severe air pollution may contribute to the differences in PBL evolution between observations and ERA-5 data. It is the large gap between PBL observations and reanalysis data under heavily polluted conditions that originally motivated this study to discover any core processes causing it. Because PBL development is directly linked with entrainment, the biases in PBL growth suggest that ERA-5 data cannot accurately estimate entrainment under polluted conditions.

Based on comprehensive field observations (~2200-hour data), Fig. 2 shows the dependence of observed entrainment rates on SH (see Methods). In our analyses, we only discussed the entrainment process for the PBL growth period during the daytime. Under the nocturnal PBL, entrainment process is insignificant, except for cloudy conditions33. These statistical results demonstrate considerable differences in the relationships among the three distinctly different pollution levels. By performing the Student’s T-test, the relationships between entrainment rates and sensible heat are found to be statistically significant at the 99% confidence level under three different pollution levels. In general, entrainment rates rise with increasing SH. However, under polluted conditions, this positive relationship is relatively flat. In particular, an increase in SH leads to a limited enhancement in entrainment rates when SH is greater than 100 W m−2. Although SH is the primary force driving the entrainment process, the entrainment rate significantly differs for the same SH. For a fixed SH, the PBL entrainment rate is considerably lower under polluted conditions, indicating that other factors may notably suppress entrainment rates when severe air pollution is present. The role of aerosol radiative forcing has not been well demonstrated in modulating entrainment fluxes.

Fig. 2: Entrainment rates for different pollution levels.
figure 2

Observed distribution of entrainment rate as a function of sensible heat fluxes. Box-and-whisker plots show the percentile values of 10th, 25th, 50th, 75th, and 90th. Blue, green, and red dots represent the mean values of PBL entrainment rates under clean (percentile: 0–33%), moderate (percentile: 33–67%), and polluted conditions (percentile: 67–100%), respectively. Hereafter, low cloud cases are excluded to avoid heat fluxes induced by boundary-layer clouds.

Full size image

The role of aerosol vertical distribution

Interactions between aerosols and PBL thermodynamics are associated with aerosol radiative forcing, which is notably affected by the aerosol vertical distribution and absorption13,34. To investigate this, we differentiate the latter into four categories according to aerosol extinction coefficient and absorption: a decreasing structure with a peak near the ground for weakly absorbing and absorbing aerosols and an inverse structure with a peak in the middle or upper PBL for weakly absorbing and absorbing aerosols (see more details in Methods). Following the previous studies35, SSA is used to differentiate between absorbing aerosols and weakly absorbing aerosols. Considerable differences in entrainment rate are revealed in Fig. 3a for the two vertical structures and different aerosol absorption levels. Relative to the decreasing aerosol structure, the entrainment rate is lower by about 80% for the inverse aerosol structure with absorbing aerosols.

Fig. 3: The PBL entrainment associated with aerosol vertical structures.
figure 3

a The black dot and whisker represent the average value and standard deviation, respectively, of the observed entrainment rates for decreasing/inverse aerosol structures. The width of the color-shaded areas represents the distribution of entrainment rates. In the violin plot60, we applied a 3.5 mm s−1 smooth window for the density distribution. The pink dots indicate the mean values of entrainment rates calculated from ERA-5 reanalysis data. The green dashed line with the shaded area shows the mean and standard deviation of SH. b Aerosol-induced changes in heat fluxes in the upper PBL for decreasing/inverse aerosol structures. The green dots represent the corresponding mean values of \(- A_i\overline {(w^\prime \theta ^\prime )_s}\), which roughly indicate the expected entrainment fluxes. In (a, b), blue and red areas indicate weakly absorbing (single-scattering albedo, SSA > 0.9) and absorbing (SSA < 0.85) cases, respectively. c Height-dependent correlation coefficients between entrainment rates and aerosol extinction coefficients. Note that entrainment rates are fixed for the atmosphere column, and aerosol extinction coefficients are height-dependent.

Full size image

Despite the drastic variations in entrainment rate for different aerosol vertical structures and absorption levels, changes in SH are much smaller than changes in entrainment rate. The differences in SH between absorbing and weakly absorbing cases are less than 10%. On the other hand, for the inverse aerosol structure, the entrainment rate for absorbing cases is only one-third of the value for weakly absorbing cases, while SH is similar. Absorbing and scattering aerosols are both able to decrease the amount of shortwave radiation reaching the surface, suppressing SH. Additionally, absorbing aerosols exert important heating effects on the atmosphere, changing thermodynamic structures and modulating the entrainment process.

For the absorbing case with the inverse structure, the entrainment rate is low (16.6 mm s−1), implying very slow growth in the PBL during the daytime, presumably because of API. It is also echoed by the systematic differences between observations and reanalysis data, which may be attributed to the effects of aerosols in suppressing PBL development. For this case, the mean value of the entrainment rate derived from ERA-5 is larger than observed by 147%.

Regarding the inverse aerosol structure, our radiative transfer calculations show strong heating of the atmosphere by absorbing aerosols (see Methods and Supplementary Fig. 2). Furthermore, the different vertical structures and aerosol absorption levels also tend to generate a potential temperature inversion within the PBL (Supplementary Fig. 3). Based on the K-closure theory25,36, we theoretically estimate the impacts of aerosols on heat fluxes in the upper PBL (see Methods and Fig. 3b). Under the inverse aerosol vertical structure, abundant aerosols located at the upper PBL can significantly suppress vertical heat fluxes, leading to a reduction in entrainment fluxes at the PBL top (Fig. 3 and Supplementary Fig. 3).

Moreover, we found a linkage between humidity profiles and aerosol vertical structures (Supplementary Fig. 4). Humidity is notably higher under the inverse aerosol structure. Although the shortwave water vapor heating effects are higher under high humid conditions, the longwave water vapor cooling effects are also notably higher. The two contrasting effects lead to insignificant differences in water vapor radiative effects between the decreasing and inverse structures. Thus, radiative forcing due to water vapor may not be a key factor contributing to the changes in entrainment rates.

Furthermore, we also examined the general relationship between entrainment rates and aerosol extinction coefficients in the vertical (Fig. 3c). Note that the vertical distribution of aerosol extinction coefficients was derived from lidar (see more details in Methods). Following previous studies37, we normalized the height coordinates according to the position of the PBL top. Entrainment rates are notably associated with aerosol loading in the vertical. In particular, negative correlations (~0.45) between the entrainment rate and aerosol loading are found in the upper PBL, suggesting an important interaction between aerosols and the entrainment process.

Sensitivity of the entrainment rate to aerosol loading

As previously mentioned, there are considerable differences in PBL evolutions between observations and ERA-5, likely associated with aerosols. One reason may be that the ECMWF model does not assimilate the daily variations in aerosol measurements. Also, cloud-free entrainment fluxes in model simulations are conventionally diagnosed from surface fluxes24, which does not consider the impacts of aerosols on vertical heat fluxes.

To further illustrate this point, we compare the responses of entrainment rates to aerosol loading (i.e., PM2.5 concentrations) from different data sources (Fig. 4). To limit the impacts of SH variations, data were divided according to three SH scenarios: high SH (percentile: 66–100%), medium SH (percentile: 33–66%), and low SH (percentile: 0–33%).

Fig. 4: Sensitivity of entrainment to PM2.5 observations based on different data sources.
figure 4

Red, green, and blue lines show how entrainment rates change as PM2.5 changes under a high SH scenario (percentile: 66–100%), the medium SH scenario (percentile: 33–66%), and the low SH scenario (percentile: 0–33%), respectively. The entrainment rates are derived from a the linear scheme using radiosonde and SH observations as input, b ERA-5 reanalysis data, c PBLH observations (\(w_e = dz_i/dt - w_i\)), and d PBLH observations removing wind shear effects. Due to the availability of radiosonde data, we only compare results during noontime (12:00–14:00 local time) in the warm season (June to September). The grey shaded areas represent the averaged linear regression. The black dash lines indicate the regressions are statistically significant at the 99% confidence level.

Full size image

With observations of SH and radiosonde, we can calculate entrainment rates using the linear scheme (see Methods). Since radiosonde data were collected during midday, calculations covering the period 12:00–14:00 local time (LT) in the warm season (June to September) were done. By assuming the linear relationship between the entrainment fluxes and surface SH fluxes, the linear scheme calculates the entrainment rate as a function of SH. As a result, the entrainment rates from the linear scheme do not change much with PM2.5 for a given SH scenario (Fig. 4a). The ECMWF also uses a linear scheme in their PBL parameterization to estimate the entrainment rate24,26. Figure 4b shows that the entrainment rate from ERA-5 slightly decreases with PM2.5 for a given SH scenario.

Despite the weak responses of entrainment to aerosols in the linear scheme and reanalysis data, Fig. 4c indicates that entrainment rates derived from PBLH observations dramatically decrease as the aerosol loading increases. With wind shear from radiosonde data and theoretical calculations (see Methods), we remove the wind shear effects in the observed variations of entrainment rate (Fig. 4d). By enhancing the turbulent heat fluxes, the wind shear is the main contributor to the observed variations in entrainment rate in the clean condition (Supplementary Fig. 5a). On the other hand, wind shear has weak effects on entrainment rates for polluted conditions.

As the key factor in dictating the entrainment, the turbulent heat fluxes within the PBL are largely dictated by three processes: (1) convections triggered by surface fluxes, (2) turbulent fluxes produced by wind shear, and (3) the temperature gradient from radiation36,38. By constraining the first two processes, we still find a notable response of entrainment to the increase in aerosol loading (Fig. 4d). Hence, we proposed that such a clear decrease pattern in the entrainment rate is significantly contributed by aerosol radiative effects, which can exert significant, nonlinear impacts on the vertical heat fluxes within the PBL (Fig. 3).

The observed sensitivity of entrainment rates to PM2.5 is much greater than theoretical calculations, suggesting that the important impact of aerosols on entrainment is not considered in the linear scheme. Since the linear scheme is widely used in PBL parameterizations, this would lead to considerable underestimations in the responses of entrainment to aerosols in model simulations carried out for various purposes.

Discussion

Following quantitative analyses of observations, theoretical calculations, and reanalysis data, we developed a conceptual scheme to demonstrate the impacts of API on the entrainment process (Fig. 5). As air masses are entrained from the free atmosphere to PBL, the entrainment heat flux is negative, and its magnitude presents a notable decrease under polluted conditions. By blocking sunlight, aerosols cool the surface and suppress SH. The lessened SH, however, cannot fully explain the observed relationships between aerosol loading and the entrainment rate (Figs. 3 and 4). Absorbing aerosols can induce strong heating effects on the lower atmosphere, which helps stabilize the atmosphere and slow PBL growth, especially under polluted conditions (Fig. 1 and Supplementary Fig. 1). When the inverse structure induced by absorbing aerosols is present (Supplementary Fig. 2), turbulence and PBL development are suppressed. This inhibits the entrainment of clean air from the free atmosphere into the PBL, facilitating the accumulation of air pollution in the upper PBL. In turn, high aerosol loading in the upper PBL further favors the formation of an inversion in potential temperature profiles (Supplementary Fig. 3), suppressing heat fluxes and the entrainment process. A positive feedback loop thus develops, which can amplify the previously discussed aerosol effect, reinforcing the interaction between the entrainment process and aerosols (Fig. 3). We refer to the strong link between entrainment and aerosols as “aerosol-entrainment coupling” (the pink double arrow in Fig. 5).

Fig. 5: A schematic diagram describing the aerosol-entrainment coupling.
figure 5

The background grey arrow sketches the vertical transport of humidity, aerosols, and heat fluxes. Orange, curved arrows represent solar radiation. The blue dash-dotted line represents the position of the PBL top (zi). Black, curved arrows indicate sensible heat fluxes. Red, curved arrows indicate entrainment at the PBL top. \(F_{H_S}\) and\(F_{H_{zi}}\) represent surface SH fluxes and the entrainment heat fluxes, respectively. The blue arrows represent the suppression effects of aerosols on the surface sensible heat and vertical heat fluxes. The blue, shaded area indicates the perturbation in heat fluxes induced by aerosols. Red arrows represent the accumulation of aerosols under the weak entrainment condition. Due to these interactions, the entrainment and aerosols are coupled here (marked as the pink double arrow).

Full size image

Aerosol absorption and its vertical distribution play key roles in modulating these impacts. Aerosol-entrainment coupling can negate the widely used linear scheme between SH fluxes and entrainment fluxes, although the radiative effect of aerosols may be partially accounted for by SH, which is proportional to incoming solar radiation at the surface. Because the linear scheme does not include aerosol-induced heating in the vertical, aerosol-entrainment coupling may be poorly represented in associated model simulations (Figs. 3 and 4).

Moreover, PBL processes cannot be accurately estimated by numerical models, especially under severe air pollution conditions39,40. As a key part of API, aerosol-entrainment coupling is a contributing factor to the uncertainties. Without explicit consideration of aerosol-entrainment coupling, we may not accurately predict or estimate PBL development in a polluted environment. This issue calls for further explicit PBL parameterization to account for such a strong response of entrainment to aerosols.

Methods

Descriptions of datasets

Beijing is one of the most urbanized and densely populated metropolises over the world, and suffers from severe air pollution with rich absorbing aerosols30. In this study, we utilized comprehensive observational data collected at a superstation (39.80°N, 116.47°E) in Beijing, China. Supplementary Fig. 6 shows the location of this station and surrounding topography. Simultaneous observations made by a micro-pulse lidar (MPL), radiosonde, surface meteorological instruments, and the eddy covariance technique were available from July 2017 to October 2019. During this period, lidar measurements were available for 592 days, and surface flux data from the eddy covariance technique (Li-7500A)41 were available for 826 days. Radiosondes were launched at 13:15 LT during the summertime (June to September), besides the operational launches at 07:15 and 19:15 LT. Observed variables include meteorological measurements and profiles of temperature, humidity, wind, and pressure. The variable vertical resolution of the radiosonde depends on height and is typically smaller than 10 m42. The lidar records the backscatter signals every 30 s with a vertical resolution of 30 m. Caused by the correlation of incomplete laser pulses, the near-ground blind zone of MPL was ∼0.20 km. Following the standard procedures, we applied the background subtraction, after-pulse, saturation, range, and overlap corrections to the lidar backscatter to derive normalized signals34. To quantify aerosol optical properties during the study period, level 1.5 single-scattering albedos (SSAs) and aerosol optical depths (AOD) from the Aerosol Robotic Network (AERONET) Beijing_RADI station (40°N, 116.38°E) under cloud-free conditions were employed43. To avoid regional differences and to get the representative value of PM2.5, we acquired averaged PM2.5 concentrations from five air quality sites located within twenty kilometers of the superstation in Beijing (Supplementary Fig. 6). In addition to the observational data, we also use ERA-5 data to obtain the PBLH and vertical velocity. ERA-5 provides the hourly measurements at a 0.25° × 0.25° longitude-latitude grid. The vertical resolution of ERA-5 data is 25 hPa in the lower atmosphere (700–1000 hPa).

PBLH derived from lidar

Lidar techniques have been widely used to retrieve aerosol, PBL, and cloud properties13,40,44,45. PBLH is retrieved from the lidar data over Beijing during the daytime (0800–1900 LT). In the method section, we use zi to represent PBLH for convenience. To estimate the zi from lidar measurements over Beijing, we utilized a well-established method46, which has been developed and evaluated by the long-term measurements obtained from the Atmospheric Radiation Measurement Southern Great Plains site. For the range between 0.25 to 4 km, we initially identified the step signals in the wavelet transform function collocated with a gradient greater than a preset threshold. Furthermore, we calculated the shot noise (ε) contributed by dark currents and background light for the backscatter profile, and then selected the threshold as 3ε. The start retrieval of zi (at 0800 LT) was constrained by the PBL depth calculated from the corresponding radiosonde. For a specific time, we collocate all step signals in lidar backscatter profiles and then select the step signal as the current zi based on the temporal continuity. The maximum variation of zi within 10-min is set to 0.2 km. During the growth period (0800–1400 LT), we first attempt to select zi higher than the previous zi position. If there is no available step signal higher than the previous zi position, zi would be selected as the closest step signal below the previous zi position. We further identified the boundary-layer clouds to diagnose the zi under the cloudy condition. Lidar-retrieved zi agree well with those derived from radiosonde, with a correlation coefficient of 0.81 (Supplementary Fig. 7). This correlation is statistically significant at the 99% confidence level. For the validation purpose, we use the Liu and Liang47 method to retrieve zi from noontime radiosonde measurements.

Classification of aerosol vertical distributions

We use the Klett method to calculate aerosol extinction profiles from the lidar backscatter signals48. As an important parameter for retrieving extinction profiles, the column-averaged lidar ratio (i.e., the ratio of extinction to backscatter) is normalized here by AERONET-derived AODs at 0.5 µm. Within the blind zone, the aerosol extinction coefficient is considered to be equal. Due to the effects of multiple scattering, the lidar ratio, the overlap function, and the background noise, the overall uncertainties of aerosol extinction fall within the range of 20–30% during the retrieval process49. By using the cloud-free aerosol extinction profiles below the PBL top, we classified aerosol vertical structures into two categories: decreasing with altitude and increasing with altitude (inverse). In general, different aerosol vertical distributions lead to drastic disparities in aerosol radiative effects in the vertical (Supplementary Fig. 2).

Calculation of entrainment rates based on the linear scheme

The linear scheme is used in nearly all PBL parameterizations to characterize the entrainment rate24,25,26. In general, the entrainment rate is positively correlated with turbulence fluxes within the PBL and negatively correlated with the capping inversion at the PBL top. The generation rate of turbulent kinetic energy is considered to be linearly proportional to the surface SH. In the linear scheme, the entrainment rate can be expressed as

$$w_e = A_i\frac{{\overline {(w^\prime \theta^{\prime} )_s} }}{{(\Delta \theta )_{z_i}}},$$
(2)

where \(\overline {(w^\prime \theta ^\prime )_s}\) represents sensible heat fluxes at the surface, which can be obtained from the surface eddy covariance measurements. The capping inversion, \((\Delta \theta )_{z_i}\), represents the jump in potential temperature at the PBL top (zi), which can be directly calculated from radiosonde-based potential temperature profiles50. The entrainment parameter Ai is considered as an empirical constant in the linear scheme. The entrainment rate can thus be determined from SH and the capping inversion using Eq. (2). For a convective PBL, Ai ranges from 0.1 to 0.321,25,29. Following observational studies51,52, we set Ai as 0.3. In the linear scheme, the entrainment rate is largely controlled by surface SH. The linkage between the capping inversion and aerosol loading is weak (Supplementary Figs. 5b and 8). Thus, the variations in the capping inversion may not explain the sensitive responses of entrainment rates to aerosols. Based on the linear scheme, we can calculate entrainment rates from radiosonde and SH measurements.

Calculations of entrainment rates based on PBL growth

The PBL has a strong diurnal cycle over land. During the growth period, air in the free atmosphere is entrained into the PBL. The entrainment process directly dictates the development of the PBL. Without a term related to cumulus, the PBL growth rate (dzi/dt) can be expressed as21:

$$\frac{{dz_i}}{{dt}} = w_e + w_i,$$
(3)

where we is the entrainment rate, representing the volume of air entrained per unit horizontal area per unit time, and wi represents the vertical velocity of the large-scale motion field at the top of the PBL. Because large-scale motion is usually downward under fair-weather conditions, wi is typically negative and relatively small compared to the entrainment rate22. The vertical velocity of large-scale motion can be acquired from ERA-5 reanalysis data. The term on the left-hand side of Eq. (3) indicates the growth of Zi, which can be calculated from lidar retrievals. Due to the blind zone of the lidar, the growth term (\(dz_i/dt\)) may not be correctly calculated in the early morning. Also, the growth term is not applicable once the PBL has grown to its daily maximum height. Therefore, the growth term was calculated for the period 1000–1400 LT. Measurements were averaged over one-hour intervals. Low-cloud cases (34% of total cases) were excluded to avoid heat fluxes associated with cumulus clouds. In this way, we used lidar observations to derive the entrainment rate. We can also use Zi from ERA-5 to calculate the entrainment rate based on Eq. (3).

Calculations of aerosol-induced changes in heat fluxes

To quantify the aerosol radiative effects in the vertical scale, we applied the Santa Barbara DISORT Atmospheric Radiative Transfer (SBDART) model53 to estimate the impacts of aerosols on atmospheric heating rate under cloud-free conditions54. Comprehensive aerosol inputs include SSAs, AODs, and asymmetry factors at multiple wavelengths (i.e., 1.02, 0.87, 0.67, 0.5, and 0.44 µm) obtained from the Beijing_RADI station and vertical profiles of the aerosol extinction coefficients retrieved from the lidar. The surface reflectances from Moderate Resolution Imaging Spectroradiometer are also used as input parameters. With the SBDART model and comprehensive aerosol observations, we can calculate aerosol radiative forcing in the vertical under different conditions (Supplementary Fig. 2). As the human-emitted fine particles are dominated in the study region, the aerosol longwave radiative effects are insignificant (Supplementary Fig. 9). Similar to the previous studies55, we only consider the aerosol radiative effects at shortwave (0.2–4 \(\mu m\)).

Based on the K-closure theory36, heat fluxes can be estimated as follows:

$$\overline {w^{\prime} \theta ^{\prime} } = - K_h\frac{{\partial \theta }}{{\partial z}},$$
(4)

where Kh represents heat eddy diffusivity, which increases with the intensity of the turbulence and varies greatly within the PBL. The K-closure theory indicates that the heat fluxes are proportionate to the potential temperature gradient22,25. By taking into account the aerosol heating effect, the change rate of heat fluxes should be proportionate to the change rate of the temperature gradient. The aerosol-induced changes in heat fluxes can be expressed as

$$\frac{{\partial \overline {w^{\prime} \theta ^{\prime} } }}{{\partial t}} = - K_h\frac{{\partial ^2\theta }}{{\partial Z\partial t}},$$
(5)

where \(\frac{{\partial ^2\theta }}{{\partial Z\partial t}}\) represents aerosol-induced changes in the potential temperature gradient, which can be calculated using the SBDART model. Following several well-established approaches25,56, the heat eddy diffusivity within the PBL is expressed as

$$K_h = Pr^{ - 1}kw_sz\left( {1 - \frac{z}{{z_i}}} \right)^2 = Pr^{ - 1}k\left( {u_ \ast ^3 + \phi _mkw_ \ast ^3z/z_i} \right)^{\frac{1}{3}}z\left( {1 - \frac{z}{{z_i}}} \right)^2,$$
(6)

where k is the von Kármán constant (equal to 0.4), z represents the height from the ground level, ws is the mixed-layer velocity scale, u* is the friction velocity, \(w_ \ast = \left[ {\left( {\frac{g}{\theta }} \right)\overline {\left( {w^\prime \theta ^\prime } \right)_s} z_i} \right]^{1/3}\) is the convective velocity scale, near-surface virtual potential temperature, and Pr is the Prandtl number57. During the PBL growth period, the dimensionless term \(\phi _m\) is given by

$$\phi _m = \left(1 - 16\frac{{0.1z_i}}{L}\right)^{ - 1/4} = \left[1 + 16\frac{{0.1kz_ig\overline {\left( {w^{\prime} \theta ^{\prime} } \right)_s} }}{{u_ \ast ^3\theta }}\right]^{ - 1/4},$$
(7)

where L represents the Monin-Obukhov length scale. In Eqs. (6, 7), u* and \(\left( {w^\prime \theta ^\prime } \right)_s\) are obtained from the eddy covariance technique, while Zi is obtained from lidar data. Therefore, we can calculate vertical profiles of \(K_h\) and further estimate \(\frac{{\partial \overline {w^\prime \theta ^\prime } }}{{\partial t}}\) due to aerosol radiative forcing in the vertical. Since the entrainment process occurs at the PBL top, we averaged the aerosol-induced changes in heat fluxes in the upper PBL (i.e., from 0.5zi to zi).

Calculations of impacts of wind shear on the entrainment rate

By using the wind shear derived from the radiosonde, we quantify the impacts of wind shear on the entrainment rate. Based on the previous studies58,59, the effects of wind shear on the entrainment rate at the PBL top can be expressed as follows:

$$w_e = w_{e0}/\left( {1 - \frac{{\left| {\Delta U} \right|^2}}{{R_iw_ \ast ^2}}} \right),$$
(8)

where \(\Delta U\) is the wind shear at the PBL top, and \(R_i\) is the Richardson number. \(w_{e0}\) and \(w_e\) represent the entrainment rate without wind shear and the entrainment rate with wind shear, respectively. In particular, \(\Delta U\) and \(R_i\) can be obtained from the radiosonde data58. Based on Eq. (8), we can remove the wind shear effects in the variations of entrainment rate (Fig. 4d).