Biomass Burning Aerosol Heating Rates from the ORACLES 2016 and 2017 Experiments

Aerosol heating due to shortwave absorption has implications for local atmospheric stability and regional dynamics. The derivation of heating rate profiles from space-based observations is challenging because it requires the vertical profile of relevant properties such as the aerosol extinction coefficient and single scattering albedo (SSA). In the southeast Atlantic, this challenge is amplified by the presence of stratocumulus clouds below the biomass burning plume advected from Africa, since the cloud properties affect the magnitude of the aerosol heating aloft, which may in turn lead to changes in the cloud properties 25 and life cycle. The combination of spaceborne lidar data with passive imagers shows promise for future derivations of heating rate profiles and curtains, but new algorithms require careful testing with data from aircraft experiments where measurements of radiation, aerosol and cloud parameters are better collocated and readily available. In this study, we derive heating rate profiles and vertical cross sections (curtains) from aircraft measurements during the NASA ObseRvations of CLouds above Aerosols and their intEractionS (ORACLES) project in the southeastern Atlantic. 30 Spectrally resolved irradiance measurements and the derived column absorption allow for the separation of total heating rates into aerosol and gas (primarily water vapor) absorption. The nine cases we analyzed capture some of the co-variability of

Lidar instruments are specifically designed to penetrate through the atmosphere to obtain vertical profiles of atmospheric backscatter and extinction, including those specific to aerosols. When positioned on satellites, lidar instruments can provide 65 aerosol extinction profiles for large spatial regions, which can be combined with additional observations that provide aerosol optical properties (e.g., Deaconu et al., 2019), or they can be incorporated into climate models (Mallet et al., 2019;Tummon et al., 2010;Gordon et al., 2018;Adebiyi et al., 2015;Wilcox 2010).
With current backscatter lidars, extinction profiles cannot be obtained directly due to the convolution of aerosol backscattering and extinction in the returned lidar signal. For example, the Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP) 70 aboard the Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation (CALIPSO) satellite measures the attenuated backscatter coefficients at 532 and 1064 nm, from which the extinction coefficient is typically retrieved based on modeled values of the extinction-to-backscatter ratio (also called the lidar ratio) inferred for each detected aerosol layer (Young and Vaughan, 2009). The depolarization ratio method (DRM) (Hu et al., 2007;Chand et al., 2009;Deaconu et al., 2017;Kacenelenbogen et al., 2019) also known as opaque water cloud (OWC) method takes advantage of the CALIOP 75 depolarization measurements and the transmission constraints provided by underlying low liquid water clouds to derive accurate measurements of the above cloud aerosol optical depth (ACAOD) and the mean lidar ratio of the aerosol layer above the cloud (Liu et al., 2015).
High Spectral Resolution Lidars (HSRL) are likely to be included in future space architectures. Among other advantages, they will provide extinction profiles directly (Hu et al., 2007;Hair et al., 2008). This would be useful with the new concept of 80 Heating Rate Efficiency (HRE) first introduced in this paper (Section 5). It represents the heating rate at any altitude as obtained from aircraft observations per aerosol extinction coefficient, which, by contrast to the heating rate itself, varies little throughout the profile. The HSRL-derived extinction profiles could be directly translated into aerosol heating rates for regions where HRE is available, bypassing, and thus directly constraining, radiative transfer calculations.
To develop the HRE concept, we use a combination of in situ and remote sensing observations from two years of the 85 ORACLES aircraft experiments (2016 and 2017). The aircraft is ideally suited for obtaining vertical information, which can be taken from varied instruments for different parameters. In situ measurements can provide profiles of water vapor while aerosol extinction profiles can be provided by either overflights with the High Spectral Resolution Lidar 2 (HSRL-2) (Burton et al., 2018) or vertical profiles with the Spectrometer for Sky-Scanning, Sun-Tracking Atmospheric Research (4STAR; Dunagan et al., 2013;LeBlanc et al., 2020) instrument. 90 For 9 cases (corresponding to those presented in Cochrane et al., 2021), we attribute the total heating of the layer to contributions from the aerosol, water vapor, and other atmospheric gases. We also examine the dependence of the aerosol heating rates and HRE on aerosol properties and cloud albedo -the same parameters that also modulate DARE (Cochrane et al., 2021). Section 2 of this paper describes the data and general methods of the heating rate calculations. Section 3 provides the results and discussion for heating rates segregated by absorber, while section 4 describes heating rate results along an 95 aircraft flight leg, known as a heating rate curtain. Section 5 introduces the new HRE concept and section 6 provides a summary.

Data
The NASA ORACLES project utilized the NASA P-3 aircraft for the 2016, 2017, and 2018 deployments, while the ER-2  100 aircraft was additionally deployed in 2016 to record high altitude measurements . Both aircraft were equipped with instrumentation to sample clouds and aerosols. The P-3 payload included both remote sensing and in situ instruments, while the ER-2 carried only remote sensing instrumentation.
To determine heating rates, we use a combination of measurements taken from the Solar Spectral Flux Radiometer (SSFR), 4STAR and HSRL-2. SSFR is a system consisting of two pairs of spectrometers that measure the upwelling (nadir) and 105 downwelling (zenith) irradiance between 350 and 2100 nm. The zenith light collector is mounted on an Active Leveling Platform (ALP), which keeps the light collector level with true horizon during flight. SSFR is radiometrically calibrated preand post-mission, with field calibrations performed throughout each deployment to keep track of instrument changes throughout the duration of the experiment (Cochrane et al., 2019). The 4STAR instrument measures aerosol optical depth (AOD) above the aircraft between 350 nm and 1650 nm, as well as provides column gas retrievals such as water vapor and 110 ozone, aerosol intensive properties and cloud properties (Segal-Rozenhaimer et al., 2014;Pistone et al., 2019;LeBlanc et al., 2021). The instrument is calibrated pre-and post-mission using the Langley method at the Mauna Loa Observatory (e.g., Schmid and Wehrli, 1995;LeBlanc et al., 2020). HSRL-2 provides vertical profiles of aerosol backscatter and depolarization at 355 nm, 532 nm, and 1064 nm wavelengths, along with aerosol extinction at 355 nm and 532 nm wavelengths, all measured below the aircraft (Hair et al., 2008;Burton et al., 2018). In 2016, HSRL-2 was mounted on the ER-2 but transitioned onto the 115 P-3 for 2017 and 2018.
A major benefit of aircraft campaigns is the flexibility to perform distinctive flight maneuvers in a way that optimizes measurements and retrievals for different instrumentation. In this work, we use two distinct flight patterns: stacked legs (socalled radiation walls) and vertical profiles (so-called square spirals) to derive a) heating rate curtains and b) heating rate profiles segregated by absorber, respectively. Figure 1 shows the locations of the radiation walls and spirals used in this study. 120 Radiation walls, the traditional maneuver for radiation science flights, consist of vertically stacked legs flown in sequence along a fixed ground track. The stacked legs bracket the aerosol layer above and below, at the bottom of the layer (BOL) and at the top of the layer (TOL). For ORACLES, the BOL leg was located just below the aerosol layer and just above the cloud layer. The column AOD (spectral), ozone, and water vapor are measured by 4STAR, and spectral scene albedo is measured by SSFR. The TOL leg is above the aerosol layer and the cloud, from which HSRL-2 measures profiles of extinction at 532 125 nm. In addition to the BOL and TOL legs, several other legs were flown within the wall, for example below and within the cloud, and within the aerosol layer. A full radiation wall often required over an hour to complete, and a square spiral was often included at one end of the radiation wall.
Since the radiation walls provide broad spatial coverage and allow us to sample varying cloud albedos and aerosol extinction profiles for the same aerosol plume, we use the radiation wall data to understand the impact of scene variability (i.e. aerosol 5 loading and cloud albedo variability) on the aerosol heating rates. We do that by calculating heating rate curtains along the linear flight path, which require spiral measurements to be used in conjunction with AOD and scene albedo measurements from radiation walls (Table 1).
The downside of the radiation wall sampling is that it does not provide measurements throughout the aerosol layer. The new square spiral maneuver, described in detail in Cochrane et al. (2019), allows for multiple measurements to be taken throughout 135 the vertical profile over a short time period (10 -20 minutes depending on the vertical extent of the aerosol layer). From the measurements in the column (the layer that encompasses both the aerosol and the water vapor layer) absorption can be derived, and ultimately heating rates segregated by absorber (primarily water vapor and aerosols, detailed in section 2.3).
Having the full column of measurements available from the spiral profiles led to the development of an aerosol optical property retrieval (SSA and the asymmetry parameter (g), as detailed in Cochrane et al., 2021) that is directly tied to the irradiance 140 measurements throughout the column. It is inherently difficult to retrieve these properties since the aerosol radiative effects can be relatively small compared to the horizontal variability of cloud albedo. The spiral sampling strategy, however, reduces cloud and aerosol inhomogeneity effects while maintaining correlation of measured irradiances throughout the spiral to the ambient cloud field (Cochrane et al., 2019).
From all useable spiral maneuvers from ORACLES 2016 and 2017 (Cochrane et al., 2021), we combine the retrieved aerosol 145 properties with the in situ measured water vapor content, 4STAR-derived spectral extinction (derived from the derivative of AOD along the vertical) and atmospheric gas retrievals (Table 2). From these inputs, we calculate the vertical heating rate profiles for these cases. Since they are constrained by the total column absorption measured by SSFR, the heating rates are directly tied to the observations, ensuring consistency between aerosol properties, the water vapor profile, and the combined radiative effects of the atmospheric constituents (radiative closure). Since this is done spectrally, we can segregate by major 150 absorbers to determine their relative contributions to the overall heating.

Heating rate calculations
To calculate heating rates for both the spiral profiles and radiation walls, we use the 1-dimensional (1D) radiative transfer model (RTM) DISORT 2.0 (Stamnes et al., 2000) with SBDART for atmospheric molecular absorption (Ricchiazzi et al., 1998) within the libRadtran library (Emde et al., 2016;libradtran. org). Table 3 lists the input parameters and their sources 155 required for the calculations. The general method is to segregate the absorption by strategically eliminating a single constituent for separate calculations (Kindel et al., 2011). First, we calculate the total heating rate from cloud top to well above the top of the aerosol layer, then remove a single atmospheric component from the profile (e.g., aerosols) and calculate the heating rate again. The difference between the two calculations provides the isolated heating rate for the individual component that was excluded in the second calculation. 160 The radiative transfer calculations are performed in this manner rather than calculating the heating rate of each component directly in order to maintain physical consistency. As incoming radiation travels through the atmosphere, less and less total radiation reaches lower altitudes. When there is a strong absorber present, such as an aerosol layer, the decrease of radiation with altitude is amplified since the layer absorbs a significant portion of the incoming radiation. Calculating absorption (and heating) for a single component directly does not consider the attenuation by other constituents, potentially leading to an 165 overestimation of the heating rates at altitudes below that of the absorbing layer.
We calculate the heating rate of a layer following the equation presented in Schmidt et al. (2021): where is the density, cp is the constant-pressure specific heat capacity of air, ∆ is the layer thickness, and Δ net is the difference of the net irradiance at the layer top and bottom, i.e., the absorbed irradiance in that layer. Since the absorbed 170 irradiance is a spectral quantity, the different wavelengths contribute varying amounts and must be integrated to find the total heating rate at each altitude. The heating rate is calculated at each layer altitude defined in the model (every 0.2 km for cloud top altitude (approximately 1 km) < altitude < 7km; every 1 km for 7 km < altitude < 12 km).
The RTM requires the spectral aerosol optical properties SSA and the asymmetry parameter (g) along with the vertical profile of aerosol extinction and the spectral albedo. SSA and g spectra are taken from the SSFR retrievals presented in Cochrane et 175 al. (2021) and assumed to be vertically homogenous since the retrieved values represent the entire layer. The spectral albedo measured by SSFR just above the cloud top defines the "surface" beneath the aerosol layer at the altitude of the cloud top.
Aerosol extinction is taken either from HSRL-2 (for radiation walls) or derived from the 4STAR AOD profile (for spirals) and detailed in section 2.3 and 2.4. Atmospheric gases are defined by the standard tropical atmosphere included in the libRadtran package (Anderson et al., 1986). 180 The water vapor and ozone profiles, however, are modified to reflect the specific conditions encountered during the time of the measurements. Specifically, due to the consistently high water vapor observed in the biomass burning plume (e.g., Pistone et al 2021), the vertical distribution of the water vapor within the aerosol layer is first determined by in situ measurements taken by the P-3 hygrometer during the spiral profiles. Beyond the top of the aerosol layer (TOL), the values revert to those found in the standard atmosphere. The full column of water vapor (the water vapor path) is then scaled such that the total 185 column (between BOL and top of atmosphere (TOA)) is equal to the 4STAR retrieval at the lowest flight altitude. One example profile is shown in Figure 2. Similarly, the ozone profile is scaled to the 4STAR ozone retrieval measured at the BOL. For the spiral cases, the water vapor and ozone profiles are consistent with those used within the aerosol retrievals of Cochrane et al.,

Heating Rate Segregation by Absorber 190
The vertical profiles measured during the spiral flight patterns allow us to separate heating rates for different constituents beyond that of just the aerosol. This separation process relies on the principle that different atmospheric gases, aerosols, and water vapor absorb radiation in different wavelength ranges.
To separate the total heating into component-specific heating rates, we strategically remove components one by one and rerun the radiative transfer calculations. For example, to determine the aerosol heating rate profile, we set the aerosol extinction 195 profile to zero and calculate the heating rate profile. The difference between the total heating rate (where all constituents are included: the aerosol layer, water vapor, and ozone (all from measurements) as well as the libRadtran standard carbon dioxide and oxygen) and the total heating rate spectrum with the aerosol layer removed is the aerosol heating rate profile. The same technique is applied to each other constituent. Figure 3 illustrates the different heating rate spectra between the different components at one example altitude. While minimal for our cases, very thick absorbing aerosol layers may induce shading 200 effects on the downwelling irradiance.
The aerosol extinction profile at each wavelength is obtained from the 4STAR measurements (derivative of AOD with respect to altitude), which provide the required spectral dependence throughout the profile. This is only possible when the profile of AOD measurements is available (i.e., from a spiral). This approach is consistent with the methodology used to retrieve the aerosol optical properties (SSA and g); Cochrane et al., 2019;Cochrane et al., 2021). Data conditioning is required such that 205 the AOD profile decreases monotonically with altitude to eliminate any unphysical (negative) extinction values. These extinction profiles are consistent with those used in the aerosol property retrieval and have a 20 m vertical resolution; one example profile is shown in Figure 2. The directly retrieved HSRL-2 aerosol extinction profile at 532 nm is close to the 4STAR derived aerosol extinction profile. Any differences are due to the different location where these were acquired (wall vs. spiral), or due to the data conditioning applied to the 4STAR data. 210 Table 2 presents the input values for the RTM heating rate calculations required input parameters for each of the spiral cases.

Heating rates from spiral profiles
Total, aerosol, and water vapor heating rate profiles separated by year are presented in Figures 4a, b, and c. For the selected cases analyzed here, the aerosol plume was lower in altitude in 2017 than 2016, which can be seen in Figures 4a (total heating rates) and 4b (aerosol heating rates). The peak aerosol heating rates were similar in both years, approximately 4-6 K day −1 , 215 with the exception of one 2017 case for which the aerosol loading was significantly higher than other cases.
In 2016, the aerosol and water vapor layer were generally collocated in altitude (Figures 4b and 4c), supporting the common assumption that aerosol and water vapor advect jointly from the continent. In 2017, by contrast, the primary aerosol loading resides below the free tropospheric water vapor loading. The maximum water vapor heating (excluding the boundary layer) was between 2-4 K day −1 at the specified solar zenith angles in both 2016 and 2017, though the 2017 water vapor heating rate 220 peak is much broader than in 2016. These differences are possibly due to the differing locations and meteorological conditions, examined further in a campaign-wide analysis by Pistone et al., 2021. Compared to Deaconu et al. (2019), who derived instantaneous heating rates from MODIS/POLDER/CALIOP satellite instrumentation and the ECMWF ERA-Interim product from June-August 2008, we generally find lower peak aerosol heating rates: 4-6 K day −1 versus 9 K day -1 . This is likely due to differences in aerosol loading and vertical structure, since Deaconu et 225 al. 2019 analyzed data much closer to the continent than any of our ORACLES cases. The instantaneous water vapor heating rate estimates, albeit at different solar zenith angles, are consistent: we calculate 2-4 K day −1 compared to 3 K day -1 .
The water vapor contributions reported in this work as well as in Deaconu et al. (2019) are significantly larger than those in Adebiyi et al. (2015), who find that the maximum of the instantaneous water vapor SW heating rate is only ~10% of the maximum aerosol heating rate for clear sky near the island of St. Helena averaged over September-October: approximately 230 0.12 K day −1 (water vapor) relative to 1.2 K day −1 (aerosol). By contrast, our averaged maxima show the water vapor heating rate to be ~60% of the aerosol heating rate: 2.8 K day −1 compared to 4.6 K day −1 . One reason is that the water vapor heating rates reported in Adebiyi et al. (2015) represent the difference between composite humidity profiles constructed from days with light and heavy aerosol loadings (AOD<0.1 and AOD>0.2 respectively), and rely on the positive correlation between aerosol and water vapor. 235 Comparing our results to other observations or model estimates previously reported in the literature (e.g., Marquardt Callow et al., 2020;Baró Pérez et al., 2021) is also challenging because heating rates depend on numerous parameters beyond spectral extinction profiles and single scattering albedo: most importantly on the albedo of the underlying scene, layer depth and vertical resolution, as well as on variations in location (aerosol type) and solar zenith angle. Also, some studies report column-averaged, other than peak heating rates. This challenge is one of the motivations for the development of the heating rate efficiency 240 (Section 5) which turns the heating rate (an extensive parameter proportional to the extinction) into an intensive quantity. If generally adopted, this new concept will allow for improved comparisons between studies.
Nevertheless, we first calculate the traditionally reported column-averaged heating rate (from cloud top to the top of the aerosol layer). Figure 5 shows the breakdown of the column-averaged total heating rate into its individual components, averaged for all cases. As expected, the largest contribution of heating stems from the aerosol (55.7% of the total), followed by that from 245 water vapor (37. 5%).

Heating Rate Dependence on Scene Parameters: Heating Rate Curtains
For the radiation walls, we calculate heating rate curtains to examine the relative dependencies of the aerosol heating rates on the AOD and the underlying albedo. This requires a) aerosol optical properties and water vapor profile (both obtained from a spiral retrieval and held constant), b) the SSFR-measured albedo spectra, the AOD spectra and column ozone retrievals from 250 4STAR, all measured from the BOL leg of the radiation wall, and c) aerosol extinction profiles at 532 nm of the full aerosol layer measured by HSRL-2 from the TOL leg (Table 2). Cases for which these criteria were not met were excluded from analysis. As mentioned above, radiation walls were typically accompanied by a spiral at one end. The spiral provides the aerosol optical properties along with the water vapor profile, while the collocated BOL and TOL legs provided the horizontal aerosol and cloud variability. Unfortunately, the spatial collocation is of little help since the time lag between sequential 255 sampling of the TOL and BOL legs is large enough that the cloud field is likely to change. In addition, some radiation walls included TOL and BOL layers that were not perfectly collocated and/or there was not an associated spiral that produced a valid aerosol retrieval.
Therefore, it was important to find a way to examine the heating rate dependence on the scene variability without knowing the albedo and AOD spectra underlying a specific aerosol extinction profile. We therefore assessed the dependence statistically 260 by considering 1) the albedo variability and 2) the AOD. For 1), we used the collection of measured spectra along the BOL leg, regardless of whether the BOL and TOL legs were collocated. From that, we created five representative albedo spectra for the heating rate calculations: the minimum, maximum, mean, and mean ± one standard deviation of the data collection. For part 2, we pursued a similar approach with the HSRL-2 extinction profiles from the TOL leg, which provide the vertical distribution of the aerosol. We paired this with the 4STAR-derived AOD from the bottom of the spiral (or from the leg-265 averaged BOL retrievals if a spiral was not available) as the basis for spectral extrapolation, which maintains spectral consistency across all heating rate calculations, as well as with the aerosol property retrievals. It should also be noted that for the two cases with no useable spiral, the aerosol optical properties were set to the ORACLES mean SSA and g spectra for the 2016 and 2017 deployments (Table 3 Cochrane et al., 2021) and the water vapor profile (and 4STAR column water vapor for scaling) are obtained from the closest spiral of the day regardless of whether there was a valid aerosol retrieval. 270 For each radiation wall, we calculate aerosol heating rates for extinction profiles from HSRL-2 along the TOL, approximately one profile every 5 seconds (interpolated from the 10 second resolution reported in the HSRL-2 data product) which is approximately every 0.005 degree of latitude or longitude. We extend the 532 nm extinction measured by HSRL-2 to the entire spectrum using the representative 4STAR AOD spectra in the following manner: 1. Integrate HSRL-2 extinction β ,-. to get AOD at 532 nm: where x is the location and z is the height.
In this approach, we are making two implicit assumptions: 1) that the 532 nm extinction can be used as a proxy to the full 285 extinction spectrum and 2) that the spectral shape of extinction does not vary with location or altitude. This may not be fully representative of the aerosol layer: For example, there could be different amounts of coarse mode aerosol in various layers of the plume, and the SSA is not necessarily constant with altitude (LeBlanc et al., 2020;Redemann et al., 2021;Dobracki et al., 2021). One possible method to refine assumption 2 would be to introduce an altitude-specific spectral extinction adjustment based on the HSRL-2 Angstrom exponent product, which could be investigated in further studies. 290 Figure 6 shows one example of the aerosol heating rate curtains calculated using the HSRL-2 extinction profiles. The high resolution of HSRL-2 allows us to resolve the minor variations within the aerosol plume while providing sufficient data to analyze the aerosol layer statistically. For each case, the albedo, AOD, and aerosol extinction profile vary across the wall, while the intensive aerosol optical properties are assumed to remain constant. For the case shown in Fig. 6, the thickest part of the aerosol layer occurs between 2km and 3km, with higher heating rates to the south. Although the AOD along the wall ranges 295 only from 0.19 to 0.22 at 532 nm, the high variability in the albedo causes large differences in the heating rate values.
The heating rate curtains provide enough data to examine the heating rate dependence on both the AOD and albedo along the wall for each of the cases. Each plot in Fig. 7 shows a thick solid profile (black) for the total heating rate and one for the aerosol heating rate (gray). These profiles represent the mean along the entire wall of HSRL-2 extinction profiles at the mean SSFR albedo value. The overlaid error bars represent the standard deviation (1-sigma) of the heating rate due to variable extinction 300 for the mean albedo value. The dashed lines represent the mean for all extinction profiles for the lowest albedo and the highest albedo encountered during the wall.
In some cases, such as 20160920 (7a) and 20160924 (7b), the variation in AOD results in larger heating variation than caused by the changing albedo below the layer. In other cases, 20170813 (7c) and 20170824 (7d) the variation in heating due to changing albedo along the wall is greater than that introduced by the AOD. In the 20170826 (7e) case, the variation in albedo 305 affects the total heating rate more than the AOD variation, while for the aerosol heating rate, the AOD and albedo variation introduce approximately equal variations.
For every case, the mean total (aerosol) heating rate ranges between 4-8 K day −1 (2-3 K day −1 ). The aerosol heating values are similar to those found for the region by Keil and Haywood (2003;2.3 K day −1 ), Gordon (2018; 1.9 K day −1 ), and Wilcox (2010; 2.0 K day −1 ). In contrast, we find slightly larger aerosol heating rate values than Tummon (2010; 1 K day −1 ) and Adebeyi et 310 al. (2015; 1.2 K day −1 ). Of course, we do not necessarily expect agreement with prior studies since there are numerous differences between them such as different observational periods, locations, approach, and aerosol optical properties.
A direct comparison (in both observation location and time) for the September 24 th , 2016 radiation wall (Fig 7, panel f) between the ALADIN-Climate (Aire Limitée Adaptation dynamique Développement InterNational) regional model (Nabat et al., 2015) and our direct calculations show a difference of approximately 1.5 K day −1 in peak heating for the average aerosol shortwave 315 heating rate profile, with the peak heating values from the ALADIN model slightly higher in altitude. The aerosol extinction profile from the ALADIN model is larger in magnitude than that of ORACLES. Likely the difference in extinction is the main driver in the differences between heating rates. The negative values between 7-10 km in the ALADIN profile may possibly come from differences in high-cloud properties between the two model simulations (Mallet et al., 2019).
Heating rate variability between cases, locations, and campaigns strongly depends on drivers such as extinction, cloud albedo, aerosol optical properties and the vertical distribution. Unlike the related quantity DARE, which is usually studied at a single level such as the top of atmosphere or the surface, heating rates extend throughout the entire column. In the past, it has been difficult to condense this into a single value to report and compare, leading to varying definitions and reported values in the literature (e.g., peak heating rate, column average heating rate.) 325 The main drivers of the heating rates are key to understanding the origin of the variability and overcoming the factors that limit our ability to understand the heating rate variability. For aerosols, the AOD is the most significant driver for the overall heating rate value. In Fig 8a, the column-averaged aerosol heating rate for each spiral case is shown as a function of the AOD at 550 nm, with some cases labeled by their 550 nm albedo value. As expected, aerosol loading, i.e., AOD, is the main driver of the heating rate, but a simple relationship is not expected because the vertical structure of the plume and its thickness vary from 330 case to case. Still, the dependence can be confirmed with radiative transfer calculations. We calculate aerosol heating rates as a function of AOD (black dashed line in Fig. 8a) for SSA=0.83, g=0.54, and an albedo of 0.60 (values listed for 550 nm, the full mean ORACLES spectra are used within the RTM). The vertical profile for these calculations was set to that of the 20160920 #2 spiral case.
For water vapor, the main controller of the heating rate is the water vapor path (the layer-integrated water vapor content), 335 shown in Fig 8b. Each case is labeled by the aerosol SSA value at 550 nm. The deviation of the individual cases from a linear relationship cannot be explained by the scene albedo at the water vapor absorption bands since the dependence (not shown) is too weak. The most plausible explanation is therefore the vertical distribution of the water vapor. From Figures 8a and 8b, we can see that layer heating by water vapor dominates over aerosol-induced heating up to a mid-visible aerosol optical thickness of about 0.25 (for a case with a water vapor path of ~1.2 g cm -2 , SSA=0. 83). 340 From Fig 8a, we can see that the variability in the AOD between cases overwhelmed any signal from the smaller heating rate drivers such as SSA and albedo. The problem is that the variability of AOD and its vertical distribution are intertwined. To isolate the vertical distribution of the aerosol extinction coefficient from its column integral, we divided the heating rate at any given altitude (z) by the aerosol extinction coefficient, rather than working with the column-averaged heating rate as in Fig. 8 and in previous studies. We named this "intensive" parameter the relative heating rate efficiency, borrowing from relative 345 aerosol forcing efficiency (Redemann et al., 2006), which is also independent of the aerosol loading. The HRE is defined as: where H6IJKJL is the aerosol heating rate profile, 67# is the aerosol extinction profile at 532 nm, and ↓ is broadband downwelling, spectrally integrated from 350nm -2100 nm. The normalization by the incident broadband irradiance at any given altitude was included in the definition to account for the 350 changing SZA and self-dimming (related to the AOD) of the layer as we step further down into it. Using HRE instead of HR provides an alternate way to examine the dependence of aerosol heating on SSA, Albedo, and SZA. Figure 9 shows an example of the vertical profile of the aerosol heating rate, extinction, and HRE, which clearly shows that the HRE for the aerosol two sub-layers (at 2 km and 4.5 km) is rather stable, and comparable between the two layers, despite the variable extinction and heating rate profiles. 355 To examine the dependencies of the heating rate efficiency on its drivers, we used the extinction profile, SSA retrieval, and measured cloud albedo for one spiral case (20160920 #2) as a reference and determined the HRE variation for changes in the SSA, cloud albedo, and SZA. The SSA spectrum changes relative to the reference case are determined through the co-albedo, where the reference case is scaled to range from 0.01 to 0.24 at 532 nm. The scaled co-albedo (1-SSA) spectra are then translated back into SSA for input into the RTM. The cloud albedo spectrum changes relative to the reference case are obtained 360 via cloud retrievals based on the reference albedo spectrum, and then varying the retrieved cloud optical thickness to modulate the spectral albedo (for more details, please see Appendix A.3.2 in Cochrane et al., 2021.) While the SSA and cloud albedo were varied in the same way for the nine cases, the extinction profiles, water vapor profile remained case-specific. Fig 10 shows the dependence of HRE on its main drivers. We can see that the HRE (in contrast to layer-averaged heating rate) varies little from case to case, and that its variability is tightly constrained by the albedo and SSA. It increases with increasing 365 cloud albedo, consistent with the finding of Adebiyi et al. (2015). The same finding would be true for any bright surface (not just a cloud), although more complicated spectral dependencies would have to be considered. The black dashed line indicates the mean fit line, which extends from albedo values 0 to 1. From the fit line, it is apparent that the HRE increases by approximately a factor of 2 across the full albedo range from 0 to 1. This can be understood intuitively, considering that the layer at an albedo of 1 is essentially illuminated "twice" -from the top and from the bottom. Of course, the heating of the layer 370 is not exactly doubled since the illumination from the bottom (the upwelling) is less than from the top (the downwelling) due to the partial attenuation of radiation.
Fig 10b shows that the HRE decreases with increasing SSA. This is expected, since an increasing SSA indicates less absorption (relative to scattering), leading to lower heating rates. The gray symbols (calculated for the widened SSA range for one case only) show that HRE goes to 0 as SSA goes to 1 (purely scattering, no absorption). Decreasing the SSA from 0.9 to 0.8 almost 375 doubles the HRE, which is also expected. It should be emphasized that the clear and easy-to-interpret dependence of HRE on the albedo and SSA could not have been achieved with the previous concepts of layer-mean or maximum heating rate.
However, the benefits of the HRE concept also have a limit in that the extinction, SSA and albedo all have a spectral dependence, which is important to consider when deriving an inherently broadband quantity such as layer heating. In addition, one needs to consider the dependence of the heating rate on the solar geometry. Fig 10c shows that HRE increases as a function 380 of SZA, with a sharp increase at the higher SZAs (lower sun elevations). This is due to the increase in path length through the aerosol layer as the sun moves towards to horizon.  Figure 10), the HRE is 0.025 K / day / km -1 / W m -2 with a standard deviation of only 0.002, approximately 8% in relative terms. This small 385 variability shows HRE could be used to translate extinction profiles in the region directly into aerosol heating rates if mid-visible cloud albedo and SSA are also known. In other words, the variability in extensive parameters (e.g., extinction) is higher than intensive parameters (e.g., SSA, g) and therefore, regionally and seasonally defined HRE are useful. If available for a specific region, the HRE concept would allow a direct translation from mid-visible extinction to heating rate. Of course, if SSA varies appreciably within the layer, that dependence may have to be made explicit. Alternatively, if in the future the 390 absorption coefficient were available at sufficient accuracy in addition to the extinction coefficient, the HRE could be redefined to normalize by the absorption coefficient, thereby accounting for the SSA vertical dependence.

Summary and Discussion
Observations from the 2016 and 2017 ORACLES experiments allow us to introduce a method of determining heating rate profiles as directly as possible by linking the heating rates to SSFR-measured irradiance profiles. The spectrally resolved 395 irradiance measurements and the derived column absorption allowed the separation of heating rates by absorber, most importantly the separation of water vapor heating from aerosol heating. We found that for many cases, the water vapor heating rate is nearly as large as the aerosol heating rate, on average 38% compared to 56% of the total heating (respectively), highlighting the large influence of the atmospheric water vapor distribution on the total heating rate distribution -even for optically thick aerosol layers. We also found that layer heating by water vapor coincident with aerosol dominates over aerosol-400 induced heating up to a mid-visible aerosol optical thickness of about 0.25 (for one case with a water vapor path of ~1.2 g cm -2 , SSA=0.83).
Analysing the dependence of the heating rate on the driving parameters (AOD, albedo, SSA) showed that the primary parameter affecting the aerosol heating rate is the AOD, and to a lesser extent the cloud albedo and aerosol SSA. For the mean SSA spectrum encountered during ORACLES (i.e., SSA=0.83 at 550 nm), the heating rate increases by ~0.5 K day -1 per 0.1 405 increase in aerosol optical depth (for albedo=0.6 at 550 nm). The heating rate variability, however, is highly case dependent because of the co-variability of the driving parameters.
The vertical distribution of the aerosol layer in relation to the underlying cloud also introduces variability from case to case that cannot easily be evaluated. The impact of the vertical distribution makes heating rates more difficult to generalize than other radiative effects, such as the direct aerosol radiative effect (DARE, Cochrane et al., 2020). The new HRE parameter, 410 however, does show potential for such a generalization because it defines the heating rate per extinction coefficient. The HRE, in contrast to the heating rate, has a clear dependence on albedo and SSA that could not have been determined through either maximum or layer-averaged heating rate concepts. By reducing the complexity of the convoluted relationship between heating rates and its drivers, the HRE parameter makes it possible to investigate the relationship between heating rates and other parameters besides the aerosol loading and extinction profile. In the future, the HRE relationships established in this work 415 could be formalized into a general parameterization applicable for the ORACLES study region. However, it will be important to consider that this type of broadband generalization will require spectral dependencies of extinction, SSA, and albedo.
For one preliminary comparison between our calculations and the ALADIN regional climate model output, we found consistent peak aerosol heating rates. The heating rate curtains, calculated along radiation walls using HSRL-2 extinction profiles, can currently only be derived from aircraft observations. To arrive at a similar product from satellite retrievals, one would have 420 to ensure that the heating rates are consistent with a radiative flux constraint (at the very least at TOA), in addition to filtering out cloud inhomogeneity effects. This will be relevant for planned space-borne missions.  Figure 2. 20160920 4STAR-derived 532 nm aerosol extinction profile (dark blue), averaged HSRL-2 532 nm aerosol extinction profile across the radiation wall (red dashed), and the water vapor content profile measured during the spiral profile (cyan). The 4STARderived extinction (532 nm) and the water vapor mixing ratio from the spiral profile is at a different location relative to the wallaveraged HSRL-2 extinction profile.    Table 3. The heating rate contribution of the category labeled "Other" is primarily gas (ozone, oxygen, and carbon dioxide) Figure 6. Heating rate curtains calculated using HSRL-2 measured extinction for the 20170813 radiation wall (shown here in separate plots: North (right) and South (left)). Peak heating of ~4.8 K/day occurs between 2 and 3 km. The underlying albedo, shown at 532 nm in green, is significantly higher on the left plot than the right (i.e., further south), contributing to higher aerosol heating rates. The AOD, shown at 532 nm in blue, does not vary significantly across the wall. Missing results between -7.08 and -6.51 are due to in-cloud sampling that replaced above-cloud albedo measurements and serve as the break point between North and South ends of the wall.  model for which these calculations are performed is set to every 0.05 km for altitude < 7km; every 1 km for 7 km < altitude < 12 km.

Author contributions
SPC collected SSFR data, performed the bulk of the analysis, and wrote the majority of the paper with input from the other authors. KSS collected SSFR data, helped with the methodology development and data analysis, and helped with developing, 650 writing, and editing the paper. HC, PP, and SK helped with the data collection of SSFR. JR was one of the PIs for the ORACLES campaign and provided 4STAR data. SL was the PI of the 4STAR instrument. KP, MK, MSR, YS, and CF provided 4STAR data. RF, SB, and CH provided HSRL data. MM provided ALADIN model output. PZ was on the leadership team for the ORACLES project and helped advise the interpretation of coincident aerosol and water vapor heating rates. All the coauthors helped in the reviewing and editing of the paper.