Pressure Drag for Shallow Cumulus Clouds: From Thermals to the Cloud Ensemble

This study takes the first step to bridge the gap between the pressure drag of a shallow cloud ensemble and that of an individual cloud composed of rising thermals. It is found that the pressure drag for a cloud ensemble is primarily controlled by the dynamical component. The dominance of dynamical pressure drag and its increased magnitude with height are independent of cloud lifetime and are common features of individual clouds except that the total drag of a single cloud over life cycle presents vertical oscillations. These oscillations are associated with successive rising thermals but are further complicated by the evaporation‐driven downdrafts outside the cloud. The horizontal vorticity associated with the vortical structure is amplified as the thermals rise to higher altitudes due to continuous baroclinic vorticity generation. This leads to the increased magnitude of local minima of dynamical pressure perturbation with height and consequently to increased dynamical pressure drag.


Introduction
Reasonable estimation of the vertical velocity within parameterized shallow cumulus clouds in climate models is critical for representing the atmospheric convective mixing that is responsible for half of the spread in model climate sensitivity (Sherwood et al., 2014). A steady-state vertical momentum equation (Simpson & Wiggert, 1969) is vertically integrated to diagnose the in-cloud vertical velocity and has been widely used in many convection parameterizations (Bretherton et al., 2004;Cheinet, 2003;de Rooy & Siebesma, 2010;Gregory, 2001;Kim & Kang, 2011;Neggers et al., 2009;Pergaud et al., 2009;Rio & Hourdin, 2008;Siebesma & Teixeira, 2000;Siebesma et al., 2007;Soares et al., 2004;Sušelj et al., 2012). It is usually expressed in a general form as 1 2 force, not entrainment, that dominates the balance with the buoyancy acceleration within an ensemble of shallow clouds. Further supporting evidence was found in deep convection by tracking finer scale thermals within clouds (Sherwood et al., 2013) or near cloud top (Romps & Charn, 2015;Romps & Oktem, 2015).
It is then natural to ask how the pressure drag can be represented in a parameterization. Theoretical studies (Jeevanjee & Romps, 2016;Morrison, 2016aMorrison, , 2016bPaulius & Garner, 2006;Weisman & Klemp, 1997) indicate that pressure perturbation effects on vertical velocity depend upon the aspect ratio H/R, where R is the updraft radius and H the updraft height. These effects can be incorporated into the vertical momentum equation through a rescaling of the buoyancy term using a virtual mass coefficient as shown by Morrison (2016b). However, their derivation only deals with the updraft center and neglects the effects from the dynamical part of the pressure perturbation (see section 4.1). Meanwhile, for parameterization purposes, it is the mean vertical velocity of a cloud that is often needed, rather than the velocity at updraft center. Furthermore, a recent study (Gu et al., 2020) found that a partitioning into the mean velocities of the "core" and "cloak," namely, the strong and weak parts of the updrafts respectively, can significantly improve the representation of vertical transport in a simple bulk transport equation. This indicates the importance of predicting vertical velocities in different parts of clouds and thus the need for a better understanding of pressure drag, not only in the cloud center but also off the central axis.
The pressure drag has been investigated from the microscopic perspective of a single updraft or thermal (Morrison, 2016a;Romps & Charn, 2015;Sherwood et al., 2013) and from the macroscopic perspective of a cloud ensemble (de Roode et al., 2012). But little is known about how the pressure drag of a cloud ensemble, which is the focus of most convection schemes, can be related to that of a single cloud or successive rising thermals. This study aims to take the first step toward bridging the gap by relating the behaviors of the cloud ensemble with those of individual thermals within clouds. We use a large eddy simulation of shallow cumulus clouds and a cloud tracking algorithm to distinguish the behavior of individual clouds within the ensemble (section 2). We construct budget equations for the vertical velocity of each cloud and for the ensemble (section 3) in an attempt to (1) quantify the contributions from dynamical and thermodynamic components to the pressure drag (section 4.1), (2) demonstrate the relationships between the pressure drag for the cloud ensemble and that for individual clouds with rising thermals embedded (sections 4.2 and 4.3), and (3) understand the physical mechanism for the vertical distribution of pressure drag (section 5).

Large Eddy Simulations
This study focuses on a large eddy simulation of shallow cumulus clouds using the Met Office-Natural Environment Research Council (NERC) Cloud model (MONC; Brown et al., 2015;Brown et al., 2018), based on the Barbados Oceanographic and Meteorological Experiment (BOMEX). The grid spacing is 25 m in all directions, and the domain size is 15 km 2 × 3 km. The Smagorinsky-Lilly scheme (Lilly, 1962;Smagorinsky, 1963) is used for the parameterization of subgrid turbulence. Further details of the model configuration can be found in Gu et al. (2020). Our analyses cover a period of equilibrium from Hour 5 to Hour 6 of the simulation, with 1 min output frequency.

Cloud Tracking Algorithm
To understand the pressure effect experienced by an individual cloud over its life cycle, we need to be able to track each cloud object. The first step is to label cloudy points for which we use the criteria that cloud liquid water q l > 10 −5 kg kg −1 . Contiguous labeled points are combined to form an individual cloud object by checking the neighboring grid points in six directions, corresponding to the six faces of the grid cell around the cloudy points, until no more cloudy points are found (Gu et al., 2020).
After the cloud objects at all output times have been identified, a matching up of cloud objects at adjacent output times is conducted. The algorithm is an extension of Muetzelfeldt (2020), which itself combines methods from Stein et al. (2014) and Plant (2009), allowing it to track 3-D cloud objects instead of 2-D cloud objects. We calculate a translation vector between one output time and the next, which produces the maximum correlation of the cloud fields. The translation vector is obtained using the convolution of Fourier transformed 2-D fields at the two adjacent output times, the fields themselves being the vertical projection of all 3-D cloud objects onto the horizontal plane. A particular cloud object is then considered to be the same cloud object from one output time to the next if, when projected forwards to the next output time, it overlaps with or is touching a cloud object at the next output time. There are some complications when building up the life history, depending on how the projected cloud objects overlap with those at next output time, which can result in merging and splitting of cloud objects (see Muetzelfeldt, 2020;Plant, 2009 for details). For convenience of budget analysis over the life cycle of each cloud object, here we retain only one cloud object if there are multiple objects identified as overlapping or touching. If one cloud merges with another cloud, we continue to track both clouds as one merged cloud until it dissipates. If a cloud splits into several clouds, we only continue to track one cloud object. The retained cloud object has the closest cloud depth to the parent cloud at previous output time. In this way, each cloud object can be appropriately tracked over its lifetime without losing much information (supporting information Figure S1). In our study, a total of 4,448 cloud objects have been tracked, and subsampling tests confirm that this is sufficient to provide robust statistical results.

Constructing the Vertical Velocity Budget Equation of the Cloud Ensemble From That of Individual Clouds
To build a bridge between an individual cloud and the cloud ensemble, we construct the vertical velocity equation for the cloud ensemble from that of individual clouds. At a given time, we assume that multiple cloud objects are located within the domain, each having a vertically coherent structure. Each cloud object is denoted with a subscript i > 0, and, for convenience, the environment is considered as an extra object denoted by i = 0. An atmospheric quantity within the object i is denoted i , the spatial average of this quantity over the object is i , and the perturbation from the average over the object is ′ i = i − i . The area fraction of each object is denoted by a i . Note that a i should be height dependent as a 3-D cloud object changes area coverage with height. Following Tan et al. (2018), conditionally averaging the continuity equation and the vertical momentum equation over the volume occupied by cloud object i and combining the two equations lead to a budget equation for the averaged vertical velocity of cloud object i where is air density, i is the fractional entrainment rate, B i represents the averaged buoyancy source that is calculated using perturbations from the domain-averaged density, ( p nh ∕ z) i represents the averaged non-hydrostatic pressure gradient force, and S i represents other averaged source terms (e.g., subgrid tendency, Coriolis force, or damping tendency near the top of domain).
The budget equation for the vertical velocity of the cloud ensemble can be obtained by first summing the conditionally averaged continuity and vertical momentum equations over all cloud objects and then combining together (see Text S1 for more details), to yield where ∑ i denotes a summation over all i > 0 and the total area fraction of the cloud ensemble is a = ∑ i a i . w c is the average vertical velocity over all cloud objects and is defined as w c = ∑ i a i w i ∕a. E i and D i denote rates of lateral entrainment and detrainment, respectively, and may be expressed in terms of fractional rates i and i , according to E i = i a i w i and D i = i a i w i . The entrainment and detrainment terms in Equation 3 will be determined together as the residual from the budget equation since all of the other terms can be explicitly calculated from the model output. For the purposes of this study, we will focus on the two major balanced terms: the buoyancy source and the pressure effect term. We intend to provide further analysis of entrainment/detrainment in a future study.  (b) only continuously tracked clouds from Hour 5 to Hour 6. "PGF," "BUOY," "VIS," "ADV," "SUB," and "ENT" represent, respectively, the pressure gradient force (blue), buoyancy source (red), subgrid tendency (green), advection (yellow), subplume transport (gray), and entrainment terms (magenta) in equation (3). The damping tendency applied near the top of the domain to prevent gravity waves reflecting is zero within the cloud layer and thus is not shown. In (c), the buoyancy source (red) and pressure gradient force (blue solid) terms are again shown (as in panel b), and the latter is decomposed into its components: the dynamical pressure gradient force ("PDGF," blue dotted) and thermodynamic pressure gradient force ("PBGF," blue dash-dotted). (d) is the same as (c) but for the buoyancy source and pressure drag of the updraft cores (defined as the cloudy points with top 0.5% upward motion). The thick black line in each plot denotes the zero tendency. z b is the depth of the well-mixed subcloud layer (≈500 m).

Pressure Drag for the Cloud Ensemble and an Individual Cloud
In section 4.1, we evaluate the vertical velocity budget for the simulated cloud ensemble. We then investigate the relationship between the pressure drag of the ensemble and the dynamics of individual clouds. We do not consider in detail the dynamics of individual thermals as this has been investigated in previous studies (Hernandez-Deckers & Sherwood, 2016Moser & Lasher-Trapp, 2017;Sherwood et al., 2013). However, the role of thermals is reflected in the budget for individual clouds and is discussed in section 5. Figure 1b shows the vertical velocity budget for all tracked cloud objects. The overall picture of the role of each term is consistent with the budget for the full cloudy region across the domain (Figure 1a). In most of the cloud layer, buoyancy is the major source term and is offset mainly by the pressure drag. Within the inversion layer (1,500-2,000 m), however, the buoyancy changes sign and the entrainment and subplume terms also become important. The pressure drag and buoyancy terms in the budget of tracked clouds are larger than those for the budget of the full cloudy region. This occurs because, during splits, the tracking keeps the cloud object whose cloud depth is closest to the depth at previous time. The neglected objects are typically less buoyant than those retained. Nonetheless, the vertical distributions of pressure drag and buoyancy in the two budgets agree pretty well.

The Cloud Ensemble
The buoyancy source increases from cloud base to a maximum slightly below 1 km above which the buoyancy decreases slowly until the inversion layer, where it acts to decelerate the flow. The pressure drag increases slowly from the cloud base upward to the inversion layer within which it decreases slightly. Therefore, it may not be appropriate to absorb the pressure drag into a reduced buoyancy term (de Roode et al., 2012;Morrison, 2016b), since these two terms do not change synchronously with height, especially in the inversion layer.
To understand the change of pressure drag with height, the total pressure perturbation is decomposed into dynamical (p D ) and thermodynamic (p B ) components as follows: where u is the wind vector. Note that "≈" in the first equation of Equation 4 is because we neglect a small Coriolis contribution. The pressure effect can then also be decomposed into its dynamical and thermodynamic components as shown in Figure 1c. The dynamical part accounts for most of the pressure drag and also dominates the vertical variations below the inversion layer. The thermodynamic part has a small increase in magnitude from cloud base to 650 m and is relatively steady with height below the inversion layer compared to the dynamical part. Within the inversion layer, the dynamical drag continues to increase with height but is largely compensated by the increasing thermodynamic drag produced by the rapid decrease in cloud buoyancy. Therefore, it may be reasonable to rescale the thermodynamic part of the pressure drag with the buoyancy source, but this does not work for the dynamical part.
The dominance of dynamical pressure drag in the vertical velocity budget of the clouds is not in contradiction with the assumption of theoretical studies (Morrison, 2016a(Morrison, , 2016b in which the thermodynamic pressure perturbation is considered to be the more important component along the updraft central axis. A budget analysis for the cloudy updraft cores (the top 0.5% of upward motions in the cloudy air) confirms that the dynamical pressure drag is smaller than its thermodynamic counterpart in such cores below the inversion layer (Figure 1d), indicating that their relative roles are largely reversed compared to the clouds as a whole. This is because the rotational and deformation contributions to ∇ 2 p D (section 5) are much smaller than the contribution of buoyancy forcing to ∇ 2 p B in an axisymmetric updraft core. The positive tendency from the dynamical pressure gradient force near cloud base is due to convergence and therefore supports upward acceleration. Only the thermodynamic pressure gradient force serves as a drag on clouds there.

Dependence on Cloud Lifetime
We wish to understand whether the increased pressure drag with height and its dominance by the dynamical pressure component, described above for the cloud ensemble, is a common feature of most clouds within the ensemble or whether it might arise from some smaller subset of clouds with particular sizes or lifetimes. In this subsection we consider the vertical velocity budget for clouds with different lifetimes. The budgets were also constructed for clouds with different maximum volumes during their life cycles, but such results are not presented here. They lead to similar conclusions, in all likelihood because the cloud maximum volume is closely related to cloud lifetime.  (Figures 2b and 2c). The contribution from the long-lived clouds (Figure 2d) is relatively small (about 20%) due to their small population but nonetheless exceeds that from the many short-lived clouds of less than 5 min. For each of the cloud categories with lifetimes longer than 5 min, the pressure drag always increases with height in the cloud layer, and its magnitude and vertical distribution are dominated by the dynamical component, as for the full ensemble. One may question whether some of the tracked clouds may have been miscategorized, with only a part of the life cycle having been captured due to the finite analysis period considered in this study. For example, if a cloud enters its decaying stage at the start of our time window for tracking, its lifetime will have been underestimated. Additional analyses (not shown) were performed with only those clouds whose entire life cycles fell within the tracking period, and the conclusions above were unchanged.

Individual Cloud Over Its Life Cycle
Given that the main features in the vertical structure of pressure drag do not depend on the lifetime of clouds, it is plausible to hypothesize that such features may hold at the level of an individual cloud. Figure 3 shows the pressure gradient force and buoyancy source for three example clouds, averaged over their own life cycles. The three clouds have different lifetimes and are chosen to be representative of various clouds examined. It is clear that the total pressure gradient force is dominated by its dynamical component within the cloud layer (below 1,800 m) and that both the total and the dynamical component have an overall tendency to increase with height in that layer. However, they also present some marked local oscillations in the vertical, in contrast to the smoother profile of the buoyancy source term. For some clouds (e.g., Figure 3a), the amplitude of oscillation is sufficiently large for the pressure gradient force to locally accelerate the upward motion.
To understand the oscillations of pressure effect, we pick a single cloud (Cloud 1, Figure 3a) to investigate the detailed structures along a west-east vertical cross section passing through the cloud liquid water centroid at 1,000 m. Figure 4 shows the distributions of buoyancy (Figure 4a), vertical velocity (Figure 4b), pressure perturbation (Figure 4c), the pressure gradient force (Figure 4d), and the dynamical pressure perturbation (Figure 4f) near the middle of its life cycle (7 min). The dynamical pressure gradient force has a similar distribution to the total pressure gradient force and thus is not shown. There are multiple local maxima of vertical velocity and buoyancy in the vertical (800, 1,600, and 1,800 m), showing that the cloud is composed of multiple rising thermals (Figures 4a and 4b). Each thermal in B and w can usually be identified with two flanking local minima in the pressure perturbation field (Figures 4c and 4f). There is weak westerly vertical shear above 1 km height in the BOMEX simulation. The downdrafts associated with the cloud top overturning circulation occur preferentially on the downshear side of the cloud (Figure 4b) and are enhanced through evaporative cooling along the cloud edge (Figure 4a). Due to the enhancement of downdrafts on the downshear side, the local minima of pressure perturbations are found to be lower but stronger than their counterparts on the upshear side (Figure 4c).
The result of these pressure perturbations is that the pressure gradient force presents couplets of negative and positive tendencies off the central axis, with a region of negative tendency above each positive tendency region (Figure 4d). Therefore, the dominant dynamical pressure drag arises due to the spatial distribution of local extremes of pressure perturbation, which is, in turn, a response to the flow structures of the constituent thermals. Multiple rising thermals stacked in the vertical lead to multiple local minima of pressure perturbation in a cloud and hence to a vertical oscillation of the pressure drag which can be further complicated by the lateral downdrafts.

Increased Pressure Drag With Height
An explanation is still required as to why the pressure drag increases with height. Our physical interpretation focuses on the increasing magnitude of the local minima in dynamical pressure perturbation since this is the direct reason for the increased pressure gradient force with height, as shown in Figures 1-3. We consider an alternative form of the dynamical pressure specification from Equation 4. This can be rewritten as (Markowski & Richardson, 2011) where e ij is the deformation tensor and is the vorticity, such that where u 1 = u, u 2 = v, u 3 = w, x 1 = x, x 2 = y, and x 3 = z.
For a well-behaved field, the dynamical pressure perturbation is a spatially smoothed form of the source terms on the right-hand side of equation (5) with p D ∝ −∇ 2 p D . Thus, a minimum in the dynamical pressure 10.1029/2020GL090460 perturbation must be associated with rotation rather than deformation, since only rotation implies the negative dynamical pressure perturbation. For an understanding of the rotation itself, we turn to the vorticity equation: where k is the unit vector in the vertical. The Boussinesq approximation has been made, and viscous effects are neglected. Production of vorticity is due to tilting/stretching and to horizontal gradients of buoyancy. Following a rising thermal, if it initially does not have vorticity, then the vorticity will only be produced through the baroclinic generation term due to the horizontal gradient of buoyancy. Note that the buoyancy always has a local maximum near the central axis from cloud base to cloud top, despite the negative buoyancy values within the inversion layer (Figure 4a). This point is further supported by a composite analysis of the normalized distribution of buoyancy within the clouds at different vertical levels ( Figure S2). Due to this distribution of buoyancy, horizontal vorticity will be continuously generated as the thermals rise through the cloud as indicated by the second term in equation (8). The enhanced horizontal vorticity couplets along with the rising thermals result in the enhanced magnitude of the dynamical pressure perturbations with height on the flanks of the updraft center. The spatial distribution of − 1 2 | | 2 (Figure 4e) successfully captures the three couplets of local minima of dynamical pressure perturbation off the central axis (800, 1,600, and 1,800 m), the increased magnitude with height, and even a lower local minima on the downshear side of each pair (Figure 4f). Note also that the dynamical pressure perturbation has a local positive maximum between 1,100 and 1,300 m. This occurs because of the deformation due to the convergence near the bottom the thermal (Figure 4f). As a result, the dynamical pressure gradient force, and hence the total pressure drag, should generally increase with height for an individual cloud over its life cycle despite the oscillations associated with rising thermals.

Discussion and Summary
Our study demonstrates that the thermodynamic component of pressure drag in the vertical momentum budget of a cloud ensemble has a similar vertical distribution to the buoyancy source but with the opposite tendency ( Figure 1c). Thus, it can reasonably be parameterized through a virtual mass coefficient to reduce the buoyancy source as in Morrison, (2016aMorrison, ( , 2016b. However, the thermodynamic pressure drag is overwhelmed by the dynamical pressure drag in the budget for the cloud ensemble. This is not in contradiction with theoretical studies such as Morrison, (2016aMorrison, ( , 2016b) that found the dynamical pressure perturbation effect to be negligible along the central axis of an updraft because we are focusing here on the whole cloud ensemble rather than the cloud center. Further analysis shows that the thermodynamic pressure drag does dominate over its dynamic counterpart within the updraft core in most of the cloud layer ( Figure 1d). This demonstrates that the dynamical and thermodynamic components of pressure drag have different impacts on different parts of shallow cumulus clouds. In the cloud cores with the strongest upward motion, the thermodynamic pressure perturbation is responsible for most of the drag while the dynamical pressure drag dominates the region outside the core with weaker updrafts.
The fact that the dynamical pressure drag increases with height means that it can neither be absorbed as a reduced buoyancy source nor scaled with the entrainment term (or the square of vertical velocity). It must therefore be parameterized with a different approach. Our findings support a treatment for pressure drag that should take as its starting point for multiple drafts with different strengths, as in the "core-cloak" conceptual model of Gu et al. (2016a) for example. A recent theoretical study (Yano, 2020) also suggests that the inhomogeneity of vertical velocity must be introduced in the plume model; otherwise, the buoyancy and entrainment effect will perfectly cancel out with the counterbalancing force from the pressure, leaving a pure drag force to prevent a steady state solution.
We have linked the behavior of pressure drag for a cloud ensemble to that of the individual thermals within each cloud. The dominance of dynamical pressure drag and its increasing magnitude with height hold not only for the cloud ensemble but also for most individual clouds over their life cycles, albeit with individual clouds exhibiting noticeable vertical oscillations. Detailed investigation indicates that it is the successive rising thermals within the clouds that give rise to the oscillations. In addition to the thermals, the downdrafts outside the cloud can also enhance or distort local minima of the dynamical pressure perturbation, further complicating the picture of dynamical pressure drag. Nonetheless, we are able to explain the increase in the overall pressure drag with height. Enhancement of horizontal vorticity during ascent occurs through baroclinic generation that is aided by the internal distribution of buoyancy within the cloud, and this vorticity, in turn, provides the dominant source term for the dynamical pressure and thus the pressure drag.

Data Availability Statement
The data set and scripts used for this study are publicly available from the University of Reading Research Data Archive at DOI (https://doi.org/10. 17864/1947.259).