## Abstract

Entrainment is critical to the development of the atmospheric convective boundary layer (CBL), but little is known about how entrainment is impacted by the aerosol radiative effect. An aerosol radiation transfer model is used in conjunction with large-eddy simulation (LES) to quantify the impact of aerosol shortwave radiative heating on entrainment and thermodynamics of an idealized dry CBL under aerosol-loading conditions. An entrainment equation is derived within the framework of a zero-order model (ZOM) with the aerosol radiative heating effect included; the equation is then examined against the LES outputs for varying aerosol optical depths (AODs) and free-atmosphere stratification scenarios. The results show that the heat flux profiles become more nonlinear in shape as compared to the case of the clean (no aerosol pollution) CBL, with the degree of nonlinearity being highly dependent on the AOD of the layer for the given type of radiation-absorbing aerosols. As AOD increases, less solar radiation reaches the surface and thus the surface heat flux becomes smaller, and both actual (LES) and ZOM-derived entrainment flux ratios decrease. This trend is opposite to the clean CBL where the LES-predicted flux ratios show an increasing trend with diminishing surface heat flux, while the ZOM-calculated flux ratio remains constant. The modified dimensionless entrainment rate closely follows the −1 power law with a modified Richardson number. The study suggests that including the aerosol radiative effect may improve numerical air quality predictions for heavy-air-pollution events.

## 1. Introduction

An entrainment zone or a capping inversion is a transition layer between the turbulent convective boundary layer (CBL) and the stably stratified free atmosphere. Convective entrainment, a continual turbulent exchange of energy and substances between the mixed layer and the free atmosphere, is a crucial mechanism for the development of the CBL, as well as for distribution of heat, water vapor, CO_{2}, and air pollutants in the CBL (Betts and Ball 1994; Barr and Betts 1997; Huang et al. 2011). Accurate quantification of entrainment is important for parameterization of turbulent exchange processes in regional numerical weather and air quality prediction models (Hu et al. 2010a,b; Nielsen-Gammon et al. 2010). Difficulties in entrainment modeling stem from the complexity of physical processes that control interactions between the CBL and the free atmosphere, and the scarcity of observational data on entrainment (Moeng et al. 1999; Huang et al. 2011). Thermal and dynamic properties of the CBL entrainment can be modified by the presence of atmospheric aerosols due to their shortwave radiative absorption effect (Yu et al. 2002; Barbaro et al. 2013). In particular, uncertainties of entrainment quantification under conditions of strong aerosol pollution tend to be larger as compared to the entrainment predictions for the CBL without aerosols or with a low load of aerosols.

Parameterizations of entrainment are vital to regional or global numerical models due to the unresolved scales of motions within the entrainment zone. Such parameterizations are usually constructed in terms of the so-called integral parameters of entrainment (Fedorovich et al. 2004). The commonly used integral parameters include entrainment zone thickness; entrainment heat flux ratio, defined as the negative of the entrainment to surface sensible heat fluxes; and entrainment rate (or entrainment velocity), which is the rate of change in time of the CBL depth in the absence of large-scale subsidence; here is defined as the level of the heat flux minimum within the entrainment zone. These parameters are usually considered as functions of the free-atmosphere stratification, surface fluxes of heat and momentum, and scalars in the entrainment parameterizations (e.g., Lilly 1968; Betts 1974; Deardorff 1979; Fedorovich 1995; Conzemius and Fedorovich 2006, 2007; Pino et al. 2006; Kim et al. 2006; Garcia and Mellado 2014).

Proposed entrainment parameterizations are usually evaluated with the so-called CBL bulk models. Typically, the bulk models are classified into three categories in terms of the degree of complexity in representing the entrainment zone structure: zero-order model (ZOM), first-order model (FOM), and general structure model (GSM). They show a similarity in describing the CBL properties (height-constant temperature and velocity) but are different with respect to representation of entrainment zone. The ZOM assumes an infinitesimally thin capping inversion layer with discontinuities in temperature and velocity profiles (e.g., Lilly 1968; Fedorovich 1995). The FOM, on the other hand, considers a linear change of temperature in entrainment zone with a finite depth (e.g., Betts 1974; Sullivan et al. 1998). The GSM uses a self-similar representation of the temperature profile within the entrainment layer (e.g., Deardorff 1979; Fedorovich and Mironov 1995). Although FOMs and GSMs in some respects are more realistic than ZOMs in representing the entrainment-zone structure, their added complexity has limited their applications in atmospheric modeling. In comparison, ZOMs are relatively simple but show a competitive performance in prediction of the basic entrainment relationships (Fedorovich 1995; Fedorovich et al. 2004).

Entrainment integral parameter relationships obtained within the conceptual framework of bulk models for the clear and shear-free CBL have been evaluated in Fedorovich et al. (2004) using data from laboratory experiments and numerical large-eddy simulations (LESs). Typically the integral parameters of entrainment are evaluated within the so-called equilibrium entrainment regime, when the entrainment heat (buoyancy) flux remains in a quasi-steady balance with the surface heat (buoyancy) flux and the dissipation rate of turbulence kinetic energy (TKE) integrated throughout the CBL (Liu et al. 2018). The laboratory experiments showed that, in the shear-free CBL, the entrainment rate normalized with the convective velocity scale is inversely proportional to the Richardson number , where *g*, , , and are, respectively, gravity acceleration, constant reference temperature value, zero-order potential temperature increment across the entrainment layer, and surface heat flux (Deardorff et al. 1980). The FOM provides a similar –1 entrainment law in terms of dependence on Ri in the range of 13.6 ≤ Ri ≤ 43.8 (Sullivan et al. 1998). Fedorovich et al. (2004) evaluated the ZOM-derived entrainment-rate dependencies with the LES data, and found that exponents in the power-law dependencies that involve Ri may differ for different Ri ranges: for instance, the normalized entrainment-rate dependence on Ri in the shear-free, linearly stratified atmosphere follows a −1 power law at low Richardson number values (Ri < 10) and a −3/2 power law at Ri > 10. It was found by Moeng et al. (1999), however, that the normalized entrainment rate does not follow the −1 power-law dependence on the Richardson number for the cloud-top CBL. In this case, the normalized entrainment velocity depends not only on the interfacial Richardson number, but also on the radiative flux divergence above the minimum buoyancy flux level. It is not known to what extent the discussed relationships between the entrainment rate and Richardson numbers are applicable in the aerosol-polluted CBL.

The CBL entrainment flux ratio is variable and dependent on several factors such as surface heat flux, wind shear, and stratification in the free atmosphere. In the quasi-equilibrium entrainment regime, the ZOM-retrieved entrainment flux ratio shows a universal value of about 0.2 for the varying stratification conditions in the free atmosphere (Fedorovich et al. 2004; Liu et al. 2018). This value is generally supported by the observations under the weak mechanical turbulence conditions (Stull 1976) and data from laboratory tank experiments (Deardorff et al. 1980). On the other hand, the LES-predicted entrainment flux ratios for the shear-free CBLs (with entrainment flux evaluated at the level) show a pronounced dependence on the stratification strength in the free atmosphere (Sorbjan 1996; Fedorovich et al. 2004) and on the surface heat flux. For given stratification strength in the free atmosphere, the actual flux ratios decrease with increasing surface heat flux under the free convection conditions in the clear CBL (Garcia and Mellado 2014; Liu et al. 2018). Because of the aerosol radiative effect, solar radiation reaching the surface reduces and thus surface heat flux decreases, while aerosol-induced heating redistributes the heat energy inside the CBL and modifies the vertical profiles of potential temperature and heat flux. Thus the relation of entrainment flux ratio with surface flux can be altered as aerosol concentrations increase substantially. However, it is not currently established how the entrainment flux ratios change with the surface heat flux when the aerosols are present in the CBLs.

Aerosols exert significant impact on the CBL flow structures and entrainment. They reduce surface heat flux, modify vertical profiles of potential temperature, and change stability of the CBL due to aerosol absorption of the shortwave (SW) atmospheric radiation (Yu et al. 2002; Kedia et al. 2010; Barbaro et al. 2013; Ding et al. 2016). The impact depends on the vertical distribution of aerosols (Raga et al. 2001; Barbaro et al. 2013). Following the methodology of Betts (1974) and Van Zanten et al. (1999), Barbaro et al. (2013) derived an equation for the normalized entrainment rate within the FOM framework for the aerosol-loaded CBL and found that follows a −*n* power law of the Richardson number, with the value of *n* coming from studies of Fernando (1991) and Pino and Vilà-Guerau de Arellano (2008). They pointed out that the ratios of depend not only on surface heat flux and entrainment characteristics, but also on aerosol absorption of SW radiation. However, their study was limited to the cases with small aerosol optical depth (AOD; i.e., 0.2), and the FOM-derived entrainment rate has not been evaluated with the LES data.

In the present study, an entrainment-rate equation for the CBL under aerosol-loaded conditions is derived within framework of the ZOM. A radiation transfer model is coupled with the National Center for Atmospheric Research (NCAR) LES (e.g., Moeng 1984; Sullivan et al. 1996; Huang et al. 2008, 2009, 2011) to validate the ZOM entrainment-rate equation. The specific objectives of this study are the following: 1) to investigate the response of vertical profiles of potential temperature and heat flux to the aerosol radiative heating, 2) to determine factors that control the entrainment flux ratio in presence of the radiative heating, and 3) to derive and validate, by means of LES, a ZOM entrainment-rate equation for the aerosol-polluted shear-free CBL.

## 2. Methodology

### a. Coupling LES with radiative transfer model

The LES code used in this study was originally developed by Moeng (1984) and later refined by Sullivan et al. (1996), Patton et al. (2005), and Huang et al. (2008, 2009, 2011). The code is designed for the clear CBLs and does not include the impact of the aerosol radiative heating. In this study, an aerosol radiative heating term has been added in the filtered LES heat balance equation written in terms of potential temperature:

where denotes resolved-scale (in the LES sense) potential temperature; , , and represent resolved-scale *x*, *y*, and *z* components of wind, respectively; , , and are the three components of subgrid-scale potential temperature turbulence fluxes in the *x*, *y*, and *z* directions, respectively; and (K m s^{−1}) is the kinematic radiative flux related to the radiative energy flux *F* (W m^{−2}) through the air density (assumed constant) and specific heat at constant pressure .

The Santa Barbara Discrete Ordinates Radiative Transfer (DISORT) Atmospheric Radiative Transfer (SBDART) model (Ricchiazzi et al. 1998) is utilized in our study to calculate the radiative flux and heating rate at different heights. The SBDART model has been widely used for the calculation of aerosol radiative forcing (Kim et al. 2004; Tripathi et al. 2007; Gao et al. 2008; Kedia et al. 2010; Liu et al. 2012). The input parameters of SBDART include longitude, latitude, date, time, wavelength, height above the ground, surface albedo, and aerosol optical parameters such as AOD, single scattering albedo (SSA), and asymmetric factor *g*_{f}. The definitions and specific values of these three parameters are given in section 2c. The surface albedo in our calculations is set to 0.2. It is assumed that all the aerosols are confined to the CBL interior. Given the constancy of AOD and other aerosol optical parameters (e.g., SSA and *g*_{f}) adopted in the SBDART model during the whole simulation time period, the radiative flux *R* (or *F*) in each atmospheric layer may be calculated based on the radiative transfer equation (see Liou 2002). The corresponding heating rate due to aerosol absorption is given by .

In the coupled model system, the LES provides the CBL depth to the SBDART, the SBDART-calculated surface downward SW radiation is sent to the land surface module (LSM), and the LSM-predicted surface fluxes of sensible and latent heat are then used in the LES to drive the development of the CBL (Lee et al. 2012; Fig. 1). Meanwhile, the SBDART-predicted shortwave radiative heating rate is horizontally uniformly prescribed at each model level. Note that in this study the aerosol longwave (LW) radiative cooling is not included because it is assumed that LW cooling and the gaseous SW radiative heating nearly offset each other during the daytime CBL development (Barbaro et al. 2014).

### b. ZOM entrainment-rate equation under conditions of radiative heating

Similar to the cloud-topped CBL case (e.g., Lilly 1968, 2002; Moeng et al. 1999), the CBL heat balance equation in the presence of aerosols and associated radiative heating, written in terms of horizontally averaged potential temperature , appears as

where *Q* is the kinematic turbulent heat flux obtained by the horizontal averaging and *R* is the net SW radiative flux considered in section 2a. In the case of quasi-steady CBL development, when the mean vertical gradient of *θ* remains unchanged with time (Moeng et al. 1999); that is,

the total heat flux (*Q* + *R*) is linear with height.

By integrating (2) with respect to *z* from 0 to *z*_{i}, the height of the minimum heat flux taken as CBL top (see Fig. 2), using the ZOM representation of the CBL structure, we obtain the following CBL integral heat budget equation:

where is the potential temperature vertically averaged over the CBL, is the potential temperature jump at the top of the CBL, is the surface heat flux, and denotes the difference of *R* values at the level and at the surface, so (Fig. 2). Details of the derivation are given in the appendix.

Integrating Eq. (2) from 0 to and using Eq. (4) provide expression for the total heat flux as function of height:

where is entrainment heat flux defined within the ZOM framework. Equation (5) indicates that the total heat flux is a linear function of height *z* within the CBL under the adopted ZOM assumptions.

Following Moeng et al. (1999) and Lilly (2002), we use the closure relationship proposed by Deardorff (1976) for the entrainment flux in the radiatively heated CBL to express the entrainment heat flux as

where *A*_{h} = 0.5. As follows from Eq. (6), this value of corresponds to the established ZOM estimate of 0.2 for in the CBL without radiative heating, where is linear. The validity of Eq. (6) will be evaluated using the LES outputs for the radiatively heated CBL in section 3d.

Equation (7) indicates that in the presence of aerosol the ZOM entrainment flux depends not only on the surface heat flux, as in the ZOM of aerosol-free CBL, but also on the SW radiative flux resulting from absorption of radiation by aerosols within the boundary layer.

Let us denote . Equation (7) then becomes

A new convective velocity scale that includes effect of the radiative heat flux may then be introduced as

Accordingly, Eq. (8) then takes the form

Furthermore, the modified Richardson number may be introduced using the newly defined convective velocity:

This allows us to rewrite Eq. (8) in the form of the following entrainment relationship:

where is the dimensionless CBL entrainment rate in the presence of aerosol radiative heating. Thus, the modified normalized entrainment rate continues to follow the same inverse relationship with the modified Richardson number as its counterpart in the ZOM of CBL without aerosols (Fedorovich et al. 2004), where Eqs. (9), (11), and (12) are reduced, respectively, to

Here, , while is equal to 0.2 in the regime of equilibrium entrainment.

### c. Simulated flow configurations

We used the LES domain with size of 5 km × 5 km × 1.92 km and grid spacing of 50 m, 50 m, and 20 m in the *x*, *y*, and *z* directions, respectively. A larger domain of 10 km × 10 km × 1.92 km with finer grid mesh of 25 m × 25 m × 10 m was also used to test the effects of grid spacing on the findings. The initial potential temperature profile included a height-constant section with 290 K beneath 640 m and a section with potential temperature growing with height above 640 m at three prescribed rates: *dθ*/*dz* = 3, 6, and 9 K km^{−1}. A very dry and cloud-free boundary layer is assumed in all the simulations in order to remove the impact of water vapor and cloud effects. Geostrophic wind is set to zero (i.e., shear free) in all the numerical experiments due to weak winds observed on heavy aerosol pollution events (usually less than 5.0 m s^{−1}), which has negligible impact on CBL development and entrainment (e.g., Pino and Vilà-Guerau de Arellano 2008; Liu et al. 2018).

The radiation-absorbing aerosol is assumed to be distributed uniformly within the CBL, which has been confirmed by many observational and modeling studies (e.g., Steyn et al. 1999; Barbaro et al. 2013). AOD is a commonly used dimensionless parameter representing the degree to which the aerosols impede the transmission of light through the atmosphere. A total of 18 simulations were conducted through combining six AOD values (i.e., AOD = 0, 0.3, 0.6, 0.9, 1.2, 1.5) and three potential temperature gradients (i.e., 3, 6, and 9 K km^{−1}) in the free atmosphere. Information about simulations is summarized in Table 1, with the individual simulation parameters encoded in the case label. For example, the CTLN3 represents the control case without aerosols (AOD = 0) and with a free-atmosphere potential temperature vertical gradient of 3 K km^{−1}, while A03N6 is the case with AOD = 0.3 and a potential temperature gradient of 6 K km^{−1}. Two other optical parameters of aerosols (SSA and *g*_{f}) are taken from the observations conducted in the Yangtze River Delta region, China (Liu et al. 2012). Here SSA is the ratio of scattering to the total extinction of radiation and *g*_{f} represents the relative strength of forward scattering (Barbaro et al. 2013). Values of SSA = 0.9 and *g*_{f} = 0.6 were used in all the simulations. Latitude and longitude are set to 32.21°N and 118.70°E, respectively. The surface heat flux in LES is calculated by a land surface module with surface energy balance equation described in Huang et al. (2009). The SW radiative flux is calculated by the SBDART model (Ricchiazzi et al. 1998).

Horizontal periodic boundary conditions were used in the simulations. All the simulations represented the CBLs at the noontime on 24 January 2015, and for each case the surface heat flux was kept practically unchanged during the simulation. The 200-time-step samples of LES data were postprocessed after the turbulence reached a quasi-steady state (spinup time is about 1 h with slight difference among the cases presented in this study) as described in Patton et al. (2005) and Huang et al. (2008). The time step of integration was dynamically determined as explained in Huang et al. (2008). In this study it typically varied between 1.6 and 2.0 s.

## 3. Results and discussion

### a. Turbulence statistics

Table 2 presents hourly-averaged values of basic parameters for all simulated CBL cases. As presented data indicate, the surface heat flux diminishes greatly with increasing AOD. In fact, it drops by 63% as AOD increases from 0 to 1.5. Over the same AOD range, the boundary layer depth , diagnosed from the height of the minimum heat flux, shows a smaller reduction, of 2%–4% only. This reduction is smaller than the finding of Barbaro et al. (2013), where the CBL depth showed reduction by 8%–17% for different vertical distribution of aerosols. In other words, the contribution of the surface heat flux to the development of the CBL is largely compensated by the SW radiation absorption by the aerosol as less radiative energy reaches the ground. Meanwhile, the convective velocity , a measure of convection intensity in the CBL, shows a moderate decrease, by 28%–30%. The Richardson number Ri, a quantity representing the stability of entrainment zone, is less than 10 for all the cases and shows a decreasing trend with increasing AOD, which means that the entrainment zone becomes less stable with increasing AOD. The change in Ri results from the reduction in the potential temperature jump Δ*θ* across the entrainment zone caused by aerosol radiative heating of the air beneath the inversion. The entrainment zone depth Δ*z*, defined as the difference between the height at which the heat flux first becomes zero below *z*_{i} and the height at which the heat flux disappears above *z*_{i}, is increased by 16%–28% as AOD increases from 0 to 1.5. A similar result is also shown in Barbaro et al. (2013), where the entrainment zone was deepened by 18% and 91% for different vertical distributions of aerosols. Moreover, the entrainment-zone depth is reduced as free-atmosphere stratification becomes stronger, which is similar to the findings in the clean CBL (e.g., Fedorovich et al. 2004).

The entrainment rate shows an opposite trend with changing AOD under different stratification conditions in the free atmosphere. Specifically, displays a large increasing trend with increasing AOD under conditions of weak stratification (3 K km^{−1}), a slight decreasing trend under conditions of strong stratification (9 K km^{−1}), and little change under with moderate stratification (6 K km^{−1}). The change is dependent on the relative contributions of surface heat flux and aerosol SW radiation absorption to the development of the CBL. Aerosol radiative heating weakens the strength of capping inversion, and so its relative effect on changing the CBL depth is more significant under weak stratification conditions. As a consequence, the offset effect of aerosol on the reduction of the CBL thickness due to the decreased surface heat flux is stronger under weak stratification conditions than the offset effect under strong stratification conditions.

### b. Vertical profiles of potential temperature and heat flux

The impact of aerosol radiative effects on the CBL thermal regime is readily illustrated by changes in the potential temperature field. Figure 3 shows the profiles of potential temperature averaged over 200 time steps (i.e., over about 6–7 min) and over the horizontal planes under different aerosol loading and stratification conditions. Figures 3a–c correspond to three different stratifications. In each panel, vertical profiles of obtained with different AOD values are compared. The thermal field in the CBL displays a clear three-layered structure: a near-surface region with comparatively large vertical gradient, a mixed layer with a weakly varying (almost constant) , and an entrainment zone (capping inversion layer) with the temperature gradient larger than in the overlying free atmosphere. This thermal structure is in general terms similar to that observed in the clean atmospheric CBL (Kaimal et al. 1976; Stull 1988), indicating that the additional heat source associated with aerosol absorption of SW radiation does not essentially modify the three-layered structure of the potential temperature in the aerosol-loaded CBL, even though it has an impact on the shape of , which will be discussed in the next section. Radiation absorption by aerosols is warming the entire CBL and is larger with a heavier aerosol load. The mean potential temperature within the mixed portion of the CBL increases by 0.28 K as AOD grows from 0.3 to 1.5 for the CBL case with the free-atmosphere stratification of 3 K km^{−1} (Fig. 3a) and by 0.34 K with the free-atmosphere stratification of 9 K km^{−1} (Fig. 3c). At the same time, stronger free-atmosphere stratification attenuates the growth of the boundary layer, resulting in a shallower CBL. The maximum CBL depth (around 1100 m) in the case with the weakest stratification (Fig. 3a) is much greater than that in the case with the strongest stratification (approximately 820 m; Fig. 3c). The reduction of the CBL depth is attributed to the capping effect of the thermal stratification in the free atmosphere. In a shallower boundary layer, the aerosols induce larger heating effect, which explains the large increase of potential temperature as AOD changes from 0.3 to 1.5 under conditions of very strong stratification (9 K km^{−1}; Fig. 3c).

Another effect of the aerosol-related radiative heating is the modification of the static stability in the CBL. As indicated in Fig. 3, the height at which the potential temperature profile becomes stable (i.e., when *dθ*/*dz* becomes positive; this level is marked in the figure by a short horizontal dash) shifts downward as the AOD increases. For example, this height is 420 m in the CTL case, and it changes to 240 m in the A15 case. This effect of growing stability with increase of the aerosol loading becomes more prominent as the free-atmosphere stratification strengthens (Fig. 3c). Barbaro et al. (2013) found a similar earlier stable feature when the aerosols located at the top of the CBL. These observed tendencies are consistent with results from the previous modeling studies showing that the aerosol radiative effect enhances the stability in the upper portion of the CBL (e.g., Gao et al. 2015; Ding et al. 2016; Petäjä et al. 2016; Qiu et al. 2017; Li et al. 2017).

Figure 4 presents the profiles of simulated kinematic heat flux together with profiles of the SBDART-calculated SW radiation flux for the CBL cases with the free-atmospheric stratification of 6 K km^{−1}. The observed profiles support our previous analyses, based on Eq. (3), predicting that the turbulent heat flux in the CBL with aerosol loading departs from a linear shape due to the aerosol SW radiation absorption. However, the total heat (thermal energy) flux (the sum of sensible heat flux and SW radiation flux affected by the aerosol-related absorption) remains linear, in agreement with well-mixed potential temperature field within the CBL (see Fig. S1 in the online supplemental material). Figure 4a demonstrates linearity of the sensible heat flux with height when the CBL does not contain aerosols (CTL case; black line). The flux-profile deviation from the linear behavior becomes more prominent with increase of presence of aerosols in the CBL (the increase of the AOD) for a given SSA. The nonlinearity degree appears to be less pronounced as SSA increases (see Fig. S2). For the outer stratification of 6 K km^{−1}, the surface heat flux drops by 63% (Table 2) as the AOD increases from 0 to 1.5. At the same time, there is an overall decrease of the sensible heat flux over the entire CBL. The entrainment heat flux shows an overall decreasing trend as AOD increases. This is apparently a combined result of the reduction of the energy supply from the surface due to the weaker surface heat flux, and less free-atmospheric air engulfed into the CBL associated with the (slightly) stable stratification in the upper part of the CBL.

### c. Entrainment flux ratio

Now we turn our attention to the entrainment flux ratio , the negative of the entrainment heat flux to surface heat flux ratio. Figure 5 presents the entrainment flux ratios predicted by the LES and evaluated from the ZOM for all simulated CBL cases. The ZOM entrainment ratio is calculated by relating the entrainment flux calculated with Eq. (7) to the surface heat flux . The LES entrainment ratio is obtained by relating the minimum heat flux determined from the LES flux profiles to the surface heat flux . Notable differences between values retrieved from LES and values obtained with ZOM are evident in this figure. The LES entrainment ratio shows a decreasing trend as AOD increases. This trend is insensitive to the grid resolution and domain size (see Fig. S4). In other words, the LES ratio is reduced as the surface heat flux diminishes. This tendency is opposite to what was found in the clean CBL, in which the LES entrainment ratio shows an increasing trend as surface heat flux decreases (Liu et al. 2018). As seen in Fig. 4a, both the entrainment and the surface heat fluxes are reduced as AOD increases. However, the reduction of the former due to the aerosol radiative effect is larger than that of the latter. For instance, as AOD increases from 0 to 1.5, the surface heat flux is reduced by 63%, whereas the entrainment flux is reduced by 82%. The change of the LES-predicted flux ratio can be traced to the LES-predicted heat flux profiles shown in Fig. 4a. As AOD increases, the heat flux profile curvature becomes more prominent, as the height at which the heat flux switches from positive to negative values decreases. Meanwhile, the entrainment zone progressively extends downward, and the minimum heat flux in the entrainment zone becomes smaller (see, e.g., the cyan line in the plot). Such features are characteristic of the aerosol-loaded CBL and are not observed in the clean CBL (black line).

For given stratification, the ZOM flux ratio is reduced as AOD increases. On the other hand, the ZOM flux ratio remains constant (about 0.2) when no aerosols are present (Fedorovich et al. 2004; Liu et al. 2018). It follows from Eq. (7) that it is the aerosol radiative heating term [i.e., ] that causes the reduction of the ZOM ratio in the presence of aerosol and associated heating.

The LES and ZOM flux ratios for the same AOD are quite different. This is apparently because the LES ratio is evaluated locally, whereas the ZOM ratio is obtained from the interpolated profiles, so it is more integral in nature. The LES-predicted value is highly dependent on the surface heat flux, stratification conditions, and aerosol radiative effect, whereas the ZOM ratio shows much less variability. The difference between LES and ZOM values of tends to be smaller as the stratification in the free atmosphere becomes stronger (Fig. 5c). A similar trend was noticed for the difference between LES and ZOM entrainment ratios in the clean CBL (Fedorovich et al. 2004; Liu et al. 2018).

For a given AOD, the ZOM flux ratio is not sensitive to the free-atmospheric temperature gradient. However, the LES flux ratio is positively correlated with the stratification in the free atmosphere, in a manner that is similar to the case of the clean CBL (Fedorovich et al. 2004; Liu et al. 2018). Close examination of the LES heat flux profiles under different stratification (Fig. S3) suggests that, as the stratification becomes stronger, the negative entrainment flux is constrained to a shallower region below , with very small entrainment flux above the level, while the actual (LES) heat flux profile becomes closer to the ZOM profile.

The entrainment flux ratio dependence on the surface heat flux, free-atmosphere stratification, and aerosol loading may have important implications for numerical weather prediction and air quality modeling. Typically, such ratios are assumed to be constant in these models. For example, the Yonsei University (YSU) planetary boundary layer parameterization scheme uses an entrainment flux ratio of 0.15 when applied in the Weather Research and Forecasting (WRF) Model (Hong et al. 2006; Hu et al. 2010a). Our study suggests that the entrainment flux ratio should be parameterized as a function of all abovementioned forcings (including aerosol optical depth) rather than kept constant.

### d. Evaluation of entrainment equations

In the derivation of the ZOM entrainment-rate equation [Eq. (6)], entrainment flux is assumed to be proportional to the heat flux averaged over the entire CBL depth, following parameterization of Deardorff (1976). The proportionality coefficient *A*_{h} is taken constant and equal 0.5 for clean CBL case, in which *A*_{h} = 0.5 corresponds to the commonly used ZOM entrainment flux ratio of 0.2. It is not known a priori whether this proportionality assumption is suitable for the aerosol-loaded CBL. Figure 6 shows the calculated *A*_{h} for various AOD values and stratification strengths. To determine *A*_{h}, we first calculated the entrainment heat flux from LES as according to the ZOM methodology. Next, we calculated the average heat flux throughout the CBL from the LES data. Finally, *A*_{h} was computed as the ratio of these two fluxes. As seen from Fig. 6, the value is close to 0.5 for most of the considered cases. A mean *A*_{h} value of 0.41 with standard deviation of 0.2 is found when AOD is larger than 1.2. The scatter mainly originates from the distortion of the heat flux profile, which becomes substantial when AOD is larger than 1.2. Overall, we found that Deardorff’s (1976) closure assumption works reasonably well in the entrainment-rate equation for the radiatively heated CBL.

The entrainment flux equation [Eq. (7)] for an aerosol-loaded CBL is derived within the framework of ZOM using the closure assumption of Deardorff (1976), and with the aerosol radiative effect included. We now check the validity of this equation against the LES output. The term on the left-hand side of the equation is obtained from the LES data postprocessed with the ZOM representation of the CBL structure, while the term on the right-hand side is computed using the surface sensible heat flux and calculated radiation flux profiles assuming *A*_{h} = 0.5. The comparison of the terms is presented in Fig. 7. The overall agreement between the terms indicates that the assumptions involved in the derivation of Eq. (7) make sense, even though the uncertainties of may cause a variability of calculation on the left-hand side of the equation. We may thus conclude that obtained entrainment equation provides a robust ZOM parameterization of entrainment heat flux for the aerosol-loaded CBL. Note that this finding is also insensitive to the domain size and grid spacing (see Fig. S5).

Equation (7) allows quantification of the relative contributions of aerosol radiative heating effect and surface heat flux to the entrainment flux. For instance, with free-atmosphere stratification of 6 K km^{−1}, as AOD increases from 0 to 1.5, the entrainment flux is reduced by 63% due to the reduction in the surface heat flux, and by an additional 22% due to the aerosol radiative heating. In other words, although the diminishing surface heat flux is a main factor causing the reduction of entrainment flux under such conditions, the aerosol radiative effect is quite strong. Putting it differently, the traditional ZOM parameterization would overestimate the entrainment flux in the considered CBL case because it neglects the aerosol radiative heating contribution.

Figure 8 compares the ZOM-calculated entrainment flux with the flux predicted by LES at level . As pointed out by Fedorovich et al. (2004) and Liu et al. (2018), the LES prediction, representing a local entrainment flux, gets closer to the ZOM-calculated flux with stronger free atmospheric stratification or weaker surface heat flux. Similar trends are found in Fig. 8 for CBL cases with aerosol heating. The entrainment flux obtained using ZOM overestimates actual flux value under weak free atmospheric stratification (Fig. 8a), but shows better agreement with the LES-predicted flux as AOD increases.

The modified Richardson number Ri_{R} in Fig. 9 shows a decreasing trend with increasing AOD; a similar decreasing trend is also seen in the traditional Ri (Table 2). This trend is primarily due to the decreasing temperature jump Δ*θ* (Table 2) induced by the aerosol absorption of SW radiation. The Richardson number is a measure of static stability within the entrainment zone. Thus, one may conclude that aerosol loading weakens this stability. Similar results were obtained by Barbaro et al. (2013). In response to the less stable entrainment zone, the CBL is expected to grow faster, leading to a deeper CBL. However, the mean CBL depth is slightly decreased as AOD increases (Table 2). A similar contradiction is also found in Barbaro et al. (2013). Such attenuation effect of aerosol-related heating on the CBL growth may be associated with the fact that aerosol heating, besides reducing the temperature jump, stabilizes the flow, in terms of potential temperature gradient, within the upper portion of the CBL, beneath the entrainment zone (Fig. 3). On the other hand, the modified Richardson number Ri_{R} shows an increasing trend with the free-atmospheric stratification since is closely related to the potential temperature gradient in the free atmosphere. Moreover, it may be noticed from Fig. 9 that Ri_{R} has values less than 10 for the conditions investigated in this study.

According to the study of Fedorovich et al. (2004), the normalized entrainment rate *E* in the CBL without aerosol loading follows the −1 power law as function of Ri at low Richardson numbers (Ri < 10). Figure 10 shows the relationship of the modified normalized entrainment rate with the modified Richardson number Ri_{R}. Dependence of *E* on Ri is also included for comparison. It is obvious that *E* does not follow a −1 power law with Ri in the aerosol-polluted CBL. The larger the AOD, the larger is the deviation of *E* from the power-law line (Fig. 10b). However, the modified entrainment rate *E*_{R} follows closely the −1 power law with Ri_{R} (Fig. 10a). In other words, the traditional relationship between the dimensionless entrainment rate and the Richardson number based on the capping inversion temperature jump is no longer appropriate for heavy-aerosol-polluted conditions.

## 4. Summary and conclusions

In this study, an entrainment equation has been derived within the framework of a zero-order model (ZOM) to account for the radiative effect of aerosols in the atmospheric convective boundary layer (CBL). A radiation transfer model (SBDART) has been coupled with LES to account for the impact of aerosol heating on thermodynamics of the CBL and entrainment. The new ZOM entrainment equation has been evaluated with the LES results for the dry CBL with different aerosol loadings and free-atmospheric stratifications.

In the studied CBL, the solar radiation reaching the surface is reduced due to the scattering and absorption effects of aerosols, whereas the CBL air is warmed by the aerosol radiative heating. The resulting redistribution of heat energy inside the CBL modifies the vertical profiles of potential temperature and heat flux. Because of aerosol radiative heating, a stably stratified flow region is formed in the upper portion of the CBL beneath the capping inversion. The profiles of heat flux in the aerosol-loaded CBL become more nonlinear in shape as compared to the case of the clean CBL, and the degree of nonlinearity depends on the AOD of the layer.

Both actual (LES) and ZOM-derived entrainment ratios show a decreasing trend with increasing AOD. The study hints that the entrainment flux ratio in atmospheric models needs to be parameterized as a function of aerosol concentration or aerosol optical depth in order to improve weather and air quality predictions associated with the heavy-air-pollution events.

The LES results show that the closure assumption for the entrainment heat flux proposed by Deardorff (1976) is appropriate for a ZOM of the CBL with the aerosol radiative heating being accounted for. The proposed ZOM entrainment equation quantitatively agrees with the LES output processed in terms of the ZOM. This ZOM equation works especially well for CBL cases with strong free-atmospheric stratification. The traditional ZOM entrainment-rate relationship between the dimensionless entrainment and Richardson number has been found no longer valid for heavy aerosol pollution conditions. However, the modified dimensionless entrainment rate closely follows the −1 power law with modified Richardson number.

In the study, uniform aerosol optical properties have been assumed throughout the entire CBL. In real life, these properties may vary within multiscale three-dimensional turbulent flow characteristic of the atmospheric CBL. In the future studies, the turbulent transport of aerosol will be incorporated in the LES along with aerosol sources and sinks to better represent the interactions of aerosol radiative effects with dynamics and thermodynamics of the aerosol-polluted CBL. Moreover, the aerosol radiative heating in current study is incorporated in the conducted LES through the bulk (slab) representation of the aerosol-loaded layer, which is in conformity with zero-order representation of the entrainment zone (i.e., a discontinuity interface; see Fig. 2). Thus we derived an equation within the ZOM framework and demonstrated that it could be successfully applied to parameterize entrainment rate. Nevertheless, once a feasible and physically consistent method is developed, derivation of a first- (or higher-) order model of entrainment in aerosol-loaded CBL accounting for entrainment-zone depth and local (as opposite to bulk) heating in LES is required to further provide new insights into entrainment affected by aerosol radiative heating. The impact of single scattering properties and vertical distribution of aerosols on the CBL dynamics and entrainment will be further investigated through numerical experiments.

## Acknowledgments

The research was supported jointly by the National Natural Science Foundation of China (Grant 41575009), the National Key Research and Development Program of China (Grant 2017YFC0210102), and the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD). The first author also acknowledges the support Postgraduate Research and Practice Innovation Program of Jiangsu Province (Grant KYLX16_0945). The computing for this project was performed at the OU Supercomputing Center for Education and Research (OSCER) at the University of Oklahoma (OU). We thank three anonymous reviewers for their helpful comments, which greatly improved the manuscript.

### APPENDIX

#### Derivation of Eq. (4) from Eq. (2)

Integrating Eq. (2) from *z =* 0 to the upper edge of the zero-order discontinuity interface provides

where is the surface heat flux, and denotes the difference of *R* between the level of and the surface.

According to the Leibniz rule, we have

Using the ZOM representation of the potential temperature profile in CBL, we obtain

where is the height-constant potential temperature value within in the CBL. Thus,

where . Here is the potential temperature gradient in the free atmosphere, and is the surface temperature value obtained by the linear extrapolation from the level *z*_{i} to *z* = 0. Substituting Eq. (A5) into to Eq. (A1), we come to the ZOM equation [Eq. (4)] of the integral budget of heat in the CBL with radiative heating:

## REFERENCES

*An Introduction to Atmospheric Radiation*. Academic Press, 583 pp.

*An Introduction to Boundary Layer Meteorology*. Kluwer Academic Publishers, 666 pp.

## Footnotes

Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/JAS-D-18-0107.s1.

^{e}Additional affiliation: Jiangxi Province Key Laboratory of the Causes and Control of Atmospheric Pollution/School of Water Resources and Environmental Engineering, East China University of Technology, Nanchang, China.

© 2019 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).