Unified Parameterization of Convective Boundary Layer Transport and Clouds With the Thermal Plume Model

The representation of stratocumulus clouds, and of the stratocumulus to cumulus transitions which are ubiquitous features of marine boundary layer clouds, remains a challenge for climate models. We show how a mass flux representation of boundary layer convective structures combined with an eddy diffusivity scheme, the “thermal plume model,” first developed to represent cumulus clouds, can also adequately simulate stratocumulus and the stratocumulus to cumulus transition in a climate model. To achieve this, the detrainment formulation, in which detrainment increases for increasing negative buoyancy, has to be slightly modified: the buoyancy of a thermal plume parcel of air is computed by comparing the virtual potential temperature θv,th of the parcel with that of the surrounding environment θv,env at a given distance above instead of at the same level. This is consistent with the picture of detrained air parcels that experience some overshoot and reach a final destination at a level lower than the one at which they effectively leave the cloud or organized convective plume. The impacts of this modification are documented both for selected cases of stratocumulus, in comparison with large‐eddy simulations, and in full 3‐D climate simulations, in comparison with satellite observations of cloud cover. The modified scheme provides a uniform treatment of the dry convective boundary layer, of cumulus clouds, of stratocumulus, and of the transition from stratocumulus to cumulus. It is included in the most recent version of the LMDZ atmospheric general circulation model.


Introduction
Stratocumulus clouds cover on average 20% of the Earth's surface (Wood, 2012). They are prevalent in the eastern subtropical oceans, which are characterized by strong large-scale subsidence. They strongly impact the incoming solar radiation as they reflect a large part of the shortwave (SW) radiation. The underestimation of stratocumulus shadowing effect is in part responsible for systematic warm biases obtained in coupled atmosphere-ocean models over the east tropical oceans (Hourdin et al., 2015). Stratocumulus clouds have a smaller effect on the outgoing longwave radiation at the top of the atmosphere because of the similar temperatures of stratocumulus cloud top and sea surface. As a result, stratocumulus have a negative net radiative effect on the Earth's radiative balance (Hartmann et al., 1992) and small variations of the coverage and thickness of stratocumulus clouds can produce a radiative effect comparable to those associated to the increase of greenhouse gases. Part of the dispersion of climate sensitivity across models can be attributed to differences in the modeled responses in regions of stratocumulus and in transition between stratocumulus and cumulus clouds over the subtropical oceans in particular (Bony & Dufresne, 2005;Nam et al., 2012;Klein et al., 2017;Zelinka et al., 2016;Zhang et al., 2013). While the structure of the boundary layer and the conditions that favor the formation of stratocumulus clouds are now well understood, the realistic simulation of stratocumulus and transition from stratocumulus to cumulus clouds in global climate models still remains problematic (Neggers et al., 2017).
The review by Wood (2012) provides a very complete description of the processes at work in stratocumulus. Like fair weather or trade wind cumulus, stratocumulus clouds are fed by convection of water evaporated at the surface, but they are typically only one to a few hundred meters thick and can cover extensive regions for months. They are capped indeed by a strong inversion, created usually by the large-scale subsidence associated with the descending branch of the Hadley-Walker circulation. The inversion is characterized by a sharp transition in most variables. Convection in stratocumulus-topped boundary layers is also driven for 10.1029/2019MS001666 a large part by the infrared radiative cooling that occurs near the cloud top. The characterization and role of organized convective structures in the maintenance and diurnal cycle of stratocumulus clouds have been investigated in many studies, in particular, using large-eddy simulations (LES; Brient et al., 2019;Davini et al., 2017;de Roode et al., 2016;Dussen et al., 2013;. The air feeding the cloud layer is mainly transported by positively buoyant updrafts rising from the surface. When updrafts reach the (very strong) inversion layer, they are rapidly stopped (over a distance of a few tens of meters typically) and create flat domes that spread horizontally at the inversion level. Evaporative cooling of cloudy air, when mixed with the dry troposphere, can also create negatively buoyant downdrafts which are often strong enough to penetrate down to the surface and to contribute, together with thermal updrafts, to an efficient mixing throughout the boundary layer (Brient et al., 2019). The coupling between convection, turbulence, and radiation in this very shallow transition layer is particularly challenging for parameterizations.
Stratocumulus clouds show a marked diurnal cycle. During daytime, the stratocumulus layer becomes slightly warmer than the subcloud layer due to the absorption of solar radiation and often thus becomes partially decoupled from the subcloud layer. As a result, turbulent mixing and the upward transport of moisture from the surface is reduced and the cloud thickness reduces.
The transition of stratocumulus to cumulus has been more specifically investigated in the past decade using a test case built on observations of the ASTEX field campaign in the northeast Atlantic (Dussen et al., 2013) and three composite transition cases (Sandu & Stevens, 2011) built from reanalyses in the northeast Pacific, based on a Lagrangian setup in which an air mass characterized by a stratocumulus-topped boundary layer is advected over warmer sea surface temperature and weaker subsidence. Intercomparisons have shown that explicit LES, with grid cell size of a few tens of meters, are sufficiently in agreement with each other to give a reliable picture of the processes at work (Dussen et al., 2013;de Roode et al., 2016) and to be used as a guide for the development of parameterizations and as a reference for their evaluation. As explained by Neggers et al. (2017; see Figure 1 in their paper), the conceptual sequence of the transition is well established: "After an initial period of gradual deepening, a thermodynamic decoupling takes place within the originally well-mixed boundary layer, after which a shallow cumulus cloud base emerges below the capping StCu layer. Subsequently the boundary layer deepening continues, with the capping cloud layer thinning and eventually breaking up. The transitional situation, consisting of shallow cumuli rising into a capping cloud layer, is sometimes recognized as a separate cloud regime." Studies of stratocumulus physics insist on top entrainment which characterizes the mixing occurring at cloud top and results from the interaction of various processes: radiation, turbulence, microphysics, and large-scale subsidence. Errors in the representation/parameterization of cloud top-entrainment mixing, and its coupling with boundary layer humidity and temperature, have been pinpointed as one key driver for low-cloud biases in climate models (Brient et al., 2019a). Grenier and Bretherton (2001) or Duynkerke et al. (2004) already identified top-entrainment parameterization as key for the representation of stratocumulus confirming previous works by, for example, Randall (1980) and Deardorff (1980), and several models include an explicit entrainment closure, in which a flux at the boundary layer top is computed as an entrainment velocity times the jump in the transported quantity w ′ q ′ = − w e Δq. Aircraft observations or LES have been used to constrain the entrainment velocity and evaluate different parametrizations (Caldwell & Bretherton, 2009;Caldwell et al., 2005). Some of those parameterizations are based on convective velocity, computed from the buoyancy of air parcel lifted from the surface and take into account cloud top evaporation (Lilly, 2002;Turton & Nicholls, 1987). Some parameterizations directly account for the transport by the downdrafts of negatively buoyant mixtures that sink away from the cloud top (Duynkerke et al., 1995). Attempts have been done as well to account directly for a radiative effect in this entrainment closure (Lock et al., 2000). Note that the way the term q is computed in entrainment closures is an issue by itself (Grenier & Bretherton, 2001;Lock, 2001).
Parameterizations differ also by the way they treat the boundary layer turbulent transport. While the role of nonlocal countergradient transport by convective structures has been recognized for a while (since at least the pioneering work of Deardorff, 1966), some schemes account for local down gradient turbulent transport only (Bretherton & Park, 2009), assuming that nonlocal transport does not matter for well-mixed quantities. Other schemes make use of the Troen and Mahrt (1986) approach, popularized by Holtslag and Boville (1993), in which a countergradient term is added to the local gradient of the transported quantity so as to

10.1029/2019MS001666
provide the correct convective flux, even in a slightly stable atmosphere. This approach is used in particular by Lock et al. (2000).
Here, we use a now classical approach for the representation of the nonlocal turbulent transport in the convective boundary layer which combines a mass flux representation (inherited from cumulus parameterizations; see, e.g., Arakawa & Schubert, 1974;Siebesma & Cuypers, 1995) of the organized ascending plumes, cells, or rolls that can transport directly conserved variables across several vertical model layers, with a classical K diffusion representation of small scale turbulence (see, e.g., Hourdin et al., 2002;Pergaud et al., 2009;Rio & Hourdin, 2008;Siebesma et al., 2007;Soares et al., 2004). This approach, often called "eddy diffusion mass flux" (EDMF) was first proposed by Chatfield and Brost (1987). It favors a continuous representation of dry and cloudy convective boundary layers, cumulus clouds being formed at the top of convective updrafts. The particular "thermal plume model" used as a basis for the present study relies on a single updraft aiming to represent the average impact of an ensemble of plumes (Hourdin et al., 2002;Rio & Hourdin, 2008). Such a mass flux parameterization can be further coupled with a statistical representation of the subgrid-scale distribution of water, in order to represent in particular the cumulus clouds that are created at the top of convective plumes or rolls Pergaud et al., 2009;Rio et al., 2010). The coupling of the thermal plume model with a bi-Gaussian distribution of the subgrid-scale saturation deficit (the difference of the total water with its saturation value) was shown to improve significantly the representation of cumulus clouds, both in 1-D simulations  and in 3-D global climate simulation with the LMDZ model .
However, in its original version, this scheme was unable to maintain stratocumulus clouds, in particular over the east tropical oceans. The problem arises from a tendency of the mass flux scheme to overshoot into the inversion layer and to bring too much tropospheric dry air into the boundary layer, breaking the fine equilibrium that controls the maintenance of the very thin stratocumulus. In order to overcome this difficulty, a trick was introduced in the previous version LMDZ5B of the model: The mass flux scheme was intentionally deactivated if a strong inversion was detected at cloud top. That way, the mass flux scheme was systematically deactivated under strong subsidence, and the stratocumulus were thus represented with a local Mellor and Yamada turbulent scheme coupled to a single-mode statistical cloud scheme .
The activation of different schemes, depending on environmental conditions, is a classical feature of many models. The parameterization by Lock et al. (2000) clearly advocates for the need of "two discrete mixing schemes [which] is to some extent undesirable [but] has been found essential" to avoid that cumulus rapidly evolve into stratocumulus. Their scheme is based on the identification of six regimes that depend on the co-occurrence in the same atmospheric column of stable and convective layers and presence of cumulus or stratocumulus clouds. In a similar way, Bretherton and Park (2009) identify several layers of different stability in the column, activating different parameterizations depending on the conditions encountered with, in particular, various entrainment parameterizations at the upper and lower interfaces of each turbulent layer. Köhler et al. (2011) propose to unify the treatment of dry convective and stratocumulus-topped boundary layer using an eddy-diffusion mass-flux approach, which corresponds to the situation sketched in Figure 1c of Lock et al. (2000). This situations is contrasted with the cases of cumulus clouds or coexistence of cumulus and decoupled stratocumulus (sketched in Figures 1e and 1f, respectively of Lock et al., 2000), where an inhibition exists at the base of the cumulus layer, reducing the cloud fraction (to typically less than 10%). Köhler et al. (2011) activate the mass flux component only if the contrast in potential temperature between 700 hPa and the surface is greater than 20 K. For other conditions, the mass flux scheme is deactivated and the Tiedtke scheme is used to represent cumulus convection.
Unification of the treatment of boundary layer clouds was a major driver of new approaches for the boundary layer parameterization, such as the Cloud Layers Unified by Binormals scheme based on a partial third-order turbulence closure coupled to joint probability distribution functions of vertical velocity, temperature, and moisture Guo et al., 2014).
Here we show how a slight modification of the thermal plume model proposed by Jam (2012) leads to a reasonable representation of boundary layer convection, cumulus clouds, stratocumulus, and transition from stratocumulus to cumulus clouds with a unified scheme. It is shown in particular that the compensating subsidence of the thermal rising plumes at the top of stratocumulus acts as a parameterization of boundary layer top entrainment and that this parameterization seems to reasonably account well for the sensitivity of this top entrainment to the cloud top radiative cooling and subcloud layer convective activity. This modification of the thermal plume model is one of the main improvements of the Institut Pierre-Simon Laplace (IPSL) climate model IPSL-CM for the coming sixth phase of the Coupled Model Intercomparison Project (CMIP; see e.g., Taylor et al., 2012), and the present paper also aims to serve as a reference of this model version for this particular aspect.
The paper is organized as follows: Section 2 presents a description of the boundary layer convection scheme and the model setup used to evaluate the parameterization. Section 3 documents the deficiencies of the original parameterization and presents the modification of the detrainment formulation. Finally, we present in section 4 the impacts of this modification in 1-D and 3-D simulations and evaluate the improvements of such a parameterization via comparisons with LES on different case studies and satellite observations.

The Thermal Plume Model
In the mass flux approach for the parameterization of boundary layer convection considered here, the collective behavior of a population of thermal plumes (or cells, or rolls) is represented through a unique thermal plume. thermal plume model (Hourdin et al., 2002;Rio & Hourdin, 2008;Rio et al., 2010), each column is divided into a mean ascending thermal plume of mass flux f = w th (where is the air density, and and w th the fractional cover and the vertical velocity of the plume), and a compensating subsidence in the environment of mass flux −f. The value of a model state variable within the thermal plume th is computed using the stationary plume conservation equation: where e and d are the lateral entrainment and detrainment of air toward and away from the plumes. The assumption is made that the environmental value of variable entering the plume equals its large-scale value. For variables conserved by the convective transport, such as liquid potential temperature l or total water q t , the source term is set to S ≡ 0. The continuity equation for the plume, which relates f to e and d, is also given by equation (1) for a source term equal to 0 and with a tracer equal to unity (S ≡ 0, ≡ 1 and th ≡ 1 in equation (1)). The vertical velocity is obtained by equation (1) as well with a source term S = a 1 B − a 2 ew 2 th ∕ where B = g( v,th − v )∕ v is the buoyancy ( v being the virtual potential temperature) that accelerates the plume and the second term a drag effect, with a 1 = 2∕3 and a 2 = 0.002. The total boundary layer vertical transport of reads where the eddy diffusivity K z = l mix S(Ri) √ TKE, l mix being a turbulent mixing length and S(Ri) a stability function that depends on the local gradient Richardson number Ri. The turbulent kinetic energy TKE is integrated in time from a local prognostic equation, following Yamada (1983) with technical implementation details given by Vignon et al. (2017). The fractional cover is then deduced as = f∕( w th ).
Parameterizations of the entrainment and detrainment rates = e and = d were derived and tuned so as to fit the collective behavior of the population of updrafts resolved in a LES using an original tracer based conditional sampling of organized updrafts Rio et al., 2010). The entrainment rate depends on the plume buoyancy and vertical velocity: where 1 = 0.9, values consistent with previous studies (Gregory, 2001;Siebert & Frank, 2003). The plume is mainly entraining in regions of positive buoyancy. It is the opposite for the detrainment rate which is favored in regions where buoyancy is negative, as suggested by observations (Bretherton & Smolarkiewicz, 1989). A satisfactorily correlation is obtained between LES results and parameterization with the following definition of : with c = 0.012 s −1 and d = 0.5. The first term corresponds to the buoyancy contribution to detrainment rate, while the second term tends to account for the fact that evaporation around the clouds can reinforce the negative buoyancy of extracted air parcels, a mechanism enhanced when the moisture contrast q t between the clouds and the environment increases.

The Cloud Scheme for Boundary Layer Clouds
For cumulus clouds, the distribution F(s) of the saturation deficit s within a mesh, at any height, can be approximated by a double Gaussian distribution. Thanks to a tracer based sampling of LES results,  demonstrated that one mode corresponds to thermal plumes and the second one to their environment. Based on these findings, a statistical cloud scheme was derived, using five parameters: the plumes fraction , the mean saturation deficits within environment s env and plumes s th , which are directly given by the thermal plume model, and their associated standard deviations, env and th , for which a parameterization was proposed. Considering that the major source of standard deviation for s is the exchange of air between the plume and its environment and that dispersion of s values is enhanced when the contrast s th − s env increases, standard deviations are parameterized as follows: and where b, c env , c th , 1 , and 2 are tunable parameters, and the last term, bq t , is a minimum width of the environment distribution introduced for a value of ≈ 0.
The values of b = 2 × 10 −3 , c env = 0.92, c th = 0.09, 1 = 0.4, and 2 = 0.6 were chosen using LES results by fitting independently the in-thermal and environment Gaussian distributions. In case of cumulus clouds, it is essentially the mode associated with the thermal plume which is saturated so that the cloud fraction is comparable to the fractional cover by thermal plumes. In the case of stratocumulus clouds, as will be shown hereafter, a 100% cloud cover can be reached for a plume fractional cover of a few percent only, the Gaussian modes associated with the plume and with its environment being both saturated.
The thermal plume model is activated before the cloud scheme. The condensation is taken into account in the computation of liquid potential temperature (considered as conserved variable in equation (1)) and virtual potential temperature involved in the buoyancy computation. Once e, d, and f are determined, equations (1) and (2) are applied to the total water and liquid potential temperature to compute tendencies associated with the boundary layer transport. From the thermal plume model computation as well, the parameters of the bi-Gaussian subgrid-scale distribution F(s) of the saturation deficit can be estimated as explained above. From this distribution, the cloud fraction ∫ ∞ 0 F(s)ds and cloud liquid content ∫ ∞ 0 sF(s)ds are finally computed. The computation of the conversion from cloud water to rainfall follows Sundqvist (1978): Rainfall starts to precipitate above a critical value (here 0.65 g/kg) with a time constant of half an hour. Following Sundqvist (1988), a fraction of the precipitation is reevaporated in the layer below and added to the total water of this layer before the statistical cloud scheme is applied. Note that the same cloud scheme is applied with a single mode of width s,env = b q t env when the thermal plume model is not activated (for stratiform clouds, for instance) while a different scheme is used for deep convection. Equations and details on the cloud scheme are given in Hourdin et al. (2013). The relevance of this bi-Gaussian parameterization for stratocumulus is discussed later on in the paper.

LES Case Studies
This study makes use of a series of LES simulations of classical case studies of boundary layer clouds.
The modification of the thermal plume model presented here was inspired and evaluated using results of the composite stratocumulus to cumulus transition case discussed by Sandu and Stevens (2011). This case, referred to as REF hereafter, was aimed at analyzing the transition from unbroken sheets of stratocumulus to fields of scattered cumulus and the processes controlling them. It was built by compositing the large-scale conditions encountered along a set of individual Lagrangian 3-day trajectories performed for the northeastern Pacific during the summer months of 2006 and 2007. The stratocumulus deck presents a pronounced diurnal cycle and begins to break up during the second day while the boundary layer deepens. Shallow cumulus appear after 24 hr under the gradually thinning and dissipating stratocumulus. The LES reference simulation used in this study was performed with the UCLA model, with a 35-m horizontal resolution and a 5-m vertical resolution up to 2,500 m. Two variations of this REF case, corresponding to a slower and a faster transition in cloud fraction, were derived in a similar manner by compositing over the trajectories experiencing the fastest, respectively, the slowest decrease in cloud fraction over the first 2 days (FAST and SLOW hereafter). The setup of the REF, FAST, and SLOW cases and the LES simulations are described in more detail in Sandu and Stevens (2011).
The scheme is further evaluated independently on a series of test cases that concern both stratocumulus clouds (FIRE), marine (RICO), and continental (ARM CU) fair weather cumulus, all performed with Meso-NH nonhydrostatic model (Lac et al., 2018;Lafore et al., 1998).
The FIRE I experiment, conducted off the coast of California, provides a comprehensive observational set of data on marine stratocumulus during July 1987. From this field campaign, a steady state case of stratocumulus with a well-observed diurnal cycle was derived and used for intercomparison purpose (Duynkerke et al., 2004). The LES reference simulation used here was performed with a time step of 1 s, a 50-m horizontal resolution on a domain of 200 × 200 grid points, and a 10-m vertical resolution.
The ARM case is derived from observations collected on 21 June 1997 at the Atmospheric Radiation Measurement site in Oklahoma, USA (Brown et al., 2002). This idealized case is typical of the diurnal cycle of shallow convection over land. The LES was performed using a horizontal domain of 128 × 128 points with a resolution of 50 m, a vertical domain of 100 levels, a resolution of 40 m, and a time step of 1 s.
The RICO experiment (RICO for Rain In Cumulus over the Ocean) focuses on precipitation processes at play in trade wind shallow cumulus. During RICO, significant precipitation was observed frequently, offering a unique opportunity to study the dynamics of shallow cumuli and precipitation. The LES simulation is performed with a horizontal resolution of 100 m × 100 m on a domain of 64 × 64 grid points, and a vertical resolution of 40 m in the first 4 km. The time step is 1 s and the duration 24 hr.
In order to evaluate the parameterized mass flux or detrainment against LES results, a thermal plume sampling is used to identify the coherent ascending structures in the LES : A grid point is supposed to belong to a thermal plume if its vertical velocity w is positive and if the difference between the total specific humidity q t and its mean q t exceeds q (where q is the standard deviation of q t at a given vertical level). This sampling was shown to be similar to a surface-based tracer sampling in Couvreux et al. (2010).

The LMDZ Climate Model
The parameterization development work presented here was conducted and tested in the framework of the global climate model of Laboratoire de Météorologie Dynamique LMDZ. LMDZ is the atmospheric component of the IPSL coupled atmosphere-ocean model, IPSL-CM, used in particular for the CMIP exercises. The model is a flexible tool which includes, besides the full 3-D coupled climate configurations, many configurations like a single column model (SCM) version extensively used for parameterization development or nudged configurations in which the synoptic situation is constrained. Nudged simulations were used extensively to demonstrate, by comparison with day-to-day site observations, the improvement of the representation of the meterology above continental surfaces induced by the boundary layer parameterizations presented above Diallo et al., 2017). All those configurations use exactly the same package of Fortran sources, an essential point to make a full use of the combination of SCM and 3-D simulations, as done in the present study.
IPSL participated in the fifth version of the CMIP exercise, CMIP5, with two versions of the model 5A and 5B which differ by the way turbulence, convection, and clouds are parameterized. The 5B version  capitalizes over 10-year research development on physics parameterization in the LMDZ team and on an extensive use of LES/SCM comparisons. It was the first climate version in which the thermal plume model described above was used to represent the vertical transport by the organized structures of the boundary layer convection, as well as a new parameterization of the pools of cold air created below cumulonimbus by reevaporation of convective rainfall (Grandpeix & Lafore, 2010;Rio et al., 2009). The thermal plume model, coupled to the bi-Gaussian cloud scheme, improves the representation of cumulus The layer thickness is shown on the x axis against altitude on the y axis for the whole atmospheric column (left) and focused on the first 3 km (right). clouds, over both the continents and oceans. However, as already explained in section 1, a trick had to be introduced to deactivate thermal plumes specifically over the regions of stratocumulus clouds due to a too strong mixing across the inversion as shown later in this paper ).
The LMDZ6A model used here is an improved version built on LMDZ5B which was developed recently for CMIP6. This 6A version includes improvements in many aspects of the model, which will be described in a forthcoming publication. In particular, the modification of the detrainment presented here, which satisfactorily reproduces the stratocumulus and cumulus to stratocumulus transition, is included in this new model. The 6A version is also the result of a long phase of tuning of the radiative balance for the CMIP6 exercise.
When defining a 3-D configuration of a climate model, compromises must be made between resolution and model accuracy. The vertical resolution is critical when designing a model configuration for long-term coupled climate simulations as those coordinated in the frame of the CMIP exercises. The 79-layer vertical discretization (L79) retained for the CMIP6 configuration of the IPSL climate model was designed in part to improve the representation of the stratospheric circulation, with an upper layer located at about 80 km above surface. The resolution was also improved in the first kilometers compared to the previous CMIP5 L39 configuration in order to better represent the boundary layer transport and associated clouds. In the L79 configuration, the first layer is centered at 10 m above surface. In the first 3 km above surface, the layer thickness z varies almost linearly with altitude with z ≃ 0.11z, so that 25 model layers are dedicated to the first 2 km.
In order to distinguish between formulation and resolution effects, two other vertical grids are used in the present study. The first one, based on 130 layers, essentially covers the first 4 km of the atmosphere. With 71 layers below 2 km, this resolution is comparable to that of LES (20 m). This version was designed specifically for LES/SCM comparisons and is used in particular in section 3 to present the scheme modification. An intermediate L95 grid was designed for both SCM and 3-D climate simulations starting from the L79 standard grid and adding layers in the first 5 km only. In the first 3 km, z ≃ 0.05z leading to 45 layers below 2 km. The three grids are shown in Figure 1. The SCM simulations are run with a time step of 10 min for the three-grid configurations. The 3-D global climate model (GCM) is run with a time step of 15 min for the boundary layer and clouds parameterizations and 1 hr for radiation. Figure 2 compares the LES simulation for the ARM cumulus, the FIRE stratocumulus, and the REF transition cases with the LMDZ SCM simulations performed using the thermal plume model combined with the bi-Gaussian cloud scheme. The version used on the graphics corresponds to the LMDZ6A version, except for the modification of the detrainment formulation which is intentionally deactivated on this first illustration. As expected, the ARM simulation is satisfactorily well reproduced with LMDZ, since the cloud scheme was tuned using this case . Nevertheless, the thermal plume model coupled with the bi-Gaussian cloud scheme is unable to reproduce the uniform 100% cloud cover observed in the FIRE stratocumulus case. The behavior is even worse for the transition case. During the first 3 hr, a significant cloud cover is simulated. It is, however, already underestimated and decreases before reaching the cloud top seen in LES. After the first 3 hr, the cloud cover strongly decreases reaching values similar to those for marine cumulus. The analysis below mainly focuses on this more challenging transition case.

Initial Formulation: Plumes Overshooting Too Far Above Inversion
The statistical cloud scheme could be questioned for the stratocumulus clouds which often show a trimodal rather than bimodal distribution of subgrid-scale water distribution, the third mode being associated with negatively buoyant downdrafts initiated at cloud top by mixing of cloud air with the very dry tropospheric air aloft. However, as discussed in previous studies Lewellen & Yoh, 1993;Perraud et al., 2011), using a bi-Gaussian distribution of the saturation deficit s leads to a satisfactory representation of clouds except for a slightly underestimated fraction in the lower part of the cloud layer. This is confirmed in Figure 3 in which we compare the effective cloudiness of the LES (contours, fraction of points with condensed water larger than 0.001 g/kg) with the cloudiness obtained applying the bi-Gaussian scheme in a diagnostic mode (color shading) as follows: By sampling thermal plumes in the LES, the mean humidity and temperature and hence the saturation deficit are computed separately for the thermal plumes s th and their environment s env , as well as the thermal plume vertical velocity w th and fractional cover , from which s,th and s,env are computed with equations (5) and (6) and hence the cloud fraction. This illustrates that the bi-Gaussian distribution, whose modes are attributed to thermals and their environment, leads to a satisfactorily representation of cloud fraction, both for cumulus and stratocumulus, without having to consider a third mode associated with downdrafts.
How do we explain, then, that the statistical parameterization, which fits the LES clouds when used in a diagnostic mode in LES, does not work when implemented in a single column model?
In fact, this erroneous behavior of the original thermal plume model for stratocumulus is caused by an overestimated top entrainment of dry air from the free troposphere as illustrated in Figure 4. The (idealized) initial water vapor profile (thin curve in Figure 4a) shows a very sharp inversion at the boundary layer top at 940 m above surface, typical of stratocumulus environment. After only 2 hr, the inversion rises by about 40 m in the SCM simulation. While the vertical humidity change in the first half of the boundary layer above surface is reasonably well captured by the SCM, the excessive moistening above the inversion is compensated by a drying by about 1 g/kg in the cloud layer. This 10% reduction decreases accordingly the relative humidity and leads to the rapid disappearance of the cloud deck (Figure 4b). In the model, it is the mass flux parameterization which is responsible for this exchange at the inversion level, the vertical diffusion being negligible at this level at the beginning of the simulation. Note that Zhang et al. (2013), in an intercomparison study of various boundary layer schemes on the CGILS (CFMIP/GASS Intercomparison of Large eddy and Single column models) case studies of cumulus and stratocumulus clouds, found that the models with active shallow convection schemes also tend to systematically reduce the cloud cover by stratocumulus clouds, resulting in a positive rather than negative radiative feedback for stratocumulus conditions. In our  In the upper part of the cloud, the lateral entrainment to the plume is almost 0. The vertical evolution of the thermal mass flux is thus controlled by lateral detrainment only ( f∕ z ≃ − f ). Since the mass flux is somewhat underestimated in the middle of the stratocumulus, the overestimated overshoot above cloud is most probably due to an insufficient reduction of the mass flux, due to an underestimated detrainment below cloud top. We show in Figure 4d the time evolution of the vertical profile of the detrainment rate superimposed with the cloud fraction in the LES. Here, the detrainment is computed from the thermal plume sampling, using the conservation equation where is deduced by inverting the conservation equation (1) for total water = q t (with S = 0), after decomposition of the derivative of f and use of the continuity equation leading to All along the simulation, the SCM results (Figure 4e), like the LES, show two maxima of the detrainment rate, one in the bottom of the cloud layer where the plumes go through a layer of negative buoyancy before acquiring additional buoyancy due to water condensation and the second at the top of the thermal plumes, starting typically 200 m below cloud top in the LES. During the first 3 hr, consistently with the mass flux profiles shown in Figure 4, the detrainment increases too close to the cloud top and extends too far above it. This leads to an overestimated moistening of the first layers just above cloud top by the thermal plume parameterization and, by compensation, to a drying of the upper part of the cloud layer, leading to the rapid disappearance of the cloud deck. As a consequence of the disappearance of the cloud deck, the radiative cooling at boundary layer top in the SCM is greatly reduced and so is the convection. Finally, the growth of the boundary layer depth is about twice slower than in the LES. Note that despite the bad behavior of the parameterizations, the model captures a diurnal cycle in phase with the LES.

The Modified Detrainment Formulation
The misrepresentation of detrainment may come from the fact that the thermal plume model does not account for the existence of downdrafts that carry dry and negatively buoyant air from the top of the boundary layer into the cloud (top entrainment or downdrafts). This process may contribute to increase lateral detrainment by mixing moist parcels from the plume with drier parcels coming from the free troposphere. Accounting more explicitly for those dry downdrafts would also, however, probably contribute to a faster drying of the cloud layer.
Another aspect to consider is that the plumes overshoot over the cloud top before spreading horizontally and downward. A large part of the overshooted air thus comes back below the averaged cloud top, in a reversible

at Hours 3 (daytime) and 15 (nighttime) for the REF test case and
Hour 10 for the ARM cumulus case, as computed directly from the LES sampling using equations (7) and (8) or with the modified parameterization (equation (4), the buoyancy B being computed from equation (9)) applied to the values of temperature and humidity sampled from the LES. The y axis is the altitude z divided by cloud top height Z top . The right column corresponds to the vertical profile of potential temperature and cloud fraction at the same time in the LES. LES = large-eddy simulation.
transformation. Thus, the deceleration of the plume over the cloud discontinuity is sensitive to the environment above the discontinuity but the air is effectively detrained for a large part below the inversion. Note that analyzing air parcel trajectories in LES simulations of cumulus clouds, Heus et al. (2008) ended up with the conclusion that such a detrainment through overshoot is at work all along the frontier of cumulus clouds, as illustrated in their Figure 16. In fact, both the overshoot effect and mixing with dry air may contribute to this tendency of the rising plumes to give their air back to the environment at a somewhat lower level than the level at which air parcels are extracted from the plumes.
This consideration may thus explain why the effective detrainment is aware of the environmental conditions at a level above.
Based on such considerations, we tested a very simple modification of the scheme proposed by Jam (2012) to make the detrainment aware of the environmental conditions at a level above. It consists in computing the buoyancy from the difference of the virtual potential temperature within the thermals at an altitude z with the potential temperature in the environment at a somewhat higher level z + z, so that buoyancy reads In the formulation retained here, we assume that z is scaling with z. The results of Heus et al. (2008; Figure  14c of their studies) suggest typically z ≃ 0.1z. In the following, we use z = Az, A being considered as a new adjustable parameter. Note that B ′ has been used in the formulation for the entrainment as well, although there is probably not much justification for it. Its effect on entrainment is weak at inversion because entrainment is negligible or null at this level and below inversion due to the well-mixed nature of the boundary layer. This point may be reconsidered in the future.

Evaluation of the Modified Detrainment in LES
Before showing SCM and 3-D tests with the modified detrainment in the next section, we first compare the detrainment obtained directly by LES sampling from equations (8) and (7) (black dots in Figure 5, same values as shown in Figure 4d), with the new detrainment formulation (equation (4), the buoyancy B being computed from equation (9) The parameterization generally underestimates the detrainment in the subcloud layer; this point will deserve further investigation and is not analyzed here. At Hour 3 (first row), at noon local time, we clearly see the two maxima in the detrainment. The lower one, at 0.6-0.75Z top , is well captured by the parameterization, but the second one at 0.9 hr is located much too close to the cloud top. The situation is even worse during the first night (second row). In both cases, the modified parameterization leads to a detrainment estimation active in the whole cloud layer. The best agreement with the LES sampling is obtained for A = 0.05 (while the other tests strongly overestimate the detrainment) for Hour 3 and A = 0.10 for Hour 15, even if the modified parameterization tends to overestimate the detrainment in the upper part of the cloud.
Note that the modified parameterization seems to have an even more positive impact on the representation of the detrainment for the ARM cumulus case as illustrated in the last row of Figure 5 for Hour 10 of the simulation. Further investigations should be conducted, however, to confirm that it is the signature of an improvement of the physics of the scheme rather than an error compensation.

Best Value of the A Parameter for the REF Transition Case
In Figure 6, we explore the sensitivity to the A parameter value for the REF transition case. Compared to the LES simulation (top panel), the simulation with the unmodified detrainment parameterization (second panel) completely misses the stratocumulus deck (the first two panels are identical to the right panels in Figure 2). At the opposite, if the thermal plume model is deactivated (in which case all the boundary layer transport is done by the Mellor and Yamada turbulent scheme as was the case in the LMDZ5B version of the model in the regions of stratocumulus), the cloud deck is strong but the increase in cloud top height is not captured and the cloud evolves into a stratus lowering fog. An almost similar behavior is obtained if too large a value of the A parameter is used. Finally, with a value of A = 0.07, the time evolution of the cloud cover is rather well captured by the model, including its diurnal modulation. There seems to be an optimal value of the A coefficient, since the cloud top height is decreasing both when increasing or decreasing it from its optimal value of 0.07. This value is the one which was selected for the CMIP6 configuration of the LMDZ model.
To illustrate how the modified detrainment enables the maintenance of the cloud deck, we show in Figures 7a and 7c, for the REF transition case and FIRE case respectively, and for the first 2 hr (12 model time steps), the vertical profile of the internal variables of the thermal plume model and eddy diffusion parameterizations responsible for the vertical transport. For both cases, the profiles are very similar for the A = 0 (blue, unmodified detrainment) on the one hand and A = 0.07 simulation (orange) on the other hand. The mass flux f (solid curves) overshoots about 50 m higher than the A = 0.07 simulation (80 m for FIRE). This overshoot induces a moistening above the cloud layer and a drying of the cloud layer by the compensating subsidence which takes its air in the very dry free troposphere (solid blue curve in Figures 7b and 7d). The effect of vertical diffusion is similar in both cases: It transports upward in the first hundred meters the Results are shown both for A = 0.00 (blue, without detrainment modification) and A = 0.07 (orange). The vertical axis is the altitude normalized by the top cloud altitude defined with a 50% threshold on cloud cover. The light gray and gray shaded areas correspond to the altitude range of clouds, corresponding to a cloud cover larger than 2% and 50%, respectively, as obtained in the A = 0.07 simulations.
air evaporated at surface (then transported to the cloud layer by thermals); within the cloud layer, in the A = 0.07 simulation, eddy diffusion contributes to a slight drying by mixing with the dry tropospheric air. Eddy diffusion is reactivated in the cloud because of the instability of the temperature profile, created by the cloud top radiative cooling to space. Because of the disappearance of the cloud for A = 0 in the REF case, the eddy diffusion disappears as well. Note that for similar reasons as the one mentioned for plume detrainment, it is a critical point (not illustrated here) for the maintenance of a cloud deck to use a parameterization for eddy diffusivity that drops to 0 at the inversion in order to limit mixing between the cloud layer and the dry free troposphere. It is the case of the Yamada (1983) scheme used here.

Boundary Layer Top Entrainment
In order to interpret the results presented in relation with other modeling studies, a direct link can be made between the mass flux model and classical "entrainment parameterizations" or "entrainment closures" (see, e.g., Lilly, 2002;Lock et al., 2000).
The numerical computation of the flux of a scalar quantity (F = w ′ ′ ) by the thermal plume model is based on upstream schemes and reads where the mass flux F k+1/2 and scalar flux F k+1/2 are computed at intermediate level k + 1∕2, between layer k where the air rising in the plume is coming from, with scalar concentration th,k , and layer k + 1 from which the air is advected downward in the environment, with scalar concentration k+1 .
At the inversion level, th,k has typical values of the cloud layer and k+1 of the free troposphere so that the flux may be expressed approximately as w ′ ′ t ≃ − w e Δ , being the jump in scalar at the inversion and w e = f∕ . Note that, since ≪ 1 in the cloud, w e is in fact almost the vertical velocity in the compensating subsidence w subs (counted positively downward), related to the mass flux by f = w th = (1 − )w subs . Figure 8a shows, for Hour 13 of the REF transition simulation, the averaged vertical flux of liquid potential temperature F l , compared to the profile of l . In practice, the maximum of the flux can be somewhat shifted downward from the maximum gradient of l . However, an estimate of this top-entrainment velocity can be obtained by dividing the maximum value of the thermal plume model flux by the maximum value of the vertical difference between two adjacent grid cells 1 l = l,k+1 − l,k . Such an estimation is compared in Figure 8b (red curve) with the entrainment velocity diagnosed as a residual between the deepening of the boundary layer and the vertical large-scale subsidence w e = Z top ∕ z + W lsc (where the large-scale subsidence velocity W lsc is counted positively downward). This diagnosed entrainment velocity is satisfactorily well reproduced in the SCM (black dashed curve) compared to the LES (black solid curve). The estimate (a) Vertical profile of the liquid potential temperature turbulent transport F l by the thermal plumes (black curve) and by eddy diffusion (red) is compared with the vertical profile of l (orange), showing the collocation of the maximum thermal plume downward transport with the maximum vertical variation of l ; (b) boundary layer top-entrainment velocity as diagnosed from the residual of the deepening of the boundary layer and imposed subsidence in the LES (black solid curve) and SCM simulation (black dashed curve). The red curves show the ratio of the maximum downward flux of liquid potential temperature to the maximum temperature variations between two adjacent layers ( 1 l = l,k+1 − l,k ) apart from a flux interface level (k + 1∕2) or between the next layers above and below ( 2 l = l,k+2 − l,k−1 ). The green curves correspond to the downward velocity f∕ in the three upper levels where the thermal flux model is active. SCM = single column model; LES = large-eddy simulation.
from the ratio of the maximum value of F l to the maximum value of 1 l (red curve) shows somewhat larger values, due to the fact that the jump in temperature is done in practice in more than one layer. If computing the temperature jump over a larger distance as 2 l = l,k+2 − l,k−1 (about 80 m), the agreement with the diagnosed entrainment velocity is much better. On the same graph, we superimpose the downward vertical velocity computed from the thermal plume mass flux ≃ f∕ in the upper three layers affected by the thermal plume model. At the beginning, the entrainment velocity is close to the velocity in the top layer. Progressively during the simulation, the effective entrainment velocity becomes closer to downward velocities in layers farther from the thermal plume top, since thermal plumes overshoot farther above the inversion.
Since the rate of lateral entrainment to the plume is 0 in the upper part of the cloud, the specification of the detrainment rate in the upper levels directly controls the vertical mass flux ( f∕ z = − f) and in turn the effective top entrainment as demonstrated above. So the requirement that the plumes do not detrain too far above inversion can be reworded as a requirement for the top entrainment by compensating subsidence not to be active too far above the inversion.  The top entrainment simulated with the thermal plume model is clearly able to capture the strong diurnal cycle observed in LES. This is illustrated further in Figure 9a which shows the time evolution of the mass flux (black curve) and eddy diffusivity (red curve) averaged over the depth of the cloud layer. The mass flux shows a strong diurnal cycle with a maximum during nighttime. During daytime, the additional heating by solar radiation within the cloud reduces the convective mass flux and, subsequently, the boundary layer top entrainment. This modulation is controlled in the thermal plume model by the plume fractional cover as seen both from the green dashed curve in Figure 9a and from the vertical profiles of the mass flux, vertical velocity, and fractional cover in Figure 9b for the first night of the simulation (dark blue) and for the following daytime conditions (yellow). The reduction of the plume fraction during daytime is due to the larger detrainment rate at cloud base (due to the larger negative buoyancy of the plumes there during daytime). At the opposite, the local turbulence systematically reinforces during daytime, probably due to the local destabilization of thermal profile by the increased absorption of solar radiation below cloud top, while the cooling to space is confined to a very thin layer at cloud top. Finally, with a smaller convective flux to feed the cloud and a more intense in-cloud turbulence, the decoupling of the stratocumulus layer is larger during daytime.
Note that the averaged thermal plume velocity in the clouds increases continuously during the course of the simulation, following the evolution of the imposed sea surface temperature. The FIRE case shows a very similar evolution of the turbulent quantities (not shown here) and a thermal plume velocity that does not evolve, consistently with the constant SST forcing of this case.
Finally, we see that the thermal plume model provides a physical parameterization of the boundary layer top entrainment and transport that compares reasonably well with LES and is sensitive both to radiation and surface instability.

Sensitivity to Forcing and Vertical Resolution
The simulation of the cloud cover is illustrated as well for the FAST and SLOW transition cases in Figure 10. If the scheme fails at dissipating the cloud top deck at the end of the FAST simulation, the sensitivity to the various forcing of cloud height as well as the diurnal cycle are reasonably well captured. Figure 11 compares, for the REF transition case, the time evolution of the cloud fraction vertical profile with the L130 (top panels) L95 (middle) and L79 (bottom) vertical grids. Not surprisingly, degrading the resolution has a negative effect on the representation of the stratocumulus. For the L79 simulation in particular, some periods appear in which the cloud deck disappears during several hours. However, the qualitative behavior is still captured with this version, despite vertical resolutions of typically 200 m at the altitude of the cloud deck, which is typically the thickness of the cloud deck. The collocation of the cloud top and detrainment maximum is preserved even with the coarsest grid.
To confirm this, we compare results obtained with the three vertical resolutions available for the three 1-D cases, using an effective height of the cloud deck computed by weighting the altitude of the cloud top by the fourth power of the cloud fraction (cf. the caption of Figure 12). On the left-hand side of Figure 12, the time evolution is shown for the optimal value A = 0.07. This confirms the good behavior with a L130 discretization, when comparing the effective height from the LES (thick curves) and that of the SCM (full circles). The cloud effective height tends to be underestimated at lower resolution compared to LES (thick Figure 11. Sensitivity test to the vertical resolution of the simulated cloud cover (left column; %) and detrainment rate (right; m −1 ) for the REF transition case with the L130 (top, already shown), L95 (top), and L79 (bottom) vertical grids. The L79 configuration is that retained for CMIP6. The 2%, 5%, 10%, 20%, and 50% (black) isolines superimposed on the color shades correspond to LES cloud fraction for the left and SCM cloud fraction for the right column. LES = large-eddy simulation; SCM = single column model. curves). However, in most cases, the difference is of the order of the thickness of model layers at most, which can be seen if computing the effective height by shifting the cloud fraction value one layer above ("+" signs). In the SLOW case, however, the error can exceed 2 z. The dissipation of the cloud deck at the end of the FAST case already mentioned is seen through the sudden decrease of the effective height on the last 3 hr of the simulation. Note also that for the L79 grid, the cloud deck tends to be more intermittent as already seen in Figure 11, explaining that the effective height suddenly drops in the simulation. It is noteworthy, however, that the sensitivity of the effective height to the value of the parameter A seems quite consistent across resolutions, as seen by comparing the value of the effective height in the middle of the third day, for the three cases and the various values of A (right column in Figure 12). Results are averaged between 04 and 08 a.m. on July 18 to avoid the disappearance of the cloud deck. A value between 0.05 and 0.10 seems relevant in most cases.
We finally compare the vertical profiles of cloud cover for three other cases of cumulus and stratocumulus ( Figure 13). Results are shown for two resolutions either without modification of the detrainment formulation (A = 0) or with the value retained in the standard LMDZ6A version (A = 0.07). With A = 0, the stratocumulus deck is maintained for the FIRE case, opposite to what is obtained for the transition case. The modification of the detrainment formulation, with A = 0.07, improves, however, the representation of the stratocumulus FIRE case, with a more stable cloud top and better representation of the dirunal cycle. Note, however, that whatever the resolution or tuning, the SCM strongly underestimates the cloud fractional cover above cloud base for the cumulus cases over ocean (RICO). This point would probably deserve further investigation. But the good representation of the maximum cloud cover at cloud base and height of this cloud base is nevertheless reasonably well simulated. It is noticeable that the results are a bit closer to LES for RICO for the L79 grid with the modified detrainment. For the ARM continental cumulus case as well, the same configuration is the closest to LES. This result may be due in part to the fact that the L79 version with A = 0.07 is the version for which the cloud free parameters were tuned for the reference LMDZ6A configuration. This tuning in particular includes, besides the detrainment/entrainment formulation or parameters controlling the bi-Gaussian scheme, parameters that control the cloud microphysics (in particular the reevaporation of condensed water which has a significant impact on cloud profiles). This point will be documented in an other study.

Three-Dimensional GCM Simulations
Three-dimensional tests are now presented for the L79 and L95 versions with the horizontal grid of the configuration retained for IPSL-CM6, based on 144 by 142 points equally distributed in longitude and latitude, respectively. The simulations are run over 4 years starting from 1 January, with an imposed mean seasonal Analysis concerns averages over the last 3 years of the simulations. Figure 14 shows vertical cross sections of cloud cover in longitude average between 25 • S and 15 • S. The cross sections covers the Pacific (on the left of the graphs), South America, and South Atlantic ocean (on the right). The transition from low stratocumulus on the east side of the two oceanic basins to deeper cumulus clouds is clearly visible in this cross section in the observations (upper panel of the left column). The behavior of the 3-D model when varying A is consistent with what is seen in SCM experiments. For A = 0, the stratocumulus essentially disappear. When increasing too much its value, the stratocumulus are getting too thick. Again, it seems that A = 0.07 is an optimum value for the A parameter. The results are essentially similar between the two configurations despite a resolution almost twice as high in the L95 than in the L79 vertical grid.
The major deficiency of the modified scheme with A = 0.07 is probably its tendency to create two distinct maximum of cloud fraction in regions of trade winds cumulus, the first one at 700 m above surface as in the observations, and the second one near 2 km. This two-maximum vertical structure is visible as well in the observations, but much less marked and concerns a much shorter transition region in longitude, near 135 • W over the Pacific ocean and 35 • W over the Atlantic. The vertical structure is upside down in the simulation westward of this transition region, with a cloud cover strongly underestimated below the upper maximum. The tendency of the modified scheme to amplify the cloud cover at 2 km is in fact visible already in the RICO SCM simulation shown in Figure 13. This deficiency may be related as well to the difficulty to dissipate the upper stratocumulus cloud at the end of the FAST transition case (Figure 10). It may be due to the microphysical scheme which would produce insufficient precipitation from the upper to the lower levels in the cloud or to some missing physics at the inversion.
In Figure 15 we show maps of the annual mean coverage by low clouds. The low clouds are compared to the GOCCP climatology constructed from the Calipso Lidar observations (upper panel on the left). The model results are postprocessed running the Calipso part of the COSP satellite simulator to be directly comparable with observations (Bodas-Salcedo et al., 2011;Chepfer et al., 2008). It means, however, that part of the difference between the observed and simulated low clouds may come from differences in their overlap by high clouds. In all the simulations, the coverage by cumulus clouds is somewhat underestimated, a deficiency which could be partly overcome by accounting for vertical cloud heterogeneities (Jouhaud et al., 2018). The cloud cover over the Pacific Intertropical Convergence Zone is too strong over the east side of the oceanic basin and too weak over the warm pool. This bias is clearly correlated to the rainfall bias of the LMDZ6A version (not shown). The comparison of the extension of the zone covered by clouds at more than 50% is best fitted for A = 0.07. The root-mean-square error (shown in the top of each panel) reaches a minimum value of ≃11% for A = 0.07 and A = 0.1, to be compared with 22% for the A = 0 case and 14% for A = 0.15. The mean bias also goes through 0 between A = 0.07 and A = 0.1.
On the right-hand side, the reflected SW radiation at top of atmosphere is compared with Ceres-Ebaf observations (Loeb et al., 2009). The amplitude of contrasts simulated when A = 0.07 is consistent with observations. However, the maximum of reflected radiation in the stratocumulus region is shifted east, a bit too far from the coast. This remaining deficiency may be related to the tendency of the scheme to maintain a stratocumulus 100% covering deck at the end of the FAST transition case, which may explain why the reflected radiation is too large far from the coast. The underestimation of the radiative forcing close to the coast is probably due to an underestimation of the liquid water path, since it corresponds to a region of 100% coverage.
The minimum value of the root-mean-square error is reached for A = 0.05 and A = 0.07 for the reflected radiation, with a value of 13 W/m 2 to be compared with 15 W/m 2 for A = 0 and 20 W/m 2 for A = 0.15. All Figure 15. Maps of low-cloud cover (%; left, clouds with cloud top below 680 hPa) and reflected SW radiation at top of atmosphere (right; W/m 2 ). Observations (top panels, GOCCP for clouds and Ceres-Ebaf for radiation) and simulations (with COSP simulator for the left column) with the L79 vertical grid. The mean biais, root-mean-square error (RMSE), and correlation (CORR) with observations is displayed at the top of each panel. GOCCP = GCM Oriented Calipso Cloud Product. together, the mean bias and root-mean-square error are at minimum values for A = 0.07. The correlation is also over 80% for A = 0.07 both for cloud cover and reflected radiation, while it is of 50% only for cloud cover and 70% for the reflected radiation for A = 0. If the reduction is not spectacular, it is significant. Moreover, it is obtained thanks to an improvement of the model physical content; the various cloud boundary layer regimes are now handled without switching from one scheme to another depending on the environmental conditions.
The previous considerations are confirmed in Figure 16 that compares the low level clouds and reflected SW radiation averaged between 25 • S and 15 • S, for both the L95 and L79 vertical resolutions. Once again, the amplitude of the longitudinal variations are well simulated with A = 0.07. While the longitude structure is well phased for the cloud cover, the radiative forcing is shifted westward.

Conclusion
A simple modification of the detrainment formulation in the thermal plume model proposed by Jam (2012) allows a reasonable simulation of both cumulus and stratocumulus clouds. The modification consists in computing the buoyancy of the plume that enters in the formulation of detrainment by comparing the buoyancy of a parcel with respect to the environmental air typically 10% above. Whether this formulation should be justified by some overshoot of air parcels when detraining from the updrafts or by the role of subsiding shells or coherent downdrafts, what seems important is to make the thermal plume start to detrain air laterally before reaching cloud top, or detrain efficiently enough at cloud top.
With this modification, the parameterizations of LMDZ represent boundary layer transport and clouds in a unified way, without switching from one parameterization to the other depending on the environmental conditions. The thermal plume model seems to provide a reasonable parameterization of the cloud top entrainment in stratocumulus clouds and its modulation by surface fluxes and radiation. Due to the coarse vertical resolution of climate models, papers dealing with local explicit parameterization of top entrainment like Grenier et al 2001or Lock et al. 2000 underline the need for a careful treatment of the local longwave radiation divergence at cloud top. It does not seem to be that an issue with our model. Indeed, as soon as the clouds cover 100% of the surface and the optical thickness is large enough, the cooling to space at cloud top which is just the black body emission will be well represented whatever the resolution of the model. It is this cooling to space which is a first-order driver of the stratocumulus-topped boundary layer convection together with surface fluxes. Radiative flux divergence of course also affects locally the stability within the cloud layer, which may in turn influence the cloud physics, and it is what happens indeed for the turbulent kinetic energy in our model. The good representation of the top entrainment is directly related to the representation of the lateral detrainment from the thermal plume which controls the shape of the mass flux in the upper layers and, in turn, the velocity in the compensating subsidence. We have demonstrated that the modified thermal plume model indirectly provides a reasonable parameterization of top entrainment. Not only the cloud cover is reasonably well reproduce but also the deepening of the boundary layer, a key for the control of surface fluxes and sea surface temperature biases (Hourdin et al., 2015). Note also the remarkable consistency between SCM and 3-D results. Both in particular suggest an optimal value between 0.05 and 0.1 for the A coefficient.
The study suggests that the L79 grid designed for the LMDZ6A standard configuration, with a layer thickness of 110 m at 1 km above surface, is too coarse in the boundary layer, which is not really a surprise in view of the real cloud thickness and importance of the discontinuity at cloud top. The proposed formulation appears, however, to be quite robust to vertical resolution, and to work satisfactory well for both cumulus and stratocumulus. The study thus reinforces the relevance of the use of mass flux parameterizations of organized structures for a good representation of the convective boundary layer. It also confirms once again the efficiency of the LES/SCM framework to improve and adjust new parameterizations for the boundary layer in global models.
The multicase agreement between LES and SCM is, however, obviously not perfect, nor is the comparison of 3-D GCM with Calipso observations. This may come both from aspects of the formulations in the parameterization that would deserve further improvement or from the choice of model free parameters. The version of the model presented here is a result of a long phase of tuning, involving iterations between 1-D test cases and 3-D climate simulations, and results from compromises between the simulations of various metrics which were chosen as targets for the tuning. It may not be by chance that the results of the reference configuration, for A = 0.07, behaves better for the coarser grid, L79, for the cumulus cases (Figure 13), since it is the grid with which the tuning was done.
The new set of parameterizations seems to capture a large part of the processes necessary to represent the cloud cover and associated radiation in a climate model. However, significant biases remain and shortcomings persist: an underestimation of cloud cover above cloud base for trade winds cumulus, an upside-down vertical structure of the cloud cover in this region, with a maximum at 2 km above surface, an underestimation of stratocumulus, and albedo close to the coast. These deficiencies may be due in part to processes not accounted for in the current parameterization, such as the negatively buoyant downdraft forced by evaporative cooling of air parcel at the upper bound of stratocumulus clouds. Also, the current microphysics parameterization in LMDZ is very crude. All those aspects would deserve further investigation. In order to reach effective improvements, those studies will require a more detailed investigation of LES cases and more systematic exploration of the sensitivity to free parameters using automatic tools, such as those developed in the Uncertainty Quantification community. A more physical parameterization of the z value in the computation of B ′ could probably be developed as well.