The impact of cloudiness and cloud type on the atmospheric heating rate of black and brown carbon in the Po Valley

We experimentally quantified the impact of cloud fraction and cloud type on the heating rate (HR) of black and brown carbon (HR_BC and HR_BrC). In particular, we examined in more detail the cloud effect on the HR detected in a previous study (Ferrero et al., 2018). High-time-resolution measurements of the aerosol absorption coefficient at multiple wavelengths were coupled with spectral measurements of the direct, diffuse and surface reflected irradiance and with lidar–ceilometer data during a field campaign in Milan, Po Valley (Italy). The experimental set-up allowed for a direct determination of the total HR (and its speciation: HR_BC and HR_BrC) in all-sky conditions (from clear-sky conditions to cloudy). The highest total HR values were found in the middle of winter (1.43 ± 0.05 K d⁻¹), and the lowest were in spring (0.54 ± 0.02 K d⁻¹). Overall, the HR_BrC accounted for 13.7 ± 0.2 % of the total HR, with the BrC being characterized by an absorption Ångström exponent (AAE) of 3.49 ± 0.01. To investigate the role of clouds, sky conditions were classified in terms of cloudiness (fraction of the sky covered by clouds: oktas) and cloud type (stratus, St; cumulus, Cu; stratocumulus, Sc; altostratus, As; altocumulus, Ac; cirrus, Ci; and cirrocumulus–cirrostratus, Cc–Cs). During the campaign, clear-sky conditions were present 23 % of the time, with the remaining time (77 %) being characterized by cloudy conditions. The average cloudiness was 3.58 ± 0.04 oktas (highest in February at 4.56 ± 0.07 oktas and lowest in November at 2.91 ± 0.06 oktas). St clouds were mostly responsible for overcast conditions (7–8 oktas, frequency of 87 % and 96 %); Sc clouds dominated the intermediate cloudiness conditions (5–6 oktas, frequency of 47 % and 66 %); and the transition from Cc–Cs to Sc determined moderate cloudiness (3–4 oktas); finally, low cloudiness (1–2 oktas) was mostly dominated by Ci and Cu (frequency of 59 % and 40 %, respectively).

HR measurements showed a constant decrease with increasing cloudiness of the atmosphere, enabling us to quantify for the first time the bias (in %) of the aerosol HR introduced by the simplified assumption of clear-sky conditions in radiative-transfer model calculations. Our results showed that the HR of light-absorbing aerosol was ∼ 20 %–30 % lower in low cloudiness (1–2 oktas) and up to 80 % lower in completely overcast conditions (i.e. 7–8 oktas) compared to clear-sky ones. This means that, in the simplified assumption of clear-sky conditions, the HR of light-absorbing aerosol can be largely overestimated (by 50 % in low cloudiness, 1–2 oktas, and up to 500 % in completely overcast conditions, 7–8 oktas).

The impact of different cloud types on the HR was also investigated. Cirrus clouds were found to have a modest impact, decreasing the HR_BC and HR_BrC by −5 % at most. Cumulus clouds decreased the HR_BC and HR_BrC by −31 ± 12 % and −26 ± 7 %, respectively; cirrocumulus–cirrostratus clouds decreased the HR_BC and HR_BrC by −60 ± 8 % and −54 ± 4 %, which was comparable to the impact of altocumulus (−60 ± 6 % and −46 ± 4 %). A higher impact on the HR_BC and HR_BrC suppression was found for stratocumulus (−63 ± 6 % and −58 ± 4 %, respectively) and altostratus (−78 ± 5 % and −73 ± 4 %, respectively). The highest impact was associated with stratus, suppressing the HR_BC and HR_BrC by −85 ± 5 % and −83 ± 3 %, respectively. The presence of clouds caused a decrease of both the HR_BC and HR_BrC (normalized to the absorption coefficient of the respective species) of −11.8 ± 1.2 % and −12.6 ± 1.4 % per okta. This study highlights the need to take into account the role of both cloudiness and different cloud types when estimating the HR caused by both BC and BrC and in turn decrease the uncertainties associated with the quantification of their impact on the climate.

1 Introduction

The impact of aerosols on the climate is traditionally investigated with a focus on their direct, indirect and semi-direct effects (Bond et al., 2013; IPCC, 2013; Ferrero et al., 2018, 2014; Ramanathan and Feng, 2009; Koren et al., 2008, 2004; Kaufman et al., 2002). Direct effects are related to the sunlight interaction with aerosols through absorption and scattering; indirect effects are related to the ability of aerosol to act as cloud condensation nuclei affecting the cloud formation and properties; and semi-direct effects are those related to a feedback on cloud evolution affecting other atmospheric parameters (e.g. the thermal structure of the atmosphere) (IPCC, 2013; Ramanathan and Feng, 2009; Koren et al., 2008, 2004; Kaufman et al., 2002). Both direct and indirect radiative effects of anthropogenic and natural aerosols are still the major sources of uncertainties on climate (IPCC, 2013). Recent studies show, for example, that the aerosol direct radiative effect (on a global scale) may switch from positive to negative forcing on short (e.g. daily) timescales (Lolli et al., 2018; Tosca et al., 2017; Campbell et al., 2016). This is due to the fact that aerosol is a heterogeneous complex mixture of particles characterized by different size, chemistry and shape (e.g. Costabile et al., 2013), greatly varying in time and space both in the horizontal and vertical dimension (e.g. Ferrero et al., 2012). On a global scale, most of the values reported for the aerosol direct radiative effect were derived from models (Bond et al., 2013; Koch and Del Genio, 2010). This has the advantage of providing fields of continuous direct radiative effect in space and time. However, inaccuracies related to simplified model assumptions on chemistry, shape and the mixing state of particles can affect the results (Nordmann et al., 2014; Koch et al., 2009); this amplifies the uncertainties on the related global and regional aerosol effects on the climate (Andreae and Ramanathan, 2013). The aerosol direct radiative effect has been usually determined in clear-sky conditions both in model simulations and measurements. The clear-sky approximation is useful when comparing measurements to radiative-transfer modelling outcomes during experimental campaigns performed in fair-weather conditions (e.g. Ferrero et al., 2014; Ramana et al., 2007); however, in general this simplification cannot capture the complexity of the phenomenon in the majority of weather conditions (Myhre et al., 2013). In fact, clouds are one of the most important factors influencing the solar radiation reaching the ground. By scattering and absorbing the radiation, clouds can affect the radiation magnitude and modify its spectrum, especially in the ultraviolet (UV) region (López et al., 2009; Thiel et al., 2008; Calbó et al., 2005). During cloudy conditions the global irradiance is usually reduced; however, the presence of clouds sometimes results in short-term enhancement of global irradiance (Duchon and O'Malley, 1999). In some specific cases, the scattering of radiation from the sides of the cloud may enhance global irradiance in the UV to the levels higher than those in clear-sky conditions (Mims and Frederick, 1994; Feister et al., 2015). Mims and Frederick (1994) determined that the scattering from the sides of cumulus clouds can enhance the total (global) UV-B solar irradiance by 20 % or more over the maximum solar-noon value when cumulus clouds were close to (but not when blocking) the solar disk. In a similar way, Feister et al. (2015) concluded that the scattering of solar radiation by clouds can enhance UV irradiance at the surface – for example, cumulonimbus clouds, with top heights close to the tropical tropopause layer, have the potential to significantly enhance diffuse UV-B radiance over its clear-sky value. UV radiation also interacts with aerosols and particularly with those featuring significant absorption values in this spectral region. UV represents an important region for brown carbon (BrC) absorption with respect to other light-absorbing aerosol (LAA) components (e.g. black carbon, BC). Thus, the presence of clouds could influence the impact of different LAA species on the climate in a different way.

Up to now, the role of cloudiness and cloud type on the aerosol direct radiative effect was poorly investigated. Matus et al. (2015) recently used a complex combination of the CloudSat's satellite multi-sensor radiative flux and heating-rate (HR) products to infer both the direct radiative effect at the top of the atmosphere and HR profiles of aerosols that lie above the clouds. The study showed how results were affected by the cloudiness (e.g. cloud fraction) and, for the southeastern Atlantic, reported a direct radiative effect ranging from −3.1 to −0.6 W m⁻², going from clear-sky to cloudy conditions.

A further investigation by Myhre et al. (2013) reported results of modelling simulations during the AeroCom (Phase II) project: in all-sky conditions (thus including the effect of clouds) they estimated an all-sky direct radiative effect of −0.27 W m⁻² (range of −0.58 to −0.02 W m⁻²) for total anthropogenic aerosols, with this being about half of the clear-sky one. The most important factors responsible for the observed difference were the amount of aerosol absorption and the location of aerosol layers in relation to clouds (above or below). In fact, the presence of LAA (mainly BC, BrC and mineral dust) might have important effects on the radiative balance. It is estimated that, due to its absorption of sunlight, BC is the second most important positive anthropogenic climate forcer after CO₂ (Bond et al., 2013; Ramanathan and Carmichael, 2008); BrC contributes ∼ 10 %–30 % to the total absorption on a global scale (Ferrero et al., 2018; Kumar et al., 2018; Shamjad et al., 2015; Chung et al., 2012). As a main difference compared to CO₂, LAA species are short-lived climate forcers, thus representing a potential global warming mitigation target. However, the real potential benefit of any mitigation strategy should also be based on observational measurements, possibly carried out in all-sky conditions.

It also noteworthy that the HR induced by LAA can trigger different atmospheric feedbacks. BC and mineral dust can alter the atmospheric thermal structure, thus affecting the atmospheric stability, the cloud distribution and even the synoptic winds such as the monsoons (IPCC, 2013; Bond et al., 2013; Ramanathan and Feng, 2009; Koch et al., 2009; Ramanathan and Carmichael, 2008; Koren et al., 2008, 2004; Kaufman et al., 2002). These feedbacks should be quantified on the basis of HR measurements in all-sky conditions. In agreement with this, both Andreae and Ramanathan (2013) and Chung et al. (2012) called for model-independent, observation-based determination of the absorptive direct radiative effect of aerosols. Since cloudiness and cloud type change on short timescales similarly to aerosols, long-term, highly time-resolved measurements (covering different sky conditions) are necessary to unravel the impact of LAA on the HR.

Satellite-based studies investigated the role of cloudiness and cloud type on the HR of aerosol layers above clouds (Matus et al., 2015). To our knowledge, there has been no experimental investigation of cloudiness and cloud type impact on the HR of aerosol layers below clouds, where most of the aerosol pollution typically resides. Cloud–aerosol feedbacks can strongly depend on the HR magnitude in cloudy conditions. As a matter of fact, the atmospheric heating induced by absorbing aerosol is traditionally related to a decrease of atmospheric relative humidity and less cloud cover (semi-direct effect). This effect can further increase the amount of the incoming solar radiation that reaches Earth's surface (and any close-to-surface LAA layers), leading to a positive feedback characterized by additional warming and a further decrease in the cloud amount (e.g. Koren et al., 2004). However, Perlwitz and Miller (2010) reported a counterintuitive feedback: the atmospheric heating induced by tropospheric absorbing aerosol could lead to a cloud cover increase (especially low-level clouds) due to a delicate interplay between relative humidity and temperature. The study concluded that high absorption by aerosols was responsible for two counteracting processes: a large diabatic heating of the atmospheric column (thus decreasing relative humidity) and a corresponding increase in the specific humidity able to exceed the temperature effect on relative humidity, with the net result of increasing low cloud cover with increasing aerosol absorption. This is an important result that underlines the importance of measuring the atmospheric HR in cloudy conditions as a constraint and/or input for more comprehensive climate models to shed light on the sign and magnitude of the related feedbacks on cloud dynamics.

This study attempts to experimentally measure for the first time the impact of different cloudiness levels and cloud types on the HR exerted by near-surface LAA species. The study was performed in Milan (Italy), located in the middle of the Po Valley (Sect. 2), which is an air pollution hotspot in Europe; its meteorological conditions are similar to those of a multitude of basin valleys (surrounded by hills or mountains) in which low wind speeds and stable atmospheric conditions promote the accumulation of aerosol (Zotter et al., 2017; Moroni et al., 2013, 2012; Ferrero et al., 2013, 2011a; Barnaba et al., 2010; Carbone et al., 2010; Rodriguez et al., 2007). Cloud presence cannot be neglected over the investigated area considering that, in the last 50 years, the annual average cloudiness, expressed in oktas, was estimated to be ∼ 5.5 over Europe (Stjern et al., 2009) and ∼ 4 over Italy (Maugeri et al., 2001). This feature is similar to 80 years of data of cloud cover in the United States (Crocke et al., 1999). To determine the LAA HR, we used a methodology previously developed in Ferrero et al. (2018) and further extended here to explore the effects of cloudiness and different cloud types on the HR of BC and BrC. More specifically, this work introduces the following novelties: (1) it describes the interaction between cloudiness and light-absorbing aerosol, presenting the aerosol HR as a function of cloudiness, and in turn estimates the systematic bias introduced by incorrectly assuming clear-sky conditions in radiative-transfer models; (2) it introduces a cloud type classification and investigates the impact of both cloudiness and cloud types on the total HR; and (3) it separates BC and BrC contributions and investigates their relative impact on the total HR as a function of sky conditions. The results presented in this study add an important piece of information in the general context of cloud–aerosol interactions and their influence on the HR.

Figure 1 (a) Location of the Milan sampling site in the Po Valley, Italy; (b) the U9 sampling site on the rooftop (10 $\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$) of the University of Milano-Bicocca. The copyright holder of Fig. 1 is Google Maps (© Google Maps).

2 Methods

Aerosol, cloud and spectral irradiance were measured in Milan (Italy) on the rooftop (10 $\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$) of the U9 building of the University of Milano-Bicocca (45^∘ 30^′38^′′ N, 9^∘12^′42^′′ E; Italy; Fig. 1). The site is located in the midst of the Po Valley, one of the most industrialized and heavily populated areas in Europe. In the Po Valley, stable atmospheric conditions often occur, causing a marked seasonal variation of aerosol concentrations within the mixing layer (Barnaba et al., 2010), well visible even from satellites (Ferrero et al., 2019; Di Nicolantonio et al., 2007, 2009; Barnaba and Gobbi 2004). A full description of the aerosol behaviour in Milan, at the University of Milano-Bicocca, and of the related properties (vertical profiles, chemistry, hygroscopicity, sources and toxicity) is reported in previous studies (Diemoz et al., 2019a; Lorelei et al., 2019; D'Angelo et al., 2016; Curci et al., 2015; Ferrero et al., 2015, 2010; Sangiorgi et al., 2011, 2014; Sandrini et al., 2014). In the framework of the present work it is important to underline that the U9 experimental site is particularly well suited for atmospheric radiation measurements: it is characterized by a full hemispherical sky view and equipped with the instruments described in Sect. 2.1. The measurement set-up allowed for the experimental determination of the instantaneous aerosol HR (K d⁻¹) induced by absorbing aerosol as detailed in Sect. 2.2. The methodological approach used to quantify the cloud fraction and to classify the cloud type is reported in Sect. 2.3. Finally, Appendix A resumes the nomenclature used in the present work.

2.1 Instruments

The aerosol, cloud and radiation instrumentations have been installed at the U9 sampling site in Milan since 2015. The site location is shown in Fig. 1. The complete instrumental set-up (Fig. S1 in the Supplement) is described hereafter.

2.1.1 Light-absorbing aerosol measurements and apportionment

Measurements of the wavelength-dependent aerosol absorption coefficient b_abs(λ) were obtained using a Magee Scientific AE-31 aethalometer. This allowed for multi-spectral measurements (7-λ: 370, 470, 520, 590, 660, 880 and 950 nm) in the wide UV–VIS–NIR (ultraviolet–visible–near-infrared) region, not available from other instruments (e.g. multi-angle absorption photometer, MAAP; particle soot absorption photometer, PSAP; and photoacoustic) (Virkkula et al., 2010; Petzold et al., 2005). This spectral range is needed for the HR determination (Sect. 2.2). The use of aethalometers also presents the advantage of global long-term data series (Ferrero et al., 2016; Eleftheriadis et al., 2009; Collaud Coen et al., 2010; Junker et al., 2006) that could allow for deriving historical data of the HR in the future.

To account for both the multiple-scattering enhancement (the elongation of the optical path induced by the filter fibres) and the loading effects (the non-linear optical path reduction induced by absorbing particles accumulating in the filter), the AE-31 data were corrected by applying the Weingartner et al. (2003) procedure (Ferrero et al., 2018, 2014, 2011b; Collaud Coen et al. 2010). As detailed by Collaud Coen et al. (2010), the Weingartner et al. (2003) procedure compensates for all the aethalometer artefacts (the backscattering is indirectly included within the multiple-scattering correction), showing a good robustness (negative values are not generated, and the results are in good agreement with other filter photometers), and, most importantly, it does not affect the derived aerosol absorption Ångström exponent (AAE) (fundamental for HR determination, Sect. 2.2). Overall, the multiple-scattering parameter C was 3.24 ± 0.03, as obtained by comparing the AE-31 data at 660 nm with an MAAP at a very similar wavelength (637 nm, Müller et al., 2011) (regression between AE-31 and MAAP in Fig. S2 in the Supplement). This value lies very close to that suggested by the Global Atmosphere Watch (GAW) programme (WMO/GAW, 2016), i.e. C = 3.5. The physical meaning of the similarity between the obtained C value (3.24 ± 0.03) and the GAW one implies that Milan (in the middle of the Po Valley) is characterized by continental-type aerosols (e.g. Carbone et al., 2010) and consistent with the global average. To verify the reliability of the obtained C value, it was also computed following the Collaud Coen et al. (2010) procedure. They defined the reference value of C (C_ref = 2.81 ± 0.11) for the AE-31 tape based on data from pristine environments (Jungfraujoch and Hohenpeissenberg sites, where aerosol has a single-scattering albedo of ∼ 1); at the same time, Collaud Coen et al. (2010) defined C for any type of aerosol as follows:

$$\begin{array}{}\text{(1)}& C=C_{\text{ref}}+\mathit{\alpha }{\displaystyle \frac{{\mathit{\omega }}_{\mathrm{0}}}{\mathrm{1}-{\mathit{\omega }}_{\mathrm{0}}}},\end{array}$$

where α is the parameter for the Arnott et al. (2005) scattering correction (0.0713 at 660 nm) and ω₀ the single-scattering albedo. In wintertime in Milan, within the mixing layer, the single-scattering albedo was found to be 0.846 ± 0.011 at 675 nm by Ferrero et al. (2014). From Eq. (1), it follows that the computed C in Milan is 3.20 ± 0.35, in keeping with the experimental one (3.24 ± 0.03). Details concerning wavelength differences are discussed in the Supplement (“Measured and computed C factor”). The loading effects were dynamically determined following the Sandradewi et al. (2008b) approach, while the final equivalent BC (eBC) concentrations were obtained applying the AE-31 apparent mass attenuation cross section (16.6 m² g⁻¹ at 880 nm).

The abovementioned compensation procedures introduce an uncertainty in the absorption coefficient measurements. Collaud Coen et al. (2010) tested these procedures in different locations and estimated the global accuracy of the Weingartner et al. (2003) correction (applied in the present work) to be ∼ 23 %. Moreover, Drinovec et al. (2015) showed a good agreement between AE-31 aethalometer data (corrected using Weingartner et al., 2003) and those of the new version, AE-33, with a slope close to unity and R² > 0.90. Thus, the Collaud Coen et al. (2010) accuracy estimation is considered as the worst scenario.

As the spectral signature of the aerosol absorption coefficient b_abs(λ) reflects the different nature of absorbing aerosol (BC and BrC), once b_abs(λ) is obtained, it can be apportioned to determine the contributions of BC and BrC, respectively. This result can be achieved considering that BC aerosol absorption is characterized by an absorption Ångström exponent, AAE ≈ 1 (Massabò et al., 2015; Sandradewi et al., 2008a; Bond and Bengstrom, 2006). Conversely, BrC absorption is spectrally more variable, with an AAE from 3 to 10 (Ferrero et al., 2018; Shamjad et al., 2015; Massabò et al., 2015; Srinivas and Sarin, 2013; Yang et al., 2009; Kirchstetter et al., 2004). The wavelength dependence of the absorption coefficient of BrC can be described by the simple harmonic oscillator reported in Moosmüller et al. (2011): the much lower absorption in the IR (infrared) region (compared to UV) is a consequence of the resonances in the UV from which the IR region is far removed. This calculation also yields to decreasing AAE values with increasing wavelengths. This is equivalent to the band-gap model with the Urbach tail as detailed in Sun et al. (2007) and references in Moosmüller et al. (2011), where the key factor is the difference between the highest occupied and lowest unoccupied energy state of the molecules included in the BrC ensemble. In this study we determined AAE_BrC following the innovative apportionment method proposed by Massabò et al. (2015). This allows for apportioning b_abs(λ) to BC and BrC and for determining, at the same time, the AAE_BrC assuming that all BrC results from biomass burning. The method by Massabò et al. (2015) was previously applied to the Milan U9 measurements leading to an annual average of AAE_BrC = 3.66 ± 0.03 (Ferrero et al., 2018).

2.1.2 Radiative, meteorological and lidar measurements

Spectral irradiance measurements were collected using a multiplexer–radiometer–irradiometer (MRI; Fig. S1; details in Cogliati et al., 2015) which resolves the UV–VIS–NIR spectrum (350–1000 nm) in 3648 spectral bands (3648-element linear CCD array detector; charge-coupled device; Toshiba TCD1304AP, Japan) for both the downwelling and the upwelling radiation fluxes. The instrument was developed at the University of Milano-Bicocca using an optical switch (MPM-2000-2x8-VIS, Ocean Optics Inc., USA) to sequentially select between different input fibres fixed to up- and down-facing entrance fore-optics. The configuration used in the present work connects each spectrometer to three input ports: (1) the CC-3 cosine-corrected irradiance probes to collect the downwelling irradiance, (2) the bare fibre optics with a 25^∘ field of view to measure the upwelling radiance from the terrestrial surface and (3) the blind port that is used to record the instrument dark current. A 5 m long optical fibre with a bundle core with a diameter of 1 mm is used to connect the entrance fore-optics to the multiplexer input, while the connection between the multiplexer output ports and the spectrometers is obtained with 0.3 m long optical fibres. The set-up allows for sequentially measuring dark current and both up- and downwelling spectra simultaneously with the two spectrometers. The two spectrometers used are high-resolution HR4000 holographic grating spectrometers (Ocean Optics Inc., USA). Finally, the multiplexer–radiometer–irradiometer was equipped with a rotating shadow band to measure separately the spectra of the direct, diffuse and reflected irradiance (F_dir(λ), F_dif(λ), F_ref(λ)). The reflected irradiance originated from a Lambertian concrete surface (due to its flat and homogeneous characteristics which represents the average spectral reflectance of the Milan urban area well; Ferrero et al, 2018).

Broadband (300–3000 nm) downwelling (global and diffuse) and upwelling (reflected) irradiance measurements were also collected using LSI Lastem radiometers (DPA154 and C201R, class 1, ISO 9060, 3 % accuracy). Diffuse broadband irradiance was measured using the DPA154 global radiometer equipped with a shadow band whose effect was corrected (Ferrero et al., 2018) to determine the true amount of both diffuse and direct (obtained after subtraction from the global) irradiance. Next, MRI spectra were normalized and completed with normalized literature spectra (Ferrero et al., 2018) to cover the broadband range (300–3000 nm) and irradiance intensity measured by standard LSI Lastem pyranometers, allowing for the HR to be evaluated over the whole short-wave range (b_abs(λ) was estimated outside the AE-31 range using its AAE). The approach was previously validated (Ferrero et al., 2018): the HR in the strict UV–VIS–NIR range (350–950 nm of the AE31 and the MRI) accounted on average for 86.4 ± 0.4 % of the total broadband values.

In addition to radiation measurements, temperature, relative humidity, pressure and wind parameters were measured using the following LSI Lastem sensors: DMA580 and DMA570 for thermo-hygrometric measurements (for T and RH: range of −30– +70 ^∘C and 10 %–98 %, accuracy of ± 0.1 ^∘C, and ± 2.5 % sensibility of 0.025 ^∘C and 0.2 %), the CX110P barometer model for pressure (range of 800–1100 hPa and accuracy of 1 hPa) and the CombiSD anemometer (range of 0–60 m s⁻¹ and 0–360^∘) for wind measurements.

The experimental station U9 is also equipped with an automatic lidar–ceilometer operated by ISAC-CNR in the framework of the Italian Automated LIdar-CEilometer network (ALICEnet, http://www.alice-net.eu, last access: 22 March 2021) and contributes to the EUMETNET (European meteorology network) E-Profile network (https://www.eumetnet.eu/, last access: 22 March 2021). This is a Jenoptik Nimbus 15k biaxial lidar–ceilometer operating 24 h d⁻¹, 7 d per week. It is equipped with an Nd:YAG (neodymium-doped yttrium aluminium garnet) laser that emits light pulses at 1064 nm with an energy of 8 µJ per pulse and a repetition rate of 5 kHz. The backscattered light is detected by an avalanche photodiode in the photon-counting mode (Wiegner and Geiß, 2012; Madonna et al., 2015). The vertical and temporal resolution of the raw signals are 15 m and 30 s, respectively. In order to improve the signal-to-noise ratio of the backscatter signal, the signal is processed with temporal averages of 2 min. The full overlap is obtained at an altitude of some hundred metres above the observation site, and overlap correction functions are applied in the first layers. The Nimbus 15k lidar–ceilometer is able to determine cloud base height (CBH), penetration depth, and with specific processing mixing layer height and vertical profiles of aerosol optical and physical properties (e.g. Diemoz et al., 2019a, b; Dionisi et al., 2018). We used the U9 ceilometer data for cloud layering and relevant cloud base height, as the system can reliably detect multiple cloud layers and cirrus clouds (Wiegner et al., 2014; Boers et al., 2010; Martucci et al., 2010) within its operating vertical range (up to 15 km). Given the vertical resolution of the instrument, expected uncertainty of the cloud base height derived by the lidar–ceilometer is less than ± 30 m.

Global and diffuse irradiance measurements, coupled with the ceilometer data, were used to determine the sky cloud fraction and to classify the cloud types by following the methodology presented in the Sect. 2.3.

2.2 Heating-rate measurements

The instantaneous aerosol HR (K d⁻¹) induced by LAA is experimentally obtained using the methodology reported and validated in Ferrero et al. (2018), where the reader is referred to for the details of the approach. Here we briefly summarize the method.

The heating rate is determined from the air density (ρ, kg m⁻³); the isobaric specific heat of dry air (C_p, 1005 $\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$); and the radiative power absorbed by aerosol per unit volume of air (W m⁻³), which describes the interaction between the radiation (either direct from the sun, diffused by atmosphere and clouds, and reflected from the ground) and the LAA (BC and BrC in Milan). The HR is determined as follows (Ferrero et al., 2018):

$$\begin{array}{}\text{(2)}& \begin{array}{rl}\text{HR}=& {\displaystyle \frac{\mathrm{1}}{\mathit{\rho }{C}_p}}\cdot \sum _{\text{dir,dif,ref}}{\int }_{\mathit{\theta }=\mathrm{0}}^{\mathit{\theta }=\mathit{\pi }/\mathrm{2}}{\int }_{\mathit{\lambda }=\mathrm{300}}^{\mathit{\lambda }=\mathrm{3000}}{\displaystyle \frac{F_{\text{dir,dif,ref}}(\mathit{\lambda },\mathit{\theta })}{\text{cos}\left(\mathit{\theta }\right)}}\\ & \cdot b_{\text{abs}}\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }\mathrm{d}\mathit{\theta },\end{array}\end{array}$$

where the subscripts dir, dif and ref refer to the direct, diffuse and reflected components of the spectral irradiance F of wavelength λ impinging on LAA with a zenith angle θ (from any azimuth).

Under the isotropic and Lambertian assumptions (as used in Ferrero et al., 2018), Eq. (2) can be solved, becoming

$$\begin{array}{}\text{(3)}& \begin{array}{rl}\text{HR}=& {\text{HR}}_{\text{dir}}+{\text{HR}}_{\text{dif}}+{\text{HR}}_{\text{ref}}\\ =& {\displaystyle \frac{\mathrm{1}}{\mathit{\rho }{C}_p}}\cdot \left[{\displaystyle \frac{\mathrm{1}}{\text{cos}\left({\mathit{\theta }}_z\right)}}\underset{\mathit{\lambda }}{\int }{F}_{\text{dir}}\right(\mathit{\lambda }\left)b_{\text{abs}}\right(\mathit{\lambda })\mathrm{d}\mathit{\lambda }\\ & +\mathrm{2}\underset{\mathit{\lambda }}{\int }{F}_{\text{dif}}\left(\mathit{\lambda }\right)b_{\text{abs}}\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }+\mathrm{2}\underset{\mathit{\lambda }}{\int }{F}_{\text{ref}}\left(\mathit{\lambda }\right)b_{\text{abs}}\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }],\end{array}\end{array}$$

where θ_z refers to the solar zenith angle, while F_dir(λ), F_dif(λ) and F_ref(λ) are the spectral direct, diffuse and reflected irradiances. Equations (2) and (3) are related to the concept of actinic flux (Tian et al., 2020; Gao et al., 2008; Liou, 2007); an extended description, as well as its demonstration, is detailed in the Supplement.

As the intensity of the irradiance components is a function of cloudiness and cloud type (Sect. 2.3), Eq. (3) enables assessing the impact of the latter components on the aerosol absorption of short-wave radiation and thus on the corresponding HR (Sects. 3.2 and 3.3).

The most important advantages and limitations of this measurement-based approach to derive the LAA HR are as follows. The advantages are as follows:

• no radiative-transfer assumptions needed (i.e. no assumption of clear-sky conditions), as the parameters input to Eq. (3) are all measured quantities;

• possibility to follow the rapid HR dynamic to investigate the HR temporal evolution, as measurements of spectral irradiance and absorption coefficient are carried out with high temporal resolution; and

• possibility to derive the HR in all-sky conditions, as measurements of spectral irradiance and the absorption coefficient are independent from atmospheric conditions enabling us to investigate the impact induced by the clouds.

The limitation is as follows:

• The HR is independent of the thickness of the investigated atmospheric layer and refers to the vertical location of the atmospheric layer in which it is experimentally determined. In the present work the HR was determined into the near-surface atmospheric layer.

With respect to this limitation, it should be mentioned that BC and HR vertical profile data previously collected at the same site and in other valley basins revealed that the HR was constant inside the mixing layer (Ferrero et al., 2014). In fact, above our observational site, vertical profile measurements with a tethered balloon and a lidar–ceilometer have been performed since 2005, mostly showing homogeneous concentrations of aerosol (and related extinction coefficient) within the mixing layer, particularly in daytime (Ferrero et al., 2019). The same condition was verified by the lidar–ceilometer data collected during the present campaign (Fig. S3 in the Supplement). The methodology is therefore believed to be also representative for the whole mixing layer if the aerosol vertical dispersion is homogeneous within this layer. This might not be the case for other regions, where the upper troposphere is impacted by high levels of BrC from biomass burning (Zhang et al., 2020), but Ferrero et al. (2019) showed that in Milan 87.0 % of aerosol optical depth signal was built up within the mixing layer, with 8.2 % being in the residual layer and 4.9 % being in the free troposphere.

2.3 Cloudiness and cloud classification

2.3.1 Cloudiness

The cloudiness was determined following the approach reported in Ehnberg and Bollen (2005) that enables calculating the fraction of the sky covered by cloud in terms of oktas (N), overall leading to nine classes, corresponding to the values of N ranging from 0 (clear-sky conditions) to 8 (completely overcast situation). As reported by Ehnberg and Bollen (2005), the amount of global irradiance (F_glo) is related to the solar elevation angle ($\mathit{\pi }/\mathrm{2}-{\mathit{\theta }}_z$) and to the cloudiness following the Nielsen et al. (1981) equation:

$$\begin{array}{}\text{(4)}& F_{\text{glo}}\left(N\right)=\left[{\displaystyle \frac{\begin{array}{l}{a}_{\mathrm{0}}\left(N\right)+a_{\mathrm{1}}\left(N\right)\sin \left(\frac{\mathit{\pi }}{\mathrm{2}}-{\mathit{\theta }}_z\right)\\ +a_{\mathrm{3}}\left(N\right){\sin }^{\mathrm{3}}\left(\frac{\mathit{\pi }}{\mathrm{2}}-{\mathit{\theta }}_z\right)-L\left(N\right)\end{array}}{a\left(N\right)}}\right],\end{array}$$

where N represents one of the nine possible classes of sky conditions expressed in oktas (from 0 for clear-sky conditions to 8 for completely overcast) and a, a₀, a₁, a₃ and L are empirical coefficients that enable computing the expected global irradiance for each okta class (F_glo(N)), at a fixed solar elevation angle ($\mathit{\pi }/\mathrm{2}-{\mathit{\theta }}_z$). Their values, extracted from the original work of Ehnberg and Bollen (2005), are summarized in Table S1 in the Supplement. Overall, Eq. (4) allows for determining the unique okta value N by comparing the measured global irradiance (F_glo) with F_glo(N) at any given time.

With this approach, the cloudiness can be used to evaluate the interaction between incoming radiation and LAA in cloudy conditions but does not provide the opportunity to discriminate between cloud type. The following section describes the methods applied to overcome this limitation by implementing a cloud classification scheme.

2.3.2 Cloud classification

The identification of clouds classes is by common practice still largely performed by human observations based on the reference standard defined by the World Meteorological Organization (WMO; https://cloudatlas.wmo.int/en/home.html, last access: 22 March 2021). However, these observations lack the required time resolution which was needed in the present work to couple highly time-resolved HR data with cloud type. Cloud classification literature reports a huge quantity of papers and reviews aimed at classifying clouds by means of different techniques and their integration to avoid the limits of a simple human inspection. Most of these rely on different ensemble of instruments: (1) ground-based, (2) remote-sensing- or satellite-based, or (3) installed on meteorological balloons (Tapakis and Charalambides, 2013). Some examples are reported in Singh and Glennen (2005), Ricciardelli et al. (2008), Calbó and Sabburg (2008), and Tapakis and Charalambides (2013).

To exploit the full potential of our measurements, we needed a cloud type classification method able to follow the high temporal resolution of the observations including the high spatial and temporal variability of clouds.

Among the abovementioned instrumental ensembles, ground-based instruments provide measurement of the incident solar irradiance for detecting the effect of clouds (Calbò et al., 2001). The concept of using irradiance measurements to estimate cloud types was first introduced in the work of Duchon and O'Malley (1999), which is based on the fact that clouds with different velocities and optical depths cross the slowly changing path of the solar beam over different time durations. Given the available irradiance data (Sect. 2.1), in the present work, the cloud classification starts from the Duchon and O'Malley (1999) method which was successfully applied in the geographical context of the Po Valley (Galli et al., 2004). In particular, we used irradiance measurements (F_glo) to compute two parameters, R_t and SD_t, as follows:

$$\begin{array}{}\text{(5)}& {\displaystyle }& {\displaystyle }{R}_t={\displaystyle \frac{\mathrm{1}}{\mathrm{20}}}\sum _{i=t-\mathrm{10}}^{i=t+\mathrm{10}}{\displaystyle \frac{F_{\text{glo}\left(i\right)}}{F_{\text{glo\_CS}\left(i\right)}}},\text{(6)}& {\displaystyle }& {\displaystyle }{\text{SD}}_t={\mathit{\sigma }}_{t\pm \mathrm{10}}(F_{\text{glo}(t\pm \mathrm{10})}\cdot {\text{Sf}}_{t\pm \mathrm{10}}),\end{array}$$

where R_t is the 20 min moving average ratio between the observed global irradiance (F_glo) and the modelled clear-sky irradiance (Robledo and Soler, 2000) expected at the same place (F_{glo_CS}) at time t. R_t describes the time-dependent cloud efficiency in reducing the incoming solar radiation (R_t = 1 in perfect clear-sky conditions, while R_t ∼ 0 in completely overcast conditions). SD_t represents the 20 min SD (standard deviation) of the scaled global irradiance (F_glo⋅Sf) centred at the time t and describes the temporal stability of clouds in the atmosphere (e.g. persistent stratus clouds are characterized by SD_t ∼ 0, while cumulus clouds in good weather are characterized by higher values of SD_t).

The scaling factor Sf_t (Duchon and O'Malley, 1999) is given by

$$\begin{array}{}\text{(7)}& {\text{Sf}}_t={\displaystyle \frac{\mathrm{1400}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}}{F_{\text{glo\_CS}\left(t\right)}}}.\end{array}$$

Figure 2 Cloud classification based on broadband solar radiation following Duchon and O'Malley (1999). Each row represents a different cloud type on a specific day as a case study. The left column represents the time series of global and diffuse measured solar irradiance (F_glo and F_dif) and modelled clear-sky irradiance (F_{glo_CS}), while the right column contains the scatter SD–R plot of the observed SD of irradiance (SD) vs. the fraction of modelled clear-sky irradiance (R). In panel (h) different colours are related to different times (hours) of the day as reported in the legend.

Figure 3 SD–R plot of the whole dataset concerning the cloud base altitude grouped into three levels, namely low-level clouds (< 2 km), mid-altitude clouds (2–7 km) and high-altitude clouds (> 7 km).

Visualization of the SD vs. R (SD–R plot) results thus represents a first tool in distinguishing different cloud categories as a function of their efficiency in reducing the incoming solar radiation (R) and their persistency (SD). The potential of the SD–R plot is presented in Fig. 2a–h; it shows four examples of the temporal evolution of the observed F_glo, F_{glo_CS} andF_dif (left column) and the corresponding SD–R diagrams (right column). Explored more in detail are the following:

1. The first case (Fig. 2a) shows F_glo following F_{glo_CS} without any significant temporal deviation, thus leading to a cluster of data in the SD–R diagram (Fig. 2b) characterized by R ∼ 1 and SD ∼ 0 W m⁻². These conditions are those associated with clear-sky (CS) conditions by Duchon and O'Malley (1999).

2. The second case (Fig. 2c) shows F_glo completely dominated by the diffuse irradiance (F_dif) throughout the day (note that in Fig. 2c F_dif is superimposed on F_glo); this condition differs completely from the CS case, as both R and SD approach 0 (Fig. 2d). Duchon and O'Malley (1999) associate these conditions with the presence of persistent stratiform clouds.

3. The third case (Fig. 2e) reports F_glo approaching F_{glo_CS} and being at the same time characterized by small amplitude oscillations. In this case R ranges between 0.75 and 1, and SD ranges from 0 to ∼ 100 W m⁻² (Fig. 2f). The cluster of data is thus more dispersed than that of the CS case featuring a larger variation in R and SD. Duchon and O'Malley (1999) attributed this situation to the presence of cirrus (Ci), underlining that in some borderline cases a misclassification between CS and Ci (just based on SD–R plot) could be possible.

4. The last case (Fig. 2g) represents a transition from a CS situation (before noon) to cloudy conditions (after midday) characterized by a significant scatter of F_glo. Figure 2h clearly shows that the sky condition evolves from the CS toward cloudy sky, shifting the R data from ∼ 1 down to ∼ 0.25 and increasing SD from ∼ 100 to ∼ 500 W m⁻². According to Duchon and O'Malley (1999), the arrival of cumulus during a “good-weather” day could be the reason for such behaviour (Cu cloud movement in the sky results in fast sun–shadow transitions). Also, in this case, the SD–R plot alone cannot exclude the presence of other cloud types responsible for a similar behaviour (e.g. altocumulus, Ac; cirrocumulus, Cc; and cirrostratus, Cs). Note that in order to show the variation of data in the SD–R diagram (Fig. 2h) as a function of time, an hourly resolved colour code was assigned to the data points; the corresponding regions in Fig. 2g were delimited by dashed lines with the same colour code.

Overall, Fig. 2a–h shows the potential (and limits) of the SD–R plots for a preliminary broad sky–cloud classification. As mentioned, the SD–R diagram alone leaves margins of misclassification, especially because it is impossible to retrieve the required information when different cloud types at different levels are present simultaneously.

In the present work, we attempted a further refinement of cloud classification, including the information of the cloud base height (CBH) and the number of cloud layers obtained from the automated lidar–ceilometer measurements. The cloud base height is a key parameter in the characterization of clouds (Hirsch et al., 2011), since its estimation limits the number of potential cloud classes (that the SD–R classifier has to discriminate between), thus maximizing the efficiency of the Duchon and O'Malley (1999) classification algorithm. In fact, ceilometer instruments were developed and are commonly used in airports to operationally detect cloud layers, and their use for aerosol-related studies is more recent. Furthermore, the use of ceilometer data for cloud classification and cloud study purposes does not represent an absolute novelty in the scientific literature as demonstrated by recent works by Huertas-Tato et al. (2017) and Costa-Surós et al. (2013). The availability of CBH information allows for dividing cloud types in three fundamental categories (Tapakis and Charalambides, 2013): low-level clouds ( < 2 km), mid-altitude clouds (2–7 km) and high-altitude clouds ( > 7 km). From a general perspective the high-altitude cloud category includes cirrus (Ci), cirrocumulus (Cc) and cirrostratus (Cs); mid-altitude clouds include altocumulus (Ac), altostratus (As) and nimbostratus (Ns); low-level clouds include cumulus (Cu), stratocumulus (Sc), stratus (St) and cumulonimbus (Cb) (Tapakis and Charalambides, 2013; Ahrens, 2009; Cotton et al., 2011).

We colour-coded the SD–R diagram in Fig. 3 using the ceilometer-based information on cloud altitude. The plot shows that, on average, low-level clouds are located on the left side of the SD–R diagram (stratiform clouds), while high-altitude clouds are conversely on the opposite side (Ci and Cu clouds); finally, mid-altitude clouds mostly cover the central part, describing all the possible transitions and combinations from St to Cu and Ci, e.g. altostratus (As) and altocumulus (Ac).

Overall, adding the CBH information to the SD–R plot enabled us to identify eight cloud types: St (stratus), Cu (cumulus) and Sc (stratocumulus) as low-level clouds; As (altostratus) and Ac (altocumulus) as mid-altitude clouds; Ci (cirrus) and Cc–Cs (cirrocumulus and cirrostratus merged in one single class) as high-altitude clouds.

Table 1 Final criteria adopted for cloud classification. SD represents the SD of the measured global irradiance with respect to the theoretical behaviour in clear-sky conditions; R represents the ratio between observed global irradiance (F_glo) and the modelled irradiance (F_{glo_CS}) in clear-sky conditions; and finally the cloud layer is the number of cloud layers detected by the lidar.

A summary of the threshold values of R, SD and cloud level used here to the final cloud classification is given in Table 1, with the R and SD limits being based on the works of Duchon and O'Malley (1999) and Harrison et al. (2008) and those of the CBH being derived considering the cloud properties at midlatitudes.

Figure 4 Cloud classification based on the improved broadband solar radiation following Duchon and O'Malley (1999) and Harrison et al. (2008) coupled with lidar data of cloud base height. From left to right: stratus (St), altostratus (As), stratocumulus (Sc), altocumulus (Ac), cirrocumulus and cirrostratus (Cc–Cs), cumulus (Cu), cirrus (Ci), and finally clear-sky (CS) conditions. The SD–R plot reports in grey the single data of the whole dataset, while centroids and the 99 % confidence interval of each cloud type are plotted in a colour scale related to the cloud base level.

Finally, to avoid misclassification due to the presence of multiple cloud layers, the analysis was limited to those cases where only one cloud layer was detected by the ceilometer (8405 single layer cases, representing 61 % of all measurements). Another reason for limiting the analyses to one cloud layer is due to the main aim of this work: to quantify the effects of different cloudiness and cloud types on the LAA HR. We wanted to avoid conditions with multiple-layer clouds, as this would result in confounding information for the purpose of the present study.

Figure 4 shows the SD–R diagram of all data (grey) with superimposed R and SD mean values and a 99 % confidence interval for each of the eight identified cloud classes, plus clear-sky (CS) conditions. The final cloud classification was obtained for the period from November 2015 to March 2016, during which all necessary parameters were available (Sect. 3).

Since this methodology is applied for the first time in the Po Valley, a complete validation of the aforementioned approach is reported in Appendix B (“Cloud type validation”). It includes two validation exercises: the first was carried out comparing the present automatized cloud classification with a visual cloud classification based on sky images collected during 1 month of the wintertime field campaign; the second was carried out comparing the present automatized cloud classification with the one discussed by Ylivinkka et al. (2020). In fact, simultaneously to the submission of our work, Ylivinkka et al. (2020) proposed a classification based on the coupling of irradiance and CBH measurements. Overall, based on these comparisons, agreement with our classification is 80 % with the visual approach and 90 % with the Ylivinkka et al. (2020) methodology, with these results further demonstrating the reliability of the cloud classification algorithm used in our study.

3 Results and discussion

Data measured over Milan from November 2015 to March 2016 are presented in Sect. 3.1, with this period covering the simultaneous presence of radiation, lidar–ceilometer and absorption information necessary for the analysis. The role of cloudiness and cloud type on the total HR is discussed in Sect. 3.2; the impact of clouds on the HR is discussed with respect to the light-absorbing aerosol species, BC and BrC, in Sect. 3.3. All data are reported as the mean ± 95 % confidence interval.

Figure 5 High-time-resolution data (5 min) for eBC, global irradiance (F_glo, yellow line) cloud base height (CBH), cloudiness (oktas) and the related heating rate (HR) from 1 November 2015 to 1 April 2016.

Figure 6 Monthly averaged values of (a) eBC and HR values and their direct, diffuse and reflected components (HR_dir, HR_dif and HR_ref); (b) global radiation values (F_glo) and their direct, diffuse and reflected components (F_dir,F_dif and F_ref).

Table 2 Monthly averaged data and the confidence interval at 95 % of temperature (T), pressure (P), equivalent black carbon (eBC), absorption coefficient (b_abs) and heating rate (HR) divided into their direct (dir), diffuse (dif) and reflected (ref) components and, finally, global (F_glo), direct (F_dir), diffuse (F_dif) and reflected (F_ref) irradiances.

^* denotes aethalometer data referring to λ = 880 nm.

3.1 eBC, irradiance, HR and cloud data presentation

Highly time-resolved data (5 min) of eBC, F_glo, CBH, cloudiness (oktas) and the resulting HR are shown in Fig. 5; their monthly average values are presented in Fig. 6a and summarized in Table 2.

The lower eBC and b_abs (880 nm) values (monthly averages of 1.54 ± 0.04 µg m⁻³ and 7.6 ± 0.2 Mm⁻¹) were recorded in March, while their higher values were found in December (6.29 ± 0.09 µg m⁻³ and 31.1 ± 0.5 Mm⁻¹, respectively) with a maximum value of 27.44 µg m⁻³ (135.7 Mm⁻¹). In December, the average PM₁₀ and PM_2.5were also at their maximum, with 73.1 ± 0.6 and 69.3 ± 0.6 µg m⁻³, respectively (source: Milan Environmental Protection Agency, ARPA Lombardia, https://www.arpalombardia.it/Pages/Aria/Richiesta-Dati.aspx, last access: 25 March 2021), and the eBC accounted for ∼ 10 % of PM mass concentration. These high values of eBC and PM₁₀ and PM_2.5 agree with those observed previously in wintertime in the Po Valley, when strong emissions in the Po Valley are released into a stable boundary layer (Sandrini et al., 2014; Ferrero et al., 2011b, 2014, 2018; Barnaba et al., 2010). During the investigated period, the lower monthly irradiance value was observed in December (F_glo of 141 ± 4 W m⁻²; Table 2), while the higher value was in March (F_glo of 310 ± 7 W m⁻²). The higher monthly average HR was recorded in December (1.43 ± 0.05 K d⁻¹), while the lower one was in March (0.54 ± 0.02 K d⁻¹; see Fig. 6a and Table 2). Even though the HR monthly behaviour is correlated with eBC (Table 2; R² = 0.82, not shown), it is also useful to compare the maximum-to-minimum ratio of the eBC monthly mean (December to March, eBC ratio of 4.10 ± 0.12) to the same for the HR (2.65 ± 0.16). This ratio is higher for eBC because the incoming irradiance was lower in December (F_glo of 141 ± 4 W m⁻²; Fig. 6b) with respect to March (F_glo of 310 ± 7 W m⁻², ratio of 0.45 ± 0.02), partially compensating the marked wintertime increase of eBC. This is due to the interaction of LAA with F_dir. In fact, once F_dir is scaled by cos (θ_z) (Eq. 3, Sect. 2.2, Fig. S4 in the Supplement) it is quite constant throughout the year (and perfectly constant only in clear-sky conditions). Conversely, the diffuse and reflected irradiance, under the isotropic and Lambertian assumptions (Eq. 3), remain seasonally modulated (Fig. S4).

These observations illustrate the importance of both the amount and the type (direct, diffuse and reflected) of radiation that interacts with LAA. In brief, any process able to influence the total amount and the type of impinging irradiance (e.g. presence or absence of clouds, cloudiness, and cloud type) will result in a different HR, even at constant LAA concentrations (and their absorption). The investigation of this aspect is the main focus and added value of this study. High-resolution data (Figs. 5 and S4) provided a first hint to the importance of cloud presence on the HR; a sharp global irradiance decrease was observed in cloudy conditions, especially in the presence of low-level clouds (low CBH) and high cloud cover (7–8 oktas).

Thus, both cloudiness and cloud type were carefully determined as detailed in Sect. 2.3.1 and 2.3.2. Overall, during the whole campaign, the average cloudiness was 3.58 ± 0.04 oktas with the higher monthly value in February (4.56 ± 0.07 oktas) and the lower one in November (2.91 ± 0.06 oktas). These data are in line with the mean cloudiness over Europe (∼ 5.5 oktas; Stjern et al., 2009) and over Italy (∼ 4 oktas; Maugeri et al., 2001). Moreover, during the campaign, clear-sky (CS) conditions were only present 23 % of the time, with the remaining time (77 %) being characterized by partially cloudy (35 %, 1–6 oktas) to totally cloudy (42 %, 7–8 oktas) conditions.

Figure 7 (a) Time frequency (%) of the cloud type classified over the U9 site (CS means clear-sky conditions); (b) contribution (%) of each cloud type to the okta values measured over the U9 site.

Figure 8 Cloudiness associated with each cloud type.

Cloudy conditions are therefore frequent. The frequency of specific cloud type occurrence is given in Fig. 7a. The dominating cloud type was St (42 %), followed by Sc (13 %), Ci and Cc–Cs (7 % and 5 %, respectively). The contribution of each cloud type to the cloudiness is reported in Fig. 7b. While St clouds were mostly responsible for overcast situations (7–8 oktas, frequency of 87 % and 96 %), Sc clouds dominated the intermediate cloudiness conditions (5–6 oktas, frequency of 47 % and 66 %), and the transition from Cc–Cs to Sc determined moderate cloudiness (3–4 oktas). Finally, low cloudiness values (1–2 oktas) were mostly dominated by Ci and Cu (frequency of 59 % and 40 %, respectively). As mentioned (Sect. 2.3.2 and Fig. 4), low-level clouds ( < 2 km) include stratus (St), cumulus (Cu) and stratocumulus (Sc); mid-altitude clouds (2–7 km) include altostratus (As) and altocumulus (Ac); and high-altitude clouds ( > 7 km) include cirrus (Ci), cirrocumulus and cirrostratus (Cc–Cs). Thus, it is clear that the higher cloud cover (higher okta value) is due to a higher frequency of low–mid-altitude clouds. This is evident in Fig. 7b, which reports the average CBH for each okta. The CBH was related to oktas (Fig. S5a in the Supplement), underling the linkage (together with Fig. 7b) between the fraction of the sky covered by clouds and the cloud type responsible for it, at least at the measuring site. Indeed, the cloudiness is a non-linear function of the cloud type, as cloud types are related to the meteorological patterns; e.g. highly persistent stratiform clouds generate cloudy weather in conditions with lower wind (see the Supplement for further details). Figure 8 summarizes the average cloudiness associated with different cloud types showing an okta rise from conditions dominated by cirrus clouds (0.51 ± 0.05 oktas) to stratus clouds (7.20 ± 0.04 oktas). This is in agreement with the recent work of Bartoszek et al. (2020), who associated a higher cloudiness level with the presence of stratiform clouds. The possible role of wind on cloud type is explored in Fig. S6 in the Supplement (“Wind speed, cloudiness and clouds”).

Figure 9 Monthly averaged values of (a) HR values and their direct, diffuse and reflected components (HR_dir, HR_dif and HR_ref) during winter and spring both in clear-sky (CS; 0 oktas) and cloudy (CLD; 7–8 oktas) conditions. (b) $\mathrm{HR}/\mathrm{eBC}$ values together with their direct, diffuse and reflected components (${\text{HR}}_{\text{dir}}/\text{eBC}$, ${\text{HR}}_{\text{dif}}/\text{eBC}$ and ${\text{HR}}_{\text{ref}}/\text{eBC}$); the direct, diffuse and reflected irradiance (F_dir, F_dir and F_dif); and the global irradiance (F_glo).

3.2 Cloud impact on the heating rate

3.2.1 The role of cloudiness

Figure 6a already provided the first indication of the important influence of clouds on the total HR. In fact, it shows the magnitude of the absolute (and relative) contribution of the diffuse component (HR_dif) with respect to the total HR revealing that, on a monthly basis, the diffuse contribution accounted on average for 40 ± 1 % (of the total HR). In most cases this was comparable or even higher than the HR_dir. The only exception was in November 2015 when the lower HR_dif (Fig. 6a) and F_dif (Fig. 6b) fractions in the total HR and F_glo were measured (30.4 ± 1.4 % and 34.3 ± 2.6 % of the total, respectively), with this also being the month with the lowest average cloudiness (2.91 ± 0.06 oktas). The aforementioned data demonstrate the importance of the diffuse component of radiation. Therefore, the absolute values of the HR and its components were firstly investigated as a function of cloudiness (clear-sky and completely overcast situations, seasonal averages, Fig. 9a). In the wintertime clear sky, the direct component of the HR (HR_dir) was higher than the HR_dif and HR_ref, accounting for 1.35 ± 0.04 K d⁻¹ and explaining on average 60 ± 5 % of the total HR. Similarly, in the springtime clear sky, HR_dir was 0.47 ± 0.01 K d⁻¹, again higher than the HR_dif and HR_ref. Conversely, in a completely overcast condition (7–8 oktas), the HR_dif dominated (84 ± 1 % of the total HR) and accounted for 0.33 ± 0.01 and for 0.19 ± 0.01 K d⁻¹ during winter and spring, respectively.

In order to further investigate the role of cloudiness, we decoupled the variability of the HR induced by radiation from that due to LAA concentrations. Thus, the HR values and those of its components (HR_dir, HR_dif and HR_ref) were normalized to the unit mass of eBC ($\mathrm{K}{\mathrm{m}}^{\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}\mathrm{\mu }{\mathrm{g}}^{-\mathrm{1}}$) and reported as a function of cloudiness in Fig. 9b together with the measured irradiance (F_glo, F_dir, F_dif and F_ref): this parameter ($\mathrm{HR}/\mathrm{eBC}$) reports the efficiency of warming per mass concentration of eBC at different cloudiness levels. Overall, Fig. 9b shows the general decease of $\mathrm{HR}/\mathrm{eBC}$ for increasing cloud cover, a pattern also observed for both ${\text{HR}}_{\text{dir}}/\text{eBC}$ and ${\text{HR}}_{\text{ref}}/\text{eBC}$, which follow the respective decrease of direct and reflected irradiance. Note that at okta values of 7–8, ${\text{HR}}_{\text{dir}}/\text{eBC}$ reached values close to 0 (due to the suppression of F_dir by clouds), while ${\text{HR}}_{\text{ref}}/\text{eBC}$ was 0.03 ± 3 × 10⁻⁴ $\mathrm{K}{\mathrm{m}}^{\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}\mathrm{\mu }{\mathrm{g}}^{-\mathrm{1}}$ due to the presence of surficial albedo effect on the diffuse irradiance (F_dif). ${\text{HR}}_{\text{dif}}/\text{eBC}$ increased with increasing cloudiness up to intermediate cloudiness conditions (5–6 oktas), reaching a maximum (0.16 ± 0.01 $\mathrm{K}{\mathrm{m}}^{\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}\mathrm{\mu }{\mathrm{g}}^{-\mathrm{1}}$). This is in line with the behaviour of the diffuse irradiance: maximum of 147 ± 6 W m⁻² (at 5–6 oktas), doubling the value in overcast conditions (74 ± 3 W m⁻²; 7–8 oktas) and exceeding 150 % of that for clear-sky conditions (91 ± 2 W m⁻²). In the overcast condition (7–8 oktas) both ${\text{HR}}_{\text{dif}}/\text{eBC}$ and the diffuse irradiance reached their minimum due to the capability of clouds to effectively attenuate the incoming radiation. However, in these conditions, ${\text{HR}}_{\text{dif}}/\text{eBC}$ was still not null (0.08 ± 0.01 $\mathrm{K}{\mathrm{m}}^{\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}\mathrm{\mu }{\mathrm{g}}^{-\mathrm{1}}$), dominating the total atmospheric HR, with a contribution of 84 ± 1 %.

$\mathrm{HR}/\mathrm{eBC}$ and cloudiness data were linearly related, showing a high level of correlation (R² = 0.935, Fig. S5b). Cloudiness could thus be used as good predictor (in modelling activity) for $\mathrm{HR}/\mathrm{eBC}$.

As from Fig. S5a (Sect. 3.1), the CBH appeared related to the cloudiness, and an additional linear correlation was tested between $\mathrm{HR}/\mathrm{eBC}$ and the CBH (Fig. S5c; R² = 0.857); this relationship is weaker than that between $\mathrm{HR}/\mathrm{eBC}$ and cloudiness. Cloudiness, describing the fraction of sky covered by clouds, is a better predictor of the capability to suppress the incoming radiation (and thus the HR promoted by LAA). The relationship between the CBH and cloudiness should be also investigated in other monitoring sites around the world to explore the possibility of using the CBH (together with cloudiness) as a promising prognostic variable for the HR of LAA in future studies.

Figure 10 Diurnal pattern of eBC (a), wind speed (b), global irradiance (F_glo) (c) and HR (d). Data are averaged for clear-sky conditions (CS; 0 oktas) and cloudy conditions (CLD; 7–8 oktas).

Overall, our experimental HR data enabled us to estimate the degree of error introduced by improperly assuming clear-sky conditions in radiative-transfer calculations. Particularly, we found that the simplified assumption of clear-sky conditions leads to an overestimation of the LAA-induced HR by a factor ranging from 50 % to 470 % (50 % in low cloudiness: 1–2 oktas; 109 % in moderate cloudiness: 3–4 oktas; 148 % in intermediate cloudiness: 5–6 oktas; and 470 % in cloudy conditions: 7–8 oktas). These results clearly highlight that clouds are responsible for an important feedback on the aerosol HR that needs to be carefully quantified, pointing to the need to correctly include and model cloudy conditions in radiative-transfer calculations aimed at evaluating the real contribution of aerosol forcing on the atmospheric HR on a global scale.

3.2.2 Cloudiness and diurnal pattern of the HR

The presence of clouds can also alter the HR diurnal pattern. Figure 10a–d show the mean diurnal pattern of eBC, wind speed, F_glo and the HR in both clear-sky (0 oktas) and cloudy conditions (7–8 oktas). In clear-sky conditions, the eBC peaked at 08:00 LST (6.41 ± 0.31 µg m⁻³) during the rush hour (Fig. 10a); then eBC decreased until its minimum in the early afternoon (1.07 ± 0.10 µg m⁻³) when the wind speed reached its maximum (1.5 ± 0.1 m s⁻¹, Fig. 10b). The incoming F_glo in clear-sky conditions peaked as expected at midday with 497 ± 10 W m⁻² (Fig. 10c). This caused an asymmetric HR diurnal pattern, being characterized by a fast increase to the maximum at 10:00 LST (3.60 ± 0.18 K d⁻¹) and a subsequent slower decrease by sunset (Fig. 10d). This pattern was not present in cloudy conditions (Fig. 10d). First, eBC showed a moderate peak at 10:00 LST (4.09 ± 0.20 µg m⁻³) being quite stable during the afternoon – remaining above 3 µg m⁻³ until 16:00 LST (Fig. 10a). The eBC behaviour was consistent with that of wind speed, which only slightly rose during the day but was however always below 1 m s⁻¹ (on average 0.64 ± 0.03 m s⁻¹, Fig. 10b). The incoming F_glo in cloudy conditions peaked again as expected at midday with 103 ± 4 W m⁻² with a much slower increase during the day (Fig. 10c). The Supplement (“Wind speed, cloudiness and clouds”) and Fig. 7b show that cloudy conditions were mostly associated with stratus and very low windy conditions (0.64 ± 0.02 m s⁻¹), explaining the flat diurnal behaviour of eBC differing from the clear-sky case. Moreover, the absence of any direct irradiance in cloudy conditions (Fig. 9b; Sect. 3.1) determines that F_glo was essentially due to the diffuse irradiance whose symmetrical bell-shaped curve drove the HR behaviour (Fig. 10d), peaking at midday with a value of 0.74 ± 0.01 K d⁻¹ (much lower than in CS).

Figure 11 Impact of each cloud type on the heating rate normalized to black carbon concentration: (a) $\mathrm{HR}/\mathrm{eBC}$ and F_glo, (b) ${\text{HR}}_{\text{dir}}/\text{eBC}$ and F_dir, (c) ${\text{HR}}_{\text{dif}}/\text{eBC}$ and F_dif, and (d) ${\text{HR}}_{\text{ref}}/\text{eBC}$ and F_ref.

As a conclusion, in different cloudiness conditions, not only the absolute magnitude of the HR but also its diurnal pattern are different. This also changes the related atmospheric feedbacks, such as the influence on the liquid water content (Jacobson, 2002), planetary boundary layer dynamics (Wang et al., 2018; Ferrero et al., 2014), regional circulation systems (Ramanathan and Feng, 2009; Ramanathan and Carmichael, 2008), and finally the cloud dynamic and evolution itself (Koren et al., 2008; Bond et al., 2013). Thus, an inappropriate use of the clear-sky assumption in models will also reflect on the modelled HR-triggered feedbacks. These results also acquire relevance in the context of the counterintuitive semi-direct effect proposed by Perlwitz and Miller (2010) and referred to in Sect. 1: the atmospheric heating induced by tropospheric absorbing aerosol could lead to a cloud cover increase (especially low-level clouds). Such a feedback stresses the need for a proper inclusion of sky conditions into radiative-transfer calculations.

3.2.3 The role of cloud type

The previous sections showed the effect of cloudiness on the total LAA HR. The impact of each cloud type on the HR is addressed here, as not all clouds have the same effect on irradiance (Tapakis and Charalambides, 2013). As previously done, we refer to HR values normalized to eBC unit mass ($\mathrm{HR}/\mathrm{eBC}$) to decouple radiation and aerosol effects. Figure 11a–d show the total $\mathrm{HR}/\mathrm{eBC}$ and F_glo together with the corresponding components (${\text{HR}}_{\text{dir}}/\text{eBC}$ and F_dir, ${\text{HR}}_{\text{dif}}/\text{eBC}$ and F_dif, and ${\text{HR}}_{\text{ref}}/\text{eBC}$ and F_ref; Fig. 11b–d). The figure shows a prefect agreement between cloud type, irradiance and the corresponding $\mathrm{HR}/\mathrm{eBC}$ component (R² > 0.93; not shown). It also highlights how critical it is, for radiative-transfer calculations and HR determination, to take into account the role of each cloud type. Only the cloud influence on the ${\text{HR}}_{\text{dif}}/\text{eBC}$ is markedly different from the other components.

Figure 12 Average values of the total HR, HR_dir, HR_dif and HR_ref as a function of the cloud type.

In terms of absolute values (not normalized for eBC), Fig. 12 reveals that the HR_dir was only dominant during periods of CS and Ci clouds (HR_dir of 1.11 ± 0.04 and 0.92 ± 0.05 K d⁻¹, respectively), explaining 66 ± 3 % and 57 ± 4 % of the total atmospheric HR. In the cases of other clouds (St, As and Sc) HR_dif dominates, reaching the highest absolute contribution of 84.4 ± 3.8 %, 83.0 ± 10.7 % and 76 ± 4 % (HR_dif of 0.25 ± 0.01, 0.34 ± 0.03 and 0.66 ± 0.02 K d⁻¹), respectively.

Figure 13 Percentage decrease of the HR with respect to clear-sky conditions as a function of the cloudiness (oktas) averaged for each cloud type.

Given this impact of cloud type, the ability of cloudiness to be a good predictor for the HR (as detailed in Sect. 3.2.1), and the relationship (over the investigated site) between cloudiness and cloud type (Sect. 3.1, Fig. 7b), the synergic impact of cloudiness and cloud type on the HR was investigated and presented in Fig. 13. In the figure, we summarize the HR results in terms of percent difference from the clear-sky (CS) case by averaging the cloudiness (in oktas) for each cloud type (as detected in Sect. 3.3). Overall, the derived linear regression indicates an HR decrease of −11.9 ± 1.2 % per okta. The regression R² (0.963) was slightly higher than that reported in Fig. S5b (R² = 0.935; relationship with the cloudiness only) suggesting the need (for precise calculations) to account for the cloud types responsible for any sky coverage in agreement with a recent work of Bartoszek et al. (2020). Figure 13 also allowed us to associate the HR decrease with each specific cloud type over Milan. Particularly, Ci clouds produced a modest impact on cloudiness (0.50 ± 0.05 oktas), decreasing the HR by ∼ 3 %, while Cu clouds (1.76 ± 0.09 oktas) decreased the HR by −26 ± 8 %. Cc–Cs clouds (3.56 ± 0.14 oktas) were responsible for a −49 ± 6 % decrease of the HR. Their impact was comparable to that of Sc clouds (4.68 ± 0.10 oktas, −48 ± 4 % of the HR). Ac clouds (4.11 ± 0.18 oktas) had a higher impact, decreasing the HR by −59 ± 6 %. The highest impact was due to As (6.57 ± 0.15 oktas; −76 ± 4 % of the HR) and by St (7.19 ± 0.04 oktas) that suppressed the HR by a factor of −83 ± 4 %.

Figure 14 Monthly averaged data for the HR of both BC and BrC.

3.3 The impact of clouds on the BC and BrC heating rates

In this last part of the work we focus on the HR of the two main absorbing aerosol species: BC and BrC (obtained as detailed in Sect. 2.1.1). The monthly averaged values of the HR of BC and BrC (HR_BC and HR_BrC) are reported in Fig. 14. The highest HR_BC and HR_BrC values were recorded in December (1.24 ± 0.03 K d⁻¹ and 0.19 ± 0.01 K d⁻¹), while the lowest were recorded in March (0.46 ± 0.01 K d⁻¹ and 0.07 ± 0.01 K d⁻¹). Overall, the HR_BrC accounted for 13.7 ± 0.2 % of the total HR.

Figure 15 The HR of BC and BrC as a function of the cloudiness (oktas): (a) total HR_BC and HR_BrC, (b) direct component of both the HR_BC and HR_BrC (HR_BC,dir and HR_BrC,dir), (c) diffuse component of both the HR_BC and HR_BrC (HR_BC,dif and HR_BrC,dif), and (d) reflected component of both the HR_BC and HR_BrC (HR_BC,ref and HR_BrC,ref). Note that, due to the different magnitude of the HR_BC and HR_BrC, the y axis of the HR_BrC in the four panels was chosen as $\mathrm{1}/\mathrm{10}$ of that of the HR_BC.

The variability of the total HR_BC and HR_BrC as a function of cloudiness is reported in Fig. 15a, with panels b–d showing their direct (HR_BC,dir and HR_BrC,dir), diffuse (HR_BC,dif and HR_BrC,dif) and reflected (HR_BC,ref and HR_BrC,ref) components. Figure 15a shows that both the HR_BC and HR_BrC decreased with increasing cloudiness, going from the CS maxima (HR_BC and HR_BrC of 1.14 ± 0.03 and 0.20 ± 0.01 K d⁻¹) to the completely overcast condition minima of 0.16 ± 0.01 and 0.02 ± 10⁻³ K d⁻¹ (8 oktas; mainly due to St and As clouds; see Fig. 7b). As shown in Fig. 9b, the change of irradiance magnitude with cloudiness was different for direct, diffuse and reflected components affecting the corresponding direct, diffuse and reflected components of the HR_BC and HR_BrC (Fig. 15b–d). The HR_BC,dir and HR_BrC,dir (Fig. 15b) decreased as a function of cloudiness from 0.74 ± 0.03 and 0.11 ± 0.01 K d⁻¹ (0 oktas) to negligible levels (HR < 10⁻⁴ K d⁻¹) in completely overcast conditions. The HR_BC,dif and HR_BrC,dif (Fig. 15c) increased with cloudiness, reaching their maximum in partially cloudy conditions (at 6 oktas: 0.51 ± 0.01 and 0.09 ± 0.01 K d⁻¹). Further increasing cloudiness reduced their values to minimum values (0.13 ± 0.01 and 0.02 ± 0.01 K d⁻¹). The HR_BC,ref and HR_BrC,ref (Fig. 15d) behave similarly to the total HR_BC and HR_BrC, since the reflected irradiance is dominated by the global irradiance impinging on the ground (see Fig. 9b for a comparison); the HR_BC,ref and HR_BrC,ref decreased with increasing okta values from maximum values in clear-sky conditions (HR_BC,ref and HR_BrC,ref of 0.17 ± 4 × 10⁻³ and 0.03 ± 1 × 10⁻³ K d⁻¹) down to the overcast minimum (HR_BC,ref and HR_BrC,ref of 0.02 ± 10⁻³ and 3 × 10⁻³ ± 10⁻³ K d⁻¹). Figure 15a–d also show that the HR_BC was always greater (in absolute values) than the HR_BrC, as expected. The relative decrease of the HR_BrC from CS to completely overcast conditions was 12 ± 6 % larger with respect to that of the HR_BC. At a first glance, Fig. 15a–d could give the impression that BrC is more efficient in heating the surrounding atmosphere (with respect to BC) in CS conditions. However, any change of both BC and BrC b_abs(λ) in different sky conditions has to be taken into account to avoid any misinterpretation of the results. While the variability of BC b_abs(λ) with cloudiness was limited (with the exception of 1 okta, Fig. S7a in the Supplement), this was not the case for BrC. In fact, b_abs(λ) BrC values in high cloudiness were statistically lower than the ones in CS conditions (at 8 oktas, b_abs(λ) of BrC was −23 ± 3 % lower than in CS conditions, Fig. S7b). The relative decrease of the HR_BrC with cloudiness was therefore higher compared to that of the HR_BC. Understanding of the reason behind the observation of higher b_abs(λ) values for BrC in CS conditions is beyond the aim of the present paper (we can speculate that it could be related to the formation of secondary BrC at high radiation levels; e.g. Kumar et al., 2018).

Figure 16 Percentage decrease of the HR_BC and HR_BrC with respect to clear-sky conditions as a function of cloudiness (oktas) averaged for each cloud type.

Here we focus on the fact that the magnitude of b_abs(λ) of BC and BrC changed differently with cloudiness. Thus, in order to decouple the variability of the HR induced by the varying incoming irradiance from that due to changes in b_abs(λ), both the HR_BC and HR_BrC were normalized to the dimensionless integral of b_abs(λ) over the whole aethalometer spectrum. In this way, the magnitude of b_abs(λ) is accounted for along the whole spectrum, avoiding the choice of an arbitrary wavelength as a reference for the normalization. Similarly to Sect. 3.2.2 for the total of the LAA HR, the variability of the normalized HR_BC and HR_BrC was investigated with respect to cloudiness and cloud type; in this respect, both the HR_BC and HR_BrC were normalized to the dimensionless integral of b_abs(λ) for each cloud type. Figure 16a shows the decrease of the normalized HR_BC and HR_BrC as a function of average cloudiness for each cloud type. We found a strong linear relationship between the decrease of both the normalized HR_BC and HR_BrC (relative to CS conditions) and the mean cloudiness (in oktas) for each cloud type. Focusing on the cloud type, Ci clouds were found to produce a statistically negligible impact on cloudiness (0.50 ± 0.05 oktas), decreasing the HR_BC and HR_BrC by ∼ 1 %–6 %, respectively. Cu clouds (1.76 ± 0.09 oktas) decreased the HR_BC and HR_BrC by −31 ± 12 % and −26 ± 7 %, respectively. Cc–Cc clouds featured 3.56 ± 0.14 oktas and were responsible for a −60 ± 8 % and −54 ± 4 % decrease of the HR_BC and HR_BrC. Their impact was comparable to that of Ac (4.11 ± 0.18 oktas): −60 ± 6 % and −46 ± 4 % decrease of the HR_BC and HR_BrC. Sc clouds (4.68 ± 0.10 oktas) had a higher impact, decreasing the HR_BC and HR_BrC of −63 ± 6 % and −58 ± 4 %. The highest impact was given by As (6.57 ± 0.15 oktas; −78 ± 5 % and −73 ± 4 % of the HR_BC and HR_BrC) and by St (7.19 ± 0.04 oktas), suppressing the HR_BC and HR_BrC by −85 ± 5 % and −83 ± 3 %, respectively.

Overall, the derived linear regressions indicate a decrease of ∼ 12 % per okta for both the HR_BC and HR_BrC (with high R² of 0.958 and 0.963, respectively). In detail, the respective decreases of the HR_BC and HR_BrC were −11.8 ± 1.2 % and −12.6 ± 1.4 % per okta, with these values not being statistically different. We show that, while BC and BrC have different optical properties and wavelength dependence of absorption, their HR normalized to absorption changed without any statistical difference as a function of cloudiness and cloud type. This simplifies the models and reduces the number of details needed to be considered: once the HR_BC and HR_BrC are determined in clear-sky conditions, their dependence on the cloudiness can be determined from the simple reduction of the HR normalized to the absorption coefficient (about 12 % for both species, once dominant cloud type is known).

However, it noteworthy that the normalized HR_BrC values in Fig. 16 were always greater than or equal to the corresponding ones of BC (even if 95 % confidence interval bands overlapped). A possible explanation can be the synergic effect between the different spectral absorption of BC and BrC and the influence of clouds on the energy of the impinging radiation; this is detailed in the Supplement (“The role of average photon energy on the HR of BC and BrC”). This feature needs further investigation in other seasons and elsewhere in the world where the prevailing cloud types and the light absorption by BrC might be different.

4 Summary and conclusions

The heating rates (HRs) associated with the two major LAA species, i.e. black carbon (BC) and brown carbon (BrC) (HR_BC and HR_BrC), were experimentally determined based on radiation and aerosol measurements (at high time resolution) in the Po Valley. We determined the impact of cloud–aerosol–radiation interactions on the atmospheric heating by examining the total HR in different sky conditions. Results showed a constant decrease of the LAA HR with increasing cloudiness of the atmosphere (∼ 12 %). Our real-atmosphere, all-sky, measurement-based results suggest that using a simplified assumption of clear-sky conditions in radiative-transfer calculations might overestimate the HR by over 400 %. The effect of different cloud types on the HR was also investigated. While cirrus clouds were characterized by a modest impact, cumulus, cirrocumulus–cirrostratus and altocumulus suppressed the HR of both BC and BrC by a factor of ∼ 2. Stratocumulus, altostratus and stratus clouds suppressed the HR_BC and HR_BrC up to 80 %. The cloudiness also changed the diurnal pattern of the HR with possible feedbacks on planetary boundary layer dynamics and/or regional circulation systems.

The total HR, HR_BC and HR_BrC are affected by both cloudiness and cloud type so that inaccurate HR_BC and HR_BrC estimations can be derived from simulations if presence of clouds is ignored and cloud type is not taken into account. Most importantly, the coupling between the cloud impact on the solar radiation spectrum (and its direct, diffuse and reflected components) and the spectral-absorption properties of BC and BrC showed that the absolute HR_BC and HR_BrC vary differently with cloudiness (especially the diffuse component) but feature a very similar normalized (to the absorption coefficient) dependence on the cloudiness. This simplifies the models and reduces the number of details that need to be considered: once the HR_BC and HR_BrC are determined in clear-sky conditions, their dependence on the cloudiness can be determined from the simple reduction of the HR normalized to the absorption coefficient (about 12 % per okta for both species). These data acquire importance when discussed in the context of the counterintuitive semi-direct effect proposed by Perlwitz and Miller (2010): the atmospheric heating induced by tropospheric absorbing aerosol could lead to a cloud cover increase stressing the need for a proper determination and simulation of sky conditions during radiative-transfer calculations.

Appendix A: Nomenclature

Aerosol acronyms
AAE	Absorption Ångström exponent
AAEBC	Absorption Ångström exponent of black carbon
AAEBrC	Absorption Ångström exponent of brown carbon
babs(λ)	Wavelength-dependent aerosol absorption coefficient (Mm−1)
BC	Black carbon
BrC	Brown carbon
eBC	Equivalent black carbon concentration (µg m−3)
LAA	Light-absorbing aerosol
HR	Heating rate (K d−1)
HRBC	Heating rate of black carbon (K d−1)
HRBrC	Heating rate of brown carbon (K d−1)
Cloud and sky acronyms
As	Altostratus
Ac	Altocumulus
Ci	Cirrus
Cc–Cs	Cirrocumulus–cirrostratus
Cu	Cumulus
CS	Clear-sky conditions
St	Stratus
Sc	Stratocumulus
CBH	Cloud base height (km)
N	Classes of sky conditions in oktas (from 0 for clear-sky conditions to 8 for completely overcast)
Rt	Ratio (R) between observed global irradiance (Fglo) and the modelled clear-sky irradiance (Fglo_CS)
SDt	SD of the measured Fglo in 20 min time intervals (W m−2)
Sft	Scaling factor Sft (Duchon and O'Malley, 1999)
Other symbols and acronyms
φ	Azimuth angle (rad)
Φλ	Photon flux density at wavelength λ ($\mathrm{number}\mathrm{of}\mathrm{photons}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{nm}}^{-\mathrm{1}}$)
λ	Wavelength (nm)
ρ	Air density (kg m−3)
θ	Zenith angle (rad)
θz	Solar zenith angle (rad)
a	Empirical coefficient from Ehnberg and Bollen (2005); Table S1
a0	Empirical coefficient from Ehnberg and Bollen (2005); Table S1
a1	Empirical coefficient from Ehnberg and Bollen (2005); Table S1
a3	Empirical coefficient from Ehnberg and Bollen (2005); Table S1
AF(λ)	Actinic flux for wavelength λ ($\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{nm}}^{-\mathrm{1}}$)
APE	Average photon energy (eV)
APEdif	Average photon energy for diffuse radiation (eV)
APEdir	Average photon energy for direct radiation (eV)
APEref	Average photon energy for reflected radiation (eV)
c	Speed of light (m s−1)
Cp	Isobaric specific heat of dry air (1005 $\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$)
dif	Diffuse
dir	Direct
Fglo	Global broadband irradiance; $F_{\text{glo}}=F_{\text{dir}}+F_{\text{dif}}$ (W m−2)
Fdif	Diffuse broadband irradiance (W m−2)
Fdir	Direct broadband irradiance (W m−2)
Fref	Reflected broadband irradiance (W m−2)
Fdir,dif,ref(λ)	Spectral irradiance as a function of λ ($\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{nm}}^{-\mathrm{1}}$)
h	Plank constant (J s)
ref	Reflected
L	Empirical coefficient from Ehnberg and Bollen (2005); Table S1
q	Electron charge (C)
$R(\mathit{\lambda },\mathit{\theta },\mathit{\phi })$	Radiance at wavelength λ from zenith and azimuth angles θ and φ ($\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{nm}}^{-\mathrm{1}}$)

Appendix B: Cloud type validation

The validation was conducted in two subsequent steps. In the first step the automatized cloud classification (based on Duchon and O'Malley, 1999; including lidar cloud base height) was compared to the visual cloud classification based on sky images collected during 1 month of a field campaign.

The second validation step involved the recently published method discussed by Ylivinkka et al. (2020), which is based on the same methodological approach used in this study: the application of the Duchon and O'Malley (1999) classification improved by the knowledge of the CBH. Thus, the aim of the second step was to determine the degree of consistency between the two approaches that were developed simultaneously and independently in two different regions of the globe.

Both the two validations were evaluated by means of a confusion matrix, a special kind of contingency table, with two dimensions and identical sets of “classes” in both of them. From the confusion matrix the balanced accuracy was computed as follows:

$$\begin{array}{}\text{(B1)}& \text{Balanced\hspace{0.17em}accuracy}={\displaystyle \frac{\text{Sensitivity}+\text{Specificity}}{\mathrm{2}}},\end{array}$$

where the Sensitivity describes the true positive rate (the number of correct positive predictions divided by the total number of positives) and the Specificity describes the true negative rate (the number of correct negative predictions divided by the total number of negatives). The balanced accuracy is especially useful when the investigated classes are imbalanced; i.e. one of the classes appears a lot more often than the other, a condition useful for cloud classification (García et al., 2009).

Figure B1 SD–R diagram (a–g) and the corresponding sky images for the February–March 2017 field campaign: (a) CS conditions, (b) Ci clouds, (c) Cu clouds, (d) Ac clouds, (e) Sc clouds, (f) As clouds and (g) St clouds.

B1 Visual cloud classification

Sky images were collected during 1 month (13 February–9 March 2017) using a sky view camera (GoPro Hero4 Session installed on the U9 roof), characterized by a field of view of 95^∘ × 123^∘; the camera was oriented south (each day manually) with the same declination of the shadow band applied to the DPA154 global radiometer (for diffuse broadband irradiance measurements, Sect. 2.1.2); sky images were taken with 1 min time resolution. Visual classification of sky images, based on the principles of cloud classification published in the International Cloud Atlas (WMO). Figure B1 reports an example of the SD–R diagram (Sect. 2.3.2) with the CBH for each sky or cloud condition with the corresponding image.

To test the performance, 869 sky images were analysed, and the cloud type was determined through visual inspection. From the visual classification and the automatized one (Table 1) the following confusion matrix (Table B1) was created. The highest balanced accuracy was found for St data (95 %), while the lowest (50 %) was found for mixed cloud types (Cc–Cs) whose absolute number of cases, however, was ∼ 0.6 % of the total, probably biassing the obtained accuracy; the same happened for Cu and Ac. Overall, five of eight classes were above 68 % of balanced accuracy, while the overall balanced accuracy was 80 %, underlying the reliability of the classification algorithm, allowing for studying the impact of clouds on the LAA HR with a sufficient grade of certainty.

Table B1 Confusion matrix and balanced accuracy for each cloud type classified visually and following the algorithm reported in Table 1 within the present work.

Table B2 Final criteria adopted for cloud classification in Ylivinkka et al. (2020). Ns here represents nimbostratus, while GRE stands for global radiation enhancement.

Table B3 Cloud class homogenization adopted for comparison purposes (merged cloud type) between the present study's cloud classification and the one reported in Ylivinkka et al. (2020).

Table B4 Confusion matrix and balanced accuracy for each cloud type classified using the algorithm reported in the present study and the one reported in Ylivinkka et al. (2020).

B2 Intercomparison with Ylivinkka et al. (2020)

The second validation step involved the recently published method discussed by Ylivinkka et al. (2020), which is based on the same logical approach followed in our work: the application of the Duchon and O'Malley (1999) classification improved by the knowledge of the CBH. For this purpose, the classification scheme of Ylivinkka et al. (2020) is resumed in Table B2 following the nomenclature used in the present work. It is necessary to underline that the cloud classes determined in the work of Ylivinkka et al. (2020) differ from those reported in the present work. Particularly, while both approaches enabled the Cu, St and Sc classification, some of the cloud classes were merged in the Ylivinkka et al. (2020) study: CS and Ci (CS + Ci); Ac and As (Ac + As); and a mixed situation composed by Ci, Cc and Cs (Ci + Cc + Cs). In addition they introduced the classes Cu + GRE and Ci + GRE to account for global radiation enhancement (GRE) due to this cloud types; a possible explanation for such a difference with respect to present work could be hidden in the different latitude at which the two algorithms were developed, which is a parameter able to affect the solar zenith angle and the sun light interaction with clouds. A detailed investigation of this difference is beyond the aim of the present work. However, it is necessary to account for the classification differences in order to properly merge cloud classes with similar features to finally perform a comparison between the two methods. The cloud class homogenization is summarized in Table B3, while the final intercomparison is reported in Table B4. The confusion matrix (Table B4) revealed a global balanced accuracy of 90 %, making the two methods comparable, despite the aforementioned differences. The highest accuracy (100 %) was obtained for CS, followed by Ac + As (99 %); Cu, St and Sc reached values of 94 %, 93 % and 86 %, respectively. The lowest performance was reached for Ns, whose presence cannot be detected in the present study, generating a false positive signal in the Ac + As class; however, due to the very low number of Ns cases (1.8 %), its impact on the cloud classification can be neglected. Overall, even the second validation step pointed out the reliability of the results obtained in the present work.