## Abstract

The indirect radiative effect of aerosol variability on shallow cumulus clouds is realized in nature with considerable concurrent meteorological variability. Large-eddy simulations constrained by observations at a continental site in Oklahoma are performed to represent the variability of different meteorological states on days with different aerosol conditions. The total radiative effect of this natural covariation between aerosol and other meteorological drivers of total cloud amount and albedo is quantified. The changes to these bulk quantities are used to understand the response of the cloud radiative effect to aerosol–cloud interactions (ACI) in the context of concurrent processes, as opposed to attempting to untangle the effect of individual processes on a case-by-case basis. Mutual information (MI) analysis suggests that meteorological variability masks the strength of the relationship between cloud drop number concentration and the cloud radiative effect. This is shown to be mostly due to variation in solar zenith angle and cloud field horizontal heterogeneity masking the relationship between cloud drop number and cloud albedo. By combining MI and more traditional differential analyses, a framework to identify important modes of covariation between aerosol, clouds, and meteorological conditions is developed. This shows that accounting for solar zenith angle variation and implementing an albedo bias correction increases the detectability of the radiative effects of ACI in simulations of shallow cumulus.

## 1. Introduction

The radiative effect of the interactions between atmospheric aerosol and boundary layer clouds remains a large source of uncertainty in estimates of climate sensitivity to anthropogenic forcing (e.g., Boucher et al. 2013). Improving forecasts of the fraction of the sky covered by boundary layer clouds is also a priority for solar renewable energy applications (Perez et al. 2016). The need for an improved understanding of aerosol–cloud–radiation variability is clear, but the net effect of interactions between aerosol and clouds is difficult to understand because of the many possible cloud responses to aerosol variability (Stevens and Feingold 2009).

The first aerosol indirect effect on cloud albedo, or the “Twomey effect” (Twomey 1974) quantifies how increasing cloud condensation nuclei (CCN) leads to brighter clouds, assuming the same amount of cloud liquid is condensed. Continued research has shown that feedbacks between cloud microphysics and atmospheric dynamics mean that the total radiative effect of increasing aerosol has many complex components (Lohmann et al. 2010). Besides changing cloud brightness, aerosol–cloud interactions (ACI) can affect the total lifetime of a cloud field (Albrecht 1989), or the size distribution of clouds (Xue et al. 2008), with quantification of these effects depending sensitively on spatiotemporal analysis scale (McComiskey and Feingold 2012).

Covariation between meteorological drivers of aerosol and cloud properties in boundary layer liquid clouds was investigated in stratocumulus over land by Sena et al. (2016) using a 14-yr dataset of surface-based measurements. Similar covariation was investigated recently for marine clouds by Andersen et al. (2017) and Mieslinger et al. (2019) using multiyear datasets of satellite-based measurements. These studies showed expected strong relationships between certain variables, for example, lower-tropospheric stability and cloud occurrence, but found the role of aerosol variation on cloud brightness to be much harder to discern. The analysis of Sena et al. (2016) showed that, while holding cloud liquid constant, the surface aerosol index (a proxy for CCN concentration) was inconsistently correlated with the cloud radiative effect, meaning the Twomey effect was not easily detectable. It was suggested that covariation between shallow-cloud-controlling processes and larger-scale background meteorological forcing could be responsible for variability that reduces the detectability of the radiative effects of ACI (e.g., Liu et al. 2016).

To investigate the effects of both natural meteorological forcing and aerosol concentration covariability on the radiative effect of shallow clouds, we have availed ourselves of the Large-Eddy Simulation (LES) Atmospheric Radiation Measurement (ARM) Symbiotic Simulation and Observation (LASSO) project. LASSO follows the template of routine LES initiated at an instrumented site, similar to a project at Cabauw in the Netherlands (Neggers et al. 2012). The LASSO project combines detailed observations from the ARM Southern Great Plains (SGP) atmospheric observatory with meteorological reanalysis and idealized forcing profiles and bundles this data together for selected days. LASSO has initially focused on nonprecipitating shallow cumuli at the SGP site. Thus, LASSO is an ideal resource for LES-based investigation of the covariability between shallow boundary layer cloud processes and background meteorology. In the current work the focus is on the implications of this covariation for detectability of aerosol-related cloud brightening.

The bundles distributed by LASSO, however, do not include aerosol concentration observations or a representation of varying aerosol in simulated clouds. We have combined co-occurring observations of aerosol concentration at SGP during the 48 different days of initial conditions chosen by LASSO to be representative of nonprecipitating shallow cumulus at SGP during the summer of 2015, 2016, and 2017. Because of a lack of CCN observations on 12 days, we were only able to resimulate 36 of these 48 different LASSO days with our own implementation of aerosol variability added. These simulations with aerosol variability are not a rigorous evaluation of how well the simulations reproduce observations [for which we defer to Gustafson et al. (2019)]. Instead they create a collection of simulations that are constrained by observations and reproduce realistic cloud field properties. Importantly, instead of considering independent combinations of meteorological state and aerosol concentration (e.g., Lin et al. 2016), we have created a simulation database roughly bounded by the range of meteorological forcings and aerosol concentrations likely to naturally co-occur in the shallow continental cumulus cloud regime.

In this work, we use an analysis framework that defines cloud radiative effect variability as coming from variation in the amount and brightness of clouds (Betts 2007; Liu et al. 2011). We combine this analysis with a mutual information analysis (after Dawe and Austin 2013; Glenn and Krueger 2017), which quantifies how background meteorological variability in general masks the relationship between cloud drop number and cloud radiative effect. It will be shown that significant variability in average cloud brightness can be caused by the degree of horizontal inhomogeneity of the cloud optical depth over the domain of interest. We improve upon the cloud radiative effect analysis framework by implementing a method to control for these variables and demonstrate commensurate improvement in the detectability of the radiative effect of ACI.

## 2. Modeling, case selection, and aerosol inputs

LES was performed with the System for Atmospheric Modeling (SAM; Khairoutdinov and Randall 2003), with a horizontal grid spacing of 100 m and a vertical grid spacing of 30 m below 5-km altitude, above which the vertical spacing is incrementally increased to 300 m at 10-km altitude. The LES domain size was 24 km × 24 km in the horizontal with periodic boundaries, with a 15-km-altitude domain top. Radiation was calculated using the Rapid Radiation Transfer Model simplified for more efficient computation (RRTMG), using the plane-parallel independent column approximation (ICA). For the cloud microphysics, we used the two-moment bulk scheme of Morrison et al. (2005) in SAM, which specifies the aerosol size distribution and number concentration as a constant everywhere in the simulation domain. We altered this microphysics code to represent spatial and temporal variability of the aerosol concentration based on observations (see below). The microphysics assumes saturation adjustment, but activation of the aerosol was calculated in each LES column based on updraft speed, temperature, pressure, and an assumed cloud drop size distribution shape following Abdul-Razzak et al. (1998) and Morrison et al. (2005). Collision–coalescence, which is minimal in our simulations, was represented by Khairoutdinov and Kogan (2000).

Aerosol observations at the SGP site for LASSO days were obtained from the NOAA Earth System Research Laboratory’s Federated Aerosol Network (Andrews et al. 2019). We used data from a condensation nuclei (CN) counter and CCN counter. The CCN instrument settings varied over the 3-yr range spanned by the 36 days of this study, but typically supersaturation settings were cycled in increments of about 0.1% from a minimum of 0.2% to a maximum of 0.8% and back down again to 0.2% by the end of an hour, repeating that pattern each hour for each day. Each setting of the cycle was of approximately 6–10 min in duration.

The CCN concentration as a function of supersaturation was transformed to aerosol concentration as a function of aerosol particle size assuming the chemical properties of the species, which for simplicity, and in recognition of other larger uncertainties, we represented as 100% ammonium sulfate with a soluble ratio of 0.7. We then fit a two-mode lognormal distribution to these CCN-derived values of aerosol concentration as a function of size. This method can be considered a simplified form of that explored in Marinescu et al. (2019); our intent here is to capture the first-order temporal variability of the aerosol—capturing the first-order spatial variability will follow. To this end, we set the shape of the size distribution equal to the best-fit interpolation from the CCN measurements made in the morning for that day, and allowed the total concentration to vary over each day according to the total CN count. The data are available at approximately 10-min intervals from the CN counter, which we interpolate to 1-min intervals for input to the model.

Figure 1 shows a statistical summary of each time series of aerosol number concentrations measured by the CN instrument for all 36 simulated days in this study. We see that the constant aerosol concentration value used in default LES runs released by the LASSO project (black dashed line in Fig. 1) is at the upper end of observed aerosol number concentrations. For further analysis, we selected a subset of these simulations that pass strict thresholds for shallow cumulus convection (blue dots in Fig. 1). To minimize complications from including ice-processes and multiple cloud regimes in our analyses, we remove from consideration any simulation that does not produce a state consistent with warm-phase shallow cumulus, as described below.

LASSO simulation forcings are derived from meteorological reanalysis products that represent much larger spatial scales than the model domain. All of our simulations are initialized with the 0700 central daylight time (CDT; UTC − 5 h) local sounding from the SGP central facility and use the 300-km-scale variational reanalysis (VARANAL) LASSO product for surface and advective fluxes. The temporal resolution of the LES output used for all analysis in this work is 1 min.

A summary of the domain cloud fraction over times when clouds occur in the 36 simulations is shown in Fig. 2. We see that some simulations produce large cloud fractions, inconsistent with shallow cumulus. These large-cloud-fraction days are more likely to occur toward the end or beginning of the summer season, while days with small cloud fractions are more frequent in midsummer. The simulations with large cloud fractions have other characteristics that are indicative of deep rather than shallow convection, with large-scale ascent at 500 hPa, and produce cloud ice (see appendix). The 12 simulations selected for further analysis have no cloud ice and have a cloud fraction not exceeding approximately 50%; this is consistent with fair-weather nonprecipitating shallow cumulus convection.

Each simulation was initiated at 0700 CDT with aerosol concentration represented by a numerical prognostic tracer at each grid point that is initially equal to 1 at the surface and constant in a shallow morning mixed layer, assumed to be initially 500 m deep. The initial value assigned to the tracer decays exponentially above the mixed layer top with a decay-scale height of 2 km. This method is based on the observations and modeling of Turner et al. (2001) at SGP. As the simulation evolves through the diurnal cycle, the tracer concentration in the boundary layer is reduced as low tracer concentration air is entrained from the free troposphere. The simulated tracer value is given a scaled weight so that it always equals one at the surface, that is, set to equal the value derived from observations at the surface, but still always goes to zero at the top of the LES domain. One weighting function *W* for the whole domain is defined at each level *z*_{i} at each time step,

where the overbar indicates a mean over the entire domain. Each tracer column *T*_{r} is individually transformed to $Tr\u2032$ at each time step to maintain spatial variability; that is,

This means that our LES represents changes in the surface value of aerosol concentration due to boundary layer dilution, as well as changes due to larger-scale horizontal aerosol advection at low levels, because both are convolved in the variation of the surface observations through the course of the day. As the dilution proceeds explicitly in the model, the weighting function is continually updated based on the assumption that both dilution and large-scale advection are represented by the observations, but only the dilution is being represented explicitly by the model until the weighting function is applied.

Figure 3 illustrates the spatial and temporal aerosol variability in two of the simulations as representative examples for all the cases. The evolving vertical profile of the aerosol concentration tracer is represented by black lines indicating the height range of the layer containing simulated shallow cumulus. Figure 3 is truncated at 1000 CDT to highlight the period of cloud and boundary layer growth; the period beginning from initialization at 0700 CDT simply shows model spinup and the dry boundary layer growth. The aerosol concentration tracer follows the flow of the dynamics of the model, expressing spatial variability in the horizontal due to mixing, in addition to changes in the vertical profile, meaning that our simulated clouds exist in regions with 80%–90% of the aerosol concentration near the surface on average. But as shown in Fig. 3, with maximum cloud-top heights reaching 3 to 4 km, some of the clouds mix with air having an aerosol concentration as low as 30%–40% of the surface value.

Given variability in the meteorological conditions, this means that the cloud-averaged drop concentration $N\xafc$ will not always track the domain-average available surface aerosol number concentration,$\u2009NA\xaf$. We use the notation of an overline with a “*c*” superscript to indicate that an average is performed over all points in the domain defined as cloudy, as opposed to the plain overbar, which indicates the average over all points in the domain. Figure 4 presents these variables as time series for the same two days represented in Fig. 3. We see that whereas the 14 June 2017 simulation does have a $N\xafc$ time series positively correlated with $NA\xaf$ the 17 July 2017 simulation does not. All of the 12 selected days exhibit varying degrees of apparent temporal correlation between $N\xafc$ and $NA\xaf$, indicating that they exhibit varying degrees of ACI. We have already performed case selection to these specific idealized nonprecipitating cumuliform cases, leading us to ask why the ACI over our 12 selected days still exhibits considerable variability? Instead of attempting to untangle the radiative effects of meteorology and ACI variation in each simulation on a case by case basis, we pursue an analysis of all 12 selected days, as described in the next section.

## 3. Method

### a. Cloud averaging and inhomogeneity

To quantify the radiative effect of covariation between ACI and other background meteorological conditions in our simulations, we need to link the cloud microphysics on the smallest simulated scales to the total cloud radiative effect on the larger cloud-field scale. Many improvements to our quantitative understanding of the interactions between clouds and radiation have been made through use of the two-stream approximation to radiative transfer (Schuster 1905). An early example is Neiburger (1949), who showed that a stratus cloud’s albedo should approach unity at a particular rate as the thickness of the cloud becomes large. The cloud albedo is defined as the fraction of incoming radiation that is scattered by clouds, integrated over the upward and downward hemispheres. Meador and Weaver (1980) gave a general expression for the two-stream approximation solution for the relationship between cloud albedo *α* and cloud optical thickness *τ*_{0} in a broadband solar shortwave sense:

where *γ* depends on the degree of forward scattering, and it is assumed that radiation is only scattered and not absorbed, clouds are horizontally homogenous and in a single layer, the underlying surface is black, and the sun is at the zenith.

Since the sun is typically some angle *θ* from the zenith, we replace *τ*_{0} in Eq. (2) with *τ* ≈ *τ*_{0}/cos*θ*, as in Bohren (1987). Sagan and Pollack (1967) showed that the amount of scattered light over a full hemisphere can be approximated simply by the amount of scattering at a particular angle near 55° and −55° from the normal that defines the forward and backward hemispheres, respectively, giving *γ* = 2/(1 − *g*), with the scattering asymmetry parameter approximated as a constant *g* = 0.85 for liquid cloud drops at visible wavelengths (Hu and Stamnes 1993).

In essence, our analysis uses Eq. (2) to quantitatively link a radiative descriptor (albedo) to a descriptor of cloud liquid microphysics, for which we use optical depth *τ*_{0}. This framework was developed in Liu et al. (2011), and includes the concept of the relative cloud radiative effect (rCRE) after Betts (2007). The rCRE is intended to quantify the radiative effect of clouds in Earth’s atmosphere with minimal sensitivity to radiative variability caused by noncloud effects. The general form of this rCRE-based framework has been used by Xie and Liu (2013), Seifert et al. (2015), Sena et al. (2016), and Feingold et al. (2017) to give quantitative estimates of the magnitude of radiative-microphysical relationships. But this framework, even while idealized, is subject to the large horizontal inhomogeneity albedo bias described and corrected for in Cahalan et al. (1994), but neglected in these previous works. We apply this correction to the rCRE-based analysis framework here as follows.

As in Liu et al. (2011), the shortwave rCRE for a domain of interest is defined as the fraction of incoming shortwave radiation reflected by clouds (the average cloud albedo $\alpha \xafc$) multiplied by the fraction of the domain occupied by those clouds:

where we note the overline with a “*c*” indicates that an average is performed over all points defined as cloudy (not the full domain); that is, for any variable *x*, $x\xafc=x\u2061(\tau >1)\xaf$. We take the area of this same spatial region in our LES where *τ* = *τ*_{0}/cos*θ* > 1, and divide by the domain area to define the cloud fraction *f*. The threshold used for this definition is equivalent to a constant threshold of cloud radiative effect, similar to defining cloud shadow edges using a constant cloud shadow darkness threshold. This is most consistent with the LES radiation calculations, and therefore the calculated rCRE.

To evaluate Eq. (3), we must address the assumptions used in Eq. (2) to calculate the average cloud albedo. Liu et al. (2011) gave a slightly more complex expression than Eq. (3) that relaxed the assumption of nonzero surface albedo and nonzero cloud absorption, and this yielded better comparison between surface measurements of the shortwave rCRE and satellite and model output. Liu et al. (2011) also outlined a method to treat multilayer clouds as an “effective” single-layer cloud, as in Ramanathan (1987). This effective cloud layer can account for vertical heterogeneity, that is, multilayer stratiform clouds, but does not account for horizontal inhomogeneity in a single layer.

The problem of assuming horizontal homogeneity is the same as the “plane-parallel albedo bias” problem, which was explained using a simple model of the fractal nature of stratocumulus by Cahalan et al. (1994). The horizontally averaged albedo of a cloud field cannot be uniformly related to the horizontally averaged optical depth without knowing the subaveraging-scale distribution (horizontal inhomogeneity) of the cloud layer. Cahalan et al. (1994) showed that fields of clouds are not randomly inhomogeneous, but have characteristic spatial distributions of properties that are well described by fractals.

Our simulations approximately resolve the scale at which cloud optical depth homogeneity can be assumed (LES columns of 100 m × 100 m horizontal) over a domain large enough to include many clouds with varying optical depths where cloud heterogeneity cannot be ignored (24 km × 24 km). We take the average of Eq. (2) and define an inhomogeneity correction factor *c*_{1} after Cahalan et al. (1994):

This measure of the horizontal inhomogeneity can be calculated at any particular time in the simulations as $log10\u2061(c1)=log10\u2061(\tau 0\xafc)\u2212log10\u2061(\tau 0)\xafc$. Note the cosine of the solar zenith angle is a constant over the whole domain and is factored out of the average without affecting *c*_{1}. In general, *c*_{1} will only equal 1 if the distribution of optical depth values *τ*_{0} being averaged is monodisperse.

A similar “inhomogeneity” problem is encountered when we attempt to use the idealized relationship between optical depth *τ*_{0}, cloud liquid water path *L*, and cloud drop number concentration *N* expressed by Boers and Mitchell (1994), Zhang et al. (2005), and others, which is valid at some idealized subcloud scale. Once again, to apply this relationship to the large scale would assume cloud horizontal homogeneity, so we approximate the cloud average as

where *L* is in grams per meter squared, *N* is cloud drop number per volume of the cloudy part of the column (cm^{−3}), and $cT,p$ is a function of temperature and pressure in an adiabatic case. Following Marshak et al. (2006), a reduced constant factor *c*_{0} is likely to be the best correction for this type of averaging bias, as opposed to changing the magnitude of the powers on *L* or *N*. We obtain *τ*_{0} for the lhs of Eq. (4) and the $L\xafc$ and $N\xafc$ for the rhs from the gridscale microphysics in each LES, averaging over all grid columns with *τ* > 1 at each minute from the 12 selected simulations when cloud exists (over 100 h). We find that setting *c*_{0} = 0.07 for all times in our dataset minimizes the mean absolute error between the lhs and rhs of Eq. (5) at 1.35 in units of optical depth. Given values of optical depth typically larger than 10, we used Eq. (5) to evaluate the rCRE in Eq. (4), substituting $\tau 0\xafc$ with $L\xafc$ and $N\xafc$, giving a typical mean absolute error in rCRE of a few percent due to this assumption at any particular time in our analysis. The total error in approximating rCRE using Eqs. (3)–(5) is reduced further by setting *c*_{1} to a best-fit value at each particular time in our analysis (see section 4c).

### b. Mutual information

In parallel with using the equations defined above, we also quantify covariability between meteorological processes affecting cloud amount and cloud brightness using a statistical technique. Based on the understanding of information as entropy (Shannon 1948), one can measure the mutual information (MI) between two variables. Following Dawe and Austin (2013),

where *P*(*A*), *P*(*B*), and *P*(*A*, *B*) are the marginal and joint probability density functions (PDFs) for the variables *A* and *B*. For example, if *P*(*A*) and *P*(*B*) are perfectly independent, the mutual information equals zero, and if perfectly correlated or anticorrelated, the MI is one.

An advantageous statistical property of MI is that it does not assume a functional form for the relationship between variables, so it measures linear and nonlinear relationships equally well. Another advantage is that we can calculate the MI between some variable *A* and *B*, but first condition the PDF of *B* on a third variable *C*. This is equivalent to finding the MI between *A* and *B* along surfaces of constant *C*, revealing the information *A* and *B* share with each other while the variability that is explained by *C* is effectively removed. This is called the conditional mutual information (CMI):

For example, if the MI between *A* and *B* is independent of *C*, the MI and the CMI will be equal to each other, indicating that *C* does not contain significant information about the *A*–*B* relationship. But if some variability in *C* matches the pattern of variability in *B* (either a correlation or anticorrelation), the CMI will be smaller than the MI due to *C* containing information about *A* that is effectively redundant given *B*. Finally, if the CMI is larger than the MI by some amount, this indicates the amount of unique information about *A* that has been added by *C* in addition to *B*. Similar to a traditional budget analysis, we will use this CMI analysis to find which variables are needed to reproduce the total rCRE variability, and to quantify how much information each variable contributes.

## 4. Analysis

### a. Overview of variability

Substituting Eqs. (4) and (5) into Eq. (3) shows that the rCRE is a function of *f*, $L\xafc$, and $N\xafc$, but also of the solar zenith angle *θ* and the degree of horizontal inhomogeneity *c*_{1}. The domain-scale variation in these properties (from time step to time step, and between different days) comes from the net result of variation in the domain-scale meteorological forcing that drives the simulations. The variation of $N\xafc$ relative to the variation of $NA\xaf$ throughout a day, or between days, is also the net result of variation in the domain-scale meteorological forcings. We will use the term “meteorological variability” to indicate any meteorological properties indicative of variability in this specific context of $N\xafc$ relative to $NA\xaf$ variation, that is, meteorological variability directly related to ACI. Examples of meteorological properties directly related to ACI could include the cloud-base updraft speed, the vertical gradient of the aerosol concentration, whether the cloud is entraining air with much lower aerosol concentration than at cloud base, etc. The goal of this work is not to determine the most important meteorological variables, but to develop a framework in which the radiative effects of such variability can be quantified, allowing for future work to rank ACI-related meteorological processes by their radiative effect.

The net result of other variability in meteorological processes realized in our simulations will be referred to as “cloud field property variation.” Examples of meteorological variations that are realized as cloud field property variations are the various input soundings and advective tendencies that drive the simulations that lead to different values of $L\xafc$ and *f*. The amount that the inhomogeneity *c*_{1} is less than 1 is proportional to the width of the distribution of *τ*_{0} that is averaged in $\tau 0\xafc$. Further investigation will be needed to determine if variation in the value of *c*_{1} for the same $\tau 0\xafc$ is more likely ACI-related meteorological variability (i.e., variation in *N* for the same $N\xafc$) or cloud field property variability (i.e., variation in *L* for the same $L\xafc$). Based on the form of Eq. (5), we suspect cloud field property variability contributes to most to *c*_{1} variability. In this section, we will attempt to quantify the radiative effects of ACI-related meteorological variability relative to the concurrent radiative effects of cloud field property variability. In Fig. 5 we quantify and illustrate the range of variation exhibited by rCRE, *f*, $L\xafc$, and $N\xafc$ for the 12 selected simulations.

Based on Eq. (3), the variability along a line of constant *f* in Fig. 5 shows the extent to which rCRE varies with $\alpha \xafc$. The left panel of Fig. 5 shows that the simulations span an albedo range from about 30%–60%. The cloud fraction ranges from 0% to just above 50%, although noting the probability density scale, a cloud fraction below 20% is most common. The middle panel shows the same distribution, but now colored by the average of all of the $N\xafc$ values; and the rightmost panel shows the same but colored for the average of $L\xafc$ values. It may appear from Fig. 5 that variability in cloud albedo is more strongly related to variability in $L\xafc$ than to $N\xafc$. However, we cannot conclude which pattern of variation ultimately contains more information about the rCRE variation without considering covariability between $N\xafc$, $L\xafc$, and *f*.

### b. Mutual information

Which variables ultimately contain the most information about the total rCRE? To quantify these relationships, we calculate the mutual information between the rCRE and the *f*, $L\xafc$, and $N\xafc$ variables from all 12 simulated days, at 1-min resolution when clouds exist. To obtain confidence intervals, we randomly (with replacement) subselected 62% of our dataset, bootstrapping in this way 1000 times. For a measure of the signal strength, we find the mutual information between a random permutation of the rCRE variable (effectively noise) and *f*, $L\xafc$, and $N\xafc$, and report these results as MI^{*} and CMI^{*}. This gives a measure of the baseline signal above which we have more confidence in the result. The results of the calculations are shown in Table 1. The first column of Table 1 quantifies what can be seen qualitatively in Fig. 5: that the current value of rCRE is strongly related to the current value of *f*, less so for $L\xafc$, and even less so for $N\xafc$. That is, over all simulations, 51% of the information content of the rCRE variation is shared with the *f* variability, 18% with the $L\xafc$variability, and only 11% with the $N\xafc$ variability. These information contents are not necessarily unique sets of information about rCRE.

To understand this further, we calculate the CMI in the second column of Table 1. Conditioning the MI between rCRE and *f* on $L\xafc$ results in a CMI of 70%, similar to the unconditioned MI between rCRE and *f* of 51% added to the unconditioned MI between rCRE and $L\xafc$. This indicates that knowing the covariation between $L\xafc$ and *f* gives no additional information about rCRE if $L\xafc$ and *f* are already known. In contrast, conditioning the MI between rCRE and $L\xafc$ on $N\xafc$ increases the shared information content to 50% compared to the MI of either $L\xafc$ or $N\xafc$ alone. This is greater than the 18% + 11% that might be expected as a maximum CMI for these two variables. This indicates that not only does $N\xafc$ contain unique information about the rCRE, but even more so if $L\xafc$ is known; that is, the covariation between rCRE and $L\xafc$ acts to mask the strength of the rCRE–$N\xafc$ relationship. Most notably, the maximum information about the rCRE from two variables is obtained from the pair of (*f*, $N\xafc$), with a CMI of 75%. This relationship also shows that the pattern of covariation between rCRE and *f* acts to mask the strength of rCRE–$N\xafc$ relationship.

What does it mean for one variable to mask the dependence of another? This can be understood by considering Fig. 5 to be a 3D landscape, with a given (*f*, $\alpha \xafc$) coordinate defining a particular rCRE, and movement in a given direction describing a change in the rCRE. From the differential analysis, one might have expected the $L\xafc$ variable to give the best approximation to $\alpha \xafc$. Indeed, they do have a strong relationship, but the pattern of variation on this landscape of cumulus clouds is such that much of what the $L\xafc$ tells us about $\alpha \xafc$ can be implied from knowing both *f* and $N\xafc$. Put another way, instead of holding $L\xafc$ constant, we used MI to account for $L\xafc$ variation by its covariation with *f,* and the remaining variation in rCRE was even better explained by $N\xafc$. Holding $L\xafc$ constant is not the same as holding all cloud field property variability constant; there is some other background meteorological variability that acts to mask the strength of the relationship between $N\xafc$ and rCRE.

### c. Contributions to albedo variability

To investigate which sources of cloud field property variation are affecting rCRE in this way, we normalize the rCRE variability by *f*. From Eq. (3), this means we only consider variability in rCRE due to variation in $\alpha \xafc$. In the top row of panels in Fig. 6, we plot the $\alpha \xafc$ versus $L\xafc$, and along any particular slice of $L\xafc$, colors are used to show the variability in $\alpha \xafc$, $N\xafc$, and $NA\xaf$. The top right panel of Fig. 6 illustrates that some of the $\alpha \xafc$ variability is likely due to the increased slant path of radiation through clouds at different solar zenith angles.

To understand the radiative effect of ACI over the 12 simulated days, we compare the change in albedo from different simulations having different $N\xafc$ and $NA\xaf$ while holding all else equal. We can approximately account for the radiative effect of illumination angle variability by binning our data by $L\xafc/\u2009cos\u2061\theta $ (bottom row of panels in Fig. 6). This represents how we expect the radiative effect of a long slant path through a small $L\xafc$ to be equivalent to a shorter path through a larger $L\xafc$. The bottom right panel of Fig. 6 suggests that variation in the distribution of cloud optical depths being averaged (inhomogeneity *c*_{1}) also contributes to variation in average cloud brightness. The form of Eq. (4) shows the dependence of $\alpha \xafc$ on *c*_{1}, so we can equivalently remove the effect of this variability by binning by $\alpha \xafc/c1$ (Fig. 7). In this case, the signal of drop concentration variability on cloud albedo becomes even more clear. The surface aerosol signal still appears somewhat inconsistently related to the $N\xafc$. But now this variation approximately represents variability in ACI, because Fig. 7 illustrates the variation in albedo when variations due to cloud field heterogeneity are approximately removed.

### d. Contributions to rCRE variability

We quantify the total response of the rCRE to an aerosol-related cloud drop number perturbation, similar to Platnick and Oreopoulos (2008), Quaas et al. (2009), and others, by differentiating Eq. (3) with respect to $N\xafc$. To numerically separate variation in rCRE due to $\alpha \xafc$ from fluctuations in *f*, we first take the natural log of both sides of Eq. (3), separate terms using the logarithmic product identity, and then perform differentiation:

Equation (8) says that the total change in the rCRE for a change in $N\xafc$ can be expressed mainly by three different terms. It is not exact because the approximations in Eqs. (4) and (5) are not exact. The first term on the rhs of Eq. (8) is the partial derivative of Eq. (4) and could be approximated by the increase of albedo with $N\xafc$ along lines of constant $L\xafc$ in Fig. 6. But some of the variability in albedo with $N\xafc$ at constant $L\xafc$ is due to variability in cos*θ* and *c*_{1}. Thus, the *partial* change in $\alpha \xafc$ due only to the radiative effect of $N\xafc$ variability is better approximated by Fig. 7, that is, after normalization has been performed.

Similarly, the second term on the rhs of Eq. (8) is the *partial* change in $\alpha \xafc$ due only to the radiative effect of $L\xafc$ variability, scaled by the actual concurrent change in $L\xafc$ with the $N\xafc$ variation. As before, if we attempt to compare the change in $\alpha \xafc$ for some $L\xafc$ on one day to the change in albedo for the same $L\xafc$ on a different day, we must normalize by the inhomogeneity *c*_{1} and solar zenith angle cos*θ*, which also likely vary between days with the same $L\xafc$. The final term on the rhs of Eq. (8) is the total response of the cloud fraction concurrent with the domain-sensible perturbation of the $N\xafc$. Like the $L\xafc$ contribution term, this *f* term represents cloud field property variation occurring at the same time as the ACI-related cloud drop number variation.

We are motivated to quantify the cloud field property variation to understand if it is large enough to overwhelm the signal of an aerosol indirect effect. To quantify the concurrent contributions of each variable to the rCRE budget, we perform temporal numerical differentiation to approximate the values of the terms in Eq. (8). The results of temporal numerical differentiation can depend sensitively on the method and time scale used, so care must be taken. We quantify the relative contributions of changes to cloud brightness and amount due to both aerosol variability and concurrent background meteorological variation on a time scale relevant for shallow cumulus. This implies a time scale slightly shorter than the characteristic time scale of mixing aerosol into the clouds, or about 5 min. We have used 3, 5, 10, 15, 20, and 30 min for the differentiation interval to examine the sensitivity of this method. For the shorter time scales, the results are qualitatively similar, but starting around 15 min, the time scale is too long to resolve the relatively fast changes in the cumulus layer, and temporal differentiation begins to describe the differences between different cloud-field states, instead of quantifying the magnitude of the changes that occur as the field evolves. Here we show the results for model output evaluated using a 5-min differentiation time scale. We used the simple form of “forward” numerical differentiation:

where *i* indicates a particular time and *δ* is the temporal differentiation interval. In practice we do this for each term on the rhs of Eq. (8) as well as the lhs, but we omit this from Eq. (9) for clarity. A statistical summary of the range of values found for the differential terms in Eq. (8) is shown in the left panel of Fig. 8. We also performed differentiation as in Seifert et al. (2015) and previous work, letting contributions from inhomogeneity *c*_{1} contribute to calculated changes in albedo by always setting *c*_{1} = 1 (Fig. 8, right panel).

The terms in Fig. 8 are unitless log-relative changes; that is, for the total response of the rCRE to a perturbation of $N\xafc$ (gray), a value of 100% indicates that a doubling in $N\xafc$ is accompanied by a doubling in rCRE, whereas −100% indicates that a doubling in $N\xafc$ is accompanied by a halving of rCRE. We see that for the total response of the rCRE (gray), the median is small and positive, but the collection of all simulated rCRE responses covers a wide range and can be negative as well as positive. This is consistent with Sena et al. (2016), who found for stratocumulus over SGP that correlations between rCRE and a measure of surface aerosol concentration are just as likely to be positive as negative.

By comparing the left panel with the right panel of Fig. 8, we see that ignoring the variability of inhomogeneity *c*_{1} manifests primarily as a wider range of partial albedo changes associated with changes in $L\xafc$. By accounting for *c*_{1} we have increased the accuracy of quantification of the average cloud albedo, allowing us to better isolate the effect of varying aerosol from other processes. This illustrates how detectability of the radiative effect of changes in $N\xafc$ is increased by accounting for inhomogeneity *c*_{1}.

### e. ACI and meteorological variability

Now that we have quantified the radiative effect of variation in $N\xafc$, we can investigate the effect of meteorological conditions leading to different amounts of $N\xafc$ being realized for the same prescribed aerosol $NA\xaf$, and thus we gain understanding of the radiative effect of meteorological covariation with ACI. We can express the albedo effect of aerosol concentration variation mathematically as

and use the final term on the rhs as a metric of ACI. In Fig. 9 we plot the $N\xafc$ versus $NA\xaf$ for all 12 simulations; each point is a snapshot every 1 min when clouds are present. In the top left panel, we include dashed lines to indicate potential constant values of $d\u2061lnN\xafc/d\u2061lnNA\xaf=\u2009ACI$. We color the points by a set of meteorological fields of interest: $L\xafc$, the vertical velocity averaged over cloud-base points, the inhomogeneity *c*_{1}, and the decoupling index. The decoupling index is the difference between the lifted condensation level (LCL) and the actual cloud base, normalized by the LCL altitude. A small value of the decoupling index indicates a cloud likely well connected to the surface, and should have higher correlation between $N\xafc$ and $NA\xaf$. Because aerosol activation is a function of updraft velocity, we expect stronger average cloud-base updraft speeds to be associated with higher values of $N\xafc$ for the same $NA\xaf$. Variability in the inhomogeneity *c*_{1} and $L\xafc$, which represents cloud field property variation, are included for comparison.

Clusters of points on Fig. 9 approximately represent different simulated days, although there are some regions of overlap. Different days appear to realize very different values of the ACI metric. Some points that have higher $NA\xaf$ for the same $N\xafc$ also have higher decoupling index values. This could be part of the explanation for the lack of correlation between $N\xafc$ and $NA\xaf$, though we see from the pattern of variation with inhomogeneity and cloud-base vertical velocity that these factors could contribute as well.

One might expect a maximum value of the ACI metric to be a slope of 1, but only if we were considering some idealized subcloud scale, for example, *d* ln*N*/*d* ln*N*_{A}. The appearance of slopes steeper than one or negative on Fig. 9 is representative of the complex interactions between the domain-scale quantities $N\xafc$ and $NA\xaf$. A process potentially exhibiting strong control of ACI here is the mixing of upper-tropospheric air into clouds rooted in the boundary layer. From Fig. 3, we know some clouds can entrain and mix with free tropospheric air that has an aerosol concentration much less than the value at the surface, even if the boundary layer is well mixed. This would lead to a pattern on Fig. 9 of the $NA\xaf$ changing slowly, while the $N\xafc$ changes more abruptly, often decreasing sharply, leading to steep slopes and a large value of the ACI metric.

## 5. Discussion and conclusions

We have defined a framework for quantification of the radiative effect of variation in ACI relative to other cloud field property variations, including a correction for horizontal heterogeneity of the cloud field. We applied this framework to LES of the diurnal cycle of nonprecipitating shallow cumulus over land, and by availing ourselves of the LASSO project database, performed 36 simulations constrained by observations of meteorological forcing and aerosol concentration at SGP. We restricted our analysis to simulations of ice-free shallow cumulus convection, selecting 12 of these days for further analysis. Using mutual information (MI) we investigated fluctuations in rCRE due to aerosol variability (represented by variation in $\alpha \xafc$ with respect to $N\xafc$) relative to the radiative effect of other cloud field property variability represented by variation in $L\xafc$ and *f*.

MI analysis quantified how much variation in rCRE is explained by *f*, $L\xafc$, and $N\xafc$ when considered alone. CMI analysis quantified how much variability in rCRE could be explained by these variables when considered in context with their covariation, that is, how much unique information they contain about rCRE. We showed that although $L\xafc$ alone does share more information with rCRE than $N\xafc$, it has very little unique information relative to *f*. After *f* is accounted for, the remaining variability in rCRE is even better explained by $N\xafc$ than we would expect from the MI between rCRE and $N\xafc$ alone. This may seem counterintuitive, because after accounting for *f*, the remaining variability in rCRE should be due only to $\alpha \xafc$, which is a stronger function of $L\xafc$ than $N\xafc$ (Fig. 6).

By transforming by both inhomogeneity parameter *c*_{1} and solar zenith angle (Fig. 7), we visually explain this counterintuitive result: holding $L\xafc$ constant is not the same as keeping all else equal between cases. The dependence of albedo on cos*θ* and *c*_{1} were affecting the calculated strength of the relationship between rCRE and $N\xafc$. By accounting for and normalizing by these factors, we were able to significantly increase the detectability of calculating albedo differences due to $\u2009N\xafc$ variation alone.

We then performed a differential analysis based on Eq. (8) to quantify the magnitude of the contributions to rCRE from changes in *f*, $L\xafc$, and $N\xafc$. A summary of the terms contributing to this rCRE budget (Fig. 8) illustrates that the radiative effect of aerosol-related $N\xafc$ variation on $\alpha \xafc$ is small compared to the variability in radiative effect from $L\xafc$ and *f* variation. We note that the total change in rCRE concurrent with an $N\xafc$ perturbation (gray bar in Fig. 8) appears to have a similar range, but larger magnitude, than the concurrent *f* response to the $N\xafc$ perturbation (green bar). This illustrates that variability in cloud fraction *f* cannot fully explain variability in the total rCRE response to an aerosol perturbation; cloud-field albedo variability on similar time scales must contribute significantly to make the total rCRE change larger. This differential analysis shows that the magnitude of the $L\xafc$ response (red bar) could explain most of this larger value of the rCRE response, while the albedo variation with droplet number has a small range and consistently small positive magnitude that also contributes.

The radiative effect of ACI on a spatial scale containing multiple clouds (or a stratus cloud with horizontal heterogeneity) is different than the radiative effect of ACI on the subcloud scale. This point has been overlooked in some recent quantitative analyses (e.g., Seifert et al. 2015; Sena et al. 2016; Feingold et al. 2017) where a horizontal heterogeneity correction factor for the albedo equation was not included. In this work we have demonstrated the practicality and utility of making the plane-parallel albedo bias correction of Cahalan et al. (1994) in the novel context of analysis of ACI in LES.

This albedo bias has been discussed in the context of global-scale estimates of the radiative effect of ACI by Barker (2000) and many others, but also recently in Merk et al. (2016) and Gryspeerdt et al. (2019) directly alongside discussion of corrections for cloud subadiabaticity and 3D radiative effects. We stress that this albedo bias is due to a homogeneity assumption that is important and possible to correct for, even in idealized cases, neglecting subadiabaticity or when analyzing a numerical model with non-3D radiative transfer.

By comparing the blue box and whiskers between the left and right panels in Fig. 8, we see the effect of correcting for the albedo bias when calculating the susceptibility of the albedo to a drop number perturbation: ignoring cloud field inhomogeneity results in a high bias of the Twomey albedo susceptibility of several percent. When inhomogeneity is accounted for in the partial derivative of Eq. (4) with respect to $N\xafc$, not only is the maximum albedo susceptibility lower, but the albedo value at which that maximum occurs is lower by a factor proportional to the characteristic value of *c*_{1} (about 0.75). The inhomogeneity-corrected framework presented here could potentially be used to improve global-scale estimates of the radiative effect of ACI that are made from local-scale measurements, because *c*_{1} can also be derived from local-scale measurements. Letting *c*_{1} < 1 effectively allows one to remove errors accumulated by spatially averaging data to a larger spatial scale, given any particular scale change. The spatial-scale change from order 100 m to order 10 km and the shallow cumulus cloud type likely maximize this effect; stratus cloud fields may exhibit less heterogeneity of albedo and thus less error if *c*_{1} ≈ 1 is assumed.

We note that in the context of anthropogenic aerosol perturbations on clouds, this work shows that any quantification of the radiative effect of ACI will take place against a backdrop of considerable meteorological variability. In the case of shallow cumulus, this variability ultimately makes the radiative effect of an aerosol perturbation more difficult to detect and quantify. This is because the variability in albedo from cloud to cloud in a field of clouds (the “horizontal heterogeneity” of the cloud field albedo) can vary systematically with meteorology, masking the dependence of $\alpha \xafc$ on $N\xafc$ for the larger cloud field average scale. We have illustrated how accounting for this variability via a best-fit value that changes with meteorology can aid detectability and quantification of the radiative effects of ACI.

## Acknowledgments

This work was supported by the Office of Biological and Environmental Research of the U.S. Department of Energy Atmospheric System Research Program Interagency Agreement DE-SC0016275. The authors acknowledge the NOAA Research and Development High Performance Computing Program for providing computing and storage resources that have contributed to the research results reported within this paper. The 2017 LASSO data bundles are from the Atmospheric Radiation Measurement user facility. The released days from years 2015, 2016, 2017 at 36°36′18.0″N, 97°29′6.0″W are from the Southern Great Plains Central Facility (C1) and were compiled by W. I. Gustafson, A. M. Vogelmann, X. Cheng, S. Endo, K. L. Johnson, B. Krishna, Z. Li, T. Toto, and H. Xiao. The dataset was accessed in 2018 from the ARM Data Center in Oak Ridge, Tennessee (https://doi.org/10.5439/1342961). The data from this study are available online (https://esrl.noaa.gov/csd/groups/csd2/clouds).

### APPENDIX

#### Case Selection

Here we provide more details about the 36 simulated days based on LASSO cases with added aerosol variability based on observations. Figure A1 shows the distribution of simulated cloud depths, the decoupling index, the vertical velocity at 500 hPa, and the ice water path (IWP) for all 36 simulations. The decoupling index is the difference between the LCL and the actual cloud base, normalized by the LCL altitude, and is large for some of the simulations that were not selected. Some of these cases have shallow cloud depths and large cloud fractions, indicating a stratocumulus regime.

Many of the other days not selected have high IWP, cloud depths extending more than 3 km, and large and positive vertical velocity at 500 hPa. This indicates a deep convective cloud regime. In contrast, the 12 days selected for further analysis have zero IWP, relatively shallow cloud depths, subsidence at 500 hPa, and some of the lowest values of the decoupling index, representative of shallow cumulus.

## REFERENCES

*J. Geophys. Res.*

*Science*

*Atmos. Chem. Phys.*

*Bull. Amer. Meteor. Soc.*

*J. Climate*

*J. Geophys. Res.*

*Tellus*

*Amer. J. Phys.*

*J. Atmos. Sci.*

*Atmos. Chem. Phys.*

*J. Geophys. Res. Atmos.*

*Geophys. Res. Lett.*

*Atmos. Chem. Phys.*

*J. Climate*

*Mon. Wea. Rev.*

*J. Atmos. Sci.*

*J. Atmos. Sci.*

*J. Atmos. Sci.*

*Atmos. Chem. Phys.*

*Atmos. Chem. Phys.*

*Atmos. Chem. Phys.*

*J. Geophys. Res.*

*Atmos. Chem. Phys.*

*J. Atmos. Sci.*

*Atmos. Chem. Phys.*

*J. Geophys. Res. Atmos.*

*J. Atmos. Sci.*

*Bull. Amer. Meteor. Soc.*

*J. Meteor.*

*Found. Trends Renewable Energy*

*J. Geophys. Res.*

*Atmos. Chem. Phys.*

*J. Geophys. Res.*

*J. Geophys. Res.*

*Astrophys. J.*

*J. Adv. Model. Earth Syst.*

*Atmos. Chem. Phys.*

*Bell Syst. Tech. J.*

*Nature*

*Geophys. Res. Lett.*

*Atmos. Environ.*

*Environ. Res. Lett.*

*J. Atmos. Sci.*

*Quart. J. Roy. Meteor. Soc.*

## Footnotes

Denotes content that is immediately available upon publication as open access.

Climate Change 2013: The Physical Science Basis, T. F. Stocker et al., Eds., Cambridge University Press, 571–657