Towards reliable retrievals of cloud droplet number for non-precipitating planetary boundary layer clouds and their susceptibility to aerosol

1 Introduction

Aerosol-cloud interactions constitute one of the most uncertain drivers of anthropogenic climate change (Seinfeld et al., 2016; IPCC, 2021). Central to these interactions are the variations in cloud condensation nuclei (CCN; the subset of aerosol particles which droplets in ambient clouds form upon) concentrations from anthropogenic activities, as they modulate the N_d in clouds and the cloud optical depth, thus, affect the atmospheric radiation budget and climate. CCN activate into cloud droplets in supersaturated air generated within cloudy updrafts, and the supersaturation level required for each particle to act as a CCN is well described by Köhler theory (Köhler, 1936), modified to account for the effects of organic compounds and multiple hygroscopic species (e.g., Petters and Kreidenweis, 2007). Supersaturation in clouds is mostly generated by expansion cooling in updrafts driven by turbulence, convection, or gravity waves (Quaas et al., 2020).

Initial studies by Twomey and Warner (1967) showed that elevated concentrations of CCN in clouds increases N_d, reduces their size and thus enhances the amount of solar radiation reflected back to space. This effect, which tends to cool climate, is known as the “aerosol indirect effect” (Twomey, 1991) or “cloud albedo effect” (IPCC, 2021) on climate. Later studies showed that this N_d increment exists, but varies considerably depending on the supersaturation level that develops in clouds–which in turn depends on the vertical velocity (w), and amount of aerosol/CCN present (Ghan et al., 1998; Feingold et al., 2001; Reutter et al., 2009). It is now generally accepted that droplet formation can range from “aerosol-limited”, when supersaturation is high and N_d is very sensitive to aerosol changes, to “velocity-limited”, when the supersaturation is so low that N_d is insensitive to aerosol changes (Reutter et al., 2009). However, more recent studies showed that when approaching “velocity limited” conditions, N_d tends to come near an asymptotic upper limit, the so-called limiting droplet number (N_d^lim), which primarily depends on the intensity of turbulence (Bougiatioti et al., 2020; Kacarab et al., 2020; Georgakaki et al., 2021).

Changes in N_d, apart from modifying the average size of cloud particles, also affect the shape of the droplet size spectrum and the resulting cloud radiative properties. The first generation of global climate models (GCM) with explicit consideration of N_d initially assumed that all droplets have the same diameter (i.e., monodisperse droplet population) which can directly be linked to the effective radius and cloud optical depth. Liu and Daum (2002) showed that this assumption leads to important biases in the calculations of optical depth and cloud radiative forcing, as monodisperse droplets scatter differently than a distribution with the same N_d and liquid water content (LWC). Known as “dispersion effect”, this bias can be mitigated through the use of the effective radius factor (β) so that when it is combined with the average droplet radius and cloud LWC gives the correct cloud optical depth. Several formulations have been proposed in the literature to represent the relationship between β and cloud microphysical quantities resolved in models—N_d and LWC. Expressions proposed by Martin et al. (1994), Liu and Daum (2002), Peng and Lohmann (2003) and Rotstayn and Liu (2003) were based on observations, while more recent formulations are derived from the adiabatic parcel theory (Liu et al., 2014; Chen et al., 2016; Chen et al., 2018). All together, these expressions show a positive correlation between β and N_d for aerosol- and velocity-limited regimes.

Despite decades of research, the imprint of aerosol effects on cloud albedo and climate still remains highly uncertain (Seinfeld et al., 2016; Quaas et al., 2020; IPCC, 2021), owing to the highly-buffered, nonlinear and multiscale nature of clouds and their interactions with dynamics, radiation and aerosol (Stevens and Feingold, 2009). Therefore, there is a pressing need to have extensive global datasets of N_d aerosol concentration, and cloud-scale dynamics (w especially) to understand their relative importance (e.g., Rosenfeld et al., 2016). Recently, Quaas et al. (2020) showed that satellite remote sensing is the only approach that offers the potential of obtaining global datasets with frequent coverage; current retrieval algorithms, however, carry many uncertainties and require constraints that can only be addressed with in situ and/or ground-based remote sensing observations. Particularly important is constraining w, which is seldomly evaluated in models (Sullivan et al., 2016). Furthermore, w can change along with aerosol and affect the response of clouds, as demonstrated recently by Bougiatioti et al. (2020) and Kacarab et al. (2020) for very different cloud systems, showing that a co-variability of aerosol concentration with w magnifies the inherent response of droplets to aerosol variations by up to 5 times.

It has long been known that w exhibits significant spatio-temporal variability even at the scale of single clouds. This presents a challenge about which velocity to use to link aerosol with N_d. Morales and Nenes (2010) proposed using a “characteristic” velocity, computed from the probability distribution function (PDF) of w in the airmass or cloud of interest, to provide a PDF-averaged N_d, which in turn is used to represent cloud-scale N_d. The characteristic velocity approach can be automated and easily applied to data either from in situ probes or ground-based remote sensing systems (Doppler lidars/radars). The calculation of N_d by using this characteristic velocity and in situ ground level data has been shown to agree with in situ measurements at cloud base level of cumulus and stratocumulus clouds in numerous studies to date (Meskhidze, 2005; Fountoukis and Nenes, 2007; Kacarab et al., 2020).

A more empirical approach to determine cloud base updrafts, based on cloud base height, has been proposed by Rosenfeld et al. (2016), using ground-based Doppler lidar measurements along with satellite data from the Visible Infrared Imaging Radiometer Suite (VIIRS). Their approach gave closure to within ± 25–30%, when applied to PBL convective clouds, non-obscured by upper layer clouds or including semitransparent cirrus. Although skillful, Rosenfeld et al. (2016) point out that this methodology cannot be applied globally, because distribution of CCN and cloud base updrafts are lacking.

In this study we improve an existing remote sensing approach for the retrieval of N_d in non-precipitating Planetary Boundary Layer clouds (PBLC), by considering the spectral dispersion of cloud droplets through β. The retrieval is evaluated using published dispersal relationships by comparing the N_d obtained from a combination of in situ observations, ground-based remote sensing data and parcel theory. The same data are then used to determine an optimal dispersal- N_d relationship that minimizes the closure error between retrieved and the estimated in situ N_d. During the optimization process, the cloud-scale w distributions, CCN spectra, and cloud maximum supersaturation, along with the N_d, are obtained. We apply this methodology to a highly variable semi-urban environment and discuss future directions and perspectives.

2 Methodology

Our methodology is based on the concurrent use of remote sensing, model data and expressions for β, applied to non-precipitating PBLCs in well-mixed PBLs and associated cloud-base vertical velocities. The remote sensing retrievals of N_d^sat are compared against estimations of in situ N_d obtained from the application of the Nenes and Seinfeld (2003) activation parameterization, with extensions introduced by Fountoukis and Nenes (2005), Barahona et al. (2010) and Morales Betancourt and Nenes (2014). As input data to the droplet parameterization, we use the vertical velocity variation inside the PBL provided by the HALO wind lidar (see Session 2.3.3), as well as the aerosol size distribution and chemical composition at ground level, based on a Scanning Mobility Particle Sizer (SMPS) and a Time-of-Flight Aerosol Chemical Speciation Monitor (ToF-ACSM), respectively.

The approach used to estimate in situ N_d makes the important assumption that calculations of activated N_d using observations of vertical velocity, aerosol and chemical composition are representative of non-precipitating PBLCs. Numerous studies to date have shown that this approach, when carefully applied, can provide estimates of in-cloud droplet number for non-precipitating boundary layer clouds within a constrained uncertainty. For example, Conant et al. (2004) performed a droplet closure study, using a dataset of 21 cumulus clouds by the Center for Interdisciplinary Remotely-Piloted Aircraft Studies (CIRPAS) Twin Otter aircraft during NASA’s Cirrus Regional Study of Tropical Anvils and Cirrus Layers–Florida Area Cirrus Experiment (CRYSTAL-FACE), and evaluated the performance of the Nenes and Seinfeld (2003) parameterization against in situ observations. They found that the measured droplet concentration close to cloud base was within 15% of the predicted N_d. Meskhidze (2005), evaluated the revised parameterization of Fountoukis and Nenes (2005), against data from cumuliform clouds collected during NASA’s Cirrus Regional Study of Tropical Anvils and Cirrus Layers–Florida Area Cirrus Experiment (CRYSTAL-FACE, Key West, Florida, July 2002). On average, the predicted N_d was within 20% of observations at the base of cumuliform clouds, and within 30% of observations at different altitudes throughout the stratiform clouds; at the same time the observational uncertainty of N_d was estimated to be 30%. Fountoukis and Nenes (2007) evaluated the parameterization of Nenes and Seinfeld (2003) and Fountoukis and Nenes (2005) with data from 27 cumuliform and stratiform clouds sampled during the 2004 International Consortium for Atmospheric Research on Transport and Transformation (ICARTT) experiment. Droplet closure was achieved to within 20% of the measurements. Morales et al. (2011) also showed that considering lateral entrainment effects diagnosed from the diabaticity (diagnosed from the liquid water vertical profile) provides N_d to within 31% for non-precipitating convective clouds. Georgakaki et al. (2021), using a combination of wind lidar data and in situ aerosol measurements, predicted N_d for orographic clouds in the Swiss Alps during the RACLETS campaign that agreed with in situ observations to within 25%. We therefore apply this methodology to estimate droplet number in the clouds considered here, assuming an inherent method uncertainty of 30%.

The satellite retrievals of N_d^sat of the PBLCs were performed using the algorithm of Bennartz (2007), which uses the principal cloud parameters (e.g., cloud effective radius (r_eff) and the cloud optical thickness (τ), along with a constant value of β. The explicit expression of β, which is related to N_d^sat, can be obtained by modifying the Bennartz (2007) algorithm and solving this form numerically. Then, we examine the closure of the satellite retrieval algorithm by using the β-N_d relationships, called “β-expressions”, found in literature (M94—Martin et al., 1994; RL 03—Rotstayn and Liu, 2003; PL03—Peng and Lohmann, 2003; Z06—Zhao et al., 2006; GCMs—Rotstayn and Liu 2009; F12—Freud and Rosenfeld, 2012). Additionally, we propose a new β-expression which minimizes the closure error between N_d—N_d^sat.

We obtained the parameters r_eff, τ and the cloud top pressure (P_ct) from the Optimal Cloud Analysis (OCA) product (EUMETSAT, 2015) which is based on geostationary meteorological Meteosat Second Generation (MSG) satellite data. Moreover, we obtained the atmospheric pressure-temperature profiles from the fifth-generation European Center for Medium-range Weather Forecasts (ECMWF) reanalysis model data (ERA5) (Hersbach et al., 2018), in order to calculate the cloud top temperature (T_ct) and the cloud top height (H_ct) based on the cloud-top pressure P_ct.

A detailed description of the proposed data analysis procedure is provided in Supplementary Figure S1.

2.1 Experimental site

Our observations took place at the Athens Metropolitan Area (AMA) in the Attica peninsula in Greece. The AMA is surrounded by four mountains: Egaleo (468 m), and Parnitha (1,413 m) to the northwest, Penteli (1,109 m) to the north and Hymettus (1,026 m) to the east, with a major opening to the sea on the south-western part (Saronikos Gulf). The AMA hosts a densely populated urban area of about 3.8 million inhabitants within about 2,928.7 km². The study region, characterized by a complex topography as shown in Supplementary Figure S3, is located at the National Center for Scientific Research-Demokritos (NCSRD) station (37.995^o N, 23.816^o E, at 270 m above sea level-asl.) within the AMA on the foothills of mount Hymettus (Vratolis et al., 2020), about 8 km to the north of the Athens city center. NCSRD is a member of the World Meteorological Organisation (WMO) Global Atmospheric Watch (GAW) network, the Aerosol, Clouds and Trace Gases Research Infrastructure (ACTRIS) and the PANhellenic infrastructure for Atmospheric Composition and climatE chAnge (PANACEA).

The selected experimental site offers a great variety of aerosol types and meteorological conditions (Kostenidou et al., 2015; Vasilatou et al., 2017). The prevailing north winds transport continental polluted air masses to the AMA, mixed with local emissions along the Balkans-Aegean Sea axis (Diapouli et al., 2014; Vratolis et al., 2019). South winds are usually associated with the transport of marine aerosols, and, occasionally, with mineral dust from the Saharan desert (Diapouli et al., 2017; Soupiona et al., 2020), while under stagnant atmospheric conditions, locally-emitted anthropogenic aerosols dominate aerosol concentrations in the PBL (Kalogridis et al., 2018).

2.2 Case studies

We focus on the period between March and May 2020, where 5-day of data with PBLC formation are analyzed: 18, 28 March, and 01, 02, and 07 of April. These cases are summarized in Table 1.

TABLE 1. Presence of clouds for 5 days:18, 28 March, and 01, 02, and 07 of April 2020.

2.3 Instrumentation

2.3.1 Scanning Mobility Particle Sizer

The SMPS consists of an electrostatic classifier (TSI Inc. Model 3,080), a cylindrical differential mobility analyzer column (TSI Inc., model 3,081) and a condensation particle counter (TSI Inc. model 3,772). The SMPS provided the aerosol size distributions with a temporal resolution of 5 min and was operated at an aerosol flow rate of 1 L per minute (lpm) and a sheath flow rate of 5 lpm, with the ability to measure the particle sizes from 10 to 550 nm. Sample and sheath flows were dried to a relative humidity (RH) lower than 40 % using Nafion driers before their introduction to the SMPS. Data acquisition and analysis were performed using the non-commercial TROPOS-SMPS data evaluation software (Wiedensohler et al., 2012), assuming that the particles are spherical so that their mobility and geometric diameters are equal. To achieve the highest measurement accuracy, the technical recommendations and quality control procedures proposed by Wiedensohler et al. (2012) were followed. The instrument was also calibrated against a reference SMPS system at the World Calibration Centre for Aerosol Physics (WCCAP) in 2018 (Bousiotis et al., 2021), exhibiting a counting accuracy within 10% for the size range 30–550 nm against a reference system under controlled laboratory conditions (Wiedensohler et al., 2012). The SMPS is also calibrated at the NCSRD station using polystyrene spheres with a size of 200 nm (Vratolis et al., 2020). The 5-min data were aggregated into time bins of 15 min, to be used as input to the droplet activation parameterization. The information about the accuracy of the aerosol particle size distribution can be found in the Supplement (c.f. Supplementary Figure S4), where the daily plot of the distributions is presented.

2.3.2 Time-of-flight-aerosol chemical speciation monitor

The ToF-ACSM provided by Aerodyne Inc. (Ng et al., 2011; Fröhlich et al., 2013) measures the non-refractory submicron aerosol mass and chemical composition (ammonium, sulfate, nitrate, chloride, and organics) in real-time (here, with 60 min time resolution), with a measurement uncertainty on the order of 30% (Budisulistiorini et al., 2014). The aerodynamic particle focusing lens is combined with particle flash vaporization in high vacuum and a hard electron impact ionization. The detection of the resulting ions is done by time-of-flight mass spectrometry. At ambient atmospheric conditions, the lens system focuses on the detection of submicron aerosols of ∼40–1,000 nm vacuum aerodynamic diameters. A collection efficiency factor CE = 0.5 (Zografou et al., 2022) was used to account for the fraction of the non-refractory particles that bounce off the vaporizer and are not detected. During the measurement period, the ambient air was drawn into the ToF-ACSM through a PM_2.5 aerosol inlet. The inlet aerosol flow was dried to a RH lower than 40% by a Nafion dryer. More detailed information about the daily chemical composition can be found in the Supplement (c.f. Supplementary Figure S5).

2.3.3 Wind lidar system (HALO)

A pulsed Doppler scanning lidar system (StreamLine Wind Pro model, HALO Photonics) emitting at 1.565 μm was deployed at the NCSRD site by the Finnish Meteorological Institute (FMI). The system was operated in the stare vertical azimuth display mode. The laser backscattered signals mainly come from the laser-matter interactions from the micron-sized aerosols, which act as tracers of the wind flow over the sampled area. These time-resolved signals are then collected by a receiving telescope, digitized and processed to obtain the vertical profiles of the radial wind along the laser beam direction (Henderson et al., 2005) by measuring the Doppler frequency shift of the backscattered pulses compared to the emitted ones (Newsom and Krishnamurthy, 2020). The range resolution of the measurements is 30 m, the temporal resolution is 20 s for the stare mode, while the maximum measurement range is 2–3 km depending on the atmospheric aerosol load. The useful lidar signals are those which present a signal-to-noise-ratio (SNR) greater than 1.003. The corresponding uncertainty of the retrieved wind velocity range is 0.1–0.5 m s⁻¹ for stare mode (Newsom and Krishnamurthy, 2020). From non-precipitating PBLC a rainfall filter was applied, by excluding possible rainy moments where the fall velocity was >4 m s⁻¹. Typical PBL heights over the AMA (around 12:00 UTC) are of the order of 1,617 ± 324 m asl (Kokkalis et al., 2020). Therefore, as the retrieved vertical velocity data were confined within a height of 1,020 ± 60 m asl. the aerosols were always inside the PBL during daytime, and the updrafts were sampled at this height.

The vertical velocity data were then sampled in segments within a 4-h time moving window, and the positive updrafts were fit to a half-Gaussian PDF with a zero mean and standard deviation σ_w as follows:

The value of σ_w was determined every 15 min, for which the characteristic vertical velocity w* was obtained and used for calculating the droplet number concentration by the parameterization, as explained in section 2.4.3.

Additionally, we examined how the sampling time with a 1-, 2- and 4-h window affects the calculated value of σ_w; in all cases the value of σ_w was found not to change the N_d value more than the inherent uncertainty of 30% for the in situ N_d estimation. Based on this, we chose the 4-h segments to ensure sufficient sampling of the PDF (about 400 updraft velocities per PDF).

2.4 Modelling and data preprocessing

2.4.1 ECMWF reanalysis (ERA5)

We use the ECMWF Reanalysis fifth Generation (ERA5) hourly temperature and specific humidity data on 37 pressure levels (1,000 hPa to 1 hPa) at a resolution of 0.25 × 0.25 (Hersbach et al., 2018). For this study, we estimate the T- and q-profiles at the site using a weighted average of the surrounding grid-points. Moreover, a linear interpolation for every pressure level was applied on each variable to produce a 15-min timeseries that is compatible with the satellite timeseries. To link P_ct with the H_ct and T_ct, used subsequently for cloud detection and inversion, we apply the hypsometric equation to the grid level quantities of T and q to express their vertical profiles.

2.4.2 Hygroscopicity and CCN concentration

The particle hygroscopicity parameter (κ) (Petters and Kreidenweis, 2007) characterizes the ability of particles to absorb water vapor during the formation of cloud droplets, and is required for calculating CCN concentrations. The value of κ is estimated from ACSM measurements as follows: the ion concentrations of ammonium, sulfate, nitrate, chloride, and organics are converted to inorganic salt mixtures using the simplified ion-pairing scheme of Q-AMS model (Gysel et al., 2007), and the mole number of NH₄NO₃, NH₄HSO₄ (NH₄)₂SO₄, H₂SO₄ and organic aerosols, respectively. The value of κ is calculated using the volume fraction-averaged hygroscopicity of each organic and inorganic salt in the aerosol (Petters and Kreidenweis, 2007), with κ being 0.68 for NH₄NO₃, 0.56 for NH₄HSO₄, 0.53 for (NH₄)₂SO₄, (Petters and Kreidenweis, 2007), 0.97 for H₂SO₄ (Biskos et al., 2009) and 0.1 for organic aerosol (Duplissy et al., 2011). Aerosol types often have characteristic values of κ; for example, low values between 0.1–0.2 are associated with organic-rich aerosol (such as biogenic or biomass burning aerosol), ∼0.35 with continental aerosol–which is often a ∼50/50 mixture of organic species and inorganic salts with a hygroscopicity close to that of (NH₄)₂SO₄ (∼0.6) (Petters and Kreidenweis, 2007). The κ timeseries is calculated on an hourly basis. We also examined the importance of constraining the uncertainty of organic κ; which equals δκ_org = ± 0.064 (Duplissy et al., 2011), and we found that the propagated error into the overall κ gives a total uncertainty of ± 0.035.

The CCN at a water vapor supersaturation level s is determined by integrating the particle size distribution from the smallest activated particle with a “critical dry diameter”, D_cr, up to the largest particle size measured by the SMPS, i.e., $CCN\left(s\right)={\int }_{D_{cr}}^{\infty }n\left(D_p\right)d{D}_p$, where n(D_p) is the aerosol size distribution measured by the SMPS. D_cr is calculated from κ-Kölher theory (Köhler, 1936), as $D_{cr}={\left(\frac{4A^3}{27\kappa s^2}\right)}^{\frac{1}{3}}$, where $A=\frac{4M_w\sigma }{RT{\rho }_w}$ is the Kelvin parameter, M_w (kg mol⁻¹) is the molar mass of water, σ (J m⁻²) is the surface tension of the solution droplet, R (J mol⁻¹ K⁻¹)is the universal gas constant, T (K) is the ambient temperature, and ρ_w is the density of liquid water.

2.4.3 Droplet activation parameterization

The physically-based aerosol activation parameterization (Nenes and Seinfeld, 2003), with extensions introduced by Fountoukis and Nenes (2005), Barahona et al. (2010) and Morales Betancourt and Nenes (2014) is used to calculate the droplet number that would form in PBLCs, from knowledge of the wind vertical velocity and the aerosol characteristics. The parameterization solves the equations that describe droplet formation in an ascending air parcel containing aerosols and water vapor, specifically at the point where supersaturation becomes maximum, s_max; N_d is then equal to the number of CCN with critical supersaturation less than s_max (Nenes et al., 2001). The parameterization uses as inputs the pressure and the temperature, the aerosol size distribution data from the SMPS, the hygroscopicity parameter κ derived from the ACSM data and the updraft velocities from HALO. Given that aerosol measurements are carried out at ground level, we convert the concentrations to cloud top conditions using the ideal gas law.

As the droplet number depends strongly on the vertical velocity, which in turn varies considerably in the PBL over time, we compute the N_d that characterizes PBLCs–being a weighted average over the PDF of vertical velocities, by using the characteristic velocity, w* according to the approach of Morales and Nenes (2010). This approach is valid for boundary layers that are not influenced by deep convection, hence are on average characterized by a mean velocity of a few cm s⁻¹ (the typical speed of boundary layer ascent/descent over a diurnal cycle). The PDF-averaging approach with a Gaussian PDF and the usage of the w*, has been shown to successfully predict cloud-scale values of N_d in field studies for cumulus and stratocumulus clouds (e.g., Conant et al., 2004; Meskhidze, 2005; Fountoukis and Nenes, 2007; Kacarab et al., 2020; Georgakaki et al., 2021).

Another significant parameter taken into account was the role of the lateral entrainment, meaning the mixing of dry air in the ascending cloudy air parcel. Entrainment effects can be significant for cumulus and convective clouds (much less for stratocumulus), affecting the vertical distribution of liquid water and number of droplets (Morales et al., 2011). Given that the clouds studied here are cumuliform, we consider the modest impact of entrainment effects on droplet number using the approach of Barahona and Nenes (2007) and Morales et al. (2011); this entails using an adiabatic activation parameterization, and reducing the vertical velocity that causes activation by a factor, called “entrainment parameter”, that depends on the entrainment rate (which in effect captures the reduction in maximum supersaturation from mixing of dry air, compared to adiabatic droplet formation). Barahona and Nenes (2007) and Morales et al. (2011), concluded that the best approach for the entrainment parameter is based on the adiabaticity (i.e., deviation of the vertical water profile from the adiabatic value). Morales et al. (2011) calculated the entrainment parameter for a number of clouds sampled during the CRYSTAL-FACE (Key West, Florida, July 2002) and CSTRIPE (Monterey, California, July 2003) missions–which on average gave a parameter value of 0.68 (see data in Table 2 of Morales et al., 2011). Assuming this mixing parameter is applied uniformly to the vertical velocity distribution, we therefore consider the effects of lateral entrainment in the cumuliform cloud sampled by multiplying the characteristic velocity for adiabatic activation w* by 0.68.

TABLE 2. β-expressions based on literature according to: 1) Martin et al. (1994)—M94, b) Rotstayn and Liu (2003)—PL03, c) Peng and Lohmann, (2003)—PL03, Zhao et al. (2006)—Z06, d) Freud and Rosenfeld (2012)—F12 and e) GCMs studies–Rotstayn and Liu (2009), respectively.

According to Morales and Nenes (2010), the calculation of w* (for diabatic activation) is obtained by multiplying σ_w with a parameter λ that depends on the type of aerosol found in the cloudy updraft. For continental aerosol with a concentration range between 1,000–10000 cm⁻³, λ equals to 0.67 (see Figure 2A of Morales and Nenes (2010)) with an uncertainty 0.02, which we apply to our calculations. Given the above considerations, the characteristic vertical velocity w* used for calculating droplet formation (and used to estimate in situ N_d) is given by $w^{*}=0.68\lambda {\sigma }_w=0.68\times 0.67{\sigma }_w=0.456{\sigma }_w$. Moreover, the relevant uncertainty of the characteristic velocity w^* can be estimated as $\delta w^{*}=0.456\delta {\sigma }_w$, where the error $\delta {\sigma }_w$ is estimated as the standard deviation of σ_w in the fitting process.

FIGURE 2. The τ, r_eff provided by OCA against with T_ct, which is calculated via the ERA5 T-Profile, for the clouds studied; the labeling refers to the serial number of each cloud (c.f. Supplementary Figure S1).

Another important parameter to consider is the “limiting droplet number”, N_d^lim, at cloud top, which helps to estimate the degree to which clouds are susceptible to changes in aerosol content. N_d^lim is obtained from the relationship given by Georgakaki et al. (2021):

where N_d^lim is expressed in cm⁻³, and σ_w in m s⁻¹.

Figure 1A presents the timeseries of total aerosol concentration (N_aer) integrated over the SMPS size distributions, and σ_w, with the symbols colored by κ, for the period 18 March to 07 April 2020, while Figure 1B presents the corresponding N_d values with the symbols colored by N_d/N_d^lim.

FIGURE 1. (A) Timeseries of N_aer versus κ and σ_w values, for the period 18 March to 07 2 April0,230. (B) The corresponding in situ derived N_d values colored by the ratio N_d/N_d^lim as derived by the parameterization.

The periods of weaker turbulence (low σ_w, Figure1A) forces N_d to persistently approach N_d^lim (e.g., Figure1B, 01 April 2020). As a result, the droplet formation is strongly velocity limited, and tends to be insensitive to aerosol variations. The opposite is seen when N_d/N_d^lim < 0.5, as droplet formation becomes aerosol limited and less sensitive to σ_w variations.

2.4.4 Effective radius factor β

The use of β, which is defined as the ratio of the r_eff to the volume mean radius of cloud droplets, was introduced to mitigate the bias known as “dispersion effect”, related to the width of the cloud droplet size spectrum, and combined with the average droplet radius, to give the correct cloud optical depth. According to Liu and Daum (2002), β is related to the relative dispersion of the droplet spectrum (ε) as:

Ιn the literature there are several studies attempting to quantify ε. Martin et al. (1994), using aircraft observations of the microphysical characteristics of warm marine stratocumulus clouds from the eastern Pacific, South Atlantic, subtropical regions of the North Atlantic, and the sea areas around the British Isles, proposed an expression of $\epsilon =0.000574{{\rm N}}_d+0.2714$, while, Rotstayn and Liu (2003) based on the observations of Liu and Daum (2002), suggested an expression of $\epsilon =1-0.7e^{-0.003{{\rm N}}_d}$. Using aircraft data collected from a region influenced by polluted marine aerosol conditions, Zhao et al. (2006) suggested a value of $\epsilon =0.4$. Using Eq. 3, we inserted the above-mentioned expressions of ε to obtain the β-expressions (c.f. Table 2). Peng and Lohmann (2003) suggested an expression of $\beta =1.18+4.5{10}^{-4}{{\rm N}}_d$ using aircraft measurements over the ocean, near coastal areas, and far from the interior of continents, through the coast of Nova Scotia in the Canadian and the United States Arctic.

Freud and Rosenfeld (2012), using measurements of cloud droplet size distribution in various locations, found an average value of β equal to 1.08, as indicated by the line of best fit between $r_{eff}$ and mean volume radius. Rotstayn and Liu (2009) pointed out that most GCM studies adopt a β of 1.1 (c.f. Table 2).

2.4.5 Satellite remote sensing—Optimal Cloud Analysis product and droplet number

The OCA product determined from EUMETSAT (EUMETSAT, 2015) provides τ, r_eff and P_ct accompanied with their errors δτ, δr_eff and δP_ct with a spatial resolution of 3.6 km × 4.6 km for the study region, as derived from SEVIRI onboard the geostationary meteorological satellite (METEOSAT) with a temporal resolution of 15 min.

The OCA algorithm uses a synergy of data to calculate the cloud parameters; however, there are always assumptions which have consequent limitations affecting the accuracy of r_eff and τ. According to the description of the OCA algorithm (EUMETSAT, 2015), the algorithm utilizes the land surface reflection, along with other model parameters derived from the Numerical Weather Prediction sources.

In Figure 2 we summarize the τ, T_ct, and r_eff which were used as input to the satellite algorithm (c.f. Eq. 6) to retrieve the N_d^sat data. Additionally, we exclude from the dataset cases that involve high or multi-layer clouds, and isolate the PBL clouds by excluding clouds with P_ct lower than 800 hPa. Moreover, we applied a parallax correction, according to Koenig (2020), for each detected cloud based on H_ct in order to calculate the true geolocation of each cloud and then, we isolated the clouds that are collocated over our measuring site. This step helps to correct the displacement that is created due to the satellite angle of view, along with H_ct.

According to Zhu et al. (2018), who further developed the Bennartz (2007) algorithm, N_d^sat, can be determined as:

in which $c\left(c_w\right)=\frac{5c_w}{4{\pi }^2{Q}_{ext}{\rho }_w}$ , and $Q_{ext}$ is the Mie extinction efficiency factor, which is equal to 2 for large droplets of diameters greater than 5 μm (van de Hulst, 1958). The condensation rate c_w (Brenguier et al., 2000) provides the liquid water released in an adiabatic updraft per meter of ascent, and is given by $c_w=0.0016+4.86{10}^{-5}{T}_{ct}-3.42{10}^{-7}{T_{ct}}^2$ (Zhu et al., 2018), where c_w is in g m⁻³ m⁻¹ and T_ct expressed in ^oC. The water density ρ_w is considered constant, equal to 997 kg m⁻³, while T_ct is used to characterize the in-cloud environment.

Moreover, we calculated the propagation error δN_d^sat of N_d by using the Chain rule to the Eq. 4, which corresponds to the diffusional error due to each variable to the N_d^sat retrieval, as follows:

in which the terms $\frac{\partial N_d^{sat}}{\partial c_w}$, $\frac{\partial N_d^{sat}}{\partial \tau }$ , $\frac{\partial N_d^{sat}}{\partial r_{eff}}$, and $\frac{\partial N_d^{sat}}{\partial \beta }$ are the sensitivity of N_d^sat with respect to τ, r_eff, c_w and β, respectively. The condensation rate error δc_w was estimated from the lapse rate along with the δP_ct. δc_w δτ, δr_eff were estimated for the whole dataset and on average were found to be equal to 6 10^–6 g m⁻³ m⁻¹, 1.07, 0.76 μm, respectively.

We note here that $\frac{\partial N_d^{sat}}{\partial c_w}\delta c_w$ was on average relatively small ( ± 5 cm⁻³), and contributes to ± 1.7% on the total bias of $N_d^{sat}$. Therefore, we decided to omit it from Eq. 5. $\frac{\partial N_d^{sat}}{\partial \tau }\delta \tau$ and $\frac{\partial N_d^{sat}}{\partial r_{eff}}\delta r_{eff}$ were found on average equal to ± 30 cm⁻³, ± 76 cm⁻³, respectively, contributing ± 12% and ± 27%, respectively to the error. Furthermore, we estimated the $\frac{\partial N_d^{sat}}{\partial \beta }$ and found it equal to on average 593 cm⁻³ per unit of β. Since the uncertainty δβ is not available from published literature, we used $\delta \beta$ derived from the optimization process and found it equals to 0.28. Thus, $\frac{\partial N_d^{sat}}{\partial \beta }\delta \beta$ is estimated equal to ± 184 cm⁻³, which contributes ± 57% to the droplet error. This implies that of all parameters considered in this study, optimally constraining β term is of prime importance for the N_d^sat retrieval. The relevant results of the normalized bias of $N_d^{sat}$ regarding δc_w, δτ, δr_eff, and δβ can be found in the Supplement (c.f. Supplementary Figure S13).

For expressions where β depends on N_d^sat, β(N_d^sat), the retrieval Eq. 4 can be modified as follows:

where N_d^sat is determined from the numerical solution of Eq. 6 using the β(N_d^sat) expressions in Table 2. We discard the less reliable retrievals when the droplet uncertainty is significant, which correspond to the solutions of Eq. 6 having δN_d^sat > 600 cm⁻³, δN_d^sat/N_d^sat > 0.5, N_d^sat > 2000 cm⁻³, or N_d^sat<100 cm⁻³.

Finally, we performed closure studies between the accepted solutions of N_d^sat using each literature based β-expression, against estimations of in situ derived N_d from the parameterization (section 2.4.3). By using the M94, RL03, PL03, Z06, GCMs, and F12 expressions (Table 2), the corresponding averaged mean normalized bias (MNB) between N_d^sat and estimations of in situ N_d is equal to -17.37% ± 32.66%, 51.34% ± 69.25%, 23.51% ± 56.09%, -21.25% ± 24.91%, -28.80% ± 22.52%, and -31.99% ± 21.51%, respectively (c.f. Figure 4B and Table 3).

FIGURE 4. (A) Droplet number closure between N_d and N_d^sat using the OPT compared to N_d in respect of β_opt on the colorscale; (B) The MNB of the closure of N_d^{sat _} N_d by using each β-expression: (I) M94, ii) RL03, iii) PL03, iv) F12, v) Z06, vi) GCMs, and vii) OPT, respectively.

TABLE 3. Statistics of the performance of the closure study of N_d^{sat _} N_d for each β-expression used: OPT, RL03, M94, Z06, PL03, GCMs, and F12.

Therefore, in the case of using a constant value of β, such as Z06, GCMs, and F12, the N_d^sat values tend to be underestimated, since the estimated mean bias is of the order of 28%, while the standard deviation is reduced by 23% on average. On the other hand, by using the PL03 expression, the N_d^sat is overestimated, although comparable with those values derived when expressions of constant value of β are used (Z06, GCMs, and F12), with increased standard deviation values. In case of using the M94 explicit relation, N_d^sat is underestimated, but the mean bias is reduced by almost a factor of two, but with an increase in the standard deviation. Usage of the RL03 relation provides N_d^sat values that are considerable overestimated along with their standard deviation (c.f. Supplementary Figure S12), while the MNBs presented in box plots can be found in Figure 4B.

Concluding, that the use of a constant value of β (or ε equivalently) or a linear relation between β and N_d^sat improves the closure error, we determined optimal parameters for a linear relationship between β and N_d^sat which minimizes the error with respect to the estimated in situ N_d (section 3).

3 Results and discussion

We determined the β values from Eq. 7, using each derived values of Ν_d and the corresponding values c(c_w), τ, r_eff as follows:

We then fit the β and N_d data to a linear relationship, β_opt = a + b N_d, to determine the “optimal β-expression” (OPT). According to Bevington et al. (1993) we used a weighted linear interpolation method, which considers both δβ (calculated by using the Chain rule on Eq. 7) and δN_d (δΝ_d = 0.25Ν_d from Morales and Nenes (2010)) of each β and N_d datum, to determine the OPT.

The coefficients of OPT, a and b were estimated to be equal to 1.0421 ± 0.1979, and 4.8717 10^–4 ± 6.1084 10^–4, respectively (Figure 3), while the average δβ_opt was estimated to be equal to 0.28 for the whole dataset. Additionally, we calculated the p-value and R-value of the fit and found equal to 0.089 and 0.412, respectively, while the fitting confidence R² was found equal to ∼0.17.

FIGURE 3. The linear interpolations over the β-N_d data for the diabatic derived N_d; the labeling refers to the serial number of each cloud (c.f. Supplementary Figure S1).

Then, we applied the OPT expression into Eq. 6, to calculate the solutions of N_d^sat, while we disregarded the solutions where δβ_opt > 1, δβ_opt/β_opt > 0.5, β_opt > 2, and β_opt < 1.

Finally, we validated the accepted solutions in respect of the N_d. The results of this closure is presented in Figure 4A.

Additionally, we present the MNB boxplots, using the literature based β-expressions (M94, RL03, PL03, Z06, RL09, and F12), versus our proposed expression (OPT) in Figure 4B.

Based on the results presented in Figure 4B and Table 3, we see that the proposed β-expression OPT exhibits the lowest mean MNB value (14.53%) with a standard deviation 36.33%. The performance of each β-expression can be ranked by their MNB values, as follows: OPT (−14.53%), M94 (−17.37%), Z06 (−21.25%), PL03 (23.51), GCMs (−28.80%), F12 (−31.99%), and RL03 (51.34%) (see also Table 3) along with the resulting standard deviation values (expressed as length of the box in the vertical axis) of MNB (c.f. Figure 4B).

It is important to note here, that the OPT expression coefficients are quite close to those proposed by Peng and Lohmann (2003), although the latter provide a degree of closure with twice the bias and standard deviation compared to the OPT expression. This suggests that the Peng and Lohmann (2003) and OPT parameters is appropriate for cloud types such as those sampled in our study, and, that the methodology adopted here to estimate in situ N_d is realistic.

4 Conclusion

The study presented here expands an established droplet number retrieval algorithm for non-precipitating PBLCs (Bennartz (2007) to explicitly account for the spectral dispersion of droplets and its dependence on droplet number in terms of β. The revised algorithm uses the cloud microphysical variables τ and r_eff, as derived from SEVIRI onboard the geostationary meteorological satellite (METEOSAT) with a temporal resolution of 15 min and with a spatial resolution 3.6 km × 4.6 km, along with an improved calculation of the total condensation rate (Zhu et al., 2018) with respect to cloud top height which can be obtained by using the ERA5 atmospheric pressure-temperature profiles (Hersbach et al., 2018). We found that the optimal retrieval of N_d^sat is most sensitive to biases of the β values, rather than biases in τ and r_eff, pointing to the need for an optimal β-expression for the most accurate N_d^sat retrievals.

We then calculated the retrieved N_d^sat values by using the literature-based β-expressions and we evaluated them against the in situ N_d estimations obtained by the droplet activation parameterization of the Nenes and Seinfeld (2003). We found that droplet number is captured to within ± 29% and ± 61%; based on these results we see that by using a constant value of β, or a linear relation between ε or β to N_d^sat, such as PL03, Z06, GCMs, and F12, the N_d^sat is captured to within ± 35%. Additionally, we proposed a new β-N_d expression, based on the in situ N_d estimations, that optimizes the closure between N_d^sat and N_d within ± 33%. Furthermore, the new β-expression we obtained through the optimal fit between N_d^sat and N_d is remarkably similar to the PL03 relationship. Given that, the PL03 relationship derived from observation data suggests that our method to estimate N_d is realistic. The use of either RL03 or our optimized relationship, captures droplet number to within 30%, which is comparable to the closure levels obtained from in situ observations.

Although more work needs to be done to evaluate the extent to which our approach can be applied elsewhere in the globe, the results presented here are both encouraging and may suggest ways to develop high-value products for climate models that can take advantage of the rich ground-based aerosol datasets available to the community.

Data availability statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Author contributions

Conceptualization, RF, AN, and AP; methodology, RF, AN, and PG; data analysis, RF, AN, MK, SV, VV, and PK; model simulations, RF, AN; investigation RF, AN, and PG; writing—original draft preparation, RF, AN, AP, and PG; review and editing, RF, AN, AP, PG, KE, MK, and MG; visualization, RF; supervision, RF, AN, and AP; All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded jointly by the PANhellenic infrastructure for Atmospheric Composition and climatE change (PANACEA) research project (MIS 5021516), implemented under the Action Reinforcement of the Research and Innovation Infrastructure, and the Operational Program Competitiveness, Entrepreneurship, and Innovation (NSRF 2014–2020), co-financed by Greece and the European Union (European Regional Development Fund). Additional funding was received from the European Research Council, CoG-2016 project PyroTRACH (726165) funded by H2020-EU.1.1.—Excellent Science, and from the European Union Horizon 2020 project FORCeS under grant agreement No 821205. Open access funding provided by École Polytechnique Fédérale de Lausanne.

Acknowledgments

The authors also acknowledge the NOAA Air Resources Laboratory (ARL) for the provision of the HYSPLIT transport model; the Finnish Meteorological Institute (FMI) for providing the HALO lidar; the NCSRD for hosting the HALO lidar. ERA5 data were generated using Copernicus Climate Change Service Information and Cloud Optimal Analysis (OCA) products were provided from EUMETSAT. Neither the European Commission nor the ECMWF are responsible for any use that maybe made of the Copernicus information or data in this publication (https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-pressure-levels?tab = form; last access: 13 July 2021). Part of this work was financed by the Basic Research Program, NTUA (PEVE) under contract PEVE0011/2021.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/frsen.2022.958207/full#supplementary-material

References