Kilometer‐Scale Simulations of Trade‐Wind Cumulus Capture Processes of Mesoscale Organization

The EUREC4A (Elucidating the role of clouds–circulation coupling in climate) and ATOMIC (Atlantic Tradewind Ocean‐Atmosphere Mesoscale Interaction Campaign) joint field campaign gathered observations to better understand the links between trade‐wind cumulus clouds, their organization, and larger scales, a large source of uncertainty in climate projections. A recent large‐eddy simulation study showed a cloud transition that occurred during this field campaign, where small shallow clouds developed into larger clouds with detrainment layers, was caused by an increase in mesoscale organization generated by a dynamical feedback in mesoscale vertical velocities. We show that kilometer‐scale simulations with the Met Office's Unified Model reproduce this increase in mesoscale organization and the process generating it, despite being much lower resolution. The timing of development is associated with large‐scale convergence. Sensitivity tests with a shorter spin‐up time, to reduce initial organization, still have the same timing of development and sensitivity tests with cold pools suppressed show only a small effect on mesoscale organization. Mesoscale organization and clouds are sensitive to resolution, with higher resolution simulations producing weaker organization and less cloud, which affects changes in net radiation. The clouds also have substantial differences to observations. Therefore, while kilometer‐scale simulations can be useful for understanding processes of mesoscale organization and links with large scales, including responses to climate change, simulations will still suffer from significant errors and uncertainties in radiative budgets.

The patterns of cloud organization have different cloud radiative effects and can largely be distinguished by large-scale parameters (surface winds and inversion strength) (Bony et al., 2020). So the objective of getting models to represent trade-wind cumulus well enough to be used for climate studies could be viewed as making sure models represent the different regimes of cloud organization and the variability between these regimes. However, the four patterns of trade-wind clouds represent extremes in a more continuous distribution (Janssens et al., 2021) and trade-wind clouds can be classified as a hierarchy of different regimes with distinct cloud structures and radiative effects (Denby, 2020).
To understand the physical processes generating cloud organization, high-resolution large-eddy simulations (LES) have been used to model trade-wind cumulus clouds. Bretherton and Blossey (2017) showed that organization in trade-wind cumulus in the Pacific, downstream from the transition from stratocumulus to shallow cumulus, can be generated solely by a dynamical feedback in latent-heat driven mesoscale vertical velocities: convection preferentially develops in moist regions and once convection develops, the circulation generated by the convection acts to converge moisture toward the existing convection, making moist regions moister and dry regions drier. In the absence of other limiting factors, shallow-cumulus convection is intrinsically unstable to growth from this latent-heating driven circulation (Janssens et al., 2022). While not crucial for the development of mesoscale organization, Bretherton and Blossey (2017) also showed that the interaction between clouds and radiation can speed up the initial development of mesoscale organization. It has also been shown that the interaction between clouds and radiation is important in developing the detrainment layers that are distinctive of the flowers regime (Vogel et al., 2020) even if they are not crucial for the development of organization. Narenpitak et al. (2021) simulated the 2nd February case from EUREC 4 A-ATOMIC using an LES driven by forcings following a Lagrangian trajectory (Lagrangian LES) and showed that the development of the flowers was associated with the development of mesoscale organization generated by the mesoscale vertical velocities, consistent with Bretherton and Blossey (2017). In this study, we look at the same, 2nd February, case study using high-resolution nested simulations with the Met Office's unified model (UM).
Our simulations are at a high enough resolution to allow an explicit representation of convection, but at a much lower resolution than needed to resolve cloud processes such as entrainment, and at a much lower resolution than the LES previously used to study convective aggregation. While our simulations will not be as good at representing the cloud processes as LES, the coarser resolution enables a much larger domain size at similar or reduced computational cost. A large domain size has been shown to be important for correctly capturing mesoscale organization of trade-wind cumulus. For example, Vogel et al. (2020) used simulations with a 51.2 × 51.2 km 2 horizontal domain size to produce mesoscale organization and contrasting simulations with a 12.8 × 12.8 km 2 horizontal domain size that did not support mesoscale organization. Bretherton and Blossey (2017) suggested that the larger domain size could be important for supporting the moisture aggregation and would explain why previous simulations on smaller domains by Seifert and Heus (2013) instead saw precipitation and cold pools leading to organized convection.
High-resolution LES are not currently possible at much larger domain sizes so are difficult to use to represent interactions between cloud organization and the large scale, or spatial variations in cloud organization over larger 3 of 27 scales. These larger scale processes could be well represented in kilometer-scale simulations, provided they can represent the mesoscale organization. Since kilometer-scale simulations are being suggested for climate-change projections due to their improvements in the representation of precipitation (Kendon et al., 2014(Kendon et al., , 2019Slingo et al., 2022;Stevens, Acquistapace, et al., 2020), it is important to assess whether kilometer-scale simulations can capture these processes closely linked to uncertainties in climate sensitivity. Recently, Beucher et al. (2022) showed that kilometer-scale simulations can predict the different regimes of organization in trade-wind cumulus. In this study, we aim to address whether our simulations can represent the processes generating mesoscale organization and can therefore be used to better understand interactions between cloud organization and larger scales.
The layout of the paper is as follows. In Section 2 the model simulations are described. In Section 3.1 we introduce the idea of "quasi-Lagrangian domains" extracted from the UM simulations to better follow the development of mesoscale organization and compare with LES. In Section 3.2 the mesoscale organization in the UM simulations is quantified across the different resolutions. In Section 3.3 we quantify the processes responsible for the mesoscale organization following the analysis of Bretherton and Blossey (2017) and Narenpitak et al. (2021). In the following sections we quantify the effects of spin up (Section 3.4) and cold pools (Section 3.5) on mesoscale organization in our simulations. In Section 3.6 we quantify the sensitivity of large-scale averages to resolution, spin up, and cold pools. In Section 4 we summarize the results from this study.

Model Data
We ran simulations with the Met Office's UM using the third iteration of the regional atmosphere and land configuration (RAL3). RAL3 is designed for nested models with resolutions fine enough for convection to be explicitly represented by the model dynamics and therefore has no convection parametrization. The previous version of the regional atmosphere and land configuration (RAL2) is described in Bush et al. (2022). Key differences between RAL2 and RAL3 of relevance here are that RAL3 uses the two-moment Cloud-AeroSol Interacting Microphysics (CASIM) parametrization described by Miltenberger et al. (2018) and the parametrization of cloud fraction as a function of the gridbox mean state is done by the bimodal cloud scheme described by Weverberg et al. (2021).
Kilometer-scale simulations (1.1, 2.2, and 4.4 km horizontal resolution) were initialized at 00Z 1st February and run for 48 hr. The initial conditions and boundary conditions are from ERA5 (Hersbach et al., 2020). ERA5 has an approximately 31 km horizontal resolution so clouds are not resolved in the initial conditions and are spun up in the UM simulations. Figure 1a shows the extent of the model domain. Higher-resolution simulations (300 and 500 m horizontal resolution) were then nested within the 1.1 km simulation, fixed at the region upstream of Barbados observed by the HALO (High Altitude and LOng Range Research Aircraft) and other platforms during the EUREC 4 A-ATOMIC field campaign, with the boundary conditions updated every 30 min. The box in Figure 1a shows the extent of the nested domain and the approximate circle flown by the HALO during the field campaign. Table 1 gives a summary of the parameters that vary between simulations. To account for the gray-zone, resolutions where boundary layer eddies are neither fully resolved or fully parametrized, the turbulent mixing scheme includes a resolution-dependent blending of the non-local fluxes (Boutle et al., 2014) as part of the RAL3 configuration. Otherwise, each simulation used the same configuration (RAL3) and vertical resolution (70 hybrid-height levels decreasing in resolution from the surface up to 40 km). Figure 1b shows the vertical resolution in the lowest 5 km. Figure 2 shows a comparison of the UM simulations with data from the GOES-16 satellite focused on the region of the inner domain. The model data is shown as the total outgoing longwave flux, whereas the satellite data is shown as brightness temperature from channel 11, which captures the water-vapor window in the infrared range. While these two fields will not be exactly the same, they will capture a lot of the same features. Model data processed with a satellite simulator to mimic the brightness temperature from channel 11 does show that the following conclusions are consistent (see Appendix A).

Results
Visually, all the simulations produce somewhat similar transitions in the cloud organization to those seen in the observations: initially small, scattered clouds with some hints of lines (06Z) aggregate and develop into larger, more circular, cloud patches (12Z-18Z) followed by less aggregation and more cloud free air (00Z). In the terminology of Stevens, Bony, et al. (2020), there is sugar at 06Z developing into flowers at 12Z and 18Z followed by sugar again at 00Z.
While we are interested in the development of mesoscale organization in the UM in this study, it is worth pointing out that these simulations have some strong differences compared to the observations (which are consistent in the simulated satellite imagery in Appendix A). The UM simulations produce too much cloud (too little cloud-free air) at all times. The satellite observations show there to be linear features as early as 00Z, before the flowers develop, most notably the structure in the top half of the domain near the center, but also weaker features in other parts of the domain. The UM does show some line-like structures, but at higher resolution they are fairly indistinct from the excessive amounts of very shallow clouds and at the lower    Figure 1). The bottom row shows the 11 μm brightness temperature from the GOES-16 geostationary satellite regridded to the same resolution (4.4 km) and area as the model data.
resolution they are too large, and while there are more regions of clear air at lower resolution, they still cover too much of the domain. Similar resolution sensitivity is seen as the flowers develop (12Z-18Z): the 4.4 km simulation roughly captures the size of the flower structures but there are too many, too close together. At increasing resolution, the flowers are more broken up and there are too many very shallow clouds in between.
The 1.1 km simulation has a large amount of cirrus cloud in the south of the domain at 00Z which is not seen in the GOES satellite. This cirrus cloud is also present in the 300 and 500 m simulations which are nested within the 1.1 km simulation. The 2.2 km has cirrus clouds near the inner-domain region but much less than the 1.1 km simulation, and the 4.4 km simulation has no cirrus clouds nearby. The GOES satellite doesn't show any cirrus clouds nearby at this time so they are incorrectly being simulated by the UM. However, the cirrus clouds only appear in the area of interest toward the end of the simulations so will have little effect on the cumulus clouds and organization which is the focus of this paper.

Quasi-Lagrangian Domains
The fixed domain in Figure 2 is limiting because we are not following the air mass as it develops and the air flowing into the domain has a strong influence on the cloud organization seen. This can be seen by the fact that the nested, inner-domain, simulations have very similar cloud structures to the driving 1.1 km simulation in Figure 2. In the following sections, we will focus our analysis on what we call "quasi-Lagrangian domains" to better follow the development of the clouds. Figure 3 shows an example of extracting a quasi-Lagrangian domain. A trajectory is calculated using Lagranto (Sprenger & Wernli, 2015;Wernli & Davies, 1997) with hourly wind output from the model and a fixed height of 500 m (same height as Narenpitak et al. (2021)). The trajectory is initialized at 57.5 W, 13.5 N (the center of the inner domain) at the end of the simulation (T + 48 hr/00Z 3rd February) and tracked back to the start of the simulation. The subset of the domain extracted is taken to be the same size as the inner domain. At T + 48hr, the domain is just the subset of the kilometer-scale domain that overlaps with the inner domain. At other times the location of this domain is translated to follow the trajectory and data linearly interpolated to the new grid. The dashed line in Figure 3 shows the trajectory and extracted data for the 1.1 km simulation of 2nd February.
The method is similar to how Tomassini et al. (2017) extracted averages in a box following a model trajectory to compare with LES driven by forcings along the same trajectory. The quasi-Lagrangian domains are designed to be a rough equivalent of a Lagrangian LES, such as the simulations by Narenpitak et al. (2021) of the 2nd February case study. A key first test for the quasi-Lagrangian domain is that the developing cloud features remain within the domain to show we are tracking a coherent patch of cloud development. Animations of the model data on the quasi-Lagrangian domains does show most cloud features rotating around, but remaining within, the domain (see Supporting Information S1). In situations where wind shear or divergence had a stronger impact on displacing the cloud features from the feeding boundary-layer airmass this quasi-Lagrangian domain approach may not work and we would need a more sophisticated approach to track the boundaries of the cloud development. Figure 4 shows the same satellite comparison as Figure 2 but for the quasi-Lagrangian domains extracted from the kilometer-scale simulations. For the satellite data, we have interpolated it to the quasi-Lagrangian grid of the 4.4 km simulation. The specific choice of the 4.4 km quasi-Lagrangian grid makes very little difference to the figure. Figure 5 shows the trajectories used for extracting quasi-Lagrangian domains from simulations with different resolutions (and sensitivity tests used in later sections in this paper) and shows that the displacement between different trajectories is much smaller than the size of the domain. Another potential issue is that the trajectories from the simulations may differ from the true trajectories of the atmosphere; however, animations of the satellite data on the quasi-Lagrangian domains also show most cloud features remaining within the domain (see Supporting Information S1) indicating that we are also following the motion of the observed clouds with these trajectories. Shown is an example of extracting a quasi-Lagrangian domain following a trajectory from the 1.1 km simulation. The dashed line shows a back-trajectory from the center of the inner domain at the end of the simulation (T + 48 hr/00Z 3rd February) and with a fixed height of 500 m. The domain extracted is the same size as the inner domain but is translated to follow the trajectory. The three grids show snapshots of total column water for this subdomain extracted from the 1.1 km simulation. The teal circle is the circle flown by the HALO during the EUREC 4 A-ATOMIC field campaign.

of 27
The quasi-Lagrangian view in Figure 4 is useful because it allows us to see the cloud development following the clouds, even in the observations. We see from the satellite data that the flower clouds develop later along the trajectory: the developing line-like cloud patterns are still present in the trajectory-following region at 12Z (GOES at 12Z in Figure 4) whereas there are already flowers present in the region of the inner domain at this time (GOES at 12Z in Figure 2). This means that if we were to only look at the clouds at a fixed position we would underestimate how rapidly the flowers develop and decay.
Compared with the fixed domain evaluation in Figure 2, the UM struggles to represent the earlier development of the clouds upstream: at 06Z and 12Z the satellite shows line-like cloud features that intensify during this time whereas the UM produces mostly circular patches of cloud. There is a strong resolution dependence in the development of these cloud patches with lower resolution yielding larger cloud patches. This resolution dependence also affects the development of the flowers at 18Z, with the lower resolution producing larger flowers, while the flowers in higher resolution simulations look most similar to the satellite, the opposite of what is seen at the region of the inner domain in Figure 2. Nevertheless, the simulations do still produce a transition in cloud organization similar to that observed even if the clouds themselves look unrealistic.

Mesoscale Organization
In this section we quantify the variation of mesoscale organization in the UM simulations in terms of horizontal variations in total column water. In Bretherton and Blossey (2017) and Narenpitak et al. (2021) the transition to mesoscale organization is seen by the emergence of mesoscale anomalies in total column water where mesoscale anomalies are defined as anomalies relative to the large-scale (domain) mean over 16 × 16 km horizontal blocks. Cumulus-scale anomalies are then defined as the anomalies at the grid-scale relative to the mesoscale anomalies such that a quantity, such as the total water content (q t ), could be decomposed as   Figure 6a shows the mesoscale total column water (large-scale mean plus mesoscale anomaly) averaged over quartiles for the quasi-Lagrangian domains. All the simulations show a similar behavior with a weak decrease in total column water initially followed by a strong increase and finally decreasing or leveling out. Although all quartiles follow the same pattern, the magnitude of changes are not the same. Bretherton and Blossey (2017) used the difference between the average total column water for the moistest and driest quartiles to show the strength of mesoscale organization. Figure 6b shows the difference between the average total column water for the moistest and driest quartiles in our simulations. The mesoscale organization strongly increases during the initial development of the flowers from around 07Z-12Z and levels out before dropping off from 19Z as the flowers dissipate.
The gray lines in Figure 6 shows the same averages but calculated fixed on the region of the inner domain. This demonstrates the added value of viewing the cloud development from a quasi-Lagrangian perspective. The region of the inner domain largely shows a steady increase in mesoscale organization from the start of the day before leveling out and then decreasing; however, following the cloud development shows that the increase in mesoscale organization is stronger and faster, consistent with the differences seen between Figures 2 and 4. The strength of the mesoscale organization is stronger with lower resolution. This is particularly noticeable for the 4.4 km simulation which also has a second period where the mesoscale organization increases at around 15Z. This is consistent with seeing larger flower structures with lower resolution in the satellite comparisons (Figures 2 and 4).
The increase in mesoscale organization during the development of the flowers is in agreement with Narenpitak et al. (2021); however, the timing and magnitude are different. In Figure 3 of Narenpitak et al. (2021) there is no initial contrast in the mesoscale-averaged total column water, only starting to develop after 14Z and continuing to increase into the following day. In our simulations the contrast in the mesoscale-averaged total column water prior to the development of the flowers is about as strong as at the end of the simulations in Narenpitak et al. (2021) and is decreasing by the end of the day. Instead, it looks like the UM simulations have a better representation of the mesoscale organization: in Figure 2 of Narenpitak et al. (2021) the satellite shows more structure at the start of the simulation and develops a flower structure earlier and stronger than the LES, whereas the UM simulations develop the flower structures at a similar time to the satellite observations.
The differences in mesoscale organization in our simulations and those of Narenpitak et al. (2021) can only partly be explained by the different domain. To account for uncertainty from using a smaller domain size, we calculated the range of values of mesoscale organization possible by using a 192 × 192 km 2 grid within the quasi-Lagrangian domain of the 1.1 km simulation which is shown as the shaded region in Figure 6. The smaller domain typically estimates a larger mesoscale organization. Using the same trajectory origin and (smaller) domain size as Narenpitak et al. (2021) gives us quasi-Lagrangian domains that are roughly a subsection in the top right of the quasi-Lagrangian domains in Figure 4 and slightly upstream (see Figure 5). The mesoscale organization for these domains are shown in Figure 6b by the lines labeled "RHB." The contrast in total column water between quartiles is smaller for these domains initially but still much stronger than in Narenpitak et al. (2021). The timing of the peak in mesoscale organization is later in these domains but still earlier than in Narenpitak et al. (2021) and decreasing by the end of the day.
The growth in mesoscale organization in our simulations is largely linear, although the 4.4 km does show a second period of development, whereas the simulations in Bretherton and Blossey (2017) and Narenpitak et al. (2021) show nonlinear development. One difference is that the development of mesoscale organization in our simulations is over a much smaller time period than in Bretherton and Blossey (2017) and is more comparable to the initial linear development in the first few hours of development in their simulations. The nonlinear development in Narenpitak et al. (2021) is more similar to our simulations but with a second period of development, similar to our 4.4 km simulation.

Processes Generating Mesoscale Organization
In this section we look at the processes responsible for generating mesoscale organization in the UM simulations. Since mesoscale organization can be described as the development of moist and dry regions, Bretherton and Blossey (2017) derived a budget for the mesoscale anomalies of total water content, ′′ t , and showed that it could be understood as the combination of two processes: (a) the advection of mesoscale anomalies of moisture and (b) the "column process" described by Bretherton and Blossey (2017) as "the combined moistening effect of the moist processes and diabatically induced vertical advection across the horizontal-mean moisture gradient." The equation where u is the three-dimensional wind, and where the terms in Equation 4 can be described in three parts (as in Narenpitak et al. (2021)). The first term, P, represents the non-advective fluxes of moisture from precipitation and surface fluxes, and the [⋅] m denotes a mesoscale average. The second term, cu qt , represents the cumulus-scale fluxes of moisture. Narenpitak et al. (2021) included the vertical (B v ) and horizontal (B h where w and v are the vertical and horizontal winds respectively. The inclusion of the horizontal fluxes by Narenpitak et al. (2021) is because, unlike (Bretherton & Blossey, 2017), they did not to apply a scale separation to simplify these terms because these terms are non-negligible in their simulations. See appendix D of Narenpitak et al. (2021) for full details. The final term, is the vertical advection of large-scale moisture by the mesoscale winds and was shown to be the dominant term in the column process by Narenpitak et al. (2021).
We have calculated each of these terms for the quasi-Lagrangian domains. Figures 7 and 8 show each of these terms, except the non-advective fluxes, as vertical profiles at 10Z (i.e., when the mesoscale aggregation is rapidly increasing) and vertically integrated as a function of time respectively. Only the 1.1 and 4.4 km simulations are shown for clarity since the 2.2 km simulation falls between the two. (2017) and Narenpitak et al. (2021), the vertical advection of large-scale moisture by the mesoscale winds, C, is the most important process for increasing aggregation. It is responsible for moistening the moistest regions and drying in other regions. The advection of mesoscale anomalies, A, also acts to oppose the aggregation by removing moisture from the moistest regions but is weaker than C during the middle of the day when we see aggregation increasing. The cumulus-scale fluxes (B v and B h ) have negligible contributions to the column averages.

Consistent with Bretherton and Blossey
To quantify the effect of the non-advective fluxes (P) on mesoscale organization, we use hourly accumulated rainfall and a tracer that accounts for the net effect of surface evaporation. The tracer for surface evaporation comes from a set of moisture tracers designed to represent a Lagrangian budget of specific humidity, where the rate of change of specific humidity (q) is where the sum over i represents all the non-advective processes modifying the specific humidity of an air parcel in the UM. The total specific humidity is then given by where each q i is represented by an individual tracer that accumulates the changes from a single process in the UM at each timestep and q adv represents the initial field of specific humidity which is passively advected by the UM.
The tracers are initialized at 00Z 2nd February and then tracked until the end of the simulation. The distribution of the passive tracer (q adv ) then tells us about how water vapor changes due to advection of the initial water vapor and the other tracers tell us about the net effect of individual sources and sinks of water vapor on air parcels. Most of the tracers account for changes between water vapor and other phases of water. The processes that affect the total water content are evaporation of water into the atmosphere from the surface and removal of water from the atmosphere by precipitation.
Evaporation of water from the surface is accounted for by the tracer which tracks changes from the boundary-layer parametrization. The boundary-layer parametrization only redistributes moisture vertically (Lock et al., 2000), so the increment to the tracer tells us about how subgrid-scale turbulent mixing redistributes moisture vertically. More relevant here, is that this tracer also accounts for surface fluxes of moisture, so the column average of a single timestep increment applied to this tracer tells us about the net effect of surface fluxes on column moisture. Since the tracer is advected following the flow the effects of advection on column moisture could become more important in this tracer at later times if, for example, wind shear acts to separate evaporated moisture from the column of interest. This is not important over the short integration time in our simulations as the horizontal distribution of this tracer stays reasonably uniform.
To account for precipitation, we just use the model output of hourly accumulated rainfall rather than a tracer. This has the disadvantage that the output precipitation is associated with where it fell rather than the air parcel that it fell from; however, we do see that most of the hourly accumulated rainfall is still associated with the moistest quartile, so this is not a problem here. Figure 9 shows the change from the previous output (i.e., rate of change per hour) of column-averaged quantities by quartile of total-column water. The aggregation can mostly be explained by the dynamics rather than direct effects of non-advective moisture fluxes (surface evaporation and rainfall) because the pattern of differences between the quartiles for specific humidity is mostly explained by the passive tracer (q adv ). The boundary-layer fluxes largely just moisten all quartiles evenly with a small opposition to aggregation early on and small enhancement to aggregation later. The rainfall acts to oppose aggregation after it develops since it is dominated by removing moisture from the moistest quartile, but it is also small compared to the differences seen in the passive tracer. This small contribution of non-advective moisture fluxes to aggregation is consistent with Narenpitak et al. (2021). Bretherton and Blossey (2017) found that the aggregation in their simulations slowed down when the advection of the mesoscale anomalies became strong enough, due to the mesoscale anomalies being stronger, to oppose the aggregation from the column process. However, in our simulations, the changes in aggregation largely follow the changes in the strength of moistening from mesoscale vertical velocities on the background moisture gradient (C). The advection of mesoscale anomalies of moisture acts to oppose aggregation throughout the day, with a moderate increase when the strength of these anomalies increases, but it is C that varies more with a strong increase from 09Z-12Z which is associated with the increasing variance in total column water and a strong decrease toward the end of the day associated with the variance in total column water decreasing.

Influence of Spin-Up Time on Mesoscale Organization
Our simulations have stronger and earlier development of mesoscale organization compared to the LES simulations of Narenpitak et al. (2021); however, the LES simulations of Narenpitak et al. (2021) do develop flowers later than the observations and also appear to be too uniform. This could be related to the different ways the simulations are spun up. Our simulations are initialized at 00 UTC 1st February by interpolating ERA5 data to the UM model grid, whereas the simulations of Narenpitak et al. (2021) are initialized at 02 UTC 2nd February using information from ERA5 at a single point on a trajectory. Therefore, our simulations have an extra 26 hr of spin up and also already include some horizontal variability from the ERA5 initial conditions. To test whether we could delay the development of mesoscale organization and flowers by having a less organized state in our simulations at the time when the flowers develop in the observations, we re-ran our simulations initialized 24 hr later (00Z 2nd February) so that the model has had less time to spin up any clouds and mesoscale structure from the ERA5 initial conditions. Figure 10a shows outgoing longwave fluxes for these simulations with a later start time. The comparison is shown for quasi-Lagrangian domains in the same way as Figure 4. The difference in trajectory location for the different start times, as with the different resolutions, is small compared to the size of the domain (see Figure 5), so any differences seen are due to the simulation of the cloud development.
There are strong differences between the simulations at different start times in the cloud organization. The simulations initialized later have much smaller cloud structures. This makes sense because at earlier lead times the model is still spinning up cloud structures and is noticeably smoother at 06Z. The smaller cloud structures persist even as the later initialized model runs start to produce similar cloud structures at 12Z. This leads to the flower structures at 18Z being smaller and perhaps looking more like the cloud structures in the observations. The simulations initialized later also do better at producing line-like cloud structures at 12Z, particularly at lower resolution. Both sets of simulations look more similar at the end of the day when the flower structures have broken down.
Despite the differences in the cloud structures for the different initialization times, the general development of the clouds follows the same pattern, with little organization initially (06Z), followed by development into larger cloud structures (12Z-18Z) which are mostly gone by the end of the day (00Z), although the 4.4 km simulation does retain some larger cloud structures at this time.
This consistency in the cloud development can also be seen in the total column water. Figure 11a shows the mesoscale total column water averaged over quartiles for the quasi-Lagrangian domains from the simulations with a later start time. For comparison, the lines for the earlier start time simulations (from Figure 6) are shown in gray. The strength of mesoscale organization (Figure 11c) is always weaker in the simulations initialized later, but the timing of the increase is fairly similar. The distribution of total column water only converges toward the end of the day when the mesoscale organization decreases more rapidly in the simulations initialized earlier.
The similarity in the timing suggests that the development of the flowers is related to large-scale dynamics and the effects of the diurnal cycle rather than some threshold reached in the cloud development. However, reducing the spin-up time does have a strong effect on the details of the flower structures in these simulations.

Influence of Cold Pools on Mesoscale Organization
In Section 3.3 we showed that the mesoscale organization stops increasing because the effects of mesoscale vertical velocities, through the column process, decreases. One possible reason for decreasing strength in the column process is the development of cold pools. The flowers are associated with precipitation in the moist regions which develops cold pools. If the cold pools sufficiently suppressed convection in these regions, then the circulation making the moist regions moister would be stopped. To identify cold pools in our simulation we use the tracer that accumulates changes in moisture due to the evaporation of rainfall in the microphysics parametrization (CASIM) from the tracers described in Section 3.3. As with the tracers in Section 3.3, the tracer is initialized at 00Z 2nd February. Figure 12 shows two snapshots of this tracer from the quasi-Lagrangian domain extracted from the 1.1 km simulation. Also shown is the equivalent potential temperature, which can also be used to identify cold pools. For example, Touzé-Peiffer et al. (2022) identified cold pools in soundings from the EUREC 4 A-ATOMIC field campaign by using equivalent potential temperature to diagnose the mixed-layer height. We see that at the earlier time (12Z) the tracer has a close correspondence to equivalent potential temperature when the cold pools are initially developing directly underneath the clouds. At the later time, the cold pool activity is less distinct in equivalent potential temperature, but the tracer shows that a large amount of the domain has been affected by cold pools over the course of the day.
To test whether the cold pools influence the mesoscale organization, we re-ran our simulations with evaporation of rain switched off, which stops the cold pools from forming. Note that the evaporation of rain being switched off also stops rain re-evaporating at levels above cloud base which does happen in the UM simulations. This effect is small but may explain some of the differences seen in Section 3.6 for the simulations without evaporation. Figure 11b shows the mesoscale total column water averaged over quartiles for the quasi-Lagrangian domains from the simulations with no evaporation of rain and gray lines for the original simulations (from Figure 6). As with the simulations with different resolutions and start times, the difference in trajectory location for the simulations with no evaporation of rain is small compared to the size of the domain (see Figure 5), so any differences seen are due to the simulation of the cloud development.
The simulations with cold pools suppressed show stronger development of mesoscale organization (Figure 11d). The overall difference in mesoscale organization, however, is small compared to the original simulations with the initial increase in mesoscale organization largely unaffected. The timing of the increase and decrease in mesoscale organization is also similar in all the simulations, indicating that the cold pools only modulate the mesoscale organization in this case, rather than being the main process stopping further development.
The simulations with no evaporation having stronger mesoscale organization may indicate that the cold pools do act to suppress convection in the moistest regions. Comparing the outgoing longwave flux after the flowers have developed (Figures 4 and 10) we see that the simulations without cold pools are associated with slightly larger flower structures, most noticeable at higher resolution, indicating that the cold pools do have a small effect on the development of the flowers. However, these differences could also be related to stopping the re-evaporation of rain above cloud base.
The outgoing longwave flux shows larger differences in the regions between the flowers where the simulations without cold pools have much less very shallow cloud indicating that these clouds are primarily generated by the interactions of cold pools. The original simulations have too much of this very shallow cloud compared to the satellite observations suggesting that the UM produces too many/too strong cold pools or too many clouds from the interactions of cold pools. We can say that the original UM simulations probably are producing too many cold pools because they produce too many cloud structures compared to satellite observations ( Figure 4) and most of these cloud structures are associated with precipitation and cold pools ( Figure 12).
The differences in moisture are not only in the moistest regions and the 4.4 km simulation actually has fairly similar total column water for the moistest quartile in the simulations with and without evaporation of rain. Another difference is that the other quartiles are drier in the simulations with no evaporation of rain at the time when the cold pools are initially developing in the reference simulation. This could be due to the cold pools acting to transport moisture from moister region to drier regions, which would make sense, although it could also just be a sign of weaker convergence of moisture to the moistest regions.

Domain Averages
In this section, we compare the domain averages across the different simulations, including the sensitivity tests from Sections 3.4 and 3.5. In the outgoing longwave flux (10), we saw large differences between simulations in the size of the flower structures produced and in the amount of very-shallow clouds. To quantify these differences, Figure 13 shows the cloud fraction in the 1.1 km simulation and the difference in cloud fraction, relative to the 1.1 km simulation, for the different resolution simulations and the sensitivity tests.
In the 1.1 km simulation, the clouds are initially very shallow, with most of the cloud around cloud base, and then deepen during the development of the flowers reaching maximum heights around 2.5-3 km. The cloud fraction in the quasi-Lagrangian domain shows a smoother increase in cloud-top height than the cloud fraction in the region of the inner domain (comparing Figures 13a and 13b) demonstrating the value of following the cloud development with the quasi-Lagrangian domain.
There is a strong resolution sensitivity in the cloud fraction, with lower resolution having increased cloud fraction. Although the hectometer-scale simulations nested in the 1.1 km simulation appear to have very similar cloud structures (Figure 2), the sensitivity to resolution is clearly seen in the cloud fraction in these simulations. This resolution sensitivity is consistent with previous studies using LES at kilometer-scale resolutions that have shown increased cloud fraction at lower resolution (Cheng et al., 2010;Stevens, Acquistapace, et al., 2020) and suggests that the resolution sensitivity could explain some of the difference in the strength of mesoscale organization between our simulations and those in Narenpitak et al. (2021).
Surprisingly, the simulations with no evaporation have slightly more cloud despite plots of longwave radiation showing less very shallow cloud in these simulations in the regions between the flowers (18Z in Figure 10b). This indicates that there are compensating increases in very shallow cloud in the region of the flowers. This makes sense because cold pools form shallow stable layers under the precipitating cloud, suppressing convection, but also generate new convection along their gust fronts away from the precipitating cloud. In the simulations without cold pools there is more cloud in the flowers and an absence of the very-shallow clouds generated by gust fronts away from the flowers The simulations without evaporation also stop evaporation of rain above cloud base which may also affect the amount of cloud. Either way, the difference in cloud fraction is fairly small compared to the value of cloud fraction and differences between other simulations. The main conclusion is that the cold pools have little effect on mesoscale organization and the development of flowers in this case.
In the previous sections, we showed that the initial organization and cold pools both have little effect on the timing of the development of mesoscale organization and flowers. This indicates that the large-scale circulation dictates the timing of development of the flowers by creating conditions favorable for the strengthening of the mesoscale circulation and moisture anomalies. To quantify the large-scale circulation, and the influence of the flowers on other large-scale parameters, Figure 14 shows large-scale averages of quantities over the region of the quasi-Lagrangian domains.
The total column water is similar across all the simulations, although there is a tendency for higher resolution to be drier, showing that the mean total column water is a poor predictor of the clouds and organization. The development of the flowers and increase in mesoscale organization is associated with the time when the large-scale total-column water is increasing. This is related to variations in the large-scale moisture flux convergence. Considering our quasi-Lagrangian domains, the net influx of moisture into these domains due to the large-scale dynamics will be given by the vertical integral of the horizontal moisture flux convergence, where the velocity, v, is the large-scale horizontal wind relative to the moving quasi-Lagrangian domain. Because the quasi-Lagrangian domains follow trajectories at a fixed height, there could be a substantial contribution to the influx of moisture just due to the movement of the domain, if the relative winds at other height levels transport moisture into the domain because of wind shear, as well as the contribution due to convergence around the domain. To disentangle these two contributions, the moisture flux convergence can be separated into a convergence term (−ρq∇ h ⋅ v) and an advection term (−v ⋅ ∇ h ρq). These two terms were calculated using large scale averages over a 3 × 3 grid with the quasi-Lagrangian domain at the center to calculate the gradients. Figure 14b shows the vertical-integral of the horizontal moisture flux convergence, and Figure 14c shows the contribution from the convergence term. The contribution of moisture flux convergence to the total moisture in the domain is generally negative, becoming slightly positive around 12Z in all simulations. The convergence term is typically larger than the moisture flux convergence and therefore the contribution from the advection term is generally negative. It is the surface fluxes that have the largest contribution to increasing moisture in the domain (Figure 9c). However, the contribution from surface fluxes does not vary much in time. The increase in total column water is more closely related to the increase in moisture flux convergence and the development of the flowers coincides with the time when the contribution from the convergence term is positive.
The increase in the convergence term in our simulations coincides with an increase in large-scale vertical velocity with both terms turning positive at a similar time (Figure 14d). This is consistent with Narenpitak et al. (2021) who showed the development of flowers was associated with positive large-scale vertical velocity and a weaker vertical velocity resulted in slower development of mesoscale organization; however, it is difficult to determine cause and effect from our simulations and the large-scale convergence could be a result of the development of the flowers.
The gray lines in Figure 14 show the quantities from ERA5 averaged over the 1.1 km quasi-Lagrangian domain for comparison. The estimates of the large-scale dynamics are largely consistent between the UM simulations and ERA5: the total column water shows a similar increase during the development of the flowers and is associated with an increase in moisture-flux convergence (stronger in ERA5) and positive large-scale vertical velocity. Interestingly, the simulations initialized a day later are less consistent with the ERA5 vertical velocity, indicating that the large-scale flow is still adjusting by the time the flowers develop and more time is needed to spin up these simulations. The vertical velocity turns positive earlier in the simulation initialized a day later but is also preceded by a period of stronger negative vertical velocities indicating that the convection is initially suppressed in these simulations, which will contribute to the weaker mesoscale organization. The differences in large-scale vertical velocity are consistent with the differences in the timing of the development of mesoscale organization between our simulations and those in Narenpitak et al. (2021). In both cases the mesoscale organization increases when vertical velocity is positive. In our simulations the development of mesoscale organization has stopped by the end of the day, even when using a trajectory and domain size consistent with Narenpitak et al. (2021), whereas, in Narenpitak et al. (2021, the mesoscale organization continues to develop at this time associated with positive vertical velocities in the ERA5-derived forcing, although the vertical velocity does turn negative below about 1.5 km by the end of the day. The difference could then be due to the UM vertical velocity not matching the ERA5 forcing because (a) the dynamics have diverged due to model differences, (b) the ERA5 forcing is calculated at a single point, or (c) ERA5 is different due to the assimilation of new data. Interestingly, there is an increase in the total column water in ERA5 toward the end of the day that is not captured in the UM simulations, and also not associated with changes in large-scale moisture flux convergence or vertical velocity, showing that differences are developing in the large-scale dynamics between our simulations and ERA5.
In contrast to total column water, the liquid water is strongly dependent on resolution (Figure 14e). The differences in liquid water make sense in relation to the differences in cloud fraction. Higher resolution is associated with less liquid water which makes sense since we have seen that higher resolution is also associated with less cloud. Despite the sensitivity of liquid water to resolution, the sensitivity of precipitation to resolution is less obvious (Figure 14f). While there can be large differences in total precipitation between simulations, there is no obvious link to resolution. This implies that the differences in liquid water are not linked to the clouds that form precipitation or there are compensating changes in the intensity of precipitation.
There are differences in liquid water and precipitation associated with start time with the simulations that are initialized later developing precipitation later and producing less precipitation overall which makes sense since these simulations start with less variability in moisture and therefore will not have the regions of larger moisture content that are then associated with more liquid water.
Unsurprisingly, the simulations without evaporation of rainfall have more precipitation since the precipitation all falls to the ground rather than forming cold pools. However, the simulations without evaporation of rainfall also show more liquid water after the initial development which also contributes to the increased rainfall. The increased liquid water makes sense because there are no cold pools to suppress further development of the clouds leading to these simulations having stronger mesoscale organization. However, it is unclear whether this is due to the cold pools or the evaporation of rain above cloud base in the UM which are both stopped by switching off evaporation.
As expected, the differences in clouds are associated with differences in radiation. The average outgoing longwave flux is most sensitive to the spin up ( Figure 14h). This makes sense because the simulations with the earlier start time have higher cloud tops and therefore a lower domain-average outgoing longwave flux. This sensitivity is also evident in the shortwave flux ( Figure 14i) where the simulations with the earlier start time, and increased cloud fraction, reflect more shortwave radiation.
We would expect to also see a strong sensitivity to resolution for the same reason: lower resolution has increased cloud fraction, so should have lower outgoing longwave flux and higher reflected shortwave flux. This is true for the reflected shortwave flux which largely follows the changes in cloud fraction and liquid water content (note the similarities between Figures 14e and 14i). The sensitivity of the longwave flux to resolution is more complex as there are other compensating changes. The longwave flux is sensitive to the depth of the boundary layer (shown in Figure 14g). Prior to the development of the flowers, the boundary layer gradually deepens in all simulations with lower resolution associated with a deeper boundary layer on average. Lower resolution, and a deeper boundary layer, is associated with lower outgoing longwave flux because the clouds are very shallow at this time, so the cloud tops are closely related to the boundary-layer depth. After the flowers develop, the average boundary layer depth decreases. This is largely due to the development of cold pools because the boundary layer is diagnosed as very shallow where there are cold pools, leading to a shallower average boundary-layer depth in simulations with more cold pools, although there is still a decrease in the boundary layer depth for the simulation with no evaporation of rainfall. At this time, the decrease in average boundary-layer depth is closely related to the development of precipitation which does not have a clear sensitivity to resolution. The average outgoing longwave flux also loses the sensitivity to resolution because the effect of the flowers for decreasing the longwave radiation, which is sensitive to resolution, is compensated by differences in the more shallow clouds associated with cold pools.

Conclusions
We have run high-resolution simulations of the 2nd February case study from EUREC 4 A-ATOMIC with the Met Office's UM. This case study is of interest because the cloud organization transitions from a regime of shallow, disorganized, cumulus clouds, known as "sugar," to a regime of deeper clouds with large detrainment layers, known as "flowers." The UM simulations reproduce the observed increase in mesoscale organization associated with the development of the flowers; however, the details of the clouds have some issues: the UM produces too much very shallow cloud and the size of the deeper shallow cloud structures are sensitive to resolution, with lower resolution producing larger cloud structures.
To better follow cloud development, we focused our analysis on subdomains extracted from our simulations following boundary-layer trajectories. These "quasi-Lagrangian" domains allow us to focus our analysis on the development of organization following the clouds. The main motivation behind using these quasi-Lagrangian domains was to compare our results with the Lagrangian LES results from Narenpitak et al. (2021) where the LES is driven by forcings following a trajectory. Higher resolution nested simulations provided little extra value because of the limited domain size making them too dependent on the kilometer-scale simulation which supplied the boundary conditions. This is because the domain size of the higher resolution nested simulations was too small to capture the development of clouds in the time it took for the flow to cross the domain which highlights how useful a Lagrangian approach is to modeling these clouds at high resolution.
Consistent with Narenpitak et al. (2021), we find that the development of the flowers is associated with an increase in mesoscale organization generated by the mesoscale vertical velocities and associated moisture transport. This process, first described by Bretherton and Blossey (2017), is where latent-heating driven mesoscale vertical velocities provide a positive feedback on convection by converging moisture toward the convection, making moist patches moister and dry patches drier. It is useful that we have shown that the kilometer-scale UM simulations can reproduce this process because that means these simulations can be used to better understand the sensitivity of mesoscale organization to changes at larger scales.
In the simulations of Bretherton and Blossey (2017), they found that the development of mesoscale organization stopped once the mesoscale anomalies became strong enough that the dis-aggregating effect of advection on the mesoscale anomalies balanced the aggregating effect of the circulation generated by the effect of mesoscale vertical velocities on the background humidity profile. In our simulations, the two processes do not reach a balance, instead the decrease in aggregation happens when the aggregating effect of the mesoscale vertical velocities strongly decreases.
A possible explanation for why the aggregation stops increasing is that cold pools generated from evaporation of rainfall underneath the flowers act to suppress convection in the moistest regions, stopping the circulation from converging more moisture to these regions. However, sensitivity studies with evaporation of rainfall switched off showed that the effect of cold pools on mesoscale organization is weak. The cold pools do still have important effects on radiation by generating shallower clouds in the regions between the flowers; however, these clouds are over represented in our simulations. Although cold pools only have a weak influence on the development of the flowers in our simulations, they may still be more important for other regimes of organization, such as "gravel," which has cold pools with stronger updrafts and deeper clouds, or "fish," which has the largest cold-pool fraction of the four named cloud patterns (Vogel et al., 2021). Cold pools may also be more important in other cases of flowers. In an LES simulation of the same case study, but with the diurnal cycle shifted 12 hr, Narenpitak et al. (2023) showed that cold pools acted to transport moisture away from the moistest regions but also generated new flowers on the gust fronts.
We found that our simulations differed from the LES of Narenpitak et al. (2021) in the timing and the magnitude of the development of mesoscale organization. In the simulations of Narenpitak et al. (2021), the mesoscale organization develops from an initially horizontally homogenous state (zero mesoscale organization) and continues to increase past the end of the day, whereas in our simulations, there is already mesoscale organization present and the development of the flowers is associated with an approximate doubling in the strength of the mesoscale organization. We also find that the development of mesoscale organization is more rapid in our simulations and starts to decrease by the end of the day. The differences between our simulations and the LES simulations in Narenpitak et al. (2021) can largely be explained by three factors: (a) lower model resolution leads to stronger mesoscale organization and increased cloud fraction, (b) longer spin up time results in stronger mesoscale organization and deeper clouds, and (c) the large-scale dynamics dictates the timing of the development of the flowers.
Other factors, such as the much larger domain size in our simulations, which has been shown to result in much stronger cloud organization (Narenpitak & Bretherton, 2019), as well as differences between the models used, may also be a factor in the differences but are not investigated here.
The kilometer-scale simulations exhibit strong sensitivities to resolution. Lower resolution is associated with stronger mesoscale organization as well as larger flower structures and increased cloud fraction. This is consistent with studies that run LES at kilometer scale resolutions finding increased cloud fraction at lower resolutions (Cheng et al., 2010;Stevens, Acquistapace, et al., 2020). An important aspect of the flowers, is the presence of thin stratiform cloud layers . (Kazil et al., 2017) showed that the aspect ratio of stratocumulus cloud cells can impact the entrainment rate. Since our simulations use the same vertical resolution, the grid aspect ratio is strongly tied to horizontal resolution and may also play a role in the sensitivity of the flowers to resolution.
To test whether the development of the flowers is strongly impacted by the degree of initial organization, a new set of simulations was initialized one day later (00Z 2nd February). Because the simulation is initialized from 25 km reanalysis fields, the coupled development of small-scale cumulus and mesoscale moisture perturbations will be weaker than in the Default simulation, which was initialized 1 day earlier. We found that the timing of the development of mesoscale organization and the flowers was unchanged despite the mesoscale organization always being weaker, and the cloud structures being smaller, in the simulations with a later start time.
The development of mesoscale organization in our simulations is associated with large-scale convergence leading to an increase in moisture. It is difficult to determine cause and effect from our simulations and the large-scale convergence could just be signature of the development of the flowers. This result is consistent with Narenpitak et al. (2021), who showed that the development of flowers was associated with the forcing of positive large-scale vertical velocity and a weaker forcing resulted in a slower development of mesoscale organization. This suggests that the large-scale circulation dictates the timing of development of the flowers by creating conditions favorable for the strengthening of the mesoscale circulation and moisture anomalies. The differences in timing of the development of the flowers between our simulations and those of Narenpitak et al. (2021) are consistent with differences in the times when the large-scale vertical velocity is positive. Although our simulations are initialized from ERA5, differences emerge by the end of 2nd February, making them inconsistent with the ERA5-derived forcings used by Narenpitak et al. (2021).
The main motivation in improving models representation of these clouds is their relation to climate sensitivity. In LES studies, it has been shown that lower-resolution simulations show more sensitivity to temperature increases because the base state has more cloud cover which is reduced with warming whereas high-resolution have less cloud cover resulting in a smaller reduction with warming (Blossey et al., 2009;Radtke et al., 2021). Our simulations show the same sensitivity to resolution where the lower resolutions produce more cloud; however, while they are producing too much cloud compared to observations, there are aspects of the clouds that actually look more realistic than at higher resolution, such as the larger size of the flowers and the stronger organization. This sensitivity is seen in shortwave radiation because the larger flowers reflect more shortwave radiation. However, the sensitivity is less obvious in longwave radiation, due to other compensating changes in low clouds. This presents a problem for kilometer-scale climate projections because the sensitivity of radiative fluxes to changes in mesoscale organization will still be uncertain. The poor representation of the cloud structures compared to observations would not be fixed by tuning the radiative fluxes because the changes in radiation, not just the absolute values, are sensitive to the model setup.
Given the large problems and sensitivities of these kilometer-scale models in producing trade-wind cumulus and the associated radiation, they cannot be considered as a solution for uncertainties in climate sensitivity yet. There is still more work needed on the representation of mesoscale organization and clouds in models at these scales, and to improve how the parametrizations and model dynamics interact across the gray zone (Honnert et al., 2020). The strong sensitivity of cloud fraction and organization to resolution highlights the need for scale-aware parametrization at these scales. Including a reduced and simplified parametrization of shallow convection can improve the representation of convection in the gray zone (Tomassini et al., 2022). However, given the models reproduce the variations in mesoscale organization and associated processes, they are a useful tool for understanding processes driving mesoscale organization and interactions with larger scales. Figure A1. The same as Figure 2 but with the model output shown as simulated 11 μm brightness temperature and the output is shown for a different set of times (top header). This output comes from the same model setup but run on a different machine, so may not be exactly comparable.

Data Availability Statement
The GOES satellite data was downloaded from French national center for Atmospheric data and services, AERIS, at https://observations.ipsl.fr/aeris/eurec4a-data/SATELLITES/GOES-E/2km_10min/ and ERA5 data was downloaded from the climate data store (Hersbach et al. (2018b) and Hersbach et al. (2018a)). The regridded data, trajectories, simulated satellite data, and remaining files required to reproduce the results in this study are available in a Zenodo repository (Saffin and Lock (2023), https://doi.org/10.5281/zenodo.7568386). The code used for data analysis is available at https://github.com/leosaffin/moisture_tracers. The UM code is freely available to anyone for non-commercial use from the Met Office Science Repository Service (https://code.metoffice.gov. uk/) and all simulation data used in this study are also archived at the Met Office and are available for research purposes through the JASMIN platform (www.jasmin.ac.uk). For details, please contact UM_collaboration@ metoffice.gov.uk referencing this paper.