Shallow Cumulus Cloud Fields Are Optically Thicker When They Are More Clustered

Shallow trade cumuli over subtropical oceans are a persistent source of uncertainty in climate projections. Mesoscale organization of trade cumulus clouds has been shown to influence their cloud radiative effect (CRE) through cloud cover. We investigate whether organization can explain CRE variability independently of cloud cover variability. By analyzing satellite observations and high-resolution simulations, we show that increased clustering leads to geometrically thicker clouds with larger domain-averaged liquid water paths, smaller cloud droplets, and consequently, larger cloud optical depths. The relationships between these variables are shaped by the mixture of deep cloud cores and shallower interstitial clouds or anvils that characterize cloud organization. Eliminating cloud cover effects, more clustered clouds reflect up to 20 W/m$^2$ more instantaneous shortwave radiation back to space.


Introduction
Marine shallow cumulus clouds, as the most prevalent cloud type (Johnson, Rickenbach, Rutledge, Ciesielski, & Schubert, 1999), play a vital role in the climate system by reflecting incoming solar radiation back to space (Bony, Dufresne, Le Treut, Morcrette, & Senior, 2004;Bony & Dufresne, 2005;Bony et al., 2015).The response of these clouds to changes in cloud-controlling factors remains as one of the largest sources of uncertainty of climate projections (Schneider et al., 2017;Nuijens & Siebesma, 2019;Sherwood et al., 2020), despite recent advances (e.g., Vogel et al., 2022).Shallow cloud fields in the trades exhibit a diverse range of patterns, which have subjectively, and almost poetically, been classified as Sugar, Gravel, Flowers and Fish (Stevens et al., 2020).A comprehensive analysis by Janssens et al. (2021) shows that the quantification of such patterns needs at least two effective dimensions.Cloud fraction f c as a bulk 1D measure and the organization index I org , which quantifies the level of nonrandomness in cloud spatial distribution within a cloud field (Weger, Lee, Zhu, & Welch, 1992;Tompkins & Semie, 2017) are an example of a suitable variable choice to represent these two dimensions.How relevant this mesoscale organization is for the low-cloud climate feedback remains an open question.
The shortwave (SW) and longwave (LW) radiative effect of clouds is sensitive to organization (Denby, 2020).The daily mean cloud radiative effect (CRE) varies by approximately 10 W/m 2 among subjectively classified cloud patterns, primarily due to differences in f c , with a variability of 5 to 10 W/m 2 at a fixed f c (Bony, Schulz, Vial, & Stevens, 2020).Contrary to the case of deep convective clouds (Tobin, Bony, & Roca, 2012), where outgoing LW radiation increases with clustering, Luebke, Ehrlich, Schäfer, Wolf, and Wendisch (2022) suggest a correlation between increased I org values and reduced LW warming.No influence of clustering on SW cooling is observed in their study.For stratocumulus cloud decks, McCoy, McCoy, Wood, Zuidema, and Bender (2022) demonstrate that different morphologies, indicative of differences in the horizontal organization of the cloud decks, modulate the relationship between albedo and f c .
We aim to investigate whether -independent of f c variability -the horizontal organization of shallow cumulus cloud fields has an impact on the their net CRE.To do so, we combine satellite data with a large ensemble of large-eddy simulations (section 2).After removing the confounding effect of cloud fraction (section 3.1), we show that clustered cloud fields feature optically thicker clouds (section 3.2).This stems from clustered cloud fields containing more liquid-water path and smaller retrieved cloud droplets (section 3.2).Analyzing the simulations establishes that increased liquid-water path due to increased clustering primarily results from increased cloud geometric thickness (section 3.3).Section 4 concludes.

Description of the Data
Following previous studies (Stevens et al., 2020;Bony et al., 2020;Janssens et al., 2021), we focus on clouds over the tropical Atlantic Ocean to the east of Barbados (10 • -20 • N, 48 • -58 • W), which are representative for the trades (Medeiros & Nuijens, 2016).Our analysis covers December to May of 2002 to 2020.Our satellite dataset combines data from NASA's Moderate Resolution Imaging Spectroradiometer (MODIS) aboard the Aqua satellite, with data from the Clouds and the Earth's Radiant Energy System (CERES) instrument.We compute organization metrics from MODIS' cloud masks.For each cloudy scene, we calculate two metrics -cloud fraction (f c ) and degree of organization (I org ).The preprocessing of MODIS' cloud masks follows Janssens et al. (2021): Scenes with > 20% cirrus coverage are excluded, as are cloud fields with solar zenith angles > 45 • .In contrast to Janssens et al. (2021), we use the full 10 • × 10 • domain.Following Schulz, Eastman, and Stevens (2021), we focus solely on shallow clouds by excluding scenes with cloud-top heights z t > 4 km.After preprocessing, approximately 750 cloud fields remain for analysis.CERES provides hourly top-of-the-atmosphere SW and LW radiative fluxes for all-sky and clear-sky conditions, as well as cloud optical depth (τ c ), cloud albedo (A c ), cloud-top height (z t ), liquid water path (L), and cloud-droplet effective radius (r e ).We select CERES data around 13:30 local time, which corresponds to the overpass time of the Aqua satellite.Shortwave cloud radiative effect (SWCRE) and longwave cloud radiative effect (LWCRE) are calculated as the difference between the all-sky and clear-sky radiative fluxes at the top of the atmosphere.For each cloud scene, we calculate domain-mean values of cloud properties provided by CERES.
We amend our satellite analysis with the Cloud Botany dataset (Jansson et al., 2023).This is a large ensemble (ca. 100 members) of high-resolution (100 m) large-eddy simulations (LES) of shallow cumulus clouds with a domain size of 150 km by 150 km.It was initialized with a variety of conditions derived from ERA5 reanalysis data (Hersbach et al., 2020) of trade cumuli that cover the climatological conditions of the area under consideration.We refer to the dataset paper, Jansson et al. (2023), for details.Our motivation for employing the Botany simulations is twofold.Firstly, considering that liquid-water path L and effective radius r e might be underestimated in less organized conditions of broken cloud fields containing small clouds (Zhang & Platnick, 2011;Seethala & Horváth, 2010;Painemal & Zuidema, 2011), the simulations support that our results are physical.Secondly, the simulations provide data on cloud-base height (z b ) and cloud geometric thickness (h) so that we can investigate how they are correlated to organization.We use hourly data from hours 37-43 of the simulations (574 cloud fields in total).These times are chosen because they approximately align with the daily overpass times of the Aqua satellite to match the diurnal phase.To determine the geometric thickness h of each cloudy column, we calculate the difference between the altitudes of the highest and lowest cloudy pixels where the liquid water specific humidity is larger than zero.Subsequently, for each cloud field, we compute the domain-averaged h.We further compute the mean size of cloud objects within each cloud field using: L c = ( n 1 A i )/n, where A i represents the area of each individual cloud object i, and n corresponds to the total number of cloud objects within the field.

Results & Discussions
3.1 Cloud organization impacts CRE independent of f c variability Figure 1(a-c) shows that as clouds become more organized (increasing I org ), they reflect less SW radiation towards space (smaller magnitude of SWCRE).In addition, increased organization also leads to a decrease in LWCRE.Overall, the warming effect induced by the SW component is partially compensated by the reduced LW warming.Consequently, enhanced organization of clouds results in diminished net cloud radiative cooling.
It is crucial to emphasize that the relationships illustrated in Fig. 1(a-c) are confounded by the variability of f c , as f c is correlated with I org , SWCRE, and LWCRE (Fig. S1).
To eliminate the confounding effect of f c , we employ the concept of partial correlation analysis (Baba, Shibata, & Sibuya, 2004).For any given metric (e.g., X), we eliminate the variability associated with f c using a regression analysis, where f c serves as the regressor, c represents the coefficient, and X|f c denotes the remaining variability in X that cannot be explained by f c .
Fig. 1 (d) shows that as I org |f c increases, SWCRE|f c becomes more negative, i.e., as clouds cluster, they reflect more incoming SW radiation.The resulting SW cooling amounts to up to 20 W/m 2 .This indicates that the positive correlation observed in Fig. 1(a) is due to I org and SWCRE being negatively correlated to f c (Fig. S1, a, b).Similarly, after elimination of f c variability, the response of LWCRE to cloud organization is strongly reduced to about 1 W/m 2 (Fig. 1, e), indicating that the correlation between I org and LWCRE in Fig. 1(b) is almost solely due to their mutual correlation with f c (Fig. S1, c).The variability in LWCRE due to clustering is thus similar in magnitude to the LW radiative effect (≈ 0.75 W/m 2 ) of the "cloud twilight zone" (Eytan, Koren, Altaratz, Kostinski, & Ronen, 2020).Ultimately, as Fig. 1(f) illustrates, the dependence of net CRE on I org |f c arises almost exclusively from the dependence of the SW component on cloud organization.

Clustering and cloud optical thickness are positively correlated
In the previous section, we eliminated the impact of f c on the I org -SWCRE relationship.The remaining variability in SWCRE is caused by variations in albedo and could be modulated by a sensitivity of 3D radiative effects to cloud organization.Regarding the latter, Singer, Lopez-Gomez, Zhang, and Schneider (2021) report that biases in topof-the-atmosphere SWCRE from neglecting 3D radiative effect are negligible for both, unorganized and organized shallow cumulus cloud fields (their Fig. 7).This is especially expected for averages over large domains, as analyzed here.Indeed, a bi-linear regression with f c and albedo as regressors can explain 94 % of variability in SWCRE in our dataset (Fig. S2).This confirms that the impact of 3D radiative effects is negligible in these large cloud fields.
Having 3D effects excluded, the remaining variability in SWCRE primarily corresponds to changes in albedo and equivalently cloud optical depth τ c . Figure 2 On average, τ c |f c increases with increasing I org |f c (Fig. 2, b).Despite the seemingly modest increase in τ c |f c , such changes can result in a significant increase in albedo of 0.03 to 0.04 (Figs.S3, S4), and the notable increase in SW reflection by 20 W/m 2 (Fig. 1, d).The I org |f c -τ c |f c relationship indicates that horizontal cloud field organization, as measured by I org , directly corresponds to its vertical organization, as captured by τ c : trade cumuli are optically thicker when they are more clustered.
As theoretically expected (Han, Rossow, & Lacis, 1994), τ c is proportional to L/r e in our dataset (Fig. S5). Figure 3 shows that both, L and r e , contribute to mediating the relationship between organization and optical depth.With the f c effect eliminated, there is a positive correlation between the degree of cloud clustering and the amount of liquid water present within the clouds (Fig. 3, a).Similarly, as the level of clustering increases, clouds tend to exhibit smaller radii r e (Fig. 3, b).A cloud field that is dominated by Flowers shows smaller radii r e but higher liquid-water path L compared to the field dominated by Sugar and Gravels.Our L-organization relationship seems to be in contrast to Schulz et al. (2021) who show that individual Gravel clouds have higher liquid-water path L compared to individual Flowers.To reconcile this with our results, we need to remind ourselves that the large cloud scenes analyzed here contain a mixture of different clouds.Stevens et al. (2020) report that Gravel clouds tend to coexist with Sugar.For Flowers such a coexistence is less pronounced.Instead, Flowers feature anvils.Such shallow cumulus anvils have notable geometric thickness (200 m (Wood et al., 2018) to 600 m (Dauhut et al., 2023)).This means that anvil cloudiness is optically thicker and more reflective than typical Sugar (see also Fig. 7 in Stevens et al. (2020)).When considering two cloud fields with identical f c , it therefore seems reasonable that a Flower -dominated field features a larger domainaveraged L as compared to a field dominated by Sugar and Gravels.
Similarly, the relationship between I org , L and r e might seem unexpected: Based on adiabatic parcel lifting, we would expect L and r e to be positively correlated, while Fig. 3 suggests a negative correlation.In fact, L and r e are positively correlated in our dataset, while projecting their variability on organization (I org |f c ) suggests a negative correlation: as degree of clustering increases, L increases but r e decreases (Fig. S6).Satellite snapshots of r e (Fig. S7) reveal that Flowers have substantially smaller r e in their veils compared to their core updrafts.In contrast, Gravels exhibit a more homogeneous r e with relatively larger values (compared to Flowers) associated with strong updrafts at the edge of their cold pool fronts.These observations suggest that for Flowers, the average L is primarily influenced by their cores, while the average r e is influenced by their veils.This is consistent with the fact that L is proportional to r e 6 , resulting in a more pronounced contrast between the core and veil in L compared to r e .As a result, the average L of Flowers is larger compared to that of Gravels, while their average r e is smaller in comparison to Gravels.
It is interesting to contrast the small droplets in relatively thick anvils described here to the very large droplet sizes and optically thin veil clouds that have been reported in the context of the stratocumulus-to-cumulus transition (Wood et al., 2018).While O, Wood, Tseng, et al. ( 2018) report an increase in the corresponding ultra-clean conditions with boundary layer height, this relationship is unlikely to extend to deep trade cumulus Flowers, which can be considered shallow mesoscale convective systems with complex outflow dynamics (Dauhut et al., 2023).On the microphysical process level, ultra clean conditions have been associated with strong precipitation scavenging (O, Wood, & Bretherton, 2018), while Radtke, Vogel, Ament, and Naumann (2023) discuss that the conversion efficiency to precipitation decreases with increasing organization in trade cumulus.
Overall, our discussion of the relationship between liquid-water path and effective radius to organization and the resulting effects on optical depth highlight that organized cloud fields cannot be conceptualized with a single, typical profile of cloudiness.Instead, the spatial variability requires considering two types of clouds as discussed for Gravel and Sugar versus Flower cores and anvils.This duality resonates with the effective twodimensionality of quantitative measures for mesoscale organization (Janssens et al., 2021;Shamekh, Lamb, Huang, & Gentine, 2023).

Mean cloud geometric thickness increases with clustering
For an entraining lifting parcel, liquid-water path L ∝ f ad h 2 , where f ad represents the degree of adiabaticity and h denotes geometric cloud thickness.To explain the observed I org -L relationship, we therefore investigate the relationships of I org to h and f ad in the Botany simulations.Repeating the analysis from Sects.3.1, 3.2 for the simulation data shows qualitative agreement and thus justifies this approach (Figs.S9(a-c), S10).Note that the discrepancy in the response of r e to cloud organization between simulations and satellite data (Fig. S9, d) is expected from the fixed cloud droplet number in the simulations but does not fundamentally affect our discussion of L here.
Figure 4(a) shows that the domain-averaged geometric thickness increases by more than 100 m as cloud fields become more clustered (increasing I org |f c ).Additionally, the variability in h has a significantly larger influence on the value of L than f ad (Fig. S11).Thus, our LES-based results indicate that the simulated increase in L due to enhanced clustering (Fig. S9, c) primarily stems from the geometric thickening of cloud fields.
Figure 4(b) further explores the relationship between horizontal and vertical cloud field properties and shows that the average size of cloud objects (L c ) increases with I org |f c .This positive correlation shows that cloud horizontal extent as quantified by L c , is positively correlated to cloud vertical extent as quantified by h, consistent with the findings of Feingold et al. (2017).The figure moreover illustrates that an increase in I org |f c corresponds to a rise in the domain-average cloud-base height (z b ).Note that z b is not the lifting condensation level; instead, it represents the lowest height of a cloudy pixel within each column.This means that a higher domain-mean z b is an indication of the presence of more anvils in the field.Overall,Fig. 4(b) demonstrates that enhanced clustering leads to a higher occurrence of larger cloud objects with elevated domain-mean z b , indicating a higher degree of anvilness.

Conclusions & Outlook
We have explored the impact of shallow cumulus cloud field organization on cloud radiative effects, where confounding variability of f c was removed through partial correlation analysis (Eq. 1).Based on satellite data, our analysis shows that an increased level of clustering (I org |f c ) results in up to 20 W/m 2 higher SW reflection to space (Fig. 1, d,  f).We observe that, irrespective of f c variations, more clustered cloud fields exhibit, on average, higher liquid water path (Fig. 3, a), smaller cloud droplets (Fig. 3, b), and consequently, greater optical thickness (Fig. 2).A complementing ensemble of large-eddy simulations indicates that increased clustering leads to geometrically thicker cloud fields that feature increased anvilness (Fig. 4). Figure 5 summarizes these results.Collectively, they suggest that, eliminating the effect of f c , the distribution of horizontal cloud sizes ultimately determines the vertical extent of clouds, subsequently influencing liquid-water path and cloud optical depth, and ultimately albedo and SWCRE.
What do our results mean in terms of the cloud feedback of trade cumulus?Myers et al. (2021) (their Supplementary Fig. 10) show that in addition to an increase in seasurface temperature, which is not expected to trigger a notable response in trade cumulus cloudiness (Myers et al., 2021;Cesana & Del Genio, 2021), estimated inversion strength EIS is projected to moderately increase, and surface wind to slightly decrease.According to Bony et al. (2020) such an increase in EIS would favor high-cloud-fraction Flowers over Gravel and Sugar with lower cloud fractions.In contrast, the decreasing surface wind would favor Sugar.While our results highlight the tight relationship between horizontal and vertical cloud organization, whether cloud fraction and optical depth interact positively or negatively in response to drivers of organization remains an open question.To address this interplay, we need to further explore how mesoscale processes (Janssens et al., 2023;George, Stevens, Bony, Vogel, & Naumann, 2023;Vogel, Konow, Schulz, & Zuidema, 2021) modulate cloud fraction, liquid-water path, effective radii, and anvilness.

Increased anvilness
The probability of Gravels coexisting with Sugar is much higher compared to Flowers.

Increased clustering
Gravel and Sugar Flowers

Calculation of the degree of adiabaticity
The degree of adiabaticity is computed using the following equation: where, z represents the height, ρ(z), q l (z), and Γ ad (z) denotes density, liquid water specific humidity, and the moist adiabatic lapse rate of q l at each model height.To ensure a more precise estimation, f ad is computed only for columns where the associated z b is in proximity to the lifting condensation level (<850 m) and the cloud thickness exceeds 150 m.
The calculation of Γ ad is determined by the following equation (Schmeissner et al., 2015;Eytan, Koren, Altaratz, Pinsky, & Khain, 2021): Here, c p represents the specific heat constant, g denotes the acceleration due to gravity, L v stands for the latent heat of vaporization.The moist adiabatic lapse rate for temperature, Γ m , is computed as: In the above equation, R d corresponds to the gas constant for dry air, T represents temperature, β is the constant in the Clausius-Clapeyron equation (e s = A exp(β(T −T 0 ))), and q s denotes the saturation specific humidity, which can be approximated as ≈ 0.622 es P where P represents pressure.
For simplicity, we calculate Γ ad for the middle of the cloud layer (at z ′ = z b + z t 2 ) and assume that it linearly changes with height through the cloud layer.This assumption leads to the equation below: where, a f ad value of 1 in each column indicates a fully adiabatic cloud, while values smaller than 1 and larger than 0 indicate non-adiabatic clouds.In the final step, after obtaining f ad for each column, we proceed to calculate the domain-mean f ad for the entire cloud field.

Figure 1 :
Figure 1: Dependence of domain-mean CRE on organization.In the first row (a-c), the relationships between I org and SWCRE (a), LWCRE (b), and net CRE (c) are illustrated.The second row (d-f) shows the same relationships but with removing the f c variability.The mean values of I org (for the 2 nd row, I org |f c ) in each bin are denoted by red circles, with their size proportional to the number of points in the bin.The red dots are fitted with a dashed black line.Values below the 5 th and above the 95 th percentile of I org (for the 2 nd row, I org |f c ) are excluded from the fit.
(a) displays the variability of cloud patterns in a plane spanned by τ c |f c and I org |f c .This figure shows a continuous range of patterns, in the terminology of Stevens et al. (2020) ranging from Sugar and Gravels in the lower left corner to Fish and Flowers in the upper right corner.

Figure 2 :Figure 3 :
Figure 2: Dependence of domain-mean cloud optical depth on organization.The scatter plot of MODIS cloud features (a) and the 2-dimensional histogram (b) depict the I org |f c -τ c |f c relationship.Specifically, for the scatter plot in (a), instead of displaying individual points, the entire cloud field is visualized to enhance pattern visualization.Clouds are represented in white, while the blue background represents the ocean color (MODIS true-color images).For plot (b), the mean values of I org |f c in each bin are denoted by red circles, with their size proportional to the number of points in the bin.The red dots are fitted with a dashed black line.Values below the 5 th and above the 95 th percentile of I org |f c are excluded from the fit.

Figure 4 :
Figure 4: Dependence of domain-mean geometric thickness, average size of cloud objects, and domain-mean cloud-base height on organization.(a) The figure shows the 2D histogram of the relationship between I org |f c and h|f c .(b) The plot shows the relationship between I org |f c and the mean-field cloud object size (L c ) with contour colors representing the values of domain-averaged cloud-base height (z b ).The gray shade indicates the inter-quartile range variability of L c in each bin of I org |f c .For both plots, the mean values of I org |f c in each bin are denoted by red circles, with their size proportional to the number of points in the bin.The red dots are fitted with a dashed black line.For both plots, values below the 5 th and above the 95 th percentile of I org |f c are excluded from the fit.

Figure 5 :Figure S1 :Figure S2 :
Figure5: Summary of results.Comparing two cloud fields with identical cloud cover, the Flower -dominated cloud field features a larger domain-averaged geometric thickness, a higher liquid water path, more frequent anvils with smaller cloud droplets, and consequently, brighter clouds and therefore larger SW reflection in comparison to the cloud field dominated by the presence of Sugar and Gravels.In the diagram, theory-based equations are provided for arrows where the underlying theory exists.The gray arrows illustrate apparent paradoxes which are discussed in section 3.2.

Figure S3 :
Figure S3: The relationship between I org and cloud albedo (A c ), having the effect of f c eliminated.

Figure S4 :Figure S5 :
Figure S4: The histogram plot of cloud scene albedo (A c ) in the Aqua satellite dataset.This plot shows that the range of variability in cloud field albedo is about 0.1.

Figure S6 :
Figure S6: This figure highlights the positive correlation between L and r e .However, when we examine the L-r e relationship in relation to clustering, a clear pattern emerges: as clustering increases (more Flowers), L increases while r e decreases.Each distinct color in the figure represents a specific class of organization.For instance, the purple color represents cloud fields where the value of I org |f c falls within the range of the 0th to 10th percentile of I org |f c .

Figure S7 :Figure S8 :
Figure S7: A Flowers (a) and a Gravel cloud field snapshots (from NASA Worldview) taken by Aqua satellites on 02/Feb/2020 and 17/Jan/2020, respectively.The cloud fields are colored by the value of r e which goes from 4 (yellow) to 30 (red) µm.

Figure S9 :
Figure S9: All results derived from the hourly analysis of the Cloud Botany dataset during the second day between the 37 th and 43 rd hours.Values below the 5 th and above the 95 th percentile of I org |f c are excluded from the fitting.