Momentum Transport in Shallow Cumulus Clouds and Its Parameterization by Higher‐Order Closure

It is challenging to parameterize subgrid vertical momentum fluxes in marine shallow cumulus layers that contain a jet in the profile of horizontal wind. In a large‐eddy simulation of such a layer, it is found that the momentum flux in the direction of strongest wind magnitude has a three‐layer structure. The lowest layer, from the ocean surface up to the jet maximum, has downgradient momentum flux. The middle layer, from the jet maximum up to an altitude several hundred meters above, has upgradient (i.e., countergradient) momentum flux because of transport of low‐magnitude momentum upward through the jet maximum. In the upper layer, the layer‐average momentum flux is weak. The budget of momentum flux shows that in the middle and upper layers, both the buoyancy production term and turbulent advection (i.e., third‐order flux‐of‐flux) terms are important.


Introduction
The parameterization of momentum transport in cumulus layers has an important impact on simulations of the general circulation of the atmosphere. In the tropics, an inadequate treatment of convective momentum transport may cause significant errors in the strength of the Hadley circulation and hence surface winds (Richter & Rasch, 2008;Wu et al., 2007). Such errors may be especially pernicious in coupled ocean-atmosphere simulations, in which surface winds apply stress on the ocean and the ocean response feeds back onto the atmosphere (e.g., Strobach et al., 2019). For instance, errors in surface wind stress may be partly responsible for errors in the equatorial cold tongue in the ocean (Li & Xie, 2014;Woelfle, 2018).
Much prior research has been done on momentum transport in deep convection. For instance, it has been found that when deep convection is organized into a squall line, the across-line momentum flux may sometimes be upgradient (e.g., Gallus & Johnson, 1992;LeMone, 1983;Trier et al., 1998). (This paper will define an upgradient momentum flux as one that has the same sign as the vertical gradient of the mean wind.) Less research has been done on momentum transport in shallow cumuli. Better understanding would be desirable because examination of reanalysis data sets hints that an important contribution to the momentum budget is made by shallow convective layers (Carr & Bretherton, 2001).
The main goal of this paper is to describe a new higher-order closure parameterization of momentum fluxes. However, in order to motivate our parameterization choices, the paper will first comment on the physics of momentum fluxes. In the large-eddy simulation (LES) of a marine shallow cumulus layer that we present 1. Lower layer: From the ocean surface up to the maximum of the jet near cloud base, the momentum flux is strong and downgradient. 2. Middle layer: Immediately above the jet maximum, there is a several-hundred-meter-thick layer in which the momentum flux is upgradient. The flux is upgradient even though the cloud field is not organized. 3. Upper layer: From the top of the upgradient layer at midcloud to the top of the cloud layer, the momentum flux is weak. Hence, with respect to the large-scale wind field, the clouds are nearly "invisible," even though those clouds contain strong turbulence and strong fluxes of heat and moisture.
Parameterizing these three layers poses difficulties. A low-order closure model that parameterizes momentum flux based on downgradient eddy diffusion cannot represent an upgradient flux. The same eddy diffusion model will be challenged by the fact that the "diffusion" of scalars in the upper part of the cumulus layer is large but the "diffusion" of momentum is much weaker. Mass-flux schemes confront a different challenge: They typically neglect vertical transport of momentum within the environment of cumulus clouds, even though the environmental flux may be significant (Brown, 1999;Schlemmer et al., 2017;Zhu, 2015).
(This paper will define "environment" to mean the ambient, clear air that laterally surrounds clouds.) Faced with these difficulties, we turn to a different parameterization method, namely, higher-order closure modeling (e.g., Mellor & Yamada, 1982). In particular, we will approximate the solution of a prognostic budget equation for momentum flux. In principle, higher-order closure allows the possibility of producing an upgradient flux, and its framework includes all relevant processes, even those in the environment. However, parameterizing the terms in the budget is challenging.
In section 2, we will describe the large-eddy and single-column models used and the configuration of cloud cases. Sections 3 and 4 will analyze the physics of the middle layer of upgradient flux and the upper layer of weak momentum flux, respectively. Section 5 will introduce our higher-order closure parameterization. Section 6 will compare parameterized solutions versus LESs. Section 7 will show that the effects of prognosing momentum fluxes on clouds and thermodynamic fields are small. We will make concluding remarks in section 8.

Description of Numerical Models Used and Cloud Cases Simulated
The large-eddy model that we will use in this study is the System for Atmospheric Modeling (SAM) model (Khairoutdinov & Randall, 2003). The single-column model that we will use is based on the Cloud Layer Unified By Binormals (CLUBB) parameterization .
These two models will be used to simulate five cloud cases. This paper will focus mainly on the the Barbados Oceanographic and Meteorological Experiment (BOMEX) shallow cumulus case (Siebesma et al., 2003), which has upgradient momentum fluxes, but in order to evaluate the behavior of prognosed momentum fluxes in different cloud regimes, the paper will also show some plots from four other cases. One of them is another marine shallow cumulus cases with upgradient momentum fluxes: the Rain in Cumulus over the Ocean (RICO) case (van Zanten et al., 2011). A third case is a continental shallow cumulus case located over the Atmospheric Radiation Measurement (ARM) site in Oklahoma (Brown & Co-authors, 2002). The ARM case allows us to assess the impact of an underlying surface consisting of land rather than ocean. In order to assess the effects of prognosing momentum fluxes on noncumulus clouds, we analyze two marine stratocumulus cases based on the second DYnamics and Chemistry Of Marine Stratocumulus field experiment, Research Flight 2 (DYCOMS-II RF02; Ackerman et al., 2009) and Research Flight 1 (DYCOMS-II RF01; Stevens et al., 2005). The RF02 cloud deck produces drizzle, whereas the RF01 deck does not.

LES Model and Single-Column Model
LESs will be used in order to study processes related to momentum transport and to provide reference simulations with which to compare single-column, parameterized simulations. The LES model used in this study is Version 6.10.6 of SAM. SAM has been widely used in the past to simulate boundary-layer and deep convective cloud cases, including the five cases simulated here. SAM solves a moist, anelastic system of equations using a finite-difference discretization and a third-order Adams-Bashforth time stepper. Subgrid-scale turbulence is parameterized by the Smagorinsky closure. Sponge damping begins at 70% of the distance from the model base and gradually increases in strength to the model top. The single-column model used to parameterize these same cases is built around CLUBB Larson & Golaz, 2005). CLUBB prognoses subgrid, second-order variances and fluxes of the eastward, northward, and upward components of velocity (u, v, and w), plus liquid water potential temperature ( l ) and total water mixing ratio (vapor plus cloud liquid, r t ). In addition, one third-order moment is prognosed, that of vertical velocity (w ′ 3 ). (In this paper, an overbar represents a horizontal, spatial average, and a prime represents a deviation from that average.) In these equations, the turbulent advection and buoyancy terms are closed by assuming that the subgrid distribution is a double Gaussian, that is, a mixture of normals. The pressure and dissipation terms are closed by use of classical turbulence parameterizations. CLUBB's formulation is described in detail in Larson (2017). For the simulations in this manuscript, the host model that drives CLUBB is the CLUBB single-column model ("CLUBB-SCM"), which is independent of any global atmospheric model and includes CLUBB plus various microphysics and radiative transfer options.
Radiative transfer is calculated in both SAM and CLUBB-SCM by use of an idealized, gray model Stevens et al., 2005). The microphysics option used in both SAM and CLUBB-SCM for the DYCOMS-II RF02 simulation is the Morrison microphysics (Morrison et al., 2009) with a prescribed cloud droplet number concentration of 55 cm −3 . In RICO, the Khairoutdinov-Kogan microphysics (Khairoutdinov & Kogan, 2000) was used for both SAM and CLUBB-SCM with a droplet number of 70 cm −3 . In the other cases, precipitation was shut off in both SAM and CLUBB-SCM. The plots presented below will display both cloud fraction and cloud-averaged quantities. In SAM's output, cloud is defined as the region where liquid water mixing ratio has a value greater than 0. (None of the cases we simulate contain ice.)

Configuration of Cases
The configuration of these five cases, including the calculation of lower-boundary surface fluxes, closely follows the configurations prescribed in the aforementioned references. The surface fluxes of heat and moisture are prescribed in all cases except RICO, which diagnoses those fluxes based on the surface values of wind speed, temperature, and humidity (van Zanten et al., 2011). The surface friction velocity is prescribed in all cases except RICO and ARM, which diagnose it based on mean wind. BOMEX, ARM, and RF01 turn off precipitation; RICO and RF02 are allowed to produce drizzle. In all five cases, there is no cloud droplet sedimentation in either SAM or CLUBB-SCM. In BOMEX, the Coriolis parameter is set to 3.76 × 10 −5 s −1 , and in DYCOMS-II RF02, it is set to 7.6 × 10 −5 s −1 . The configurations are summarized in Table 1.

Upgradient Momentum Fluxes in Shallow Cumulus Cloud Layers: Phenomenology and Budgets
This section discusses the physics of upgradient momentum fluxes in shallow cumulus layers and the budget terms associated with upgradient fluxes. The purpose of this discussion is to provide context and explain what motivates our later parameterization choices.

Layer of Upgradient Momentum Flux in the BOMEX Shallow Cumulus Case
The BOMEX LES exhibits a jet in the horizontal wind, with a maximum westward wind speed at an altitude of about 770 m above sea level ( Figure 1). This altitude is above cloud base (500 m) and well below cloud top (2,000 m). Below the jet maximum, the wind speed is reduced by friction associated with the turbulent subcloud layer.
Below the jet maximum, the eastward momentum flux, u ′ w ′ , is positive, and the vertical wind gradient, u∕ z, is negative (see Figure 1). Because u ′ w ′ and u∕ z have opposite sign, the momentum flux in this lower layer is downgradient. This corresponds to a positive eddy diffusivity, K, if we assume that A standard argument to rationalize downgradient fluxes heuristically, tailored to this particular wind profile, is to suppose that on average, a parcel that rises a small distance above the lower surface (and hence has w ′ > 0) carries with it low-magnitude surface momentum. Because u∕ z < 0, the parcel finds itself moving more slowly than environmental air (and hence has u ′ > 0; see Figure 1). Thus, u ′ w ′ > 0, and the flux is downgradient.
In contrast, u ′ w ′ is upgradient in a middle layer extending from the jet maximum at 770 m up to 1,070 m. In the middle layer, u ′ w ′ remains positive, but u∕ z becomes positive too. (We define the top of this layer, somewhat arbitrarily, as being the altitude where the within-cloud momentum flux, u ′ w ′ cld , first becomes negative. This is slightly above the altitude at which u cld = u env at 1,020 m. In this paper, () cld and () env denote averages over cloud and environmental air, respectively.) What is the cause of the upgradient flux above the jet maximum? The existence of an upgradient flux in this case cannot depend essentially on the organization of the cloud field into, say, lines or clusters, because our LES of the BOMEX (nonprecipitating) shallow cumulus case does not manifest any discernible form In this simulation, cloud base is at about 500 m, and cloud top is at 2,000 m. Aspects of the profile of u ′ w ′ can be rationalized by the facts that the layered-averaged (i.e., mean) wind (black line) exhibits a jet and that the within-cloud wind profile (orange dots) has a similar shape but is elevated. For example, in the middle layer, the momentum flux is upgradient because cloudy updrafts carry weaker-magnitude momentum up past the stronger jet maximum.
of organization (Domke, 2019). Instead, we rationalize the upgradient sign of the flux by inspection of the within-cloud eastward wind, u cld , shown in Figure 1. As a buoyant parcel rises from the lower surface, it carries low-magnitude momentum with it. The within-parcel momentum is driven by an across-cloud pressure gradient toward the environmental value, but only slowly. Hence, as the parcel rises above the jet maximum at ∼770 m, it still carries low-magnitude momentum (u ′ > 0), and it carries this momentum upward (w ′ > 0) into a layer of air with u∕ z of the "wrong" (positive) sign. In this layer, the momentum flux is upgradient. Almost as soon as the momentum of the parcel reaches an altitude where it matches the environmental wind (∼1,070 m), the sign of u ′ w ′ switches, and the flux is once again downgradient. The key to the behavior of momentum fluxes is the jet in the mean wind profile u (black solid line in Figure 1), and the resulting within-cloud wind profile u cld that looks like the mean profile, but shifted upward (orange dots in Figure 1).
In summary, the existence of upgradient flux in the BOMEX case appears to depend on two ingredients: (1) a mechanism, such as convective plumes, to transport momentum vertically a significant distance with relatively slow adjustment to the environmental wind, and (2) the presence of a jet, that is, a maximum of horizontal environmental wind magnitude, within the layer of vertical transport. Stated more colloquially, upgradient momentum flux in BOMEX occurs because of "nonlocal" transport in the presence of a jet.
In atmospheric boundary layers, a strong local maximum or "jet" appears more commonly in profiles of mean wind than in l or r t . For this reason, it may be more important to adequately parameterize upgradient fluxes in momentum than in l or r t .

Budget Terms That Give Rise to Upgradient Momentum Fluxes
What physical processes give rise to upgradient momentum fluxes? To address that question, we may examine the budget of u ′ w ′ , which contains a breakdown of all processes that contribute to the u ′ w ′ tendency. We first write down a prognostic budget equation for u ′ w ′ that neglects terms involving w, which is small, and uses conventional closures for the pressure terms, following André et al. (1978) and Bougeault (1981): (2) Here t denotes time, g is the acceleration due to gravity, ′ v is the perturbation virtual potential temperature, vs is a basic-state reference temperature, = (z) is the mean air density, and is a turbulence time scale described in Larson (2017). Also, C 6 and C 7 are constants (more precisely, slowly varying functions as described in Golaz et al., 2007) that govern the strength of pressure terms related to turbulence and buoyancy fluctuations, respectively. This equation makes use of the anelastic approximation.
Each of the terms on the right-hand side tends to change the value of u ′ w ′ , and each term may be interpreted physically: • The turbulence production term generates u ′ w ′ when there exist updrafts and downdrafts (w ′2 > 0) and a vertical gradient of the mean wind ( u∕ z ≠ 0). In this circumstance, a draft brings momentum to a different altitude with a different value of u, thereby generating u ′ w ′ . • Positive buoyancy (as measured by ′ v ) tends to accelerate a parcel upward (e.g., cause w is offset by the nonhydrostatic pressure force. • The return-to-isotropy term is a pressure term that relaxes the parcel horizontal velocity, u ′ , back to the environmental value (i.e., u with a time scale . This relaxation alters u ′ w ′ when u ′ and w ′ are already correlated with each other. This pressure term represents the lateral push exerted on a parcel by a horizontal pressure gradient. This term is included explicitly in some mass-flux schemes (Gregory et al., 1997) but not in others (Kain & Fritsch, 1990;Schneider & Lindzen, 1976;Tiedtke, 1989). • The turbulent advection (i.e., flux-of-flux) term advects the flux (u ′ w ′ ) by the perturbation vertical velocity, w ′ . Because this term is preceded by a vertical derivative, its vertical integral is zero if the boundary fluxes are zero. That is, the term acts as neither a sink nor a source but instead redistributes u ′ w ′ vertically.
Which of the aforementioned terms contribute to upgradient fluxes? To address this question, we simplify and rearrange equation (2) in order to find a diagnostic expression for u ′ w ′ . Specifically, we neglect the time tendency term, move the return-to-isotropy term to the left-hand side, and multiply by ∕C 6 . We find (3) The turbulence production and return-to-isotropy terms have conspired to produce an eddy diffusivity term with a positive eddy diffusivity, K, defined as If the buoyancy/pressure and turbulent advection terms are neglected, then there remains only downgradient diffusion of momentum flux. This is how momentum is treated in some low-order closure models (e.g., Bogenschutz & Krueger, 2013).
Therefore, an upgradient flux is made possible only by the existence of nonzero buoyancy production and/or turbulent advection (flux-of-flux) terms. The sign and magnitude of these terms for the BOMEX case is indicated in the upper left panel of Figure 7. In the layer of upgradient momentum flux above the jet maximum (770-1,070 m), the buoyancy production and turbulent advection terms are positive. They tend to increase u ′ w ′ , pushing it toward upgradient (positive) values.
The turbulent advection term turns out to be large in the BOMEX shallow cumulus case (Figure 7) but small in the DYCOMS-II RF02 stratocumulus case ( Figure 13). Why? By integrating over the subgrid PDF of Larson et al. (2002), the turbulent advection term may be modeled, neglecting density variations, as (Larson & Golaz, 2005) where u is a layer average, which is averaged over both clear and cloudy areas. A similar meaning is ascribed to w ′ . We define u ′ and w ′ in this way so that In the upper layer (>1,070 m), the within-cloud momentum flux is partially offset by the environmental momentum flux, leading to a weak layer-averaged momentum flux (see also Figures 1 and 2).
Here a 1 varies slowly. Hence, the term takes the form of a simple advection of u ′ w ′ by a "convective velocity" or "skewness velocity," w ′3 ∕w ′2 (Zilitinkevich et al., 1999). The term is large when skewness is large. Large, positive skewness of vertical velocity is often associated with strong plumes of small areal coverage, as would be found in cumulus layers like BOMEX, but not in stratocumulus layers like DYCOMS-II RF02.
The turbulent advection term is the place in the higher-order budget that represents the nonlocal transport that we conjecture is so helpful to produce upgradient momentum fluxes. The nonlocal character arises from the term's vertical derivative, ∕ z, which connects two different vertical levels. A similar term with a similar vertical derivative causes nonlocal transport in mass-flux schemes (Plant, 2010;Plant & Yano, 2016).

Why is the Correlation of u ′ and w ′ Small in the Upper Part of the Cumulus Layer?
In the upper (third) layer, from midcloud (1,070 m) up to the cloud top, the linear correlation coefficient between u ′ and w ′ has small magnitudes (about 0.05 or less). These correlations are lower than those found in the subcloud layer from 0 to 500 m (see Figure 2). The weakness of the momentum transport aloft is counterintuitive, given that the transport of heat and moisture is strong there. For example, the correlation between w ′ and total water mixing ratio (vapor + liquid), r ′ t , has a magnitude of 0.2 at 1,450 m in altitude, and so does the correlation between w ′ and the liquid water potential temperature, ′ l , (not shown). If we were to attempt to model the flux using eddy diffusivity, the eddy diffusivity for momentum would need to be much weaker than that for scalars. That is, the turbulent Prandtl number is much less than 1.  Each gravity wave is assumed to have a phase speed c that approximates u cld at its launch level. At an altitude of 1,020 m, upwelling gravity waves from below carry positive u ′ w ′ , according to the first theorem of Eliassen and Palm (equation (6)). Likewise, downwelling waves from above also carry positive u ′ w ′ , although they may be halted at critical levels (u ≈ c) before they reach 1,020 m. The net effect is u ′ w ′ env > 0.
Part of the reason that the momentum transport is weak is that the upper two thirds of each cloud nearly drifts with the environmental wind. It is true that there is a nonzero perturbation u ′ within cloud, as evidenced by the fact that u cld ≠ u in the upper two thirds of the cloud layer ( Figure 1). Therefore, the clouds do serve as an obstacle to the flow. However, u ′ within cloud is weak, because the mean wind profile contains a jet, rather than changing monotonically with altitude, like the profiles of water vapor or potential temperature. Because of the jet, a cloudy parcel at 500 m in altitude has about the same eastward wind speed (−7 m/s) as the environment at 1,500 m. Moreover, u cld = u at 1,020 m in altitude ( Figure 1).
Furthermore, in the upper layer, u ′ w ′ is negative within cloud, but u ′ w ′ is positive outside of cloud ( Figure 3). The within-cloud and environmental values of u ′ w ′ partially cancel each other, further weakening the layer-averaged values of u ′ w ′ .

Why Is u ′ w ′ Positive in the Environment?
We hypothesize that the positive momentum flux in the environment is due to gravity waves (Bretherton & Smolarkiewicz, 1989;Fovell et al., 1992;Kuettner et al., 1987). (The gravity waves we discuss here occur in the clear air alongside cumulus clouds; we do not discuss the more commonly studied gravity waves that originate at cloud top and propagate upward.) The motions in the clear air between the clouds are likely to be gravity waves for two reasons. First, the clear air between clouds is stably stratified. The stable stratification provides a restoring force for gravity but does not permit buoyantly driven turbulence. For example, the temperature perturbations in our simulation are ≤0.3 K. Given the stratification, a buoyant parcel could be lifted at most about 100 m. Second, the heat flux (or moist static energy flux) is small in the cloud layer ( Figure 4). This is symptomatic of nonbreaking gravity waves. In contrast, if the environmental motions consisted, for example, of turbulent horizontal rolls driven by vertical wind shear, then the heat flux would be expected to be negative, as warm air is pulled down and colder air is pulled up.
The gravity waves may have u ′ w ′ env > 0 because they are triggered by air motions that have positive correlation of u ′ and w ′ . For instance, when a cloud first grows up out of the subcloud layer, we have w ′ > 0 and u cld > u env (Figure 1), implying that u ′ > 0 and u ′ w ′ env > 0. That is, a buoyant parcel that impinges on the base of the cloud layer pushes environmental air up and to the east (Clark et al., 1986). More quantitatively, the sign of the momentum flux by gravity waves may be further analyzed with the aid of Eliassen and Palm's first theorem (Eliassen & Palm, 1961): where () env denotes a value averaged over the environment, p ′ is a pressure perturbation, and c is the phase speed of a gravity wave. We have restricted the averaging to the environment because convective turbulence contaminates the values in cloud. Equation (6) is valid for small (linear) waves and does not require that the Coriolis force be important. To rationalize why u ′ w ′ env > 0, we examine the signs of w ′ p ′ env and ( u env − c ) in equation (6). For simplicity, let us start by analyzing u ′ w ′ env at the altitude at which u cld ≈ u env (∼1,020 m; see Figure 5). The sign of w ′ p ′ is the sign of a gravity wave's vertical group velocity (Lindzen, 2005). First, consider an upwelling gravity wave that emanates from vertical motion at a source level below 1,020 m and propagates upward to 1,020 m (see Figure 5). Such a wave might be excited, for example, by a buoyant thermal impinging on the base of the cloud layer, which is stably stratified. Because that wave is traveling upward, it has w ′ p ′ env > 0. Now assume that the phase speed of the wave, c, roughly equals u cld at its source level. If c stays constant with altitude, then as the wave propagates upward, it has (u env − c) < 0 at every level up to and including 1,020 m because of the shape of the wind profile in BOMEX. Hence, by equation (6), u ′ w ′ env > 0 at 1,020 m for upwelling waves. If there is any gravity wave that originates above 1,020 m and propagates downward, then it has w ′ p ′ env < 0. If it also has c equal to u cld at its source level, then, by similar reasoning, it will initially have (u env − c) > 0, and hence, once again, we have u ′ w ′ env > 0 down to its critical level at c = u env . Near the critical level, the downwelling wave is hypothesized to be absorbed or reflected. Hence, downwelling waves are expected to be absorbed or reflected before they have a chance to contribute significantly at 1,020 m.
If we look at the momentum flux at a higher level, say, 1,300 m, but again assume that waves emanate with c ≈ u cld , then we see that some upwelling waves that are launched within the altitude range 1,020-1,300 m have u ′ w ′ env < 0, but other upwelling waves that arrive at 1,300 m from launch altitudes below this will have u ′ w ′ env > 0, and so will downwelling waves from above that arrive at 1,020 m. Therefore, once again, we expect that the net momentum flux will be positive (u ′ w ′ env > 0).
Unfortunately, testing these hypotheses would involve analyses of the environmental air that are beyond the scope of this study. Nevertheless, we have offered the hypotheses in order to spark further interest in the question of why u ′ w ′ env has the opposite sign of u ′ w ′ cld in the upper part of the cloud layer. Figure 7. Budgets of u ′ w ′ and v ′ w ′ from the Barbados Oceanographic and Meteorological Experiment (BOMEX) shallow cumulus case. Here "buoyancy" denotes (g∕ vs )u ′ i w ′ , and "pressure" denotes the sum of the C 7 (pressure-buoyancy) and C 6 (return-to-isotropy) terms in equation (2). Cloud Layer Unified By Binormals' (CLUBB) buoyancy flux has the wrong sign at cloud base, which causes errors in other CLUBB fields near that altitude. Away from that altitude, the CLUBB profiles qualitatively match those from System for Atmospheric Modeling large-eddy simulation (LES).
Although small-amplitude, linear gravity waves transport momentum vertically, they do not transport significant amounts of heat content (Holton, 2004, and Figure 4). The different fluxes of momentum and scalars in gravity waves may help explain the different fluxes of momentum and scalars in the BOMEX shallow cumulus layer. But in addition, whereas the momentum flux in cumuli is mostly offset by the momentum flux in the environmental waves, the scalar fluxes in cumuli are not so offset.

Budget Terms That Collectively Allow Weak Momentum Flux Aloft
The budget of u ′ w ′ (Figure 7) sheds further light on why u ′ w ′ is small in the third, upper layer. In that layer, turbulence production, −w ′2 u∕ z, is large in magnitude, and if the main term balancing it were return to isotropy, −C 6 u ′ w ′ ∕ , then u ′ w ′ would be large too (see equation (2)). In fact, however, turbulence production is mainly balanced by buoyancy production, (g∕ vs )u ′ ′ v , (see Figure 7). This allows the return-to-isotropy term, and hence u ′ w ′ , to be small.
Why is u ′ ′ v large and positive near cloud top? In the upper part of the cloud layer, u cld < u env ( Figure 1) and hence u ′ < 0 in and near cloud. Furthermore, near cloud top, cloudy parcels overshoot their level of neutral buoyancy. Thus, v cld < v env ( Figure 6) and ′ v < 0. Therefore, u ′ ′ v > 0.

Implementation of Prognostic Momentum Flux Equations in a Single-Column Model
The analysis of the u ′ w ′ budget terms in sections 3 and 4 has two main implications for parameterizing u ′ w ′ : 1. The parameterization of upgradient u ′ w ′ immediately above the jet maximum is aided by a strong turbulent advection term. 2. The parameterization of weak u ′ w ′ in the upper two thirds of the cloud layer is aided by a strong, positive buoyancy production term.
In order to prognose u ′ w ′ , we need to calculate each term on the right-hand side of the corresponding budget equation, equation (2). The turbulence production term has no need of closure, and the return-to-isotropy term has already been closed in equation (2) by a simple relaxation formula. The buoyancy and the turbulent advection terms are closed with the aid of CLUBB's PDF, as described presently.
This parameterization approach differs from that of Deardorff (1972) or Holtslag and Moeng (1991). Although Deardorff (1972) recognizes the importance of the buoyancy/pressure term, and Holtslag and Moeng (1991) recognizes the importance of the turbulent advection term, neither paper closes those terms Figure 9. Budgets of u ′ w ′ and v ′ w ′ , as in Figure 7, but from the Rain in Cumulus over the Ocean (RICO) shallow cumulus case. Cloud Layer Unified By Binormals' (CLUBB) budgets qualitatively resemble those of large-eddy simulation (LES) at most altitudes.
by integration over a subgrid PDF. Integration over a PDF allows for a closer connection to the underlying governing equations (although it does require an empirical assumption about the PDF's shape).

Diagnosis of the Buoyancy Term, u
Using an Assumed PDF Shape The horizontal flux of virtual potential temperature, u ′ i ′ v , where u i denotes u or v, can be written in terms of fluxes of conserved variables, u ′ i ′ l and u ′ i r ′ t , plus the flux of cloud liquid water mixing ratio, u ′ i r ′ c Moeng, 1998): Here R d and R v are the gas constants of dry air and water vapor, respectively; 0 = R d ∕R v ; L v is the latent heat of vaporization; c p is the heat capacity of air; 0 is a reference potential temperature; p 0 is the corresponding reference pressure; and p is air pressure itself. With this formula, the problem of diagnosing u ′ i ′ v reduces to the diagnosis of two fluxes of conserved variables, u ′ i ′ l and u ′ i r ′ t , plus the diagnosis of u ′ i r ′ c . Diagnosing the conserved variables does not require any assumption about PDF shape. For conciseness, let each conserved-variable flux be denoted as it is not needed for any other purpose. Hence, we choose to Figure 10. Profiles of mean winds and momentum fluxes, as in Figure 8, but from the Rain in Cumulus over the Ocean (RICO) shallow cumulus case. Displayed are profiles from System for Atmospheric Modeling (SAM) large-eddy simulations (LESs; black solid line), a version of Cloud Layer Unified By Binormals (CLUBB) that prognoses momentum fluxes (red dashed lines), and a version of CLUBB that diagnoses momentum fluxes using downgradient diffusion (green dotted lines). Prognosing momentum flux (red dashes) produces shallow layers of upgradient flux in u ′ w ′ near 1,000 m and in v ′ w ′ just below 500 m. diagnose u ′ i ′ as simply and inexpensively as possible. For that purpose, we write a simplified form of the André et al. (1978) in which the third-order turbulent advection and Coriolis terms are omitted: Here, is the same return-to-isotropy time scale introduced in equation (2). Also, C ps u and C pi u are pressure terms related to mean shear and return to isotropy for the u ′ i ′ equation. Next, we drop the time tendency term and rearrange to find In CLUBB's current implementation of this equation, we set C pi u = C 6 and C ps u = 0.
Diagnosing u ′ i r ′ c is more difficult because r c is not conserved under condensation. Here we diagnose u ′ i r ′ c by integrating it analytically over CLUBB's subgrid double Gaussian PDF. Specifically, we integrate over the Figure 11. Budgets of u ′ w ′ and v ′ w ′ , as in Figure 7, but from the Atmospheric Radiation Measurement (ARM) continental shallow cumulus case. The turbulent production term (brown dashed line) is opposite in sign to the pressure term (purple solid line) nearly everywhere, suggesting the fluxes are downgradient. Cloud Layer Unified By Binormals' (CLUBB) budgets qualitatively match System for Atmospheric Modeling's budgets. LES = large-eddy simulation.
Analytic Double Gaussian 1 (ADG1) PDF of Larson et al. (2002). The procedure and equations are analogous to those used to diagnose w ′ r ′ c , which are described in Larson et al. (2002) and Larson (2017). Because ADG1 assumes that velocity components are uncorrelated to thermodynamic variables within an individual mixture component, the formula for u ′ i r ′ c is simple: Here 0 ≤ 1 ≤ 1 is the mixture fraction, u i1 and u i2 are the values of u i within the first and second Gaussian components, and u i is the grid-box mean of u i . Also, r c1 and r c2 are the values of r c within the first and second Gaussian components, which are calculated as in Larson et al. (2002).
The component means u i1 and u i2 are calculated analogously to the equivalent quantities for scalars, such as r t1 and r t2 , as described in Larson et al. (2002). The information needed to do so is available because u ′ w ′ is predicted.

Compute Turbulent Advection Term Using PDF.
As with the buoyancy production term, the turbulent advection term is also closed by integration over the ADG1 PDF (Larson & Golaz, 2005). The result was presented already in equation (5): where Here is a tunable parameter, and c xy denotes the linear correlation coefficient between two variables x and y. The term involving correlations is added to ensure realizability (Larson & Golaz, 2005).
It was found necessary to replace a 1 by 3a 1 in equation (11) in order to produce an upgradient flux u ′ w ′ in the BOMEX simulation presented below. (Tuning has a similar effect on a 1 , but unfortunately, it degrades the cloud fields.) This is undesirable, because it breaks the consistency that is such an attractive feature of assumed PDF closures. The need to do so may be related to the fact that CLUBB's buoyancy term (g∕ vs )u ′ w ′ is underestimated in the upgradient region (see upper right panel of Figure 7). This motivates us to improve the underlying assumption about the PDF shape in future work.

Numerical Solution
Prognosing momentum fluxes is computationally expensive, but in the case of the CLUBB code, the cost is comparable to diagnosing u ′ w ′ in terms of downgradient diffusion, as in equation (1). More precisely, the total cost of CLUBB-excluding radiation, microphysics, output statistics, and any host model computations-increases only about 2% when prognostic momentum for both u and v is turned on. The extracost is small because CLUBB is able to reuse prior computational work.
CLUBB prognoses a slightly embellished form of equation (2) that includes some extra, minor terms (see Equation 4.7 of Larson, 2017). However, for simplicity of discussion, we will treat equation (2) as CLUBB's equation. Equation (2) is identical in form to the prognostic equation for the vertical turbulent flux of a scalar , such as l or r t , if one replaces u everywhere with . In CLUBB, w ′ ′ is solved simultaneously with in a single matrix equation (see Larson, 2017, section 4.2). The buoyancy production term is treated explicitly on the right-hand side of the matrix, and the other terms in equation (2) are treated implicitly on the left-hand side of the equation. Then the left-hand side matrix is LU (Lower-Upper) decomposed (Press et al., 2007). The LU decomposition is expensive, but once it is completed, finding the solution for a particular right-hand side vector is cheap. There is one right-hand side vector for each scalar.
Prognosing momentum fluxes in CLUBB is cheap because the LU decomposition that has already been computed for the scalar fluxes can be reused for the momentum fluxes. Solving for the momentum fluxes merely requires adding a vector to the right-hand side for each horizontal momentum component. The cost of solving for this "extra" right-hand side vector is about the same as diagnosing u ′ w ′ via eddy diffusivity.

Comparison of Single-Column and LESs
Now we compare simulations of momentum-related budgets and profiles calculated by SAM and CLUBB-SCM. Our goal is to illustrate that prognosing momentum fluxes produces an overall better match to LES than does diagnosing momentum fluxes using CLUBB's previous downgradient diffusion formulation. The match is significantly better in cases that exhibit upgradient fluxes in the cloud layer (i.e., the BOMEX and RICO marine shallow cumulus cases), and the match is comparable in the cases that exhibit downgradient fluxes in cloud (i.e., the ARM continental cumulus case and the two marine stratocumulus cases, DYCOMS-II RF02 and RF01). In this paper, all five cloud cases are simulated using the same set of adjustable model parameter values, regardless of whether the simulation prognoses momentum fluxes or uses downgradient diffusion.
The budgets of u ′ w ′ and v ′ w ′ as calculated by SAM and CLUBB are shown for the BOMEX case in Figure 7.
(All CLUBB budgets that we present come from the version of CLUBB that prognoses momentum fluxes; the version that uses downgradient diffusion has no momentum flux budget.) In these budget plots, the buoyancy (C 7 ) and return-to-isotropy (C 6 ) pressure terms have been combined. CLUBB's budget terms match SAM except for buoyancy production near cloud base, where CLUBB has the wrong sign. In CLUBB's budgets, the error in buoyancy production is compensated mostly by an error in the pressure term (i.e., C 7 plus C 6 terms). Figure 8 displays the profiles of mean wind and momentum fluxes for BOMEX. Two versions of CLUBB are compared. They differ only in that one version diagnoses the momentum fluxes using downgradient diffusion as in equation (1) green dotted lines) and the other prognoses momentum fluxes as in equation (2) (red dashed lines). In order to avoid producing an upgradient sign of u ′ w ′ , downgradient diffusion produces the wrong sign of u∕ z between 770 and 1,070 m. The net outcome is to spuriously elevate and thicken the jet, leading to an excessive wind magnitude at the ocean surface. In contrast, prognosing the momentum fluxes produces a region of upgradient flux and produces a more realistic jet shape and surface wind values. Apparently, the accuracy of the parameterized winds and momentum fluxes is not seriously impaired by the errors in the budget terms.
We note parenthetically that v ′ w ′ is upgradient in the SAM-LES between about 50 and 200 m in altitude (see lower panels of Figure 8). In some respects, this upgradient region resembles the one in u ′ w ′ that we have been discussing. For instance, both upgradient regions reside immediately above an extremum in wind, and both contain vertical turbulent transport. However, unlike the upgradient region in u ′ w ′ , the upgradient , and a version of CLUBB that diagnoses momentum fluxes using downgradient diffusion (green dotted lines). In this case, CLUBB's wind profiles are not well mixed enough, but prognosing momentum fluxes (red dashes) produces better-mixed profiles than does downgradient diffusion (green dots). region in v ′ w ′ exhibits opposite signs for turbulent advection and turbulence production (Figure 7), which might be expected to be associated with a downgradient flux. The opposite signs imply that return to isotropy does not dominate the pressure term in the upgradient region of v ′ w ′ . CLUBB is able to produce upgradient values of v ′ w ′ , but they are weak and confined to the altitude range of 100 to 175 m (Figure 8). This paper will not dwell on this upgradient region because it is presumably influenced by surface-layer physics that is beyond the scope of this study.
For the RICO case, the budgets of u ′ w ′ and v ′ w ′ (Figure 9) and the profiles of mean wind and momentum fluxes ( Figure 10) look similar to those of BOMEX. As for the BOMEX case, CLUBB misrepresents the budget terms near cloud base, but nevertheless, prognosing momentum fluxes improves the profiles of wind and fluxes. SAM-LES indicates that u ′ w ′ is upgradient in the region between 750 and 1,500 m, and v ′ w ′ is upgradient between 100 and 600 m ( Figure 10). Prognosing momentum fluxes in CLUBB produces upgradient fluxes through a portion of those regions.
For the ARM continental case, the SAM-LES budgets ( Figure 11) and profiles of wind and fluxes ( Figure 12) look quite different than those for the marine cumulus cases. The eastward wind profile has no jet. Therefore, u ′ w ′ is downgradient everywhere. The lack of upgradient fluxes is reflected in the budget of u ′ w ′ . Recall that in BOMEX and RICO, the turbulent advection term, − w ′2 u ′ ∕ z, is opposite in sign to the turbulence production term, −w ′2 u∕ z, just above the jet max. This is because above the jet max, turbulent patches with large w ′2 have u ′ > 0 in a layer in which u∕ z > 0 (upper left panel of Figure 10). Therefore, turbulence production can be offset by the turbulent advection term (and buoyancy term), allowing the pressure and turbulence production terms to have the same sign and the flux to be upgradient. ARM, in contrast, has no jet maximum, and hence, u ′ and u∕ z have the opposite sign everywhere (upper left panel of Figure 12).
Consequently, turbulence production must be offset by pressure (upper left panel of Figure 11). Finally, because the pressure term is dominated by return to isotropy, the flux is downgradient.
As for BOMEX and RICO, ARM contains a region of upgradient v ′ w ′ near the ground surface. ARM's budget of v ′ w ′ resembles that of BOMEX near the ground, suggesting that the physics of upgradient flux is similar there. CLUBB is able to produce an upgradient flux from 400 to 550 m in altitude, but it is only very weakly upgradient.
The budgets of u ′ w ′ and v ′ w ′ for the DYCOMS-II RF02 stratocumulus case are shown in Figure 13. Near the inversion, SAM shows a balance of buoyancy production and pressure. CLUBB shows a similar balance, except that there is a spuriously large contribution of turbulence production.
The profiles of mean wind and momentum fluxes are shown for the DYCOMS-II RF02 case in Figure 14.
Prognosing momentum fluxes (red dashed lines) is comparable in accuracy to downgradient diffusion (green dotted lines). Figure 18. Cloud and thermodynamic fields as in Figure 17, but for the Rain in Cumulus over the Ocean (RICO) shallow cumulus case.

Journal of Advances in Modeling Earth Systems
10.1029/2019MS001743 Figure 19. Cloud and thermodynamic fields as in Figure 17, but for the Atmospheric Radiation Measurement (ARM) continental shallow cumulus case.
The story is similar for DYCOMS-II RF01. Once again, CLUBB overpredicts the magnitude of turbulence production near cloud top ( Figure 15), but prognosing momentum fluxes produces comparably accurate profiles (Figure 16), except for slightly more well-mixed profiles, in better agreement with the LES.

How Does Prognosing Momentum Flux Influence Cloud and Thermodynamic Fields?
Perhaps surprisingly at first, prognosing momentum fluxes has little impact on cloud and thermodynamic fields (see Figures 17-21). The largest impact is on cloud fraction in RICO (Figure 18), but even there the effect is minor, given the large change in horizontal winds.
The small impact is less surprising if one considers how CLUBB's moments are coupled by the higher-order equations. A change in momentum flux or mean horizontal wind has a direct influence mainly on the horizontal velocity variances, u ′2 and v ′2 . Those moments, in turn, affect the rest of the model mainly through their influence on w ′2 . Namely, the vertical component of turbulence, w ′2 , is boosted when the horizontal components, u ′2 and v ′2 , are large. This is because of return to isotropy, a nonhydrostatic pressure effect that strives to equalize the magnitude of all turbulence components. However, this effect is weak because the turbulence in our cases is generated mainly by buoyant convection, not wind shear.
In coupled ocean-atmosphere simulations, changes in winds are expected to affect cloud fields through changes in surface latent and sensible heat fluxes. However, the single-column cases that we test here are too constrained to estimate such effects.

Conclusions
A LES of the BOMEX shallow cumulus case reveals a three-layer structure of momentum flux in the direction of strongest wind: 1. From the ocean surface up to the jet maximum just above cloud base, the momentum flux is strong and downgradient. 2. Starting at the jet maximum and extending to the altitude where u ′ w ′ cld = 0, the momentum flux is upgradient. The upgradient flux is the result of parcels of air that carry low-magnitude momentum from below cloud base upward while adjusting only slowly to the environmental wind. When such a parcel rises above the jet maximum, the local environmental wind gradient, u∕ z, reverses sign. Until the parcel speed adjusts to the environmental speed (u cld = u env ), the momentum flux is upgradient. The BOMEX case has no cloud organization and hence differs from those deep convective situations in which convective organization matters (LeMone, 1983). However, Grubišić and Moncrieff (2000) do note that upgradient momentum fluxes can occur in two-dimensional squall lines if the ambient flow is jetlike. 3. Above the middle, upgradient layer, the momentum flux is weak. In part, this occurs because of the presence of a jet maximum, which means that parcels can rise long distances without suffering strong turbulent deviations u ′ from the mean horizontal wind u. In addition, the momentum flux in the upper layer in BOMEX is partly balanced by an opposite-signed momentum flux, which we hypothesize is carried by gravity waves in the environment. In contrast, scalars are only weakly transported by gravity waves in the environment. The weak layer-averaged momentum flux and strong layer-averaged scalar fluxes in the BOMEX shallow cumulus case provide an interesting contrast to noncloudy stably stratified layers, where increased stability (more precisely, Richardson number) leads to an increased ratio of momentum to scalar flux (more precisely, increased turbulent Prandtl number; see He et al., 2019, and references therein).
The budget of u ′ w ′ has implications for parameterization development. First, in order to produce an upgradient flux in a layer above the jet maximum, the turbulent advection and buoyancy production terms are crucial. Second, in order to produce weak momentum flux in the upper part of the cloud layer, it is crucial to produce a strong buoyancy production term.
To parameterize momentum fluxes, we prognose them. In those prognostic equations, we parameterize the buoyancy production and turbulent advection terms by integrating them over the ADG1 subgrid PDF. Discrepancies appear between some terms in the parameterized budget of momentum flux and the corresponding terms simulated by LES.
Nevertheless, we find that prognosing the momentum fluxes yields an upgradient momentum flux where expected and accurately simulates the wind profile in the BOMEX and RICO cases, whereas CLUBB's prior downgradient diffusion method distorts the jet profile in order to ensure that the fluxes are downgradient.