On Climate Change and Trade Cumulus Organization

Abstract

We investigate the role of mesoscale organization for the response of trade cumulus (Tc) clouds to climate change. Among four recently identified states of Tc organization, the “Sugar” state has the lowest and the “Flower” state the highest cloud fraction and cloud radiative effect. Using large‐eddy simulations, we find that the Flower Tc state is more sensitive to climate change than the Sugar Tc state. In the considered case, the short‐wave cloud radiative effect weakens by 0.28 W m⁻² in the Sugar state and by 1.5 W m⁻² in the Flower state over the course of 21st century under the RCP8.5 emissions scenario. This is accompanied by a reduction of the short‐wave cloud radiative effect variance on the mesoscale. The primary mechanism is stabilization of the boundary layer by stronger long‐wave radiative heating at the inversion associated with higher greenhouse gas levels. This weakens the boundary layer mesoscale circulation that is responsible for aggregation of moisture and formation of the Flower Tc state. Thus, in the considered case, organization on the mesoscale amplifies the positive feedback of Tc clouds to climate change. Owing to the widespread occurrence of boundary layer mesoscale circulations in the Tc regime, this mechanism could modulate the Tc response to climate change in general.

Introduction

Low clouds cool the Earth by reflecting sunlight back to space. The cooling effect depends on the response of clouds to the climate system, and in turn affects Earth's climate. This cloud-climate feedback determines climate sensitivity, the Earth's temperature response to rising levels of greenhouse gases. The latest generation of global climate models predicts a high climate sensitivity due to strong low cloud feedbacks (Zelinka et al., 2020). Observations, however, indicate a moderate climate sensitivity due to low clouds (Cesana & Del Genio, 2021; Myers et al., 2021). The discrepancy is caused in part by trade cumulus (Tc) clouds, a key source of uncertainty in climate projections (Bony & Dufresne, 2005; Medeiros et al., 2015; Vial et al., 2013). “Trade cumulus feedbacks in climate models are governed by changes in cloud fraction near cloud base” (Vogel at al., 2022; see also Brient et al., 2016; Vial et al., 2016). Climate models with a high climate sensitivity suggest a strong decrease in Tc cloud base cloudiness with warming owing to increased lower-tropospheric mixing and cloud base evaporation (Brient et al., 2016; Sherwood et al., 2014; Vial et al., 2016). A mechanism in which stronger surface fluxes due to a warming sea surface temperature (SST) deepen and desiccate the boundary layer, thus reducing Tc cloud fraction, was found in large-eddy simulations (Rieck et al., 2012).

Vogel et al. (2022) observationally refuted this mixing-desiccation hypothesis. They furthermore found that “a weak trade cumulus feedback is more plausible than a strong one,” providing an important line of evidence against high climate sensitivity. Observations (Nuijens et al., 2015) and large-eddy simulations (Blossey et al., 2013; Bretherton et al., 2013; Z. Tan et al., 2017; Vogel et al., 2016) suggest that Tc are remarkably stable against climate change, and exhibit a weak positive feedback (Vial et al., 2017). The leading hypothesis for a broad resilience of Tc against climate change is the cumulus-valve mechanism (Neggers et al., 2006). The cumulus-valve mechanism postulates a negative feedback between the cloud fraction and the cloud mass flux through cloud base as the atmosphere warms and moistens. This is a possible explanation for the weak feedback of Tc clouds under climate change.

The vertical gradient of specific humidity turns more negative with warming due to the non-linearity of the Clausius-Clapeyron relationship (Brient & Bony, 2013). This could lead to more entrainment of dry free tropospheric air into the boundary layer (Bretherton et al., 1995, 2013; Eastman & Wood, 2018) and to a reduction in Tc cloudiness with warming. Other mechanisms can enhance or reduce cloud fraction with warming. First is the radiative cooling in the boundary layer. Enhanced radiative cooling in a Tc-topped boundary layer leads to increased Tc cloudiness, and vice versa. For example, for non-precipitating shallow Tc, SST warming leads to enhanced boundary layer clear-sky radiative cooling, which destabilizes the cloud layer. Under an unchanged or stronger inversion, the boundary layer destabilization drives stronger cloud mass flux and hence increases Tc cloud fraction (Narenpitak & Bretherton, 2019; Wyant et al., 2009). On the other hand, increased carbon dioxide (CO₂) levels lead to more down-welling long-wave radiation and thus less boundary layer radiative cooling. This weakens the cloud mass flux and leads to shallower Tc clouds, but barely affects the cloud-base cloud fraction (Wyant et al., 2012). In general, the boundary layer radiative heating or cooling rate is balanced by other large-scale factors that control the boundary layer energy budget, and a change that leads to enhanced boundary layer radiative cooling could result in a negative Tc cloud-climate feedback (Narenpitak & Bretherton, 2019). These additional factors include, but are not limited to, surface heat fluxes, entrainment fluxes from the free troposphere, large-scale advection, and, if any, precipitation.

Precipitation helps balance the boundary layer energy budget and regulates the depth of Tc clouds and of the boundary layer in response to a warmer SST (Blossey et al., 2016; Stevens & Seifert, 2008; Vial et al., 2017). Bretherton et al. (2013) studied the responses of Tc feedback to various cloud-controlling factors, such as increased CO₂ and SST. The large-scale conditions in their study are those of the shallow cumulus case of the Cloud Feedback Model Intercomparison Project—Global Atmospheric System Study Intercomparison of Large-Eddy Simulation and Single Column Models (CGILS, Zhang et al., 2013). Thus, their Tc clouds are deeper, have more inversion clouds, and precipitate more than those studied by Wyant et al. (2009, 2012) and Narenpitak and Bretherton (2019). They found that increased CO₂ leads to weaker long-wave radiative cooling and reduced cloud fraction, consistent with Wyant et al. (2012), but only with a slightly shallower inversion height. While the warmer SST reduced the cloud fraction, it barely changed the inversion height of the precipitating Tc clouds. Furthermore, the precipitating Tc in the study of Bretherton et al. (2013) respond to surface wind speed differently than the non-precipitating Tc in Nuijens and Stevens (2012). These studies demonstrate that precipitation is another important factor for the Tc cloud-climate feedback.

Recent studies show that organization of shallow clouds adds complexity to their characteristics and behavior, with possible ramifications for their response to climate change (Narenpitak & Bretherton, 2019; Nuijens & Siebesma, 2019; Vial et al., 2017). McCoy et al. (2023) determined, based on an analysis of satellite observations, that the response of boundary layer cloud mesoscale cellular organization to climate change will render the global cloud feedback more positive. I. Tan et al. (2024) showed in satellite data that shifts between cloud morphologies over the Southern Ocean play an important role for the interannual short-wave cloud feedback.

Mesoscale aggregation of shallow cumulus convection is driven by gross moist instability (Bretherton & Blossey, 2017). Bretherton and Blossey (2017) found that it does not require feedbacks on the mesoscale from radiation, surface fluxes, or precipitation, but that radiative feedbacks and precipitation can enhance the aggregation. Stronger precipitation accompanying mesoscale organization is thought to result in “a more stable and drier trade-wind layer” (Vial et al., 2017). Vogel et al. (2020) found in large-eddy simulations that “the deepening and organization of shallow convection plays an important role in regulating stratiform cloudiness and thus total cloud cover in the downstream trades.”

Non-precipitating Tc clouds that organize more strongly may experience stronger long-wave cooling. Simulations of non-precipitating Tc clouds with different domain sizes that result in different degrees of self-aggregation suggest that, when the enhanced radiative cooling due to warmer SST leads to increased cloud fraction, self-aggregation might amplify the negative cloud feedback further (Narenpitak & Bretherton, 2019). This is because the enhanced radiative cooling inside the clouds destabilizes the cloud layer, leading to stronger cloud updrafts. This leads to stronger mesoscale circulation that amplifies the organization, which is found in large-eddy simulations even without the effects of climate change (Bretherton & Blossey, 2017; Narenpitak et al., 2023).

Four manifestations of Tc organization called Sugar, Gravel, Fish, and Flowers have been identified by Stevens et al. (2020) and studied in a joint field effort, the U.S. Atlantic Tradewind Ocean-Atmosphere Mesoscale Interaction Campaign (ATOMIC) and the European multinational Elucidating the Role of Clouds-Circulation Coupling in Climate (EUREC⁴A) campaign (Bony et al., 2017; Pincus et al., 2021; Quinn et al., 2021; Stephan et al., 2021; Stevens et al., 2021). Narenpitak et al. (2021) showed that the Flower cloud state observed during ATOMIC and EUREC⁴A arises by mesoscale circulation and moisture aggregation as a manifestation of gross moist instability (Bretherton & Blossey, 2017). Janssens, Vilà-Guerau de Arellano, van Heerwaarden, de Roode, et al. (2023) determined that this instability and upscale transport of moisture variance is an intrinsic feature of the non-precipitating Tc regime. The organizing mechanisms of the Gravel and Fish cloud states are currently unknown. The four Tc states cool the planet to different extents: the Sugar Tc state has the lowest and the Flower Tc state the highest cloud fraction and cloud radiative effect (Bony et al., 2020). This diversity suggests an uncharted landscape of Tc responses to climate change.

We explore the response of the Flower and Sugar Tc states to climate change from present-day (PD) to end-of-21st-century (EC) conditions using large-eddy simulations (LES). Idealized Eulerian LES serve to compare the response of the Sugar and Flower cloud states to climate change, and to identify the responsible mechanisms. Lagrangian LES initialized and forced by reanalysis data are used to quantify the sensitivity of the Sugar and Flower cloud states to climate change in a realistic meteorological setting.

The methodology is described in Section 2. Section 2.1 summarizes the adopted climate change scenario. The model is described in Section 2.2, the simulations and their setup in Section 2.3. Results of the Eulerian simulations are analyzed in Section 3, results of the Lagrangian simulations in Section 4. Section 5 provides a discussion, and Section 6 summarizes our findings and conclusions. Appendix A gives further details of the adopted climate change scenario. Appendix B describes the construction of idealized large scale conditions used in the Eulerian simulations. Supplemental text and figures are given in Supporting Information S1.

Methodology

21st Century Climate Change

We define PD as the period 2016–2025 and EC as the period 2090–2099. We split climate change from PD to EC into a contribution from a change in large scale conditions (LSC) and a contribution from change in greenhouse gas (GHG) levels. LSC comprises SST, air temperature, water vapor, zonal and meridional wind speed, subsidence, and horizontal advective tendencies of heat and moisture. The GHG species considered are carbon dioxide (CO₂), methane (CH₄), nitrous oxide (N₂O), and ozone (O₃). LSC and GHG in PD and EC are constructed from results of climate simulations in the Coupled Model Intercomparison Project 5 (CMIP5, Taylor et al., 2012) with the Community Earth System Model—Whole Atmosphere Community Climate Model (CESM1(WACCM), Marsh et al., 2013), under the Representative Concentration Pathway scenario 8.5 (RCP8.5, Van Vuuren et al., 2011). The RCP8.5 scenario is “the best match out to midcentury under current and stated policies with still highly plausible levels of CO₂ emissions in 2100” (Schwalm et al., 2020). PD and EC are characterized by the mean LSC and GHG concentrations over the respective periods in a region of the Caribbean east of Barbados. From PD to EC, CO₂ approximately doubles and SST warms by 1.8 K in this region. The construction of PD and EC LSC and GHG conditions is described in detail in Appendix A.

Model

We use the System for Atmospheric Modeling (SAM, Khairoutdinov & Randall, 2003). Cloud microphysics is described using a two-moment method which represents the cloud and rain modes with log-normal functions (Feingold et al., 1998). The cloud microphysics interacts with aerosol by activation, cloud processing, release upon droplet evaporation, and wet deposition. It is described in detail in Yamaguchi et al. (2017). Radiation is calculated from temperature, gas phase constituents including CO₂, CH₄, N₂O, and O₃, and liquid water mass and effective radius with the Rapid Radiative Transfer Model (Iacono et al., 2008; Mlawer et al., 1997, RRTMG). Profiles of temperature, water vapor, and of GHG species from the forcing data are prescribed above model top for the radiative transfer calculations. The model calculates interactive surface fluxes of heat, moisture, momentum, and of sea spray aerosol. Further model details are given in Kazil et al. (2021).

Simulations

Both Eulerian and Lagrangian simulations were conducted. The Eulerian simulations use idealized LSC, constructed to produce a Sugar cloud state during a first daytime period and a Flower cloud state in the second. The purpose of simulations with a cloud state extending over at least one daytime period is to allow temporal averaging to suppress intermittent variability, caused largely by precipitation in the Flower cloud state, in order to isolate the signal of climate change. The Lagrangian simulations use realistic LSC, and result in a transition from the Sugar to the Flower cloud state during one daytime period. Their purpose is to quantify their response assuming a present-day meteorological setting.

The simulations use a horizontal domain of 192 × 192 km with periodic lateral boundary conditions. The horizontal grid spacing is 100 m, the vertical grid spacing 50 m from the surface to 5,000 m, followed by 25 geometrically increasing levels up to domain top at 10 km. The simulation time step is 2 s. Time is measured as fractional day of year (d), with d = 0 corresponding to 1 January 2020, 00h00m00s Coordinated Universal Time (UTC). Boundary layer aerosol is initialized with 400 mg⁻¹ particles in the accumulation mode and 13 mg⁻¹ particles in the coarse mode, free tropospheric aerosol with 32 mg⁻¹ particles in the accumulation mode. The aerosol initialization is based on measurements during the ATOMIC field campaign (Quinn et al., 2021), and is described in further detail in Narenpitak et al. (2021).

Unless noted otherwise, cloud fraction is diagnosed as the fraction of locations with a vertically integrated liquid water optical depth ≥1. The inversion height is diagnosed as the altitude with the maximum vertical gradient of liquid water static energy. Buoyant production of turbulence kinetic energy (TKE) is calculated as $$\frac{g}{\langle {\theta }_{v}\rangle }{w}^{\prime }{\theta }_{v}^{\prime }$$, where g is Earth's accelaration due to gravity at the surface, w vertical wind speed, and θ_v virtual potential temperature. Brackets denote the mean over a geopotential isosurface. Primes denote perturbations from this mean. TKE is calculated as $$\frac{1}{2}\sqrt{{\left({u}^{\prime }\right)}^{2}+{\left({v}^{\prime }\right)}^{2}+{\left({w}^{\prime }\right)}^{2}}$$, where u and v are the horizontal wind speed. The vertical flux of moist variables due to resolved dynamics is calculated as ρ_airwq, where q is the moisture mass mixing ratio, and ρ_air the volumetric mass density of air.

Eulerian Simulation Setup With Idealized Forcing

The Eulerian simulations are initialized and forced with idealized LSC and GHG profiles that represent PD and EC conditions (Appendix B). Apart from the warming and moistening of the atmosphere, the change from PD to EC conditions features an increase in atmospheric stability with altitude (Figure B1h). Four simulations were conducted: simulation E_PD uses PD LSC and PD GHG, simulation E* uses EC LSC but PD GHG, and simulation E_EC uses EC LSC and EC GHG. Hence the signal of 21st century change in LSC (E_PD → E*) and GHG (E* → E_EC) can be analyzed separately. The fourth simulation E** is initialized and forced with PD LSC and EC GHG. It serves to demonstrate that our findings pertaining to the role of higher GHG hold regardless of the base state of the atmosphere, within the range of base states considered.

No nudging is used except at the two top model levels, to maintain continuity between the temperature and water vapor below and above domain top. Surface production of sea spray aerosol is turned off. The simulations are located at 15.14°N, 50.44°W, the initial time is February 2, 02h00m00s UTC (d = 31.083), and the simulations cover a period of 44 hr. The simulations use a domain translation in the zonal direction with −8 m s⁻¹.

Lagrangian Simulation Setup With Realistic Forcing

The Lagrangian simulations follow a 24 hr trajectory south-east of Barbados (Figure SF1 in Supporting Information S1) on 2 February 2020. The trajectory was constructed from the ERA5 wind field at 500 m above sea level using the Hybrid Single Particle Lagrangian Integrated Trajectory Model (HYSPLIT, Stein et al., 2015). It was chosen using Geostationary Operational Environmental Satellite-16 (GOES-16) satellite imagery as it leads to a location with larger Flower clouds compared to the more north-easterly trajectory in Narenpitak et al. (2021).

The Lagrangian PD simulation L_PD is initialized and forced with ERA5 LSC along the trajectory, that is, with PD LSC. L_PD uses PD GHG (Appendix A). The Lagrangian EC simulation L_EC is initialized and forced with the ERA5 LSC along the trajectory to which the change from PD to EC conditions was applied, and uses EC GHG (Appendix A). The simulations are nudged using Newtonian relaxation: horizontal mean temperature and water vapor are nudged toward the LSC in the free troposphere, and horizontal advective tendencies of heat and moisture are set to zero at all levels. Mean horizontal wind speed is nudged toward the LSC at all altitudes. Large scale subsidence from the LSC is prescribed. Surface production of sea spray aerosol is turned on. The simulations use domain translation at the trajectory velocity. The Lagrangian simulation approach is described in detail and evaluated against observations in Kazil et al. (2021) and Narenpitak et al. (2021).

Eulerian Simulations With Idealized Forcing

Figure 1 shows selected time series in simulation E_PD, E*, and E_EC. The cloud field is shown in Figure SF2 in Supporting Information S1. Additional time series are shown in Figure SF5 in Supporting Information S1; mean temperature, water vapor, relative humidity, and liquid water static energy profiles are given in Figure SF6 in Supporting Information S1.

[IMAGE OMITTED. SEE PDF]

Mean individual cloud area (Figure 1a), liquid water path (LWP, Figure 1b), and cloud fraction (CF, Figure 1c) evolve nearly identically until the first sunset (d = 32.878), when they begin to diverge: E_PD produces the highest and E_EC the lowest values. Values in simulation E* are at or below simulation E_PD, but above simulation E_EC, not considering variability on time scales of hours. The parting of the mean individual cloud area time series (Figure 1a) indicates the transition from the Sugar to the Flower cloud state. The cloud field is in the Sugar cloud state throughout the first daytime period and until first sunset of the simulations, and in the Flower cloud state thereafter (Figure SF2 in Supporting Information S1). The largest Flower cloud, in simulation E_PD, grows to a size of 75 km (Figure SF2d in Supporting Information S1).

We first discuss the cumulus-valve mechanism (Neggers et al., 2006). We diagnose cloud fraction and the cloud mass flux at the mean base of active cloud convection (mBACC). The mBACC is diagnosed as the altitude at which the average profile of TKE production at overcast locations first reaches or exceeds 0 above the altitude where the TKE production profile has a minimum (Supporting Information S1). The response of the mBACC E_PD → E* → E_EC (Figure 1d) can be divided into two periods that roughly coincide with the initial Sugar cloud state and the subsequent Flower cloud state. In the Sugar cloud state, mBACC descends E_PD → E* → E_EC as the boundary layer warms and moistens. In the Flower cloud state, the mBACC settles at a similar altitude across the simulations. Cloud fraction and cloud mass flux at mBACC (Figures 1e and 1f) are insensitive to the warming and moistening (Figure SF6 in Supporting Information S1) of the boundary layer E_PD → E* → E_EC. This is consistent with the cumulus-valve mechanism. Yet, once the Flower cloud state takes shape, LWP (Figure 1b) and CF (Figure 1c) decrease E_PD → E* → E_EC. Therefore, if the cumulus valve mechanism indeed acts in the simulated cloud field, as indicated by the insensitivity of the cloud base cloud fraction and mass flux, it does not protect the Flower cloud state against the effect of climate change.

Table 1 gives mean values over the Sugar and Flower periods. The Sugar period is taken as d = 32.25–32.878, the Flower period as the following 24 hr (d = 32.878–33.878). The transition from the Sugar to the Flower cloud state at d = 32.878 (sunset) is based on the parting of the mean individual cloud area time series at this time (Figure 1a). The longer averaging period for the Flower cloud state increases representativeness of the mean by better suppressing variability in LWP and CF of the Flower cloud state on time scales of hours (Figures 1b and 1c). Averaging over the shorter period for the Sugar cloud state suffices to suppress its very low variability on time scales of hours, but under-represents the nighttime contribution of the response of LWP and CF to the response of the Sugar cloud state. This has no impact on the response of the short-wave cloud radiative effect (SWCRE).

Table 1 Absolute Change in Cloud Fraction, Liquid Water Path (g m⁻²), and Short-Wave Cloud Radiative Effect (W m⁻²), With Relative Change in Parentheses, for the Sugar Cloud State and the Flower Cloud State in the Eulerian Simulations With Idealized Forcings

	Sugar	Flower
ΔCF	ΔLWP	ΔSWCRE	ΔCF	ΔLWP	ΔSWCRE
EPD → EEC	−1.2 ⋅ 10−3 (−1.3%)	−0.42 (−6.1%)	0.28 (−2.9%)	−1.3 ⋅ 10−2 (−10.0%)	−13.94 (−35.0%)	1.50 (−11.0%)
EPD → E*	−2.1 ⋅ 10−5 (0.0%)	0.11 (1.6%)	−0.02 (0.24%)	1.4 ⋅ 10−3 (1.1%)	−6.68 (−16.8%)	0.44 (−3.2%)
E* → EEC	−1.2 ⋅ 10−3 (−1.3%)	−0.54 (−7.6%)	0.31 (−3.1%)	−1.4 ⋅ 10−2 (−11.0%)	−7.26 (−21.9%)	1.06 (−8.0%)
EPD → E**	−1.4 ⋅ 10−3 (−1.5%)	−0.51 (−7.4%)	0.31 (−3.2%)	−8.9 ⋅ 10−3 (−7.0%)	−7.12 (−17.9%)	1.73 (−12.6%)

Flowers are more sensitive than Sugar to climate change (E_PD → E_EC). In terms of the relative reduction from PD to EC, CF is more sensitive by a factor of about 8, liquid water path by about a factor of 6, and SWCRE by about a factor of 4 in the Flower cloud state compared to the Sugar cloud state (Table 1). When splitting the change from PD to EC into contributions from LSC and GHG, nuances emerge. The Sugar cloud state is much less sensitive to the change in LSC (E_PD → E*) than to the change in GHG (E* → E_EC), and its overall response (E_PD → E_EC) can be largely attributed to the change in GHG. The response of the Flower cloud state to the change in LSC (E_PD → E*) is also weaker compared to its response to the change in GHG (E* → E_EC), but the contrast is less pronounced than in the Sugar cloud state. In the Flower cloud state, change in LSC contributes about 1/3 and change in GHG about 2/3 to the overall response of SWCRE E_PD → E_EC. Overall, we find that the Flower cloud state is more sensitive to climate change than the Sugar cloud state, and that in both, the change in GHG plays a greater role than change in LSC for their overall response.

Mechanisms of the Sugar and Flower Response to Change in LSC

In this section we analyze the mechanisms of the response of the Sugar and Flower cloud state to the change in LSC from PD to EC, and come to the conclusion that it is primarily the result of competition between moistening and stabilization of the atmosphere in a warmer climate, with a negligible contribution from a slowing of wind speed and a weakening of subsidence in the considered case. We also show that a change in radiative heating and cooling at the inversion contributes to the stabilization.

Figures 2a, 2c, 2e, and 2g show the response in total water flux, water vapor, liquid water, and cloud fraction E_PD → E*. The transition from the Sugar to the Flower cloud period appears as an increase in intermittence starting on the first sunset (d = 32.878). The intermittence reflects drizzle in individual Flowers, which is accompanied by temporary excursions from the mean response. Time-averaged profiles for the Sugar period (d = 32.25–32.878) are shown in Figure 3, and for the Flower period (d = 32.878–33.878) in Figure 4.

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

Warmer EC conditions cause a moistening of rising air parcels by water uptake from the surface, which results in stronger upward moisture transport E_PD → E* in both the Sugar and Flower cloud state (Figures 3a and 4a). In the Sugar cloud state, the enhanced moisture transport extends close to the inversion (Figure 3a). Liquid water follows suit and increases across the boundary layer (Figure 3b). Cloud fraction remains stable, except for a subtle reduction at the cloud base peak (Figure 3c). TKE production decreases very slightly in the cloud layer and more so just above the inversion (Figure 3d), indicating a suppression of penetrating updrafts by higher stability in EC conditions (Figure B1h). The effect of higher stability in EC conditions also manifests itself in reduced vertical and horizontal TKE in the cloud layer (Figures 3e and 3f).

In the Flower cloud state, the moisture transport enhancement E_PD → E* does not extend to the inversion but stalls below it (Figure 4a). This is a consequence of the increase in stability from PD to EC: as the inversion is located higher in the Flower cloud state compared to the Sugar cloud state (Figure 1h), and because atmospheric stability increases with altitude from PD to EC (Figure B1h), more TKE is converted to potential energy around and above the inversion in the Flower cloud state (Figure 4d) compared to the Sugar cloud state (Figure 3d). Liquid water tracks the response of the total water flux, and decreases around and above the inversion (Figure 4b). It is noteworthy that while cloud fraction also decreases around and above the inversion E_PD → E*, it increases in the cloud peak (Figure 4c). We interpret this as an adjustment of boundary layer circulation: while TKE production (Figure 4d) and vertical TKE (Figure 4e) decrease around and above the inversion, horizontal TKE shows a subtle increase in the cloud peak (Figure 4f). This indicates a strengthening of cloud top outflows at this level of the Flower cloud state which distributes liquid water horizontally due to increased stability, not unlike a lid. This will be discussed in more detail in Section 3.4.

Figures 5a, 5c, 5e, and 5g show the response E_PD → E* in stability, total radiative heating, long-wave (LW) heating, and short-wave (SW) heating. The response of LW heating (Figure 5c) turns from positive to negative values at the inversion at d ≈ 33.125. This is several hours after sunset (d = 32.878), when the parting of the mean individual cloud area time series between the simulations indicates the end of the Suger period and the start of the Flower period (Figure 1a). For the analysis of the radiative heating rates, we define the radiative period R_S as d = 32.25–33.125 and the radiative period R_F as d = 33.125–33.878. The time-averaged profiles of the radiative heating rates over the R_S and R_F periods are shown in Figure SF3 in Supporting Information S1.

[IMAGE OMITTED. SEE PDF]

Figure 5a reveals that air at the inversion becomes more stable in the Sugar period E_PD → E*, on top of the general increase in stability with altitude (Figure B1h). This is caused by positive values in the response of radiative heating just below and at the inversion in the R_S period (d < 33.125, Figure 5c), which reflect weaker radiative cooling at these levels (Figure SF3a in Supporting Information S1). This is followed by negative values in the response of radiative heating just below the inversion in the R_F period (d > 33.125, Figure 5c), which reflect an increase in radiative cooling at these levels (Figure SF3d in Supporting Information S1). Hence, from E_PD → E*, the response of radiative heating increases stability just below and at the inversion in the Sugar cloud state, and decreases stability just below and at the inversion in the Flower cloud state.

The response of total radiative heating (Figure 5c) E_PD → E* comprises the contribution of LW and SW radiation (Figures 5e and 5g). Juxtaposing the total radiative heating response (Figure 5c) and the LW response (Figure 5e) reveals that just below and at the inversion, it is the LW response which shapes the total response. The LW response is hence largely responsible for modulating the stability of the boundary layer at its top. The SW response (Figure 5g) just below and at the inversion enhances the LW response during the first daytime period in the Sugar cloud state, and offsets it in the second daytime period in the Flower cloud state. The contribution of SW is specific to the considered case, as it depends on the position of the sun and hence on latitude and time of year.

We now assess the role of change in horizontal wind speed and large scale subsidence. The geostrophic wind driving the simulations slows from PD to EC (Figures B1c and B1d). In the surface layer, where the wind speed drives surface fluxes, the reduction of the geostrophic wind speed is −0.55%. Nuijens and Stevens (2012) studied the geostrophic wind speed response of an idealized non-precipitating trade-wind cumulus case. They found that the inversion height depends strongly on geostrophic wind speed. A reduction by −50% in geostrophic wind speed, from 10 to 5 m s⁻¹, reduced the inversion height by −26%, from 2,700 to 2,000 m, once the simulations approached steady state. Using linear scaling between the relative wind speed change and inversion height change as an estimate, the reduction of geostrophic wind speed in the surface layer by −0.55% in our simulations would result in a reduction of inversion height by −0.29%, or by −7 m at an inversion height of 2500 m, to which the simulations converge (Figure 1h). Large scale subsidence weakens from PD to EC (Figure B1e). Between 1,300 m and 2,500 m, the range covered by the inversion over the course of the simulations (Figure 1h), it weakens on average by 10 m d⁻¹, which would raise the inversion height by 18 m over the 44 hr duration of the simulations. In combination, we estimate that the slowing of geostrophic wind and the weakening of subsidence E_PD → E* cause a rise in inversion height by 11 m over the course of the simulations, or 0.44%, relative to the inversion height of 2,500 m. In comparison, the LWP drops from E_PD → E* by −16.8% in the Flower cloud state (Table 1). We conclude that the response of cloud properties E_PD → E* in our simulations is primarily determined by the moistening and stabilization of the atmosphere from PD to EC, as opposed to the change in geostropic wind speed and subsidence. A stronger reduction in subsidence could, however, result in a larger adjustment of inversion height, in particular in conditions of steady-state and energy balance between the ocean and the atmosphere (Z. Tan et al., 2017).

Mechanisms of the Sugar and Flower Response to Change in GHG

In this section we discuss the response of the Sugar and Flower cloud state to the change in GHG from PD to EC, and come to the conclusion that it is caused by boundary layer stabilization due to enhanced radiative heating at the inversion.

Figures 2b, 2d, 2f, and 2h show the response in total water flux, water vapor, liquid water, and cloud fraction E* → E_EC. Time-averaged profiles for the Sugar period (d = 32.25–32.878) are shown in Figure 3, and for the Flower period (d = 32.878–33.878) in Figure 4. Both in the Sugar (Figure 3a) and Flower (Figure 4a) period, the total water flux E* → E_EC is suppressed, and the suppression increases from the surface toward the top of the boundary layer. Above z/z_i ≈ 0.7, the total water flux in E_EC is suppressed even below that in E_PD. This pattern is tracked by the response of liquid water (Figures 3b and 4b), cloud fraction (Figures 3c and 4c), and TKE production (Figures 3d and 4d), which are nearly unchanged at cloud base, but decrease with altitude in the cloud layer E* → E_EC to values even below those in E_PD. The suppression of these quantities in the cloud layer is weaker in the Sugar cloud state (Figures 3a–3d) and stronger in the Flower cloud state (Figures 4a–4d).

Figures 5b, 5d, 5f, and 5h show the response E* → E_EC in stability, total radiative heating, long-wave (LW) heating, and short-wave (SW) heating. The time-averaged profiles of the radiative heating rates in the R_S and R_F periods are shown in Figure SF3 in Supporting Information S1. Figure 5b reveals that the inversion becomes more stable E* → E_EC. This stabilization is continuous in the Sugar period (d = 32.25–32.878), and interrupted in the Flower period (d = 32.878–33.878). The stabilization is caused by positive values in the response of radiative heating around the inversion (Figure 5d), which reflect a reduction in radiative cooling at these levels (Figures SF3a and SF3d in Supporting Information S1). Comparison of the response in the total, LW, and SW radiative heating rate (Figures 5d, 5f, and 5h) reveals that the stabilization is caused by the response of LW heating, which during daytime is weakly offset by the response of SW heating. In the Sugar period, the offsetting by SW during daytime (around d = 32.7) is too weak to compensate the LW stabilization at the inversion, but in the Flower period (around d = 33.7), the SW and LW response are about equal, and the LW stabilization is temporarily suspended. Stabilization at the top of the boundary layer in a 4 × CO₂ scenario has been reported by Wyant et al. (2012) in large eddy simulations. It is noteworthy that Wyant et al. (2012) found that this stabilization was absent in super-parameterized climate simulations and cloud-resolving simulations. This suggests that such simulations might struggle to represent the cloud response via mesoscale organization to higher GHG levels.

Drivers of the Radiative Heating Response to Climate Change

In this section we discuss the drivers of the radiative heating response E_PD → E* (Figures 5c, 5e, and 5g) and E* → E_EC (Figures 5d, 5f, and 5h), with the objective of separating the effects of change in SST, T, q_v, and GHG, and the effects of cloud adjustment. The focus is placed on conditions at the inversion, where changes in radiative heating modulate boundary layer stability in addition to the general increase in atmospheric stability from PD to EC.

We calculate radiative heating rates offline from the simulations E_PD, E*, and E_EC using the Santa Barbara DISORT Atmospheric Radiative Transfer (SBDART, Ricchiazzi et al., 1998) model. The offline heating rates are calculated from hourly domain mean profiles of air temperature (T), water vapor (q_v), liquid water content (q_l), liquid water effective radius (r_eff), cloud fraction (CF), and O₃. Further quantities are SST and the mean atmospheric mixing ratios of CO₂, CH₄, and N₂O. As throughout the text, CO₂, CH₄, N₂O, and O₃ are collectively termed GHG. The quantities used in the offline radiative transfer calculation and the simulation from which they originate are denoted with, for example, E_PD(q_l, r_eff, CF, SST, T, q_v, GHG). In this case, all parameters originate from simulation E_PD. Quantities from different simulations can be combined in a heating rate calculation, which we denote, for example, with E_PD(q_l, r_eff, CF) + E*(SST, T, q_v).

The offline calculation is simplified compared to the online calculation by neglecting effects of the horizontal distribution of clouds. It is thus an imperfect but useful approach to illustrate the heating rate sensitivity to individual elements of climate change. The offline calculations and our notation are described in further detail and their results compared with the results of the online calculations in the Supporting Information S1.

Radiative Heating Response to Change in LSC

In the R_S period, the change in SST, T, and q_v from E_PD → E* creates a peak of positive values in the LW heating rate response just below the inversion (curve 1 in Figure 6a). Curve 2, which additionally accounts for the change in cloud properties E_PD → E*, exhibits a smaller peak. Hence adjustment of cloud properties offsets the LW heating rate response due to the change in SST, T, and q_v. The cloud adjustment is the slight increase in LWP from E_PD → E* in the Sugar period (Table 1) which allows for more efficient LW cooling in E*. The peak just below the inversion in curve 2 (Figure 6a) corresponds to the peak in the R_S period at the inversion in the LW (Figure 5e) and total (Figure 5c) radiative heating rate response. The stabilization around the inversion from E_PD → E* in the R_S period (Figure 5a) is hence largely driven by the change in SST, T, and q_v, whose effect is partially offset by cloud adjustment, in the form of a slightly higher LWP.

[IMAGE OMITTED. SEE PDF]

In the first daytime period (Sugar cloud state), the increase in SST, T, and q_v from E_PD → E* create negative values in the SW heating rate response just below the inversion (curve 1 in Figure 6c). Curve 2, which additionally accounts for the change in cloud properties E_PD → E*, does not exhibit this trough, indicating that the slight increase in LWP from E_PD → E* in the Sugar cloud state (Table 1) increases SW absorption in E*.

In the R_F period, the change in SST, T, and q_v from E_PD → E* creates negative values in the LW heating rate response around the inversion (curve 1 in Figure 6b). Curve 2, which additionally accounts for the change in cloud properties E_PD → E*, highlights the dominant role of cloud adjustment, which creates a large trough of negative values below the inversion and a peak of positive values above it. This dipole pattern is visible around the inversion in the R_F period in Figure 5e, and is caused by the reduction above and increase below the inversion in cloudiness in that period (Figure 2g). The reduction of liquid water above the inversion weakens LW cooling, whereas the increase in liquid water below the inversion strengthens it.

In the second daytime period (Flower cloud state), the increase in SST, T, and q_v from E_PD → E* create a dipole pattern in the SW heating rate around the inversion (curve 1 in Figure 6d), which is strengthened by the adjustment in cloud properties (curve 2). This dipole pattern reflects the dipole pattern in the Flower period in Figure 5g.

The individual contributions from the change in SST, T, and q_v from E_PD → E* to the change in the radiative heating rates are documented in Figure SF12 in Supporting Information S1. Water vapor q_v plays the leading role among these quantities, and its response (Figure 2c) shapes the response of both LW and SW heating in the Sugar and Flower periods.

Radiative Heating Response to Change in GHG

Figure 7 shows the contribution to the radiative heating rates from the change E* → E_EC in GHG (curve 1), and the contribution from the change E* → E_EC in GHG as well as of the adjustment in T, q_v, and cloud properties (curve 2). In the LW, these complement each other in driving positive values of the radiative heating response (Figures 7a and 7b). In the R_S period (Figure 7a), at the inversion, the contribution of GHG dominates, with adjustments playing a secondary role. In the R_F period (Figure 7b), adjustments contribute the majority to the LW heating response at the inversion, with GHG playing a minor role. This contrast between the R_S and the R_F periods arises because cloudiness below the inversion decreases less in the R_S period and more in the R_F period (Figure 2h). As a result, the reduction in LW cooling E* → E_EC at the inversion is weaker in the Sugar cloud state and stronger in the Flower cloud state. SW heating shows a very weak response from E* → E_EC (Figures 7c and 7d). GHG make a negligible contribution (curve 1), while adjustments cause slightly negative values in the response at the inversion of the Flower cloud state (Figure 7d).

[IMAGE OMITTED. SEE PDF]

Based on these findings, we can state the following about the stabilization at the inversion E* → E_EC (Figure 5b): In the Sugar cloud state, it is mainly stronger LW radiative heating from GHG, with a smaller contribution from adjustments, which warms and stabilizes air around the inversion. In the Flower cloud state, the importance of GHG and adjustments is reversed, with the adjustments, largely a reduction of cloud top cloudiness, being mainly responsible for warming and stabilizing the air around the inversion. The response of SW heating E* → E_EC plays a negligible role in the Sugar cloud state and a minor role in the Flower cloud state.

The contributions from individual GHG to the change in the radiative heating rates E* → E_EC are documented in Figure SF13 in Supporting Information S1. The increase in CO₂ from PD to EC plays the leading role.

The Role of Organization

In this section we study the role of mesoscale organization in the response of the Sugar and Flower cloud states to the change from PD to EC, using spectral analysis. Figure 8 shows spectra of TKE production and of vertical and horizontal TKE in the boundary layer, while Figure 9 shows spectra of the vertical flux of total water across the boundary layer and of SWCRE variance. With mesoscale we refer to wavelengths λ ≳ 6 km, the wavelength at which a trough approximately separates two peaks in the moisture flux and SWCRE variance spectra (Figure 9).

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

The spectra of TKE production and TKE (Figure 8) show a mesoscale peak that reflects the mesoscale circulation which gives rise to moisture aggregation and cloud organization in the trade cumulus regime (Bretherton & Blossey, 2017; Narenpitak et al., 2021). The TKE production peak is subdued in the Sugar cloud state (Figure 8a), indicative of its lower degree of organization, and pronounced in the Flower cloud state (Figure 8b), indicative of its greater degree of organization. The mesoscale peak also rises in the spectrum of the flux of total water across the boundary layer (Figures 9a and 9b). This peak corresponds to upward motion of moist air from the surface and downward motion of dry air from the inversion on the mesoscale. Cloudiness follows and brings about a corresponding peak in the spectrum of SWCRE variance (Figures 9c and 9d).

From E_PD → E*, TKE production exhibits a subtle decrease on the mesoscale, apparent at wavelengths λ ≳ 24 km (Figures 8a and 8b). The decrease is very subtle in the Sugar cloud state and carries over to subtle reductions on the mesoscale in vertical TKE (Figure 8c), horizontal TKE (Figure 8e), vertical moisture flux (Figure 9a), and SWCRE variance (Figure 9c). In the Flower cloud state, the response is more noticeable in TKE production (Figure 8b) and vertical TKE (Figure 8d), but it is absent in horizontal TKE (Figure 8f) and in the vertical moisture flux (Figure 9b). The SWCRE variance decreases in the Flower cloud state only between λ = 24–48 km (Figure 9d).

The reduction in vertical TKE without a reduction in horizontal TKE in the Flower cloud state E_PD → E* indicates a weakening of the vertical branch of mesoscale circulation without a weakening of the horizontal branch. We hypothesize that this is associated with the weakening of TKE production (Figure 4d) and vertical TKE (Figure 4e) at the top of the boundary layer, just below the inversion, which is accompanied by a subtle strengthening of horizontal TKE at its peak just below the inversion (Figure 4f). The increase in stability with altitude E_PD → E* (Figure B1h) appears to act as a lid which partially consumes and partially redistributes vertical TKE into horizontal TKE in the Flower cloud state at the inversion. We further hypothesize that the suppression of SWCRE variance on the mesoscale (λ = 24–48 km) in the Flower cloud state (Figure 9d) is caused by the suppression of cumulus clouds that rise above the mean inversion as a result of the increase in stability with altitude. These deeper cumulus clouds would organize on the mesoscale along with moisture aggregation. Their suppression would appear in the SWCRE variance spectrum on the mesoscale. These hypotheses have no bearing on our conclusions and we entrust them to future examination.

From E* → E_EC, TKE production and vertical and horizontal TKE (Figure 8), as well as the vertical moisture flux across the boundary layer and SWCRE variance (Figure 9) exhibit a pronounced decrease on the mesoscale. Both the vertical (Figures 8c and 8d) and horizontal (Figures 8e and 8f) components of TKE, hence both branches of mesoscale circulation, are suppressed. The suppression of the vertical moisture flux across the boundary layer (Figures 9a and 9b) and SWCRE variance (Figures 9c and 9d) is confined to the mesoscale. The suppression of the smaller mesoscale peak in SWCRE variance in the Sugar cloud state (Figure 9c) accompanies the smaller decrease of SWCRE in the Sugar cloud state (Table 1), and the suppression of the larger mesoscale peak in SWCRE variance in the Flower cloud state (Figure 9d) accompanies the stronger decrease of SWCRE in the Flower cloud state (Table 1).

In summary, we find that the change in LSC from PD to EC brings about a subtle, and the change in GHG a pronounced suppression of the Sugar and Flower cloud states on the mesoscale. The suppression of SWCRE variance is entirely confined to the mesoscale. The less organized Sugar cloud state has a lower, and the more organized Flower cloud state a higher, mesoscale peak in the SWCRE variance spectrum. The suppression of this mesoscale peak from PD to EC reflects a weaker overall reduction of SWCRE in the Sugar cloud state, and a larger overall reduction of SWCRE in the Flower cloud state. Thus, mesoscale organization amplifies the response of Tc clouds to climate change in the considered case, primarily as a result of the increase in GHG from PD to EC.

Base State Dependence

We address the question whether our findings hold when applying higher GHG levels in both PD LSC and EC LSC conditions. To this end we compare the change E_PD → E** and E* → E_EC. Quantitatively, the response of cloud fraction, liquid water path, and SWCRE to higher GHG levels E_PD → E** and E* → E_EC are nearly identical in the Sugar cloud state, and similar in the Flower cloud state (Table 1). The indicates base state independence to a GHG increase in the Sugar cloud state, and a weak base case dependence in the Flower cloud state. Table 1 also shows that higher GHG levels suppress liquid water path and cloud fraction significantly more in the more organized Flower cloud state compared to the less organized Sugar cloud state both from E_PD → E** and from E* → E_EC. This confirms that the finding that organization amplifies the response of Tc clouds to climate change is independent of the base state, within the range of the considered base states. Selected time series in E** are compared against those in E_PD in Figure SF14 in Supporting Information S1, and SWCRE variance spectra in Figure SF15 in Supporting Information S1. The responses in the time series and spectra E_PD → E** are consistent with those of E* → E_EC (Figures 1, 9c and 9d).

Lagrangian Simulations With Realistic Forcing

Figure 10 shows the time evolution in simulation L_PD and L_EC. The cloud field at sunrise, noon, and sunset is shown in Figure SF4 in Supporting Information S1. Additional time series are shown in Figure SF7 in Supporting Information S1; mean temperature, water vapor, relative humidity, and liquid water static energy profiles are given in Figure SF8 in Supporting Information S1. The transition from the Sugar to the Flower cloud state takes place during daytime, in sync with boundary layer deepening on the day of the simulations, 2 February 2020 (Narenpitak et al., 2021): the cloud field evolves from clustered cumulus clouds at sunrise (Figures SF4a and SF4d in Supporting Information S1) to non-precipitating Flower clouds at noon (Figures SF4b and SF4e in Supporting Information S1), which by sunset are precipitating and forming cold pools (Figures SF4c and SF4f in Supporting Information S1). From PD to EC conditions (L_PD → L_EC), a reduction in size of the cumulus clusters at sunrise, the Flower clouds at noon, and the Flower outflows and cold pools at sunset is discernible (Figure SF4 in Supporting Information S1). This size reduction appears in the mean individual cloud area time series in L_PD and L_EC (Figure 10a): from d = 32.625, the mean individual cloud area grows more rapidly in L_PD and less rapidly in L_EC, and remains lower in L_EC for the remainder of the simulations. The onset of growth in mean cloud area at d = 32.625 marks the transition from the Sugar to the Flower cloud state. From d = 32.625, LWP (Figure 10b) and CF (Figure 10c) drop from PD to EC. This reduction in LWP and CF causes a weaker SWCRE in L_EC relative to L_PD (Figure 10d). The reduction of SWCRE variance L_PD → L_EC is confined to wavelengths ≳3 km (Figure 11).

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

Figure 12 shows the response L_PD → L_EC in (a) stability, (b) total radiative heating, (c) LW heating, (d) SW heating, and (e) cloud fraction. Stability increases in a layer at the inversion up to d ≈ 32.625, when the cloud state transitions from Sugar to Flowers (Figure 12a). As the boundary layer deepens in the course of the transition, the inversion rises higher and the general increase in stability with altitude from PD to EC takes hold (see Section 3.1). The increase in stability around the inversion arises from positive values in the response of total radiative heating L_PD → L_EC (Figure 12b). Juxtaposing the response of LW (Figure 12c) and SW (Figure 12d) heating reveals that the positive change in total radiative heating at the inversion is dominated by the LW response, which is partially offset by the SW response. Cloudiness decreases L_PD → L_EC most strongly at cloud top, less in the Sugar cloud state (d < 32.625) and more in the Flower cloud state (d ≥ 32.625) (Figure 12e).

[IMAGE OMITTED. SEE PDF]

The response to the change from PD to EC conditions in L_PD → L_EC exhibits the same features as studied in E_PD → E* → E_EC (Section 3), but with the Sugar to Flower transition compressed into one daytime period by the ERA5 forcing, as in the study by Narenpitak et al. (2021). Equipped with the insight from Section 3.4, we conclude that the cloud and SWCRE response L_PD → L_EC is primarily caused by higher GHG levels warming the inversion, stabilizing the boundary layer, and suppressing trade cumulus mesoscale organization.

Discussion

The transition from the Sugar to the Flower cloud state in simulation L_PD and L_EC takes place at a particular phase relative to the diurnal cycle. The associated SWCRE and its response to change from PD to EC is therefore specific to the considered case and L_PD and L_EC therefore represent only one instance in a continuum of possible responses of the Sugar and Flower cloud states. Keeping these caveats in mind, it is nonetheless instructive to compare the change in cloud properties and SWCRE to existing work.

The response of boundary layer clouds to climate change has been studied with LES by Blossey et al. (2013) and Bretherton et al. (2013). The study by Blossey et al. (2013) conducted simulations of idealized climate change with different models. The study by Bretherton et al. (2013), using the approach of Blossey et al. (2013), explored individual elements of climate change, including a quadrupling of CO₂. These studies examined the steady state cloud response with simulations driven with constant, diurnally averaged insolation until convergence was achieved. For trade cumulus clouds, domain sizes of 9.6 km were used. The simulated clouds reached cloud top heights between 2 and 3 km. In this respect they are comparable to the Flower clouds simulated in this work, but without mesoscale organization. The results of Blossey et al. (2013) and Bretherton et al. (2013) are therefore a baseline against which the contribution of organization on larger scales to the climate response of trade cumulus clouds can be compared.

Table 2 compares the response of cloud properties and SWCRE to the change from PD to EC conditions (simulation L_PD → L_EC) in this work with the response of cumulus clouds to the change from simulation CTL → P2S (Blossey et al., 2013) and to the change from simulation CTL → 4CO₂ (Bretherton et al., 2013). CTL → P2S represents idealized climate change with an increase in SST by 2 K and a reduction in subsidence, but without a change in GHG. CTL → 4CO₂ represents the change due to a four-fold increase in CO₂ levels. In our work, SST increases by 1.8 K from PD to EC. Table 2 therefore gives the absolute change and the relative change normalized by change in SST from PD to EC.

Table 2 Comparison of the Trade Cumulus Cloud Response to Climate Change From Different Studies: Absolute Change and Relative Change per Change in SST of Cloud Fraction (CF), Liquid Water Path (LWP), and Short-Wave Cloud Radiative Effect (SWCRE)

	ΔCF	ΔLWP (g m−2)	ΔSWCRE (W m−2)	$$\frac{{\Delta }\mathrm{C}\mathrm{F}}{\mathrm{C}\mathrm{F}{\Delta }\mathrm{S}\mathrm{S}\mathrm{T}}$$ (% K−1)	$$\frac{{\Delta }\mathrm{L}\mathrm{W}\mathrm{P}}{\mathrm{L}\mathrm{W}\mathrm{P}{\Delta }\mathrm{S}\mathrm{S}\mathrm{T}}$$ (% K−1)	$$\frac{{\Delta }\mathrm{S}\mathrm{W}\mathrm{C}\mathrm{R}\mathrm{E}}{\mathrm{S}\mathrm{W}\mathrm{C}\mathrm{R}\mathrm{E}{\Delta }\mathrm{S}\mathrm{S}\mathrm{T}}$$ (% K−1)
CTL → P2S	−0.011	0.6	1.1	−3.1	1.0	−2.5
CTL → 4CO2	0.0	0.0	0.0	0.0	0.0	0.0
LPD → LEC	−0.015	−7.2	2.0	−6.5	−17	−7.1

Cloud fraction and LWP change from CTL → P2S to produce a reduction of SWCRE by 1.1 W m⁻² (Table 2). Cloud fraction, LWP, and SWCRE are insensitive to a four-fold increase in CO₂ (CTL → 4CO₂). This insensitivity to CO₂ is consistent with our finding that an increase in GHG acts primarily by suppressing mesoscale organization of trade cumulus clouds. Hence, for the purpose of this comparison, we take the response CTL → P2S as the baseline steady state trade cumulus response in absence of mesoscale organization, with or without any change in GHG. From L_PD → L_EC, SWCRE weakens by 2 W m⁻². In terms of a relative change normalized by a change in SST, SWCRE is nearly 3 times more sensitive to climate change in L_PD → L_EC compared to CTL → P2S (Table 2). Based on the analysis in Section 3.4, we interpret this increased sensitivity as primarily the result of GHG reducing SWCRE by suppressing mesoscale organization. In light of the caveats discussed above, the three-fold increase in sensitivity is specific to the considered case, and further study is needed to quantify the higher sensitivity of trade cumulus cloud to climate change in the presence of mesoscale organization.

Summary and Conclusions

In this work we investigated, using large-eddy simulations, the response to climate change of two recently identified manifestations of trade cumulus (Tc) organization, the less organized Sugar Tc state and the more organized Flower Tc state. Climate change is represented by applying to a given meteorology the change in large scale conditions (LSC) and greenhouse gas (GHG) levels from present-day (PD) to end-of-21st-century (EC) based on CMIP5 simulations with the CESM1(WACCM) model under the RCP8.5 emissions scenario. For the purpose of this work, we defined PD as the period 2016–2025 and EC as the period 2090–2099. LSC comprises sea surface temperature, air temperature, water vapor, zonal and meridional wind speed, subsidence, and horizontal advective tendencies of heat and moisture. GHG species are carbon dioxide, methane, nitrous oxide, and ozone. In the Caribbean, the region studied in this work, sea surface temperature increases by 1.8 K and CO₂ increases approximately two-fold from PD to EC in the RCP8.5 emission scenario simulated by the CESM1(WACCM) model.

Using Eulerian simulations with idealized initial conditions and forcings, we studied the response to climate change of the short-wave cloud radiative effect (SWCRE) from PD to EC in the Sugar and the Flower cloud state, and identified the responsible mechanisms. The simulations cover one daytime period in the Sugar cloud state and one daytime period in the Flower cloud state over the 44 hr simulation period. The simulated cloud system was first exposed to the change in LSC and then to the change in LSC and GHG.

The change in LSC from PD to EC subtly strengthens the SWCRE in the Sugar cloud state, making it more negative by −0.02 W m⁻², due to a small increase in liquid water path without a cloud fraction response. The Sugar cloud state is hence very insensitive to the change in LSC. In the Flower cloud state, the change in LSC weakens SWCRE by making it more positive by 0.44 W m⁻². We find that the response to the change in LSC from PD to EC is primarily the result of competition between the increase in specific humidity in the boundary layer, which supports cloudiness, and an increase in stability with altitude, which suppresses cloudiness. The shallower Sugar cloud state is less exposed to the increase in stability with altitude, giving moistening of the atmosphere a slight advantage, while the deeper Flower cloud state is more exposed to the increase in stability with altitude, giving it the upper hand. We estimate that the slowing of geostrophic wind speed and weakening of subsidence from PD to EC in the adopted climate change scenario make a negligible contribution to the response of cloud properties and SWCRE. Adjustments of the boundary layer to the changes in geostrophic wind speed and subsidence from PD to EC on time scales much longer than our simulations may, however, make a contribution that is not quantified in this work.

The increase in GHG from PD to EC weakens the SWCRE by 0.31 W m⁻² in the Sugar cloud state and by 1.06 W m⁻² in the Flower cloud state. The cause is stabilization of the boundary layer by stronger long-wave radiative heating at the inversion due to the increase in GHG from PD to EC. This weakens the mesoscale circulation that organizes Tc clouds, and results in a reduction of SWCRE. The suppression of SWCRE variance is confined to the mesoscale. The suppression is weaker in the Sugar cloud state, which has a weaker SWCRE variance on the mesoscale due to its lesser degree of organization, and stronger in the Flower cloud state, which has a larger SWCRE variance on the mesoscale due to its greater degree of organization. Organization hence amplifies the response of Tc clouds to climate change.

In combination, the changes in LSC and GHG from PD to EC weaken the SWCRE by 0.28 W m⁻² in the Sugar cloud state and by 1.50 W m⁻² in the Flower cloud state, with GHG responsible for all of the SWCRE weakening in the Sugar cloud state, and for two thirds of the weakening in the Flower cloud state. Cloud base cloud fraction and cloud mass flux through cloud base are insensitive to the changes in LSC and GHG in both the Sugar and the Flower cloud state. This is consistent with the cumulus-valve hypothesis, which postulates that a negative feedback between the cloud fraction and the cloud mass flux through cloud base buffers the climate response of trade cumulus clouds (Neggers et al., 2006). However, the strong response of the Flower cloud state to climate change suggests that it is not strongly buffered by the cumulus-valve mechanism, and that its response instead depends mainly on changes in cloud near the trade inversion.

To quantify the sensitivity of the Sugar and Flower cloud states to climate change in a present-day realistic meteorological setting, we used Lagrangian large-eddy simulations initialized and forced by the ERA5 reanalysis. The Lagrangian simulations follow a 24 hr trajectory in the Caribbean that was constructed from the ERA5 wind field, following the study by Narenpitak et al. (2021). The ERA5 meteorology compresses the transition from the Sugar to the Flower transition into one daytime period. Climate change is represented by applying the change in LSC and GHG under the RCP8.5 emissions scenario as in the Eulerian simulations. In these simulations, the SWCRE weakens from PD to EC by 2 W m⁻². Based on the comparison of this result with those of Blossey et al. (2013) and Bretherton et al. (2013), who studied the Tc response to climate change without mesoscale organization, we estimate that mesoscale organization increases the relative sensitivity of SWCRE normalized by change in sea surface temperature by a factor of 3 in the simulated Tc cloud case.

The presented results are specific to the particular Sugar and Flower cases considered, and not necessarily representative for the Sugar and Flower cloud states in general. One factor that prevents generalization is the specific phase between cloud evolution and the diurnal cycle in the considered cases. A different phase may change the SWCRE even without a change in cloud properties, but a different phase would also change cloud properties. The effect of the diurnal cycle on the evolution of the Flower cloud state was studied by Narenpitak et al. (2023). Furthermore, variability in the meteorology that governs the properties of the Sugar and the Flower cloud states will cause variability in their SWCRE, and therefore variability of its response to climate change. The presented results therefore offer only a first and narrow view of the continuum of possible cloud and SWCRE responses in Tc clouds with mesoscale organization. The continuum could be explored with an ensemble of simulations covering different Sugar and Flower cloud evolutions, constructed from reanalysis data. This would produce a statistically more robust quantification of the sensitivity of the Sugar and Flower cloud states to climate change.

Z. Tan et al. (2016) and Z. Tan et al. (2017) demonstrated the importance of maintaining a balanced surface energy budget with a coupled ocean model and interactive sea surface temperature in studying the steady-state response of boundary layer clouds to climate change. The evolution of the boundary layer and clouds from the Sugar to the Flower cloud state in our simulations does not proceed in steady state with the ocean. The simulations therefore use a prescribed sea surface temperature. Simulations with a coupled ocean model and an interactive sea surface temperature may yield additional insight.

Mesoscale shallow circulation has been found to be widespread in the Tc regime (George et al., 2023). Thus, the mechanism identified in this work by which mesoscale organization increases the sensitivity of Tc to climate change may be common in the Tc regime. It may express itself in the Tc aggregates that prompted the discovery of gross moist instability as the organizing mechanism of Tc clouds (Bretherton & Blossey, 2017), possibly in the Gravel and Fish Tc states (Stevens et al., 2020), and even in less prominent manifestations of Tc mesoscale organization, such as the simple clustering of Tc clouds. A greater sensitivity of the Tc response to climate change due to mesoscale organization would magnify the challenge in representing Tc feedbacks in climate simulations. Climate models do not represent mesoscale organization or differentiate between individual Tc states. The issue is compounded because observational studies (Cesana & Del Genio, 2021; Denby, 2020; Myers et al., 2021; Schulz et al., 2021; Scott et al., 2020) cannot detect the effect of increased long-wave heating with increasing GHG on clouds and their spatial organization, and hence cannot constrain associated cloud-climate feedbacks.

This work highlights the need for better understanding of cloud organization for cloud feedbacks. High resolution simulations and large-eddy simulations (LES) can provide the physical understanding of the various Tc states and their behavior in response to climate change, and help formulate improved representations of Tc in climate models that account for cloud organization. In this, particular consideration should be afforded to the dependence of self-aggregation on numerics and resolution in large-eddy simulations (Janssens, Vilà-Guerau de Arellano, van Heerwaarden, Van Stratum, et al., 2023). Still, constraints of Tc feedbacks from simulations, which need to connect large scale dynamics, mesoscale organization, boundary layer dynamics, and cloud processes, may remain encumbered by limited resolution, as in the case of climate models, or limited spatial and temporal extent, as in the case of high resolution simulations and LES. The promise of an independent approach to constrain Tc feedbacks might be found in the application of machine learning methods trained on observational data and reanalysis products (Denby, 2020) to climate simulations to identify individual cloud states, in combination with an emulator-based quantification of the cloud radiative effect based on LES (Feingold et al., 2016). Thus, more than one approach to constraining the Tc feedback to climate is emerging on the horizon, not only as an opportunity to improve understanding and quantification of cloud-climate interactions, but also as checks and balances on one another.

Appendix

Appendix - A

A1

21st Century Climate Change in the Caribbean

We represent 21st century climate change in the Caribbean with the three CMIP5 (Taylor et al., 2012) RCP8.5 (Van Vuuren et al., 2011) ensemble members r2i1p1, r3i1p1, and r4i1p1 (National Science Foundation; US Department of Energy; National Center for Atmospheric Research, 2017) produced by the first version of the Community Earth System Model—Whole Atmosphere Community Climate Model (CESM1(WACCM), Marsh et al., 2013). We averaged the CESM1(WACCM) RCP8.5 simulation ensemble for January and February over the region 10–22°N, 295–310°E, the Caribbean east of Barbados, over the PD and EC periods. The temporal averaging suppresses interannual variability. SST increases in the region by 1.8 K from PD to EC. Atmospheric data from CESM1(WACCM) CMIP5 were interpolated to the vertical grid of the fifth generation of the European Centre for Medium-Range Weather Forecasts (ECMWF) atmospheric reanalysis (ERA5, Hersbach et al., 2020). The ERA5 top altitude is approximately 72 km, while the CESM1(WACCM) CMIP5 top altitude is 54 km. Between the CESM1(WACCM) CMIP5 top altitude and the ERA5 top altitude, the top-level CESM1(WACCM) values are used. The CESM1(WACCM) CMIP5 PD and EC GHG profiles in are shown Figure A1, the PD and EC LSC profiles in Figure A2.

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

To assess the representativeness of the CESM1(WACCM) ensemble within CMIP5, we compare it to other CMIP5 models. Figure A3 shows the LSC in PD and their change to EC under the RCP8.5 emissions scenario in the region 10–22°N, 295–310°E across select CMIP5 models. These models include HadGEM-2A, MIROC5, and CCSM4, which reasonably reproduce observations in the shallow cumulus cloud regime relative to other CMIP5 models (Nam et al., 2012), and CESM1-CAM5, which uses CESM1 as in CESM1-WACCM. LSC changes consistently across models from EC to PD. Water vapor increases mainly in the boundary layer in each model. Temperature increases by a larger value in the mid-to upper troposphere compared to the surface, corresponding to an increase in tropospheric stability across the models. The change in LSC from PD to EC in CESM1-WACCM is broadly representative of the CMIP5 models which reasonably reproduce observations in the shallow cumulus cloud regime.

[IMAGE OMITTED. SEE PDF]

Appendix

Appendix - B

B1

Idealized Large Scale Conditions for Eulerian Simulations

The PD and EC initial and forcing profiles used in simulation E_PD, E*, and E_EC are shown in Figures B1a–B1g. Figure B1h documents the resulting different atmospheric stability in the simulations. PD profiles were constructed from ERA5 at 15.14°N, 50.44°W on February 2, 02h00m00s UTC. Between the surface and domain top at 10 km, the ERA5 temperature, water vapor, and subsidence profiles were smoothed and modified to cause a gradual deepening of the boundary layer with time in the simulations. At domain top, the pressure, temperature and water vapor profiles continuously connect to the ERA5 profiles to enable radiative transfer calculation above domain top. The zonal geostrophic wind speed was set to −8 m s⁻¹, the meridional geostrophic wind speed to 0. The horizontal advective tendencies of heat and moisture were constructed to balance the vertical advective tendencies in the initial temperature and moisture profiles due to subsidence. EC profiles were constructed from the PD profiles by adding the change from PD to EC in the CMIP5 CESM1(WACCM) simulations under the RCP8.5 emissions scenario (Appendix A). Apparent discontinuities in the horizontal advective tendency profiles arise from the relatively coarse vertical grid in the CMIP5 CESM1(WACCM) simulations.

[IMAGE OMITTED. SEE PDF]

The SST evolution in simulation E_PD, E*, and E_EC is shown in Figure B2. The SST was constructed to drive a gradual deepening of the boundary layer in the simulations until the end of the first daytime period, and an approximately stable boundary layer depth thereafter. The SST change from PD to EC reflects the SST warming by 1.8 K between the PD and EC in the CMIP5 CESM1(WACCM) simulations under the RCP8.5 emissions scenario (Appendix A).

[IMAGE OMITTED. SEE PDF]

Acknowledgments

We thank Dr. Marat Khairoutdinov, Stony Brook University, for providing the System for Atmospheric Modeling (SAM), and Dr. Eshkol Eytan, University of Colorado, for helpful advice on using SBDART. We thank three anonymous reviewers for improving this work. This study was supported by NOAA's Climate Program Office, Climate Variability and Predictability Program (GC19-303) and by the NOAA cooperative agreement NA22OAR4320151. The NOAA Research and Development High Performance Computing Program provided computing and storage resources. The statements, findings, conclusions, and recommendations are those of the authors and do not necessarily reflect the views of NOAA or the U.S. Department of Commerce.

Data Availability Statement

The version 6.10.10 of the System for Atmospheric Modeling (SAM, Khairoutdinov & Randall, 2003) is available from . The fifth generation of the European Centre for Medium-Range Weather Forecasts (ECMWF) atmospheric reanalysis (ERA5, Hersbach et al., 2020) is available through the Copernicus Climate Change Service (). CESM1(WACCM) model output for the RCP8.5 scenario (National Science Foundation; US Department of Energy; National Center for Atmospheric Research, 2017) is available from . The version 4.2.0 of the Hybrid Single Particle Lagrangian Integrated Trajectory Model (HYSPLIT, Stein et al., 2015) is available from . The summer 2002 release of SBDART (Ricchiazzi et al., 1998) is available at . The forcings and results of the simulations are available at .

References

Blossey, P. N., Bretherton, C. S., Cheng, A., Endo, S., Heus, T., Lock, A. P., & van der Dussen, J. J. (2016). CGILS Phase 2 LES intercomparison of response of subtropical marine low cloud regimes to CO₂ quadrupling and a CMIP3 composite forcing change. Journal of Advances in Modeling Earth Systems, 8(4), 1714–1726. [DOI: https://dx.doi.org/10.1002/2016MS000765]