Atmospheric Convection

ABSTRACT

Convective parameterization is the long-lasting bottleneck of global climate modelling and one of the most difficult problems in atmospheric sciences. Uncertainty in convective parameterization is the leading cause of the widespread climate sensitivity in IPCC global warming projections. This paper reviews the observations and parameterizations of atmospheric convection with emphasis on the cloud structure, bulk effects, and closure assumption. The representative state-of-the-art convection schemes are presented, including the ECMWF convection scheme, the Grell scheme used in NCEP model and WRF model, the Zhang-MacFarlane scheme used in NCAR and DOE models, and parameterizations of shallow moist convection. The observed convection has self-suppression mechanisms caused by entrainment in convective updrafts, surface cold pool generated by unsaturated convective downdrafts, and warm and dry lower troposphere created by mesoscale downdrafts. The post-convection environment is often characterized by “diamond sounding” suggesting an over-stabilization rather than barely returning to neutral state. Then the pre-convection environment is characterized by slow moistening of lower troposphere triggered by surface moisture convergence and other mechanisms. The over-stabilization and slow moistening make the convection events episodic and decouple the middle/upper troposphere from the boundary layer, making the state-type quasi-equilibrium hypothesis invalid. Right now, unsaturated convective downdrafts and especially mesoscale downdrafts are missing in most convection schemes, while some schemes are using undiluted convective updrafts, all of which favour easily turned-on convection linked to double-ITCZ (inter-tropical convergence zone), overly weak MJO (Madden-Julian Oscillation) and precocious diurnal precipitation maximum. We propose a new strategy for convection scheme development using reanalysis-driven model experiments such as the assimilation runs in weather prediction centres and the decadal prediction runs in climate modelling centres, aided by satellite simulators evaluating key characteristics such as the lifecycle of convective cloud-top distribution and stratiform precipitation fraction.

RESUME

[Traduit par la redaction] La paramétrisation convective est le goulot d’étranglement durable de la modélisation du climat mondial et l’un des problèmes les plus difficiles des sciences de l’atmosphère. L’incertitude dans la paramétrisation de la convection est la principale cause de la sensibilité climatique étendue dans les projections de réchauffement global du GIEC. Le présent article porte sur les observations et les paramétrisations de la convection atmosphérique en mettant l’accent sur la structure des nuages, les effets de masse et l’hypothèse de fermeture. Les schémas de convection représentatifs de l’état de la technique sont présentés, y compris le schéma de convection du CEPMMT, le schéma Grell utilisé dans le modèle NCEP et le modèle WRF, le schéma Zhang-MacFarlane utilisé dans les modèles NCAR et DOE, et les paramétrisations de la convection humide peu profonde. La convection observée possède des mécanismes d’auto-suppression causés par l’entraînement dans les courants ascendants de convection, le bassin froid de surface généré par les courants descendants de convection non saturés, et la basse troposphère chaude et sèche créée par les courants descendants de mésoéchelle. L’environnement post-convection est souvent caractérisé par un “sondage en diamant” qui suggère une sur-stabilisation plutôt qu’un retour à l’état neutre. Ensuite, l’environnement pré-convection est caractérisé par une lente humidification de la basse troposphère déclenchée par la convergence de l’humidité de surface et d’autres mécanismes. La sur-stabilisation et la lenteur de l’humidification rendent les événements de convection épisodiques et découplent la moyenne/supérieure troposphère de la couche limite, ce qui rend l’hypothèse de quasi-équilibre de type état invalide. À l’heure actuelle, les courants convectifs descendants non saturés et surtout les courants descendants à méso-échelle sont absents de la plupart des schémas de convection, tandis que certains schémas utilisent des courants convectifs ascendants non dilués, qui favorisent tous une convection facilement activée liée à une double ZCIT (zone de convergence intertropicale), une OMJ (oscillation Madden-Julian) trop faible et un maximum de précipitations diurnes précoce. Nous proposons une nouvelle stratégie pour l’élaboration de schémas de convection à l’aide d’expériences de modèles pilotées par des réanalyses, telles que les séries d’assimilation dans les centres de prévision météorologique et les séries de prévisions décennales dans les centres de modélisation climatique, assistées par des simulateurs par satellite évaluant des caractéristiques clés telles que le cycle de vie de la distribution des sommets des nuages convectifs et la fraction des précipitations stratiformes.

Keywords:

• atmospheric convection
• convective parameterization
• global climate modelling
• weather prediction

1 Introduction

Atmospheric convection is the vertical movement of buoyant air parcels, often called updrafts or downdrafts, associated with thermals, clouds, thunderstorms, and mesoscale cloud systems (Fig. 1). Atmospheric convection is a fast non-local transport of mass, heat, water, momentum and vorticity, which is often associated with phase change of water and resultant release/consumption of heat. Convective heating directly drives large-scale atmospheric circulations such as the Hadley Circulation (Hadley, 1735; Simpson et al., 1988), Walker Circulation (Walker, 1923), ENSO circulation (Bjerknes, 1969) and MJO circulation (Madden & Julian, 1971), as well as extreme weather systems such as the tropical cyclones (Riehl, 1950). Atmospheric convection is also closely connected to cloud-radiation feedback (Lin et al., 2014; Slingo, 1990), surface flux feedback (Lin, 2007; Wallace, 1992), and chemical transport (Gidel, 1983; Thompson et al., 1997) and plays an important role in global climate change (Bony et al., 2015). Atmospheric convection is the leading factor controlling the climate sensitivity of climate models and can explain half of the variance among the more than 40 IPCC models (Sanderson et al., 2010; Sherwood et al., 2014; Stainforth et al., 2005; Zhao, 2014).

Fig. 1 Photos of cloud systems. (A) A mesoscale convective system over tropical continent. (B) A thunderstorm over tropical ocean. (C) Fair weather cumulus clouds over land. (D) Stratocumulus clouds over ocean. [Courtesy of NASA].

Representation of atmospheric convection is one of the most difficult problems in global climate modelling. Since the grid sizes of global climate model are generally much larger than the convective updrafts and downdrafts, a subgrid-scale physical model is needed to describe the bulk effects of atmospheric convection, which is called convective parametrization or convection scheme. Theoretical studies of convection started from idealized convection (Benard, 1900; Rayleigh, 1916; Thomson, 1882), then advanced to realistic atmospheric convection (Bjerknes, 1938; Fujiwhara, 1939; Kuo, 1960; Lilly, 1960; Petterssen, 1939), and to more complicated convective parameterizations (Arakawa & Schubert, 1974; Kuo, 1965; Manabe et al., 1965).

In the past half century, numerous convection schemes have been developed. However, modelling atmospheric convection remains one of the key bottlenecks of climate modelling (Bony et al., 2015; Randall et al., 2003). The state-of-the-art climate models still have significant difficulty even in simulating the climatological mean surface precipitation, which equals the column-integrated latent heating (Fig. 2). The ensemble mean precipitation of 23 IPCC AR5 models show ±3 mm/day biases in tropical Pacific, Indian Ocean, Amazon and India, which is about 30% of the observed precipitation, together with a strong double-ITCZ pattern in the tropical Pacific (Huang et al., 2018, their Fig.1b; Fig. 2a). When forced by observed SST, the corresponding AGCMs show slightly smaller biases but still with significant double-ITCZ pattern (Lin, 2007). Experiments using super-parameterization also show similar magnitude of biases (Randall et al., 2016, their Fig.15-5c; Fig. 2b). The non-hydrostatic high-resolution global cloud-resolving models still have large biases (Kodama et al., 2015, their Fig.1c; Fig. 2c). The most striking result is for reanalysis, which is forced not only by observed SST, but also by a vast set of observed surface and upper air states. The precipitation biases are only partly reduced in the most recent reanalysis (Hersbach et al., 2020, their Fig.24; Fig. 2d). Overall, the convection schemes in global climate models are too easy to be turned on, leading to unrealistically frequent but weak precipitation and drizzles in the models in contrast with the episodic strong precipitation events in nature. The persistent weak precipitation suppresses variability such as the MJO (Lin et al., 2006; Hung et al., 2013; Fig. 3a), and causes precocious precipitation maximum in diurnal cycle over both land and ocean (Bechtold et al., 2004; Covey et al., 2016; Dai, 2006; Tang et al., 2021; Fig. 3b) and missing stratocumulus clouds over eastern parts of ocean basins (Lin et al., 2014).

Fig. 2 Biases of climatological mean precipitation with respect to GPCP/TRMM observations for (A) Ensemble mean of 23 CMIP5 global climate models (Huang et al., 2018). (B) An AGCM with super-parameterization of convection (Randall et al., 2016). (C) A global cloud resolving model (Kodama et al., 2015); and (D) ERA5 reanalysis (Hersbach et al., 2020).

Fig. 3 (A) MJO precipitation variance in CMIP3 models (left, from Lin et al., 2006) and CMIP5 models (right, from Hung et al., 2013). (B) Phase and amplitude of diurnal cycle of precipitation over land (left) and ocean (right) for CMIP5 AMIP models and TRMM observations (adapted from Covey et al., 2016).

Table 1 summarizes the convection schemes used in global and regional climate models grouped by the types of closure assumptions. The biases in global cloud-resolving models and models with super-parameterization suggest that atmospheric convection is not only a resolution problem. In the foreseeable future, supercomputers will not be able to conduct long-term (e.g. 100 years and longer) ensemble global cloud-resolving model runs, although such runs are basic requirements for understanding global climate and climate change. Therefore, improving convective parameterization is very important at the current stage and in the foreseeable future.

Table 1 Convection schemes used in global and regional models grouped by closure assumptions.

Closure	First-Generation Scheme	Later Schemes	Models	Other Models
Moisture Convergence	Kuo (1965, 1974)	Bougeault (1985) Frank and Cohen (1987) Tiedtke (1989) Grell (1993), Grell and Devenyi (2002) Freitas et al. (2021)	Centre National de Recherches Météorologiques, France PSU/NCAR Mesoscale Model European Centre for Medium-Range Weather Forecasts Max Planck Institute for Meteorology, Germany NCAR Weather Research and Forecasting Model NOAA National Center for Environmental Prediction NASA Global Modeling and Assimilation Office NOAA Global Systems Laboratory NCAR Weather Research and Forecasting Model NCAR Regional Climate Model	Irish Centre for High-End Computing
Flux-type CQE	Arakawa and Schubert (1974)	Moorthi and Suarez (1992) Donner (1993) Randall and Pan (1993) Grell (1993), Grell and Devenyi (2002) Freitas et al. (2021) Zhang and McFarlane (1995) Chikira and Sugiyama (2010) Wu (2012) Zhao et al. (2018)	Geophysical Fluid Dynamics Laboratory Geophysical Fluid Dynamics Laboratory Model for Interdisciplinary Research on Climate, Japan Meteorological Research Institute, Japan Same as above National Center for Atmospheric Research Canadian Centre for Climate Modelling and Analysis, Canada Model for Interdisciplinary Research on Climate, Japan Beijing Climate Center, China Geophysical Fluid Dynamics Laboratory	Bjerknes Centre for Climate Research, Norway Institute of Atmospheric Physics, China Academia Sinica, Taiwan
State-type CQE	Manabe et al. (1965)	Fritsch and Chappell (1980) Betts (1986) Kain and Fritsch (1990) Gregory and Rowntree (1990) Emanuel (1991) Del Genio and Yao (1993) Grell (1993), Grell and Devenyi (2002) Freitas et al. (2021) Yoshimura (2015)	PSU/NCAR Mesoscale Model Institute of Numerical Mathematics, Russia NCAR Weather Research and Forecasting Model Met Office Hadley Center, United Kingdom Institute Pierre-Simon Laplace, France NASA Goddard Institute for Space Studies Same as above Meteorological Research Institute, Japan	Commonwealth Scientific and Industrial Research Organization, Australia

The purpose of this paper is to review the observational constraints for convective parameterization and the most widely used convection schemes. Section 2 will review the observational studies of atmospheric convection. Section 3 will review the ECMWF convection scheme, which evolved from the Tiedtke scheme with moisture convergence closure. Section 4 reviews the Grell scheme with flux-type quasi-equilibrium closure, which is used in the National Center for Environmental Prediction (NCEP) GFS model, the NOAA GSL model, and the regional Weather Research and Forecast (WRF) model. Section 5 reviews the Zhang-McFarlane scheme with flux-type quasi-equilibrium closure, which is used in the National Center for Atmospheric Research (NCAR) and the US Department of Energy climate models. Section 6 reviews the parametrization of shallow moist convection. A summary of current challenges and suggested future directions is given in Section 7.

2 Observations of atmospheric convection

a Cloud Structure

Atmospheric convection can be divided into deep convection and shallow convection. Deep convective systems are thunderstorms, which include ordinary thunderstorms, multi-cell thunderstorms and supercell thunderstorms, and mesoscale convective systems (MCSs), which include squall lines and mesoscale convective complexes (MCCs). Shallow convection refers to cumulus clouds, which are also called fair weather cumulus or trade wind cumulus, and stratocumulus clouds.

Soon after the invention of basic modern meteorological instruments, such as rain gauge, anemometer, thermometer and barometer, scientists started to observe the surface structure of thunderstorms, especially sudden increase of wind speed and sudden drop of air temperature, which are likely downdrafts and cold pools (e.g. Planer, 1782; Rosenthal, 1786; Strehlke, 1830; Symons, 1890; Toaldo, 1794). Cold pool is a cold pocket of dense air that forms when rain evaporates during intense precipitation inside downdrafts underneath a thunderstorm cloud. Cold pools propagate away from the rain event along the surface as a moving gust front. When the gust front passes, cold pools cause a sudden increase in wind speed and a sudden drop in specific humidity and in air temperature. In 1857, G. J. Symons established an organization to study the English thunderstorms (Symons, 1889). In 1887, the Royal Meteorological Society set up a committee to study thunderstorms, and published two summary reports (Abercromby, 1888; Mareiott, 1890). Using 10 years (1925-1934) of data from surface station network, Ward (1936) constructed composite surface conditions for different types of thunderstorms. Four distinct types of pressure distribution were identified as giving rise to thunderstorms on the national forests of the Pacific Northwest. Listed according to their frequency and forecasting importance, they are: Type I. Trough; Type II. Cyclonic; Type III. Transition; Type IV. Border. With data from the earliest upper air sounding network, Neiburger (1941) analyzed the potential vorticity field of a thunderstorm.

MCSs, such as squall lines, were also discovered from surface measurements in the nineteenth century (e.g. Hinrichs, 1883, 1888a, 1888b; Ley, 1878, 1883; Stewart, 1863). Stewart (1863) noted that a sudden squall occurred almost simultaneously at Oxford and Kew, UK, which are 53 miles away from each other. He documented “a very sudden increase of pressure accompanied with a violent gust of wind”, which was a cold pool associated with convective downdrafts. Hinrichs (1883) presented nice charts of the propagation of squall lines, which he later named derechos (Hinrichs, 1888a, 1888b). Ley (1883) linked these continental squall lines to their oceanic counterparts known to English seaman, who bestowed the name “arched squalls” to all squalls which are seen in perspective to rise as arches of cloud above the horizon. The upper air structure of squall lines was examined using sounding data (Giblett, 1923; Hamilton & Archbold, 1945; Newton, 1950). Observational studies also began on stratocumulus clouds (Sverdrup, 1917; Wyatt, 1923) and trade wind cumulus clouds (Ficker, 1936; Riehl et al., 1951).

World War II led to the birth of radar meteorology (Bent, 1946; Maynard, 1945; Ryde, 1946; Wexler, 1947; Wexler & Swingle, 1947), which, together with the use of aircraft reconnaissance, made it possible to study the three-dimensional dynamical structure of convective systems. Surface precipitation was estimated using Z-R relationships (Laws & Parsons, 1943; Wexler, 1948). Maynard (1945) presented radar echoes for various types of convective systems obtained during World War II. The 1946–1947 Thunderstorm Project, which was one of the largest field experiments in the history of atmospheric sciences, integrated surface stations, upper air soundings, radar, and aircraft reconnaissance to study the dynamical structures of ordinary and multi-cell thunderstorms (Byers et al., 1946; Byers & Battan, 1949; Byers & Braham, 1948; Byers & Hull, 1949; Byers & Rodebush, 1948). Soon after, radars were used to study severe supercell thunderstorms (Stout & Huff, 1953) and MCSs (Ligda, 1956).

Since then, numerous field projects have been conducted, which have significantly advanced our knowledge and understanding of atmospheric convection. The whole atmospheric convection process is initiated by convective updrafts, which are driven by positive buoyancy force, but suppressed by entrainment of dry environmental air. The basic parameters for convective updrafts are size, vertical velocity, and entrainment rate. The size and vertical velocity of convective updrafts have been measured using aircraft flight level data and vertically pointing radars (Black et al., 1996; Giangrande et al., 2013, 2016; Jorgensen et al., 1985; Jorgensen & LeMone, 1989; LeMone & Zipser, 1980; Lucas et al., 1994; May & Rajopadhyaya, 1999; Wang et al., 2020a; Zipser & LeMone, 1980). Convective updrafts over land have a larger size and stronger vertical velocity than those over ocean. The strongest 10% updrafts have an average vertical velocity of ∼4 m/s over ocean, but >8 m/s over land (Fig. 4a,b). The width is ∼2 km over ocean and 3–4 km over land (Fig. 4c,d). Entrainment of lower/middle troposphere air is the leading factor controlling cloud-top height (Brown & Zhang, 1997; Jensen & Del Genio, 2006; Stanfield et al., 2019; Wang et al., 2020b). There are few undiluted updrafts in nature. Brown and Zhang (1997) found that the cloud-top heights of deep convection in TOGA COARE are much lower than those predicted using simple undiluted updrafts, but better explained by an entraining updraft model. Wang et al. (2020) found that the cloud-top heights observed by cloud radar in six ARM tropical sites are 2–8 km lower than those predicted using undiluted updrafts. Using global cloud-top heights observed by CloudSat/CALIPSO and associated carbon monoxide measurements, Stanfield et al. (2019) estimated that the entrainment rate is between 15%/km and 50%/km for deep convection around the world (Fig. 4e). The vertical profiles of entrainment rate have been estimated from cloud-resolving model simulations (Becker et al., 2018; Becker & Hohenegger, 2021; de Rooy et al., 2013; Del Genio & Wu, 2010; Gu et al., 2020; Hannah, 2017; Lu et al., 2018; Romps, 2010; Zhang et al., 2016). Zhang et al. (2016) found that for each cloud category, the entrainment rate is high near cloud base and top, but low in the middle of clouds (Fig. 4f). Becker et al. (2018) found that in the lower free troposphere the bulk entrainment rate increases when convection aggregates, which is against the hypothesis that convective updrafts surrounded by pre-existing convection are undiluted.

Fig. 4 (A) Average vertical velocity in the strongest 10% convective updrafts and downdrafts in oceanic convection comparing with Florida Thunderstorm Project data (from Black et al., 1996). (B) Average vertical velocity in the strongest 10% convective updrafts in Amazon and Southern Great Plain (SGP) (from Wang et al. 2020). (C) Same as A but for convective core width (from Lucas et al., 1994). (D) Same as B but for convective core width. (E) Global estimated entrainment rate binned by cloud top heights (from Stanfield et al., 2019). (F) Vertical profile of entrainment rate for TWP-Ice in Australia (from Zhang et al., 2016).

Many studies have added stochastic processes to convective updrafts using Monte Carlo buoyancy sorting parcels (Emanuel, 1991; Grandpeix et al., 2004; Raymond & Blyth, 1986), stochastically perturbed parameters (Grell & Devenyi, 2002; Grell & Freitas, 2014; Leutbecher et al., 2017), stochastic mass flux distribution (Keane et al., 2014, 2016; Keane & Plant, 2012; Plant & Craig, 2008; Wang & Zhang, 2016), stochastic entrainment (Romps, 2016; Romps & Kuang, 2010; Suselj et al., 2013), stochastic size distribution (Hagos et al., 2018), stochastic cloud types (Goswami et al., 2017; Khouider et al., 2010; Peters et al., 2017), and stochastic triggers (D’Andrea et al., 2014; Rio et al., 2009, 2013; Rochetin et al., 2014a, 2014b). In general, stochastic processes enhance the sway of the convection scheme. Due to the nonlinear nature of convection, expansion of the distribution may make convection more difficult to occur (Wang et al., 2016), enhance the variability of precipitation (Goswami et al., 2017; Peters et al., 2017), and postpone the precipitation maximum in the diurnal cycle (Rio et al., 2009, 2013). There are other modifications related to convective updrafts such as explicit calculation of updraft velocity and chemical transport (Donner, 1993; Donner et al., 2001, 2016; Jeevanjee & Romps, 2015; Lee et al., 2009; Morrison, 2016a, 2016b; Peters, 2016; Simpson & Wiggert, 1969), scale-aware parameters (Freitas et al., 2020; Grell & Freitas, 2014; Han et al., 2017b), machine learning (Gentine et al., 2018; O'Gorman & Dwyer, 2018; Schneider et al., 2017), and unified treatment of boundary layer, shallow convection and deep convection (D’Andrea et al., 2014; Park, 2014a, 2014b). See Rio et al. (2019) for a review of recent studies.

Triggering is generally needed for convective updrafts to happen, since the environment for atmospheric convection is usually conditionally unstable with convective inhibition. When entrainment effect is considered, free convection is further suppressed and initial trigger is even more important. Mapes and Lin (2005) analyzed the life cycle of deep convective systems in seven tropical field experiments and found that boundary layer convergence always leads deep convection by a few hours (Fig. 5). Vertical motion associated with boundary layer convergence can push the air parcel passing the lifting condensation level and achieving positive buoyancy. Coincidently, the two IPCC AR4 models that produced the best MJO simulations were the only ones having moisture convergence closure/trigger (Lin et al., 2006). There are also subgrid-scale triggers such as the density currents from convective downdrafts and low-level gravity waves.

Fig. 5 Lag-regression of divergence profile with respect to surface precipitation for seven field experiments. Lag 0 is the time of maximum precipitation, and lag −10 (+10) hours means 10 h before (after) maximum precipitation. The locations of the field experiments are shown in the top map (from Mapes & Lin, 2005).

The convective downdrafts include negatively buoyant downdrafts and positively buoyant downdrafts. Low-level downdrafts are often unsaturated and have negative buoyancy caused by precipitation loading, evaporation and melting (Knupp & Cotton, 1985). When touching down at the ground, the convective downdrafts in deep convection significantly decrease the boundary layer entropy, and thus suppress the local development of new convection (Barnes & Garstang, 1982; Das & Subba Rao, 1972; de Szoeke et al., 2017; Engerer et al., 2008; Jabouille et al., 1996; Johnson & Nicholls, 1983; Kamburova & Ludlam, 1966; Saxen & Rutledge, 1998; Schiro & Neelin, 2018; Young et al., 1995; Zipser, 1977). It is important to note that although Zipser (1977) suggested that the convective downdrafts are saturated, most of the other studies showed that the convective downdrafts are unsaturated, which is likely because the precipitation particles do not have enough time to evaporate. LES simulations also supported that the convective downdrafts are unsaturated (Hohenegger & Bretherton, 2011; Torri & Kuang, 2016). Convective downdrafts also enhance surface fluxes in convective wakes, which contribute to boundary layer recovery and accumulation of convective instability. The cold pools may also affect future convective organization and vertical structure (Holloway et al., 2017; Tobin et al., 2012). Upper-level warm downdrafts not driven by precipitation or evaporation have been observed in both thunderstorms (Kingsmill & Wakimoto, 1991; Knupp, 1987, 1988; Raymond et al., 1991; Yuter & Houze, 1995a, 1995b, 1995c) and MCSs (Heymsfield & Schotz, 1985; Smull & Houze, 1987; Sun et al., 1993). Sun et al. (1993) conducted thermodynamic retrievals using Doppler radar data and found that the upper-level downdrafts next to convective updrafts are generally positively buoyant. Warm low-level downdrafts with positive buoyancy have also been observed by aircraft flight-level measurements (Igau et al., 1999; Jorgensen & LeMone, 1989; Lucas et al., 1994). The positively buoyant warm downdrafts likely result from the pressure gradient forces required to maintain mass continuity in the presence of adjacent buoyant updrafts (Feynman et al., 1965; Yuter & Houze, 1995b).

Atmospheric convection often occurs in the form of well-organized convective systems. The detailed structures of deep convective systems are shown in Fig. 6. The ordinary thunderstorms generally go through three stages in their lifecycle: the developing stage, mature stage, and dissipating stage (Fig. 6a; Byers & Braham, 1948). In the developing stage the storm is formed from an updraft of air which, as in the other stages, entrains air from the environment. In this stage no rain has yet reached the ground. In the mature stage, rain is occurring, and a large part of the storm consists of a downdraft which characterizes the rain area. The updraft continues in a portion of the storm in the low and intermediate levels and in all parts of the top levels. In the dissipating stage, downdrafts are present throughout, although weak upward motion still exists in the upper levels. Fig. 6a can also be viewed as a schematic of the multi-cell thunderstorm, which shows three cells lining up side by side, and downdrafts from the old mature or dissipating cells triggering new developing cells in the nearby region.

Fig. 6 (A) Circulation within an ordinary thunderstorm in (left) developing, (middle) mature, and (right) dissipating stages (adapted from Byers & Braham, 1948). (B) Visual model of the mature phase of a classic supercell thunderstorm (adapted from Bluestein & Parks, 1983). (C) Vertical cross-section of an MCC (adapted from Fortune et al., 1992). (D) Schematic cross-section of a squall line moving from right to left, and the post-convection sounding. Circled numbers are typical values of θw in °C (adapted from Zipser, 1977).

The supercell thunderstorm is characterized by a persistent, deep, 2–10 km wide, rotating updraft in strong vertical wind shear, which is associated with forward-flank and rear-flank downdrafts (Bluestein & Parks, 1983; Browning, 1964; Fujita, 1958; Lemon & Doswell, 1979; Markowski, 2002; Markowski et al., 2018; Markowski & Richardson, 2009; Markowski & Straka, 2002; Marquis et al., 2008; Wakimoto & Liu, 1998). Fujita (1958) suggested that the tornadic supercell thunderstorms resemble a miniature hurricane with a central eye and spiral echo bands. Lemon and Doswell (1979) found that the structure of supercell storms is similar to a miniature occluded extratropical cyclone. Bluestein and Parks (1983) showed the vertical structure of classic supercell thunderstorms (Fig. 6b). The storm has a penetrating convective cloud top and wide anvil clouds, all tilting downshear and generating precipitation underneath. Wall clouds and tornado generally form upshear, close to the rear-flank downdrafts.

The MCCs show a wide spectrum of structures from chaotic systems to well-organized systems (Augustine & Howard, 1988; Blanchard, 1990; Cotton et al., 1989; Fortune et al., 1992; Kane et al., 1987; Leary & Rappaport, 1987; Maddox, 1980, 1983; McAnelly & Cotton, 1989; Smull & Augustine, 1993; Wetzel et al., 1983). The structure of well-organized MCCs is similar to a small occluded extratropical cyclone with a warm conveyer belt and a cold conveyer belt (Fortune et al., 1992). An MCC is composed of the same four components as a squall line, although the stratiform precipitation exists in a broad region surrounding the convective cores rather than trailing the convective line (Fig. 6c). In the early stage, thunderstorm-scale motions and strong convective precipitation prevail, while in the mature stage, light stratiform precipitation reaches maximum horizontal extent leading to maximum precipitation amount. Then the system slowly decays, producing lighter and lighter rainfall.

The squall lines have been examined extensively by many field experiments (e.g. Bluestein & Jain, 1985; Bryan & Parker, 2010; Gallus & Johnson, 1991, 1992; Grim et al., 2009; Houze, 1977; Houze et al., 1989; Johnson & Hamilton, 1988; Roux et al., 1984; Scott & Rutledge, 1995; Weisman, 2001; Zipser, 1977). As shown by Zipser (1977; Fig. 6d), a squall line has four components: the convective updraft, convective downdraft, mesoscale updraft, and mesoscale downdraft. In general, these four components are also the building blocks of all convective systems. The convective downdraft brings low entropy air into the boundary-layer and suppresses the convective instability, while the mesoscale downdraft warms up and dries up the lower troposphere above the boundary-layer, both of which lead to a “diamond sounding” and suppress the development of future deep convection. Gamache and Houze (1983) calculated the water budget of a tropical squall line and found that the mesoscale updraft accounts for 25–40% of the stratiform precipitation, while the remaining 60–75% is supplied by horizontal transport of condensate generated in the convective updrafts. During their lifetime, squall lines tend to evolve from a symmetric configuration to an asymmetric configuration with bow echo or comma echo (Scott & Rutledge, 1995; Weisman, 2001).

Self-aggregation of convection, the spontaneous organization of initially scattered convection into isolated convective clusters under homogeneous boundary conditions and forcing, has been found in numerous modelling studies (e.g. Bretherton et al., 2005; Held et al., 1993; Muller & Held, 2012; Wing et al., 2017). The formation mechanisms in the models include longwave radiation, shortwave radiation, surface fluxes, moisture feedbacks and advective processes (Wing et al., 2017). However, there is limited observational evidence of convective self-aggregation (Holloway et al., 2017; Tobin et al., 2012; Zuidema et al., 2017), and whether it needs to be specially parameterized needs further studies (Tobin et al., 2013).

Overall, the convective systems have four components: the always-entraining convective updrafts, unsaturated convective downdrafts, mesoscale updrafts and mesoscale downdrafts (Fig. 6d). The observed convection has self-suppression mechanisms caused by entrainment in convective updrafts, surface cold pool generated by convective downdrafts, and warm and dry lower troposphere created by mesoscale downdrafts. The post-convection environment is often characterized by “diamond sounding”, which suggests an over-stabilization of the atmosphere rather than barely returning to the neutral state. Then the pre-convection environment for the future events is characterized by slow moistening of the lower troposphere forced by moisture convergence and surface fluxes. The over-stabilization and slow moistening make the convection events episodic.

As summarized by Lin et al. (2006, 2015) and Lin (2007), there are some parametrization schemes still using undiluted convective updrafts either in the cloud ensemble or in the convective trigger. The undiluted convective updrafts will ignore the suppression effect of a dry lower troposphere, skip the slow pre-conditioning process associated with the development of shallow convection, and lead to an unrealistic quasi-equilibrium state. Possible modifications include setting up a minimum entrainment rate for the cloud ensemble (Tokioka et al., 1988) and using a strong, possibly relative humidity or buoyancy dependent, entrainment (Bechtold et al., 2014; Derbyshire et al., 2004). Such modifications will suppress deep convection and may cause cold temperature biases in the upper troposphere (Gates et al., 1999; John & Soden, 2007; Tian et al., 2013). Adding mesoscale updrafts will warm up the upper troposphere as in nature.

In most convection schemes, the convective downdrafts are saturated, and need to be modified to unsaturated downdrafts. Betts and Silva Dias (1979) and Emanuel (1981) developed models for unsaturated downdrafts, while Emanuel (1991) is the only scheme using unsaturated convective downdrafts. Parameterizations of mesoscale enhancement of surface fluxes have been developed by Qian et al. (1998) and Redelsperger et al. (2000). Parameterizations have also been developed on how convective downdrafts trigger new convection and affect convective organization (Grandpeix et al., 2010; Grandpeix & Lafore, 2010; Mapes & Neale, 2011).

Only one convection scheme has explicitly considered the mesoscale updrafts and downdrafts (Donner, 1993; Donner et al., 2001, 2011; Wilcox & Donner, 2007), although various ideas have been proposed on how to parameterize the mesoscale effects (Alexander & Cotton, 1998; Khouider & Moncrieff, 2015; Mapes & Neale, 2011; Moncrieff et al., 2017; Yano & Moncrieff, 2016, 2018). There are indications that incorporating mesoscale heating structure can drastically improve the simulation of intraseasonal variability such as MJOs in climate models (Cao & Zhang, 2017). The mesoscale heating/moistening structures are quite simple and consistent around the world (Fig. 8a,b). The key unanswered question is what controls the fraction of stratiform precipitation. Observational studies suggested that wind shear and upper-troposphere moisture enhance the formation of anvil clouds and stratiform precipitation (Hogan & Illingworth, 2003; Lin & Mapes, 2004; Saxen & Rutledge, 2000). Cloud-resolving models have a long-lasting bias of underestimating stratiform precipitation and anvil cloud area (Fovell & Ogura, 1988; Franklin et al., 2016; Fridlind et al., 2017; Han et al., 2019; Lang et al., 2003; Luo et al., 2010; McCumber et al., 1991; Morrison et al., 2015; Varble et al., 2011, 2014; Wu et al., 2013).

b Bulk Effects

The amount of surface precipitation represents column-integrated latent heating. Figure 7a shows the GPCP climatological annual mean precipitation for 1981–2010. The largest precipitation is over the tropical continents, Indo-Pacific warm pool, and the tropical/subtropical convergence zones, such as the ITCZ, Mei-Yu, and SPCZ. Satellites provide excellent data for studying climatology of tropical deep convection (Houze et al., 2015; Liu et al., 2007; Liu & Zipser, 2005, 2015; Yuan & Houze, 2010; Zipser et al., 2006). Extremely deep and intense convective elements occur almost exclusively over land, while shallow isolated raining clouds are overwhelmingly an oceanic phenomenon. Continental MCSs tend to have stronger convective regions, with some of the strongest convective regions occurring near the tropical great mountains. Oceanic systems mesoscale convective systems tend to have weaker convective regions but wider stratiform regions. Yuan and Houze (2010) found that the distribution of MCSs is similar to the distribution of total rainfall. MCSs contribute 56% of tropical rainfall, which implies that isolated thunderstorms contribute less than 44% of tropical rainfall. Partitioning between the two types of MCSs have been conducted for the United States (Anderson & Arritt, 1998; Jirak et al., 2003). Using 3 years of data, Jirak et al. (2003) found that 61% of MCSs in the United States were squall lines, and 39% were MCCs. Their data showed that the total amount of rainfall produced by the squall lines double that of MCCs. Over the globe, MCCs mainly occur over land (Laing & Fritsch, 1997). During warm seasons, MCCs contribute 21–26% of total rainfall in western Africa (Laing et al., 1999), 8–18% of total rainfall in central United States (Ashley et al., 2003), 15–21% in central South America (Durkee et al., 2009) and 10–20% in eastern South Africa (Blamey & Reason, 2013). Future studies are needed to partition rainfall of isolated thunderstorms into contribution from different types of thunderstorms (ordinary, multi-cell and supercell).

Fig. 7 (A) The GPCP climatological mean precipitation (contour interval 1.5 mm/day). The coloured boxes are regions used in b and c. (B) Normalized vertical distribution of TRMM precipitation radar 20 dBZ echo top for convective precipitation from 16 years of data (1998–2013). (C) Normalized Q₁ profiles for NAME, TOGA COARE, DYNAMO and GATE.

The distribution of deep convective cloud tops exhibits three regimes (Fig. 7b): the continental deep convection, oceanic deep convection, and oceanic congestus convection. The oceanic congestus convection was first discovered by Johnson et al. (1999) over western Pacific warm pool during the pre-conditioning stage of deep convection. Here we found that congestus convection occurs throughout the tropical and subtropical convergence zones over medium sea surface temperature. For example, the tropical eastern Pacific ITCZ has the strongest climatological mean surface precipitation, but the cloud top of deep convection is much lower than over the western Pacific warm pool. The same feature occurs over the eastern Atlantic ITCZ where the GATE experiment was conducted.

The TRMM radar echo top distributions are confirmed by the vertical heating profiles observed during field experiments. The vertical heating profile of atmospheric convection can be derived from heat budgets using sounding array or reanalysis datasets (e.g. Frank et al., 1996; Frank & McBride, 1989; Johnson et al., 2007, 2015; Lin et al., 2004; Lin & Johnson, 1996; Thompson et al., 1979; Yanai et al., 2000). The methods for deriving apparent heat source Q₁ and moisture sink Q₂ were summarized by Yanai and Johnson (1993). Zhang and Lin (1997) later developed a constrained variational analysis method to make the mass, heat, moisture and momentum budgets self-consistent. Figure 7c shows the Q₁ profiles from NAME, TOGA COARE, DYNAMO and GATE, which clearly demonstrate the three regimes of tropical deep convection. It is important to note that deep convection observed during GATE represents the congestus regime, which is different from the deep convection in TOGA COARE/DYNAMO.

The stratiform precipitation is associated with a universally consistent cloud top distribution (Fig. 8a). Q₁ profile in stratiform region is characterized by heating in the upper troposphere, but cooling in the lower troposphere (Houze, 1982, 1997; Johnson & Young, 1983). Given the standard Q₁, Q₂ profiles associated with stratiform precipitation, any Q₁, Q₂ profile could be partitioned into a convective component and a stratiform component provided that the stratiform rainfall fraction and radiative heating profile are known (Fig. 8b; Johnson, 1984; Lin et al., 2004).

Fig. 8 (A) Normalized vertical distribution of TRMM precipitation radar 20 dBZ echo top for stratiform precipitation for regions shown in Fig. 7A. (B) Partitioning of GATE Q₁ profile into convective, stratiform and radiative components (from Johnson, 1984). (C) TRMM PR convective and stratiform precipitation fractions for NH summer (from Yang & Smith, 2008). (D) CMIP5 model convective precipitation fraction for high-resolution and low-resolution ensemble means for NH summer (from Huang et al., 2018).

The stratiform precipitation contributes significantly to the tropical precipitation and thus the total latent heating (Schumacher & Houze, 2003; Yang & Smith, 2008; Fig. 8c). Over the tropical oceanic convection centres such as the Indo-Pacific warm pool, ITCZ and SPCZ, the stratiform precipitation fraction is about 40–50%, and the convective precipitation fraction is also 40–50%. Over tropical continents, the stratiform precipitation fraction is about 30–40%, while the convective precipitation fraction is about 60–70%. The global climate models significantly underestimate the stratiform precipitation fraction with most precipitation being convective (Dai, 2006; Huang et al., 2018; Fig. 8d). The CMIP5 low-resolution models produce ∼90% convective precipitation fraction over tropical oceans and continents. The high-resolution models have improved simulations, but still produce >70% convective precipitation fraction over tropical oceans (Huang et al., 2018). The stratiform precipitation in climate models is produced by the large-scale condensation schemes or partly from detrained convective condensate. However, the stratiform heating profile might be different from the observed stratiform heating profile with upper-troposphere heating and lower-troposphere cooling, which is important for generating the “diamond sounding” and suppressing future convection.

The regions free of deep convection in Fig. 7a (with <1.5 mm/day climatological mean precipitation) are generally covered by shallow cumulus and stratocumulus clouds (Agee et al., 1973; Agee, 1987; Fig. 9a). In cold seasons, the Mei-Yu front region and Gulf Stream front region are also covered by shallow cumulus clouds. There is a transition from deep convection to cumulus, then to stratocumulus when we move from the ascending branch to the descending branch of Hadley circulation and/or Walker circulation. The mesoscale organization tends to be open cells in cumulus clouds, closed cells in stratocumulus clouds, and no cellularity in the polar/subpolar stratus clouds. The heat and moisture budgets of trade wind cumulus have been calculated using sounding array data from several field experiments (Augstein et al., 1973; Betts, 1975; Brummer, 1978; Esbensen, 1975; Johnson & Lin, 1997; Nitta & Esbensen, 1974; Fig. 9b). The Q₁ profile is characterized by warming in the subcloud layer where vertical eddy heat flux convergence exceeds radiative cooling, and cooling in the upper cloud layer caused by radiative cooling and evaporation of condensate. The Q₂ profile shows moistening in the subcloud layer caused by surface evaporation, and strong moistening in the upper cloud layer caused by the evaporation of cloud water. The longwave heating profile is characterized by a strong cooling at cloud top (Larson et al., 2007; Slingo et al., 1982; Fig. 9c). Vertical velocity in trade wind cumulus has been studied using aircraft data, which increases with height in both updrafts and downdrafts and has a magnitude of ∼0.5–2 m/s (Ghate et al., 2010, 2011; Kollias & Albrecht, 2010; Lamer et al., 2015).

Fig. 9 (A) Global climatology of mesoscale cellular convection depicting the most favoured regions of open and closed mesoscale cellular convection over the oceans (from Agee, 1987). (B) Left: The observed Q₁, Q₂, Q_R, and Q₁ − Q_R for the undisturbed BOMEX period 22–26 June 1969 (from Nitta & Esbensen, 1974). Right: Schematic of trade wind cumulus layer showing effects of condensation c and evaporation e on the heat and moisture budgets (from Johnson & Lin, 1997). (C) Left: Longwave heating rate of a stratocumulus-topped boundary layer (from Larson et al., 2007). Right: Schematic depiction of the large-scale forcing and physical processes for a stratocumulus-topped boundary layer. LTS is lower troposphere stability and EIS is estimated inversion strength (adapted from Lin et al., 2014).

Numerous field experiments have been conducted to study shallow cumulus clouds and stratocumulus clouds (Albrecht et al., 1985, 1988, 1995, 2019; Austin et al., 1996; Bretherton et al., 2004b; Brocks, 1972; Curry et al., 2000; Keuttner & Holland, 1969; Lenschow et al., 1988; Lu et al., 2007; Paluch, 1979; Stevens et al., 2003; Verlinde et al., 2007; Wood et al., 2011; Zuidema, 2018; Zuidema et al., 2016). For the stratocumulus-topped boundary layer (STBL, Fig. 9b), convective instability and turbulence are driven mainly by the cloud-top longwave cooling and evaporative cooling, which are partially reduced by shortwave warming and latent heating inside the cloud layer, and the resulting STBL turbulence is enhanced by latent heating in updrafts and cooling in downdrafts. Turbulent eddies and evaporative cooling drives entrainment at the top of the STBL, which tends to deepen the STBL, maintaining it against large-scale subsidence. Drizzle reduces the liquid water path and albedo and can lead to increased mesoscale variability, stratification of the STBL, and in some cases cloud breakup. For a given cloud thickness, polluted clouds tend to produce more and smaller cloud droplets, greater cloud albedo, and drizzle suppression. Feedbacks between radiative cooling, precipitation formation, turbulence, and entrainment help regulate stratocumulus. The stratocumulus cloud cover is well correlated with the lower troposphere stability (LTS) and estimated inversion strength (EIS), both of which are measures of the temperature inversion strength (Klein & Hartmann, 1993; Norris, 1998; Slingo, 1987; Wood & Bretherton, 2006). The liquid water path of the SEP stratocumulus clouds is often close to the adiabatic value, and thus determined by cloud thickness (Zuidema et al. 2005, 2012; Bretherton et al. 2004, 2010). The cloud thickness is primarily maintained by a strongly negative cloud-radiation-turbulent-entrainment feedback (Zhu et al., 2005), and the thickness could vary due to changes in turbulent driving, vertical gradient of moisture and moist static energy, large-scale subsidence, and inversion strength (Brient & Bony, 2013; Bretherton et al., 2013; Caldwell & Bretherton, 2009; Zhang & Bretherton 2008; Zhu et al., 2007).

In global climate models, shallow cumulus clouds are treated either together with deep convection by a single convection scheme, or by a separate shallow convection scheme. Stratocumulus clouds, on the other hand, are treated by the PBL scheme coupled with microphysics and radiation schemes. PBL schemes have evolved through three stages: (1) local schemes (e.g. Louis, 1979) relating the diffusivity to the local stability, which work well for stable conditions but not for unstable conditions; (2) nonlocal schemes (e.g. Holtslag & Boville, 1993) with nonlocal diffusivity whose magnitude is determined by surface forcing, which works well for unstable conditions driven by surface forcing but not for unstable conditions driven by forcing from boundary layer top, such as the stratocumulus-topped boundary layer; and (3) nonlocal schemes with consideration of cloud-top forcing (e.g. Bretherton & Park, 2009; Grenier & Bretherton, 2001; Lock et al., 2000; van Meijgaard & van Ulden, 1998). Lin et al. (2014) examined the stratocumulus clouds and associated cloud feedback in the southeast Pacific simulated by eight CMIP5/CFMIP global climate models and found that two models could capture the observed stratocumulus clouds and associated cloud feedback, which are the only ones using cloud-top radiative cooling to drive boundary layer turbulence.

Convective momentum transport has been diagnosed for various convective systems using sounding array or reanalysis datasets (Gallus & Johnson, 1992; Hsu & Li, 2011; LeMone et al., 1984; LeMone & Moncrieff, 1994; Lin et al., 2005, 2008; Stevens, 1979; Sui & Yanai, 1986; Tung & Yanai, 2002a, 2002b; Wu & Yanai, 1994; Zhang & Lin, 1997). Wu and Yanai (1994) analyzed the momentum budget of MCSs during SEASAME and PRE-STORM experiments and found that convective momentum transport is downgradient in MCCs, but upgradient in the upper level of squall lines for momentum normal to the squall line. Stevens (1979) analyzed the momentum budget of easterly waves during GATE and found that convective momentum transport is an important term in the meridional momentum budget. Lin et al. (2005) calculated the zonal momentum budget for the MJO, but found that convective momentum transport is not a leading term in MJO’s momentum budget. However, Lin et al. (2008) examined the zonal momentum budget of the Walker Circulation (Fig. 10), and discovered that convective momentum transport, together with pressure gradient force, are the two leading forces driving the Walker Circulation. Convective momentum transport has not been included in many convection schemes. The inclusion of convective momentum transport (Zhang & Cho, 1991) in a GCM showed that it had significant effects on the Hadley circulation (Zhang & McFarlane, 1995b). Therefore, including convective momentum transport in convection schemes is very important for simulating a realistic tropical mean climate.

Fig. 10 Zonal momentum budget of the Walker Circulation, as shown by climatological annual mean (a) pressure gradient force, (b) Coriolis force, (c) advective tendency, and (d) convective eddy momentum flux convergence along the equator (5N-5S) derived from 15 years (1979–1993) of NCEP reanalysis data. Unit is m/s/day (from Lin et al., 2008).

c Closure Assumption

A convection scheme generally has two aspects: the cloud model and the closure assumption. Because the closure assumption determines when the convection will happen and how strong the convective fluxes will be, it is generally considered as a more fundamental characteristic of a convection scheme. The closure assumptions of existing convection schemes can be categorized into three groups (Table 1): (1) moisture convergence (Anthes, 1977; Bougeault, 1985; Frank & Cohen, 1987; Krishnamurti et al., 1976; Kuo, 1965, 1974; Molinari, 1985; Tiedtke, 1989); (2) flux-type convective quasi-equilibrium (CQE) (Arakawa & Schubert, 1974; Bechtold et al., 2014; Chikira & Sugiyama, 2010; Donner, 1993; Emanuel, 1995; Grell, 1993; Grell & Devenyi, 2002; Moorthi & Suarez, 1992; Randall & Pan, 1993; Raymond, 1995; Wu, 2012; Zhang, 2002; Zhang & McFarlane, 1995; Zhang & Wang, 2006; Zhao et al., 2018). (3) state-type CQE (Betts, 1986; Emanuel, 1991; Emanuel et al., 1994; Fritsch & Chappell, 1980; Gregory & Rowntree, 1990; Kain, 2004; Kain & Fritsch, 1990, 1992; Khouider & Majda, 2006; Kuang, 2008; Majda & Shefter, 2001; Manabe et al., 1965; Mapes, 2000; Raymond et al., 2007). The three types of closure assumptions were all proposed in the first generation of convection schemes. Then in the late 1980s and early 1990s, the moisture convergence closures were criticized seriously and started to fade from the global climate models. Most of the remaining convection schemes are using either flux-type CQE closure or state-type CQE closure.

There are fundamental differences between the flux-type CQE closure and state-type CQE closure. They are two different ways to decompose and constrain the change of convective available potential energy (CAPE) or the cloud work function. The flux-type CQE decomposes the CAPE change into its large-scale component and convective component and requires that the CAPE change is much smaller than any of the two flux terms. It was first proposed for the full troposphere (Arakawa & Schubert, 1974; Moorthi & Suarez, 1992; Randall & Pan, 1993; Zhang & McFarlane, 1995) and later also applied to the boundary layer (Emanuel, 1995; Raymond, 1995). There is also a variant of the flux-type CQE called free tropospheric CQE by Zhang (2002) and environmental CQE by Bechtold et al. (2014), which is applied only to the free troposphere and tends to decouple the free troposphere from the boundary layer. Observational budget analysis showed that the flux-type CQE generally is not valid at hourly time scales, but becomes valid at daily and longer time scales (Arakawa & Schubert, 1974; Donner & Phillips, 2003; Zhang, 2003). It is important to note that in climate model implementations of the flux-type CQE, a relaxation time is often introduced for convective adjustment (e.g. Moorthi & Suarez, 1992; Zhang & McFarlane, 1995). In this way, the convective instability is not removed instantly, which tends to make the thermodynamic structure of the model atmosphere shift away from the CQE.

The state-type CQE, on the other hand, provides a stricter constraint on the CAPE change by decomposing it into its boundary layer component and free troposphere component, and requires that the CAPE change is much smaller than any of the two state change terms. It was first proposed for the full troposphere (Betts, 1986; Emanuel et al., 1994; Manabe et al., 1965) and later also applied to only the lower troposphere (Khouider & Majda, 2006; Kuang, 2008; Majda & Shefter, 2001; Mapes, 2000; Raymond et al., 2007). The CQE assumption is very attractive for theoretical modelling because it leads to a very simple picture for global atmospheric circulation and climate variability, ranging from the Hadley and Walker circulations to the Madden-Julian oscillation (Emanuel et al., 1994; Emanuel, 2007). However, because the state change terms are generally much smaller than the flux terms, the validity of flux-type CQE does not guarantee the validity of state-type CQE. Brown and Bretherton (1997) examined the co-variability of ship-observed surface state and satellite-derived troposphere-mean temperature. They found that the constants of proportionality between boundary layer MSE and troposphere-mean temperature were only half of the CQE-predicted value even when the data were subject to a strict precipitation window and averaged over a large region for a long time period. Analysis of soundings from several field experiments also showed that the CAPE change is dominated by its boundary layer component (Donner & Phillips, 2003; Yano et al., 2001; Zhang, 2003).

Lin et al. (2015) examined the validity of the state-type CQE hypothesis at different vertical levels using long-term sounding data from tropical heating centres. The results show that the tropical atmosphere is far away from the CQE, with much weaker warming in the middle and upper troposphere associated with the increase of boundary layer moist static energy. This is true for all the time scales resolved by the observational data, ranging from hourly to interannual and decadal variability. It is likely caused by the ubiquitous existence of cumulus congestus and stratiform precipitation, both leading to sign reversal of heating from lower troposphere to upper troposphere and decoupling of the upper troposphere from the boundary layer. The cold pool generated by convective downdrafts and the warm/dry lower troposphere created by mesoscale downdrafts lead to an over-stabilized post-convection environment. Therefore, the oversimplified 2-phase view of the state-type CQE, which leads to instantaneous occurrence of deep convection, should be replaced by the observed 4-phase structure including the cumulus congestus and stratiform precipitation, which is associated with prolonged timescales and makes convection episodic (Fig. 11).

Fig. 11 Schematic depiction of the vertical structure of tropical atmosphere for (upper) CQE’s 2-phase view, and (lower) observed 4-phase structure. The types of convection are represented by the clouds, while the corresponding profiles of saturation moist static energy anomaly are plotted underneath them (adapted from Lin et al., 2015).

3 The ECMWF convection scheme: evolution and challenges

The convection scheme of the European Centre for Medium range Weather Forecast (ECMWF) Integrated Forecast System (IFS) has undergone significant upgrades and improvements since its original formulation and implementation by Tiedtke (1989). The scheme is also operational, and jointly developed, in the ICON model of Deutsche Wetterdienst since 2014 and is prepared for operational implementation in the Arpège global model of Météo France in 2021.

In 2012 Peter Bechtold was asked in an interview with the then new Director general, Alan Thorpe, “why do we still use a convection scheme that is now more than 20 years old”? Peter replied “it is because the basic equations and physical principles are correct and we should keep that”. This settled the issue and now in 2021 we are still using the basic Tiedtke scheme but with many corrections and extensions to better represent important processes like tropical variability, night-time convection, mesoscale convective systems, the diurnal cycle, ice processes and not to forget numerical stability. Importantly, in this decade high-resolution forecasting including at convection permitting resolutions is becoming more and more prominent. We are currently preparing the IFS and the convection scheme for the next resolution upgrade from currently ∼ 9 km for the high-resolution forecast and ∼16 km for the 50-member ensemble to a O (5 km) global ensemble prediction system in 2025–2027. Our aim has always been more accurate and extended predictions of the coupled atmosphere and ocean system. Convection plays a big art in it and hopefully tracing the past evolution of the IFS scheme and discussing the future challenges gives a reasonable idea to what is important and might be possible in global atmospheric prediction.

a Basic Characteristics

Nowadays all convection schemes used by numerical weather prediction centres are mass flux schemes and so is the IFS scheme which uses the full flux form of the mass flux equations. Such a scheme can be considered as a simplified, but reasonable approximation of the full three-dimensional subgrid convective motions under the assumption of stationarity, the neglect of sub-plume variations and non-hydrostatic pressure forces (Yano et al., 2010; Thuburn et al., 2018). The basic features of the IFS scheme are illustrated in Fig. 12. Shallow or deep convection are represented depending on the cloud depth. If neither of the types can be detected “mid-level”, i.e. local elevated convection is activated when the relative humidity at cloud base exceeds 80% and dynamic lifting is present. The convective drafts consist of a single updraft and downdraft couple that mixes with the environment through detrainment and entrainment. The microphysical processes include mixed phase processes and the generation and fallout of precipitation. Evaporation occurs implicitly through saturated downdrafts and explicitly below cloud base. The convective fluxes are obtained through a rescaling (closure) that for deep convection is based on the convective available energy (CAPE) as introduced by Gregory et al. (2000), while the closure for shallow convection is deduced from the budget of the moist static energy in the sub-cloud layer. The mass flux equations are solved implicitly and the scheme provides convective tendencies for the dry static energy, specific humidity, cloud condensate, rain/snow, momentum and chemical tracers to the IFS. Finally, the link to the prognostic cloud scheme is provided through the detrainment of cloud water/ice that is an important source term and the mass flux subsidence that leads to cloud evaporation. While these basic characteristics are likely rather similar between the different convection schemes used in global models, “details” can have important consequences.

Fig. 12 Schematic of the IFS mass flux convection scheme.

b The Evolution of the IFS Scheme

The evolution of the IFS convection scheme has been driven by the identification of forecast errors at different forecast ranges (e.g. short-range forecast errors in the analysis cycle versus observations or medium-range or seasonal forecast errors versus reanalysis and observations) that can be traced back to the convection scheme.

1 2002/2003 Revised “Trigger Function”

The term trigger function designates a simple procedure to decide on the occurrence of convection and to determine the cloud base properties. The Tiedtke scheme only considered surface based deep convection. Jakob and Siebesma (2003) developed a more realistic sub-cloud model with a strongly entraining parcel departing from the lowest model level. Convection is activated based on the kinetic energy of the parcel at cloud base and the distinction between deep and shallow convection is made depending on a cloud depth threshold of 200 hPa. However, the IFS still strongly underestimated night-time convection over land and instead produced strong grid-scale precipitation, notably during the spring convective season over the continental USA, when the forecasts for Europe were also badly affected by downstream propagating errors. Bechtold et al. (2004) then revised the convective trigger computing convective ascents departing from all model levels below 350 hPa and retaining the first parcel ascent that produces deep convection. As illustrated in Fig. 13 for May 2002 this strongly increased the (night-time) convective precipitation that occurs over the central Great Plains and strongly decreased the excessive grid-scale precipitation, total precipitation is also decreased and in better agreement with observations. Rodwell et al. (2013) concluded that an improved representation of convection over the USA results in improved forecast performance over Europe by strongly reducing the number of bad forecasts. These forecast busts predominantly occur, when there is a strong interaction of convective outflows with the jet stream.

Fig. 13 24–48 h convective and stratiform rainfall (mm day⁻¹) over North America for May 2002 with the operational IFS in 2002 (left column) and with the revised convective initiation allowing convection to depart from any model layer below 350 hPa (right column). This version became operational in 2003.

2 2007 Entrainment, Closure, Numerics

The physics changes in 2007 had probably the largest impact on the tropical forecast performance of the IFS in terms of variability, rainfall distribution and climatology in the last decade. Changes to the shortwave radiation scheme (Morcrette et al., 2008) led to increased convective precipitation over tropical land, while the convection scheme was largely revised as described in Bechtold et al. (2008): the entrainment formulation, consisting of a weak turbulent entrainment and a contribution from moisture convergence was replaced by a relative humidity dependent strong entrainment profile decreasing with height from cloud base values of O (1 km⁻¹) in agreement with large eddy simulations, the detrainment also became relative humidity dependent. Furthermore, the convective adjustment time-scale was no longer constant but computed as the convective turnover time. Finally, an implicit formulation was used to solve for the convective tendencies allowing for large mass fluxes.

The overall impact of these changes is summarized in Fig. 14. The changes in tropical variability are highlighted by the classical wavenumber frequency spectra of the outgoing longwave radiation (OLR) in the tropical band. While the satellite data display the characteristic equatorially symmetric Kelvin, Rossby and MJO spectral signatures in Fig. 14a, these are very weak in the IFS operational model until November 2007 (Fig. 14b) and in particular the Kelvin wave mode is absent. In contrast, with the revised convection (Fig. 14c) the IFS realistically predicts the main modes of tropical variability. The large-scale tropical precipitation pattern also improved, implying a weaker Hadley cell (less precipitation at the equator) and a more intense Walker cell. Furthermore, the prediction-range of the Madden-Julian oscillation is also strongly increased (Fig. 14d) from about 18 days in 2006 to 24 days in 2008, now in 2020 it is around 30 days (not shown). More on the prediction of the MJO and its teleconnections can be found in Vitart and Molteni (2010). Finally, the prediction of tropical cyclones also strongly improved as the background flow improved, with a reduction in cyclone track error (green bars in Fig. 14e) that comparably is larger than the error reduction obtained by the resolution increase from 40 km in 2005 to 25 km in 2007 (blue and red bars in Fig. 14e).

Fig. 14 Wavenumber frequency spectra of the outgoing longwave radiation from NOAA data (a) and from multi-year integrations with the IFS using the operational cycle in 2006 (b) and with the version that became operational in 2008 (c); the MJO spectral band is highlighted by the black rectangle. (d)–(e) measure the gain in prediction skill: (d) evolution of the prediction skill of the IFS for the MJO between 2002 and 2013 as given by the bivariate correlation with the observed empirical orthogonal functions for wind and outgoing longwave radiation, a value of 0.6 (red line) delimits skillfull forecasts (Vitart & Molteni, 2010), (e) statistics of cyclone positions errors (km) as a function of forecast lead time from the 40 km resolution forecasts in 2005/6 (blue), the 25 km forecasts in 2006/7 (red) and the 25 km forecasts in 2008 (green).

3 2013 Diurnal Cycle

Representing the diurnal cycle of convection over land has been and still is an important and difficult challenge for convection parametrization and even convection permitting models. Resolving the diurnal cycle likely requires convection resolving simulations at 1 km resolution (e.g. Lean et al., 2008). A lot of efforts on this subject have also gone into the IFS. Bechtold et al. (2014) discussed the limits of the convective parcel method and proposed a revised convective closure, where only the fraction of the surface heating/CAPE is released to the free troposphere that does not contribute to boundary-layer mixing. As displayed in Fig. 15 for summer 2012, the revised closure produces a more realistic diurnal cycle of precipitation as a function of local solar time (red line in Fig. 15) compared to radar observations (black line) than the default IFS model before 2013 (sky blue lines) which peaks around noon. The revised scheme was able to shift the convection (CAPE) from a maximum at local noon, coinciding with the maximum in the surface heat fluxes, to a maximum in the late afternoon. However, as also evident in Fig. 15, even with the revised closure there is still a substantial underestimation of convection during night-time in the IFS which is related to difficulties in representing propagating mesoscale convective systems.

Fig. 15 Composite diurnal cycle of precipitation (mm day⁻¹) during JJA 2013 over Europe and continental United States from radar observations (black) and from 24 to 48 h reforecasts with the IFS operational cycle in 2012 (sky blue) and the operational cycle in 2013 (red).

In summary, improvements in the diurnal cycle resulted in better forecasts in regions where the convection and the mean flow are diurnal cycle driven, like the Sahel region of Africa. The quality of the 4-dimensional variational analysis and therefore the initial conditions also benefit from improvements in the diurnal cycle of the model through a better assimilation of time-dependent satellite and conventional data. However, the overall medium-range range forecast performance is only moderately improved by shifts in the diurnal cycle as long as the amplitude, i.e. the total convective heating profile is not altered.

4 2017/2018 Mixed Phase Microphysics

Finally, the last major operational upgrade so far targeted the improvement in convective heating/cooling rates in the upper troposphere and near the melting levels, as well as a reduction in the shortwave radiation errors. Biases in the shortwave radiation penalize both the assimilation of satellite data and the coupled seasonal predictions via the feedback from the sea surface temperatures.

Further improvements in the convective heating rates and the distribution of liquid/ice clouds became possible through: revised mixed phase microphysics, including the glaciation of rain and cloud water in the updraft throughout a revised temperature interval, melting that occurs at the wet bulb temperature, detrainment of the liquid condensate phase only for shallow convection and adding the detrainment of convective rain and snow to the prognostic cloud scheme. The impact of all these changes in terms of annual mean shortwave radiation errors at the top of the atmosphere versus the CERES-EBAF satellite product is shown in Fig. 16. Here we compare annual mean bias from multi-year integrations with the model cycle operational in 2019 (Fig. 16a) and with the same model version, but all ice phase microphysical changes to the convection that have been added between 2017 and 2019 reverted (Fig. 16b). Regions that are particularly affected by the model upgrade are the southern hemisphere storm tracks and to a lesser extent the northern hemisphere storm tracks which all become more reflective through an increase/decrease of the liquid/ice phase in the upper part of the shallow and congestus clouds. Temperature/wind in the upper troposphere and temperatures near the melting level now also better match radiosonde observations (not shown).

Fig. 16 Cloud and radiation evaluation from multi-annual coupled integrations with the IFS Cy47r1. (a)–(b) difference in top of atmosphere net shortwave radiation (W/m²) between the model and the Earth’s Radiant Energy System (CERES) Energy Balanced and Filled (EBAF) product for (a) the operational model version in 2018 and (b) as (a) but with all the changes relating to the mixed phase microphysics added during 2016–2018 removed.

c Mesoscale Convective Systems and Challenges at High Resolution

The lack of night-time convection over land has already been discussed in the context of the diurnal cycle of convection. Today we consider this as the major error in the IFS forecasts of convective activity. That this error is related to the representation of mesoscale convective systems is shown in Fig. 17 which displays the evolution of convection on 12 August 2017 at 15, 18 and 21 UTC over Central Africa and the Sahel region as observed by the 10.8 μm infrared channel of Metsosat-10 (Fig. 17a–c) and the 3-hourly rain accumulations from the TRMM radar product 3B42 (Fig. 17d–f). Consistently, these observations show mesoscale convective systems, notably those near 15°N that intensify during the afternoon and early night-time hours and propagate westward.

Fig. 17 Evolution of continental convective systems over tropical Africa during 12 September 2017 in 3-hourly slots from 15 to 21 UTC as seen by Meteosat-11 infrared image at 10.9 μ wavelength (a,b,c), as well as 3 hourly accumulated rainfall (mm) from 12 to 15, 15 to 18 and 18 to 21 UTC from the TRMM 3B42 product (d,e,f), from the 4 km IFS reforecasts with (g,h,i) (operational version) and without (j,k,l) the deep convection scheme, and with the revised deep convective closure (m,n,o). The IFS reforecasts start at 11 September 2017 at 00 UTC and use the model cycle operational in 2019. There is no TRMM 3B42 data East of 25°E at 21 UTC.

To explore the potential of the IFS at future higher resolutions we have rerun this case with the operational 2019 cycle but at 4 km horizontal resolution with (Fig. 17g–i) and without (Fig. 17j–l) the deep convection parametrisation. As developed through a collaboration with G. Zängl at the DWD (Offenbach) the deep convection scheme includes a smooth reduction of the parametrized convective fluxes (Malardel & Bechtold, 2019), and therefore a transition to resolved convection with increasing resolution (higher than 8 km).

However, with the deep convection parametrisation the rainfall patterns in Fig. 17g–i are too broad scale and the night-time propagating systems at 15°N are absent; similar results are obtained with the operational 9 km horizontal resolution (not shown). In contrast, without the deep convection parametrisation the IFS better simulates the intense westward propagating mesoscale systems. Unfortunately, the amplitude of these systems is too strong as is evident from the comparison with the TRMM data and the global precipitation is overestimated by more than 10%. Also, the root mean square error of precipitation and upper-air forecast skill are significantly degraded with this version of the model.

We therefore further explored the coupling between the convection and the dynamics which is particularly delicate in the case of mesoscale convective systems that propagate and regenerate by producing their own horizontal convergence. Together with Tobias Becker we analysed output from the explicit convection runs over Africa for the whole month of August. It was found that the lack of intense continental convection in the parametrisation can be corrected for by including the vertically integrated advective moisture tendency in the convective instability closure. The results with the revised closure at 4 km are displayed in Fig. 17m–o. Indeed, the convection is now more intense than the current scheme and realistic propagating features develop when compared with the observations in Fig. 17d–f. The revised convective closure now also closely reproduces the satellite observed rainfall distribution (not shown), overall the results are now somewhere in between the current operational scheme with a CAPE closure and the simulations without the deep convection scheme.

Evaluations are ongoing and we plan to implement operationally the above CAPE closure with a moisture convergence term in 2021. As stated in the introduction, we are aiming for a O(5 km) ensemble in 2025–2027 which should also include a revised stochastic physics scheme, namely stochastically perturbed parameters (Leutbecher et al., 2017), where among many parameters from the physical parametrizations, six important parameters of the convection scheme are perturbed. We hope that the ensemble will then be able to largely explore the uncertainties in the predictions and the uncertainties in convection in particular. Ideally, one would aim for a fully prognostic description of convection as has been implemented in a regional model by Gerard (2015). However, such a scheme requires many additional prognostic variables, its closure is not straightforward as is its application in a global model with a 4D-Var data assimilation cycle. There is currently a rapid increase in regional and global applications with explicit deep convection and therefore a fully prognostic convection scheme might eventually never be used in a global model.

IFS operational predictions with explicit deep convection are not likely before 2030, but colleagues (Wedi et al., 2020) are already exploring 1.4 km explicit convection runs of the IFS at the most powerful computer SUMMIT. Preliminary results indicate that the current 9 km model with parametrized convection (effective resolution 30–60 km) and the 1.4 km model give very similar results in terms of spectral energy budget, diurnal cycle etc., with the 1 km model mainly improving scales <100 km. The 9 or 4 km run without deep convection produce however significantly worse results in most aspects including a delayed onset of convection, an overestimation of global precipitation by 5–10% and have difficulties in representing the MJO. Deep convection parametrization seems to remain competitive!

4 The Grell convection scheme

The development of the original Grell convective parameterization began as doctoral work focused on research and science. Since then, the scheme has undergone many transitions, motivated by both research interests and operational constraints. In 1974, Arakawa and Schubert (1974, hereafter AS) and Lord (1978) proposed using a spectrum of updrafts to represent clouds in numerical weather prediction models. Modification of this idealized model served as the basis for the original work that resulted in the Grell scheme. In the first phase of development, the AS scheme was modified to include additional physical processes, most importantly buoyant downdrafts (Grell, 1988). This development was also related to transport of tracers in chemical transport models, where the impact of downdrafts is substantial. Buoyant plumes (positive or negative buoyancy) are very effective transport mechanisms. In a second phase of development (1992–2012), a large effort focused on the simplest and most computationally efficient implementation of convective parameterizations to improve the operational forecasts (Grell, 1993; Grell & Devenyi, 2002; Pan & Wu, 1995). The third phase then focused again on improving the representation of physical processes. This third phase is still ongoing with a resurgent interest from young scientists (Freitas et al., 2021; Grell & Freitas, 2014; Han et al., 2017).

a Modifying and Expanding the Arakawa-Schubert Scheme

As a first step, a parameterization of downdrafts was developed (Grell, 1988) and added to the spectral AS scheme. The downdraft model followed that of an inverted plume (Houze, 1977; Johnson, 1976; Nitta, 1977) driven by the generation of negative buoyant energy, and, therefore, was restricted only to the lower troposphere. The downdraft originating level was assumed to be at the level of minimum moist static energy – in the lowest 300 mb of the atmosphere in the original implementation – to allow the most negatively buoyant downdraft. The original AS parameterization uses a closure that demands a quasi-equilibrium between the “larger scale” and the convection. In AS (1) $\left({\displaystyle \frac{dA}{dt}}\right)={\left({\displaystyle \frac{dA}{dt}}\right)}_{LS}+{\left({\displaystyle \frac{dA}{dt}}\right)}_{CU}\approx 0$(1) where A was defined as the cloud work function, a measure of available buoyancy for a particular cloud type. Subscript $LS$ stands for changes due to the larger scale, subscript $CU$ for changes due to convection. In discretized form this equation became an elegant look (Lord, 1978), but was still complex to solve, since no exact solution existed. Instead, the total change of the cloud work function was minimized using a linear programming method to come to a “quasi” equilibrium. In commonly used applications of the AS scheme, the cloud work function was moved back to climatological values, derived from large observational field campaigns. In Grell-type schemes, Eq. (1) was modified to (2) $\left({\displaystyle \frac{dA}{dt}}\right)={\left({\displaystyle \frac{dA}{dt}}\right)}_{Forc}+{\left({\displaystyle \frac{dA}{dt}}\right)}_{CU}\approx 0.$(2)

In other words, the first term was never perceived as “large” scale, but always as changes to the cloud work function due to impacts from larger scales as well as from other physical processes and the subscript is here replaced by $Forc.$ The most obvious of these physical processes are changes from boundary layer or radiation parameterizations. Later, even changes from convective activity would become part of “Forc”. Additionally, all Grell type schemes – in contrast to the original AS schemes – used forcing functions from the NWP model (large scale advection, PBL parameterization, radiation parameterization) directly to calculate the forcing term. A detailed description of the implementation of the scheme can be found in Grell (1988) and Grell (1993). This version of the parameterization was tested extensively and compared to other parameterizations (Kreitzberg & Perkey, 1976; Krishnamurthi et al., 1980) semi-prognostically using data from the severe midlatitude convective storms field experiment observed during SESAME 1979 (Grell et al., 1991). Additionally, the scheme was tested fully prognostically, again in the mid-latitude severe storm environment, this time with data from the 1985 Pre-STORM experiment (Grell, 1993). The importance of downdrafts was shown for the tests presented in the results. The AS scheme with implemented downdrafts and this application of the “quasi-equilibrium” assumption worked well in the mid-latitude severe storm environment. This included the simulation of a squall line and a mesoscale convective system. A version of the Pennsylvania State University-National Center for Atmospheric Research hydrostatic mesoscale model (MM4; Anthes et al., 1987) was used for this study with 25 km horizontal resolution and 19 vertical levels. To compare to previous studies with a version of the Fritsch and Chappell (1980) scheme (Zhang et al., 1989) the setup was identical to the one chosen by Zhang et al. (1989).

b Improvements Focused on Operational Implementations

Most significant during this time was probably the development of an extremely simple version of the scheme that later was widely used in operations and research and was known as the Grell scheme (hereafter G1). For implementation in operational applications, computational efficiency is essential. In the early 1990s we were looking for a convective parameterization in the Rapid Update Cycle (RUC), which was implemented into operations at the National Center for Environmental Prediction (NCEP) in 1994. The earlier developed version of the AS scheme was too expensive to be used in operations. Hence, the data sets and simulations shown in Grell (1993) were used to simplify the original scheme. The linear programming method solution in the spectral implementation appeared to prefer only a few cloud types. So instead of requiring the quasi-equilibrium, we assumed an equilibrium, and assumed that one cloud type was sufficient to characterize deep convection. Hence (3) ${\left({\displaystyle \frac{dA}{dt}}\right)}_{Forc}=-{\left({\displaystyle \frac{dA}{dt}}\right)}_{CU}$(3) where (4) ${\left({\displaystyle \frac{dA}{dt}}\right)}_{Forc}={\displaystyle \frac{A\left(t+dt\right)-A\left(t\right)}{dt}}$(4) and (5) ${\left({\displaystyle \frac{dA}{dt}}\right)}_{CU}=f\left(CU\right)\ast m_{b\left(CU\right)}$(5) A and f are fairly simple functions in dependence of environmental conditions and assumptions made by the cumulus parameterization. The feedback equations are all normalized by $m_{b\left(CU\right)}$ and can be expressed as (6) ${\displaystyle \frac{\mathrm{\partial }\phi \left(k\right)}{\mathrm{\partial }t}=\varrho \left[\phi \left(z\right)\right]\ast m_{b\left(CU\right)}}$(6) where $\varrho$ is the change of variable $\phi$ per unit of mass. With this simple equation and only one type of cloud allowed, it was realized that the closure is easily changed to other approaches. Stability closures would simply be expressed with (7) ${\displaystyle \frac{A}{\mathrm{\Delta }\tau }=-f\left(CU\right)\ast m_{b\left(CU\right)}}$(7) where $\mathrm{\Delta }\tau$ is some specified time interval over which the instability is removed. Moisture convergence closures were also easy implemented, allowing for direct calculation of the cloud base mass flux, since the total rainfall was related to $m_{b\left(CU\right)}$ with similarly simple equations. Even though the original version of the scheme was oversimplified (no entrainment/detrainment with height), it served as a stepping stone for further developments. At first, different closures were tested and employed for different applications. Giorgi et al. (2012) used the simple scheme in his version of the REGional climate modelling system (REGcm) with different closures.

G1 was also taken by NCEP and modified and implemented as the Simplified Arakawa Schubert (SAS, Pan & Wu, 1995) scheme, of which an improved version is still in operations in various modelling systems at NCEP today. Original modifications by Pan and Wu included simplifying the closure again by assuming that the cloud work function is moved back to a climatological value (as in the original AS approach), but also adding lateral mixing. The SAS scheme went into operation in NCEP’s Global Forecast System (GFS) in 1993. In 2010 the SAS scheme was largely revised by Han and Pan (2011). Some of these revisions were similar to those that were implemented in parallel in Boulder for use in the regional forecast systems (RUC, later the RAP). This included parameterizing the effects of the convection-induced pressure gradient force on convective momentum transport, and the detrainment of cloud condensate into the grid-scale condensate from upper cloud layers above the downdraft initiating level. To deplete more instability in the atmospheric column and result in the suppression of excessive grid-scale precipitation, the scheme was modified to make cumulus convection stronger and deeper by enhancing the maximum allowable cloud base mass flux and by having convective overshooting, respectively.

Another major upgrade for the GFS SAS scheme was made in 2017 (Han et al., 2017). A scale-aware parameterization where the cloud mass flux decreases with increasing grid resolution was developed similar to GF. A simple aerosol-aware parameterization where rain conversion and cloud condensate detrainment in the convective updraft is modified by aerosol concentration number was included. The cloud base mass-flux computation was modified to use convective turnover time as the convective adjustment time scale. The rain conversion rate was modified to decrease with decreasing air temperature above the freezing level, giving rise to more detrainment of cloud condensate in the upper updraft layers and consequently more high clouds. To suppress unrealistically spotty precipitation especially over high terrain during summer time, a convective inhibition in the sub-cloud layer was employed as an additional trigger condition.

In parallel, at NOAA’s Forecast Systems Laboratory (FSL, which is now the Global Systems Laboratory – GSL) the focus for development shifted to the use of different closure assumptions. While all convective parameterizations attempt to statistically represent the interaction between convection and the environment, the closures and assumptions used in individual parameterizations can fundamentally differ. Additionally, some of the parameters used in convective parameterizations are highly sensitive to their assumed value, which can create variability within different iterations of the same convective parameterization (Grell & Devenyi, 2002, GD). While convective parameterizations can be criticized for this variability, GD takes advantage of this diversity. As described above the simple Grell scheme could be used with many different closures, and feedback assumptions were easily changed or modified. Rather than using one set of closures and parameters, GD uses different closures and parameters to create an ensemble. This ensemble is then used to determine the value of cloud base mass flux as well as other feedback properties. The most straightforward and cost-effective approach uses the simple mean of the ensemble. This is the technique that is currently implemented in GD in the Weather and Research Forecast Model (WRF) (Skamarock et al., 2008). GD was implemented into operations in the RUC in 2002.

There are several ways in which GD can be modified so that a more optimal value of cloud base mass flux is selected, since the simple mean is not necessarily the “best” value to choose. For example, a weighted mean of the ensemble members could be calculated. Not only could this weighted mean have different weights assigned to each ensemble member, but those weights could also vary across the domain of the model. Given observational data, techniques such as the linear least squares fitting method, Bayesian data assimilation using probability density function fitting (Grell & Devenyi, 2002), and artificial intelligence could be used to select a better value of cloud base mass flux. Neural networks were also tried to determine better weights with results similar to the least square fitting approach. In a unique approach, the Firefly weighting method proposed by Yang (2008, 2010) was implemented by Dos Santos et al. (2013), based on the bioluminescence process which characterizes firefly mating behaviour. The important aspects that are used in this method by Dos Santos et al. (2013) in determining weights of ensemble members are the light intensity variation and the attractiveness formulation. For simplicity, the attractiveness is determined by the brightness, which in turn is associated with the objective function. For details the reader is referred to Dos Santos et al. (2013).

c Implementing Advanced Physical Processes – for Research and Operations

In parallel to the development of the GD scheme, discussion started about scale-awareness, since horizontal resolution continued to increase in research and operations every year. This led to the development of a three-dimensional application of the GD scheme, which became the G3d parameterization (described also in Grell & Freitas, 2014). The assumption used was that convection, with increased resolution, should impact more than one grid point. Its effects should not be limited to just one grid box, so the parameterization was applied over several grid points (not just one). The ensembles were kept as in GD. G3d became operational in 2012 with the implementation of the first version of the RAPid refresh model (RAPv1).

While the applications of the scheme over neighbouring grid points appeared as the physically most realistic approach, it introduced significant complications, since G3d had to be applied in three dimensions, requiring data communication in parallel applications. Additionally, it was not clear how many grid points would need to be impacted. In parallel, Arakawa et al. (2011) developed an elegant method reobtaining the equation for the vertical eddy transport in terms of the fractional area covered by the active cloud draft and the vertical eddy transport term given by conventional (non-scale aware) tendencies. In that way, the convection parameterization tendencies are simply scaled by a factor that depends on area coverage. Neighbouring grid points are however not impacted. This became one of the foundations of the Grell and Freitas (GF) parameterization, but with many new additional changes and features. This scaling approach is now used in many other convective parameterizations to parameterize scale-awareness, although not all developers claim to follow Arakawa et al. (2011).

GF and G3d scale-aware performance was evaluated many times. Examples are given by Grell and Freitas (2014), Fowler et al. (2016) for variable resolutions, Freitas et al. (2017, 2020, 2021). The paper by Grell and Freitas (2014) showed that the performance of G3d was very similar to GF with respect to scale awareness, but no scaling is applied in G3d. GF used horizontal grid spacings of 20, 10 and 5 km over an area of South America. The full GF scheme (GF-A) was tested on all resolutions and compared to G3d as well as a run without any convection parameterization (NO-CP) and a run with no scale awareness (GF-NS). Each experiment included 15 runs from January 1 to January 15 for 36-hour forecasts, all starting at 00UTC.

When averaged over the domain, a necessary requirement for successful scale-aware applications can usually be seen when comparing resolved versus non-resolved precipitation. This is shown in Fig. 18. With an increase in resolution the precipitation transits from mostly generated by non-resolved to mostly generated by resolved processes (explicit microphysics). In the experiment GF-NS, the parameterized precipitation is much larger than the resolved one, even for a grid spacing of 5 km. Also, the G3d scheme shows a similar ratio to the GF-A formulation.

Fig. 18 Fraction of the resolved precipitation compared to the total precipitation. The 6-hourly model areal mean precipitation rates are averaged for each experiment over the 15 runs.

The root mean square error (RMSE) and mean error (Bias) of the precipitation forecasts for the same experiments are seen in Fig. 19. The observed 6-hourly rainfall data were collected at 861 rain-gauge stations distributed over South America. For GF-A on 20, 10, and 5 km horizontal resolutions the daily mean values for RMSE and Bias (in mm/6hr) are (1.85, 0.45), (1.80, 0.32) and (1.81, 0.12), respectively. G3d has very similar performance with daily mean RMSE and Bias of 1.83 and 0.09 mm/6hr, respectively. Turning off the convective parameterization on 5 km (NO CP) leads to a mean Bias of −0.15 mm/6hr and increases the RMSE to 1.98 mm/6hr; larger than the errors of the 10 and 20 km experiments. The GF-NS showed to the worst performance, with a much higher RMSE and Bias.

Fig. 19 As in Fig. 18 except for Root Mean Square Error (RMSE) and mean error (Bias). Units are mm/6hr.

While many other authors also showed impressive performance of this Arakawa et al. (2011) approach as well as other scaling approaches – it was also adapted in the SAS scheme (described above in Section 4.2) – GF never was implemented into the operational storm-scale modelling system, since not only the success but also the flaws of the Arakawa approach became much more visible for storm scale evaluation metrics. Although its application led to improved bias and other commonly used precipitation scores in evaluations, storm-scale verification also depends on realism of precipitation distribution and simulated radar reflectivities. Physically the scaling remains somewhat unrealistic, and this is seen particularly in areas with relatively weak forcing. An example is given in Fig. 20 which shows observed and simulated precipitation over the South East US, using the operational HRRR model with and without the GF scheme, as it was implemented in December of 2020. The scaling leads to wide spread lighter precipitation over Georgia and Alabama. Using no parameterization on the other hand – will lead to a low bias for light precipitation and an over-forecast for large thresholds (Fig. 21). Through tuning (GF employs a threshold when the parameterized clouds are assumed to shrink in size) the bias can be brought in much closer agreement to observation, but the hourly comparison to radar become very unrealistic. Work is now ongoing (see the last section) to bring some of the three dimensionality of G3d back into the GF scheme.

Fig. 20 6hr precipitation forecasts from the HRRR (middle and right panel) for runs without convective parameterizations (middle) and runs with the GF scheme (right panel), compared to observations (right panel) for the same period.

Fig. 21 Frequency BIAS ration for 12hr accumulated precip, August 8–9, a HRRR run without convective parameterizations (blue) and with a version of the GF scheme (red) in dependence of threshold precipitation amounts over the 12 h period ending on August 9, 00z.

Even more foundational in GF were 2 other changes. (1) The uni-modal approach was replaced with a tri-modal formulation, which allows the three convective modes to exist (Johnson et al., 1999): shallow, congestus, and deep convection. In this unique approach, each of the modes is distinguished by an assumed average size (with a characteristic initial entrainment rate) that strongly controls its vertical depth and profile. (2) While many convective parameterizations use assumptions on entrainment and detrainment to derive a vertical mass flux profile, this was changed in GF by introducing a PDF approach. A characteristic PDF for deep convection is assumed to characterize the vertical mass flux profiles by deep convection in the grid box, and then entrainment and detrainment rates are derived from the vertical mass flux profiles. For deep convection, the maximum of the PDF is given by the stability minimum. For shallow convection, the maximum is just at or above cloud base. Figure 22 shows an example of deep convective mass flux in comparison to observations during the TWP-ICE field experiment. While the closure for deep convection was still using the original ensemble idea from GD, a set of different approaches were introduced to adequately account for the diverse shallow and congestus regimes of convection in a given grid cell. The three modes all transport momentum, tracers, water, and moist static energy. GF also changed the microphysics though inclusion of mixed-phase impacts on buoyancy in the parameterized clouds.

Fig. 22 A Comparison of deep, congestus (mid), and shallow convection massflux profiles for Single Column Model simulations (left) and radar observations (right) during the TWP-ICE field experiment. Observations are from Kumar et al. (2015) showing mass flux profiles derived from windprofilers (black) and CPOL data (red).

Another new feature is the inclusion of a diurnal cycle closure adapted from Bechtold et al. (2014), which notably improved the simulation of the diurnal cycle of convection and precipitation in Single Column Approaches as well as NASA’s GEOS modelling system (Freitas et al., 2018). Figure 23 compares the net grid-scale vertical moistening (left panels) and heating (right panels) tendencies associated with the three convective modes with (B and D) and without (A and C) the Bechtold closure and for a time dependent SCM run from January to 25 February 1999. Without applying the diurnal cycle closure, the three convective modes coexist and are triggered just a few hours after the sunrise, with the deep convection occurring three to six hours too early, producing a maximum of precipitation about 15 UTC (approximately ∼ 11h local time).

Fig. 23 The time average of the diurnal cycle of the grid-scale vertical moistening (panels A and B) and heating (panels C and D) tendencies associated with the three parameterized convective modes (shaded). The total precipitation and the GF parameterized precipitation from the deep and congestus plumes are shown by the graphic lines: black, green, and purple, respectively. The upper rows show model results without Bechtold’s closure, while in the lower row, this closure is applied (see text for further details).

Including the closure (Panels B and D), a smooth and realistic transition is simulated with a late morning and early afternoon low/mid tropospheric moistening by shallow and congestus plumes, followed by a late afternoon and early evening tropospheric drying by the rainfall from the deep cumulus. In this case, the precipitation from the deep penetrative convection is delayed (green contour) with the maximum rate taking place between 18 and 21 UTC, more consistent with the observations of the diurnal cycle over the Amazon region.

d Ongoing and Future Work

Aspects of GF that continue to be actively developed include the treatment of aerosols, inclusion of memory effects, improving microphysics in the parameterized clouds, cold pool generation and movement, and stochastic approaches, as well as further improvements of the scale-awareness (described above). As one implementation of memory effects the parameterization for aerosol wet scavenging has been revised and made proportional to the precipitation efficiency based on Wang et al. (2013). Aerosols in GF can also now be initialized using aerosol optical depth data that varies geographically based on data such as NASA’s Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2). In the future, for operational applications at NCEP, aerosol optical depth data from NOAA’s global aerosol modelling system could also be used. The coupling of the particle distribution in GF and double moment microphysics parameterizations has also been improved by adding diagnostic estimates of cloud water and cloud ice number concentrations to GF. This modification is designed to diminish the artificial modification of the particle size distribution that occurs when single-moment microphysics from a convective parameterization is used with a double-moment explicit microphysics scheme.

In parallel the SAS scheme is also under continued development. Recently, the scheme was modified to interact with the vertical turbulent mixing scheme in the GFS. Originated from the NCEP SAS scheme in 2010, continuous efforts to improve its performance have been given at the Korea Institute of Atmospheric Prediction Systems (KIAPS). The Korean Integrated Model (KIM) has been operational at the Korea Meteorological Administration (KMA) since April 2020, and a description of the revised SAS scheme in KIM is given in Han et al. (2020).

5 The Zhang-McFarlane convection scheme

The Zhang-McFarlane convection scheme was developed in the early 1990s to replace the moist convective adjustment scheme used in the Canadian Climate Centre General Circulation Model at the time. One of the objectives was to represent convection in a mass flux form so that convective transport could be represented. The paper documenting the scheme and its performance in the model was published a few years later (Zhang & McFarlane, 1995, hereafter ZM scheme). Shortly after, Zhang collaborated with scientists at the National Center for Atmospheric Research (NCAR) and implemented the ZM scheme into the NCAR Community Climate Model version 3 (CCM3) to represent deep convection in the model (Zhang et al., 1998). It has been used since then in subsequent generations of the NCAR model. Over the past quarter century, several modelling centres have also adopted the ZM scheme in their models (Bentsen et al., 2013; Lee et al., 2020; Li et al., 2013; Lin et al., 2020; Xie et al., 2018). Various modifications have also been made to address model simulation deficiencies deemed to be rooted in the parameterization of convection. This section will describe the ZM scheme and its updates.

a The Original Zhang-McFarlane Scheme

The ZM scheme was based on the quasi-equilibrium hypothesis proposed in the classic work of Arakawa and Schubert (1974). It has three main components, the trigger function, a 1-D simplified cloud model and the closure. The trigger function answers the question of whether there will be convection, given the atmospheric conditions. The 1-D cloud model determines the vertical distribution of cloud updraft and downdraft properties and convective heating and drying. The closure determines, in a statistical sense, the amount of convection, often measured in terms of cloud-base convective mass flux.

1 Trigger Conditions

The triggering of convection can depend on many factors. In current convection parameterization schemes, many different convection trigger functions have been used. In the Arakawa and Schubert (1974) scheme a threshold value of cloud work function is used for deep convection trigger. The Tiedtke scheme (Tiedtke, 1989) assumes that convection will be triggered if the atmospheric column has a net moisture convergence and the surface air is buoyant if lifted to the lifting condensation level. In the Kain-Fritsch scheme (Kain, 2004), the trigger function is determined by the vertical velocity and the height of lifting condensation level of parcels lifted from the boundary layer. The convection trigger function in the Donner scheme (Donner, 1993; Donner et al., 2001) requires that the large-scale vertical velocity at the convection initiation level integrated over a time span be able to lift the parcel to the level of free convection (LFC). The Zhang-McFarlane scheme uses a CAPE threshold for convection trigger function. CAPE is defined by (8) $CAPE={\int }_{p_t}^{p_b}{R}_d\left(T_{pv}-T_v\right)d\ln p$(8) where $T_{pv}=T_p\left(1+0.608q_p-q_l\right)$ and ${\overline{T}}_v=\overline{T}\left(1+0.608\overline{q}\right)$ are virtual temperatures of the air parcel and the large-scale mean. p_b and p_t are pressure values at the parcel’s initial level and the neutral buoyancy level, respectively. R_d is gas constant for dry air and q_l is the liquid water condensed following the reversible moist adiabat of the air parcel. The parcel originates at the most unstable layer within the planetary boundary layer. Subscript p stands for the parcel’s properties.

2 The Bulk Cloud Model

The bulk cloud model is based on the idea of entraining plume ensembles, that is, cloud updrafts consist of an ensemble of entraining plumes, each of which has a constant entrainment rate. The plume top is reached when the plume becomes neutrally buoyant. However, instead of determining the cloud properties of each plume individually, the budget equations for the ensemble-mean properties are solved. It is assumed that the plumes are in a steady state. Thus, the budget equations for the updraft ensemble mean are given by (9) ${\displaystyle \frac{\mathrm{\partial }{\eta }_u}{\mathrm{\partial }z}=e_u-d_u}$(9) (10) ${\displaystyle \frac{\mathrm{\partial }{\eta }_u{s}_u}{\mathrm{\partial }z}=e_u\overline{s}-d_u\hat{s}+L\rho c_u}$(10) (11) ${\displaystyle \frac{\mathrm{\partial }{\eta }_u{q}_u}{\mathrm{\partial }z}=e_u\overline{q}-d_u\hat{q}-\rho c_u}$(11) (12) ${\displaystyle \frac{\mathrm{\partial }{\eta }_u{l}_u}{\mathrm{\partial }z}=-d_u{l}_u+\rho c_u-R_u}$(12) where $\eta$ is the cloud mass flux normalized by its value at the cloud base. $e$ and $d$ are normalized entrainment and detrainment rates, respectively. $s=c_{p}T+gz$ is the dry static energy, $q$ is specific humidity, and $l$ is cloud liquid water content. Subscript u denotes updraft. The overbar denotes large-scale grid mean properties in GCMs. $\hat{s}$ and $\hat{q}$ are detrained s and q from updrafts. $c$ is the rate of condensation and $R$ is the rate of conversion of cloud water to rainwater given by (13) $R_u=c_0{\eta }_u{l}_u$(13) where $c_0$ is the conversion coefficient, set to $5.9\times {10}^{-3}{\mathrm{m}}^{-1}$ for land and $4.5\times {10}^{-2}{\mathrm{m}}^{-1}$ for ocean in the NCAR Community Atmosphere Model (CAM). s and q in Eqs. (10) and (11) are often combined to form the moist static energy equation (14) ${\displaystyle \frac{\mathrm{\partial }{\eta }_u{h}_u}{\mathrm{\partial }z}=e_u\overline{h}-d_u\hat{h}}$(14) where $h=s+Lq$ is moist static energy.

The cloud base is treated the same as the convection initiation level, which is at the level of highest moist static energy in the PBL. The cloud top is set at the neutral buoyancy level for each plume. The detrained air is assumed to have the same virtual temperature as the environment and is saturated. The detrained liquid water is assumed to have the same value as that in the updraft.

The vertical profile of the updraft mass flux is specified through fractional entrainment and detrainment rates. For a cloud type with a constant fractional entrainment rate λ its mass flux varies exponentially with height. Integrating over all possible λ’s that contribute to the mass flux at height z gives: (15) $M_u\left(z\right)={\int }_0^{{\lambda }_D\left(z\right)}{m}_b\left(\lambda \right)\exp \left[\left(\lambda \left(z-z_b\right)\right)\right]d\lambda$(15) where $M_u\left(z\right)$ is the total mass flux at height z and $m_b\left(\lambda \right)$ is the mass flux at the cloud base of the plume with entrainment rate $\lambda$. ${\lambda }_D\left(z\right)$ is the fractional entrainment rate of the updrafts that detrain at height z. It is assumed that ${\lambda }_D\left(z\right)$ decreases monotonically with height. Thus, clouds with $\lambda >{\lambda }_D\left(z\right)$ have no contribution to mass flux at height z. In the Arakawa and Schubert (1974) scheme,$m_b\left(\lambda \right)$ was obtained by applying a closure condition to each cloud type. The ZM introduced a simplifying assumption. Note that the bulk cloud base mass flux, denoted by $M_b$, is given by (16) $M_b={\int }_0^{{\lambda }_0}{m}_b\left(\lambda \right)d\lambda$(16) where λ₀ is the maximum fractional entrainment rate corresponding to the shallowest updraft plume of the bulk cloud. The ZM scheme assumes that the cloud base mass flux for each cloud type is independent of the cloud type. Thus, (17) $m_b\left(\lambda \right)={\displaystyle \frac{M_b}{{\lambda }_0}}$(17)

Substituting Eq. (17) into Eq. (16) and carrying out the integration yield (note that ${\eta }_u={\displaystyle \frac{M_u}{M_b}}$): (18) ${\eta }_u\left(z\right)={\displaystyle \frac{1}{{\lambda }_0\left(z-z_b\right)}\exp \left\{{\lambda }_D\left(z\right)\left(z-z_b\right)-1\right\}}$(18)

To determine ${\lambda }_D\left(z\right)$, consider a cloud type with entrainment rate $\lambda$. The moist static energy h_u of the steady-state updraft satisfies: (19) ${\displaystyle \frac{\mathrm{\partial }{h}_u\left(\lambda \right)}{\mathrm{\partial }z}+\lambda \left(h_u\left(\lambda \right)-\overline{h}\right)=0}$(19) This equation can be solved by iteration to obtain ${\lambda }_D\left(z\right)$ such that the air parcel at height z is neutrally buoyant and saturated. At the cloud base, $h_u$ is the same as the large-scale value ${\overline{h}}_b$. The entrainment and detrainment rates can be determined from ${\eta }_u$ and ${\lambda }_D\left(z\right)$.

$s_u$ and $q_u$ can be determined from $h_u$: (20) $s_u=\overline{s}+{\displaystyle \frac{1}{1+\gamma }\left(h_u-{\overline{h}}^{\ast }\right)}$(20) (21) $q_u={\overline{q}}^{\ast }+{\displaystyle \frac{\gamma }{L\left(1+\gamma \right)}\left(h_u-{\overline{h}}^{\ast }\right)}$(21) where $\gamma \equiv \frac{L}{c_p}\frac{\mathrm{\partial }{q}^{\ast }}{\mathrm{\partial }\overline{T}}$. Once $s_u$ and $q_u$ are determined, $c_u$ can be calculated from (10) and then $l_u$ and $R_u$ can be calculated from (12) and (13).

Saturated downdraft produced by precipitation evaporation is assumed to exist when there is precipitation production in the updraft. It starts at the height of minimum saturation moist static energy. The budget equations for bulk downdrafts are given by (22) ${\displaystyle \frac{\mathrm{\partial }{\eta }_d}{\mathrm{\partial }z}=e_d-d_d}$(22) (23) ${\displaystyle \frac{\mathrm{\partial }{\eta }_d{s}_d}{\mathrm{\partial }z}=e_d\overline{s}-d_d{s}_d+L{e}_d}$(23) (24) ${\displaystyle \frac{\mathrm{\partial }{\eta }_d{q}_d}{\mathrm{\partial }z}=e_d\overline{q}-d_d{q}_d-\rho e_d}$(24) (25) ${\displaystyle \frac{\mathrm{\partial }{\eta }_d{h}_d}{\mathrm{\partial }z}=e_d\overline{h}-d_d{h}_d}$(25) where subscript d denotes downdraft. The downdraft mass flux at the downdraft initiation level is assumed to be proportional to mass flux of updraft.

3 The Closure

The closure assumption is a hypothesis that relates convection to the large-scale or model-resolved-scale processes. The original ZM scheme uses a CAPE-based closure. Since the large-scale temperature and moisture changes in both the cloud layer and the subcloud layer due to convective activity are linearly related to the cloud base mass flux, CAPE change due to convection can be written as (26) ${\left({\displaystyle \frac{\mathrm{\partial }CAPE}{\mathrm{\partial }t}}\right)}_{cu}=-F{M}_b,$(26) where F is the CAPE consumption rate by convection per unit cloud base updraft mass flux. To calculate F, one can prescribe a small value of $M_b$ and compute the CAPE change due to this amount of convection. Then F is obtained by dividing the CAPE change by the prescribed value of $M_b$. Alternatively, note that the rate of CAPE change due to convection can be represented in terms of convective tendencies of temperature and moisture in both the free troposphere and the PBL, which are linearly proportional to $M_b$. One can compute convective CAPE change per unit cloud base mass flux, and thus F and it depends on the large-scale thermodynamic profiles and the cloud model. In the ZM scheme, the latter approach is used.

The closure condition is that the CAPE is removed at an exponential rate by convection with a characteristic adjustment time scale τ: (27) ${\left({\displaystyle \frac{\mathrm{\partial }CAPE}{\mathrm{\partial }t}}\right)}_{cu}=-{\displaystyle \frac{CAPE-CAP{E}_0}{\tau }.}$(27)

Thus (28) $M_b={\displaystyle \frac{CAPE-CAP{E}_0}{\tau F},}$(28) where $\tau$ is the relaxation timescale and is set to 2 h in CAM5. $CAP{E}_0$ is the threshold value of CAPE for triggering convection, and is set to 70 J/kg.

b Revisions to the ZM Scheme

There have been a number of updates to the original ZM scheme. Some are revisions to elements in the original scheme, others are additions to enhance the scheme.

1 Diluted CAPE

The CAPE calculation in the original ZM scheme assumes that the air parcel, when lifted, first follows a dry adiabat, and after reaching the lifting condensation level it follows a pseudo-moist adiabat, with no entrainment dilution. This decouples the tropospheric humidity, particularly in the lower troposphere, from convection. As such, it overestimates the frequency of convection, especially when the lower troposphere is dry. To remedy this, Neale et al. (2008) introduced entrainment dilution into the CAPE calculation. It is assumed that the parcel’s air mass is doubled every kilometre from the entrainment of the environmental air as it rises to the neutral buoyancy level. Because of the dilution, the CAPE values can be very small in dry environment even when the undilute CAPE is large (Zhang, 2009), thereby reducing the frequency of fictitious occurrence of convection.

2 dCAPE-Based Closure

The CAPE-based closure often produces too frequent convection, even with diluted CAPE. In addition, since CAPE is tied too closely to the boundary layer temperature and moisture, which are strongly influenced by surface turbulent fluxes, convection in the GCM often occurs too early during the day, leading to the wrong diurnal cycle (Dai, 2006). Zhang (2002, 2003) examined the relationships between convection and the large-scale thermodynamic fields using field observations from the U.S. DOE Atmospheric Radiation Measurement (ARM) programme at the Southern Great Plains site and the Tropical Ocean–Global Atmosphere Coupled Ocean-Atmosphere Response Experiment (TOGA COARE) intensive observation period in the tropical western Pacific. He found that under both tropical maritime and midlatitude continental conditions, convection is highly correlated with the large-scale generation of CAPE in the free troposphere (referred to as dCAPE), and proposed to use dCAPE as closure. dCAPE was first used by Xie and Zhang (2000) as a trigger condition for convection.

The dCAPE closure is based on the fact that CAPE as defined in Eq. (1) is determined by the difference between the parcel’s virtual temperature and that of its environment. Thus, the time rate of change of CAPE is given by the difference of the vertical integral of the time rate of change of the parcel’s virtual temperature and that of the environment. Observations show that CAPE change with time is dominantly (∼90% or more) governed by contributions of the parcel’s virtual temperature, which for undilute parcels is entirely determined by temperature and moisture changes in the PBL. In other words, the contribution to CAPE change from the free troposphere is negligible compared to that from the PBL, that is, (29) ${\displaystyle \frac{\mathrm{\partial }CAP{E}_{ft}}{\mathrm{\partial }t}=-{\int }_{p_t}^{p_b}{R}_d{\displaystyle \frac{\mathrm{\partial }{\overline{T}}_v}{\mathrm{\partial }t}d\ln p\approx 0}}$(29)

This forms the basis of the dCAPE-based closure. Since the net CAPE change can be represented as the sum of changes due to convection and those due to non-convective (i.e. large-scale) processes, formally Eq. (21) becomes (30) ${\left({\displaystyle \frac{\mathrm{\partial }CAP{E}_{ft}}{\mathrm{\partial }t}}\right)}_{cu}=-{\int }_{p_t}^{p_b}{R}_d{\left({\displaystyle \frac{\mathrm{\partial }{\overline{T}}_v}{\mathrm{\partial }t}}\right)}_{cu}d\ln p=-{\left({\displaystyle \frac{\mathrm{\partial }CAP{E}_{ft}}{\mathrm{\partial }t}}\right)}_{ls}={\int }_{p_t}^{p_b}{R}_d{\left({\displaystyle \frac{\mathrm{\partial }{\overline{T}}_v}{\mathrm{\partial }t}}\right)}_{ls}d\ln p$(30)

This is well supported by ARM observations (Fig. 24) during its summer 1997 field campaign. The x-axis represents ${\left(\frac{\mathrm{\partial }CAP{E}_{ft}}{\mathrm{\partial }t}\right)}_{ls}$ and the y-axis represents ${\displaystyle \frac{\mathrm{\partial }CAP{E}_{ft}}{\mathrm{\partial }t}}$.

Fig. 24 Scatter plot of the CAPE change due to the ambient virtual temperature change vs. CAPE change due to large-scale forcing from advection and radiative cooling during the ARM summer 1997 IOP. Triangles are for convective periods, crosses and dots are for non-convective periods, the latter of which are for CIN < −100 J/kg. From Zhang (2002).

Since ${\left(\frac{\mathrm{\partial }CAP{E}_{ft}}{\mathrm{\partial }t}\right)}_{cu}$ is proportional to cloud base mass flux, we have the final expression for the dCAPE-based closure: (31) $M_b=-{\displaystyle \frac{1}{K}{\int }_{p_t}^{p_b}{R}_d{\left({\displaystyle \frac{\mathrm{\partial }{\overline{T}}_v}{\mathrm{\partial }t}}\right)}_{ls}d\ln p}$(31) where (32) $K={\int }_{p_t}^{p_b}\left[\left(1+0.608\overline{q}\right)\left(-\eta {\displaystyle \frac{\mathrm{\partial }s}{\mathrm{\partial }p}}\right)+0.608\overline{T}\left(-\eta {\displaystyle \frac{\mathrm{\partial }\overline{q}}{\mathrm{\partial }p}+d_u\left(q_s-\overline{q}\right)}\right)\right]d\ln p$(32) is similar to F in Eq. (18) except it is the CAPE consumption rate in the free troposphere per unit cloud base mass flux. The integrand in the equation is nothing but virtual temperature change of the convection environment due to convective heating and drying.

This closure significantly improved the simulation of precipitation in the single column model of the NCAR CCM3 compared to the original ZM scheme (Fig. 25). It also produced a much better simulation of Madden Julian oscillation (Zhang & Mu, 2005) in the NCAR CCM3, as well as reduced the double ITCZ bias in the NCAR CCSM3 (Song & Zhang, 2009; Zhang & Song, 2010; Zhang & Wang, 2006).

Fig. 25 Precipitation time series observed (blue line) during the ARM 1997 IOP at the SGP site and simulated by a single column model using the original ZM scheme (dashed line, top) and the revised ZM scheme with dCAPE closure (dashed line, bottom).

c Addition of New Functionalities

1 Convective Cloud Microphysics

In the ZM scheme the formation of precipitation is treated through a tunable conversion coefficient c₀ in the bulk cloud model (Eq. (13)). To more accurately represent the microphysical processes that are involved in the conversion from cloud water to rainwater, Song and Zhang (2011) developed a 2-moment convective microphysics scheme, which was later incorporated in the ZM scheme in the NCAR CAM3 (Song et al., 2012). Instead of folding all microphysical processes into a single tuning parameter, the scheme explicitly treats mass mixing ratio and number concentration of four hydrometeor species (cloud water, cloud ice, rainwater and snow). The microphysical processes considered include cloud droplet activation, freezing, cloud ice nucleation, autoconversion, self-collection, collection between hydrometeor species, and sedimentation (Fig. 26).

Fig. 26 Schematic showing the microphysical processes represented in the convective microphysics parameterization of the ZM scheme.

For given convective updraft mass flux (${\eta }_u$) and detrainment rate ($d_u$) from the ZM scheme, the diagnostic budget equations for mass mixing ratio $q_x$ (in units of kg kg⁻¹), and number concentration $N_x$ (in units of kg⁻¹), where the subscript x denotes cloud water, cloud ice, rainwater, or snow in saturated updrafts, can be written as: (33) ${\displaystyle \frac{\mathrm{\partial }{\eta }_u{q}_x}{\mathrm{\partial }z}=-d_u{q}_x+{\sigma }_u{S}_x^q}$(33) (34) ${\displaystyle \frac{\mathrm{\partial }{\eta }_u{N}_x}{\mathrm{\partial }z}=-d_u{N}_x+{\sigma }_u{S}_x^N}$(34) where ${\sigma }_u$ is the fractional area occupied by convective updrafts. $S_x^q$ and $S_x^N$ are the source/sink terms for $q_x$ and $N_x$, respectively, per unit cloud area. They include autoconversion of cloud water/ice to rain/snow; accretion of cloud water by rain; accretion of cloud water, cloud ice, and rain by snow; homogeneous and heterogeneous freezing of cloud water/rain to form ice/snow; Bergeron-Findeisen process; fallout of rain and snow; condensation/deposition; self-collection of rain drops; self-aggregation of snow; activation of cloud condensation nuclei or ice nucleation. Contributions from each of these source/sink terms are parameterized in Song and Zhang (2011). The vertical profiles of hydrometeors within the updrafts are obtained by integrating (33) and (34) from cloud base upward. The detrainment of cloud ice and liquid water (both mass mixing ratio and number concentration) is then added to the grid-scale cloud microphysics scheme.

2 Stochastic Convection

In the ZM scheme, the convective mass flux represents an ensemble mean from a population of convective clouds in a statistical sense. As the model resolution increases, the stochasticity kicks in. To represent the stochasticity in convection, the Plant and Craig (2008) stochastic convection scheme was incorporated into the ZM scheme, such that the probability of the cloud base mass flux is determined by the stochastic scheme while the expected mean of the probability distribution of mass flux in the stochastic scheme is determined by the closure of the ZM scheme. In the stochastic convection scheme, the mass flux of a cloud is assumed to follow an exponential distribution, and thus the average number of clouds having a mass flux between m and m + dm is given by (35) $d\overline{n}\left(m\right)=〈N〉p\left(m\right)dm={\displaystyle \frac{〈N〉}{〈m〉}{e}^{-{\displaystyle \frac{m}{〈m〉}}}dm.}$(35) where $〈N〉$ is the ensemble mean number of clouds and is given by $〈N〉=〈M〉/〈m〉$. $〈M〉$ is the ensemble mean cloud base mass flux and is given by the closure of the ZM scheme, and $〈m〉$ is the mean mass flux of a cloud, which is prescribed as a tuning parameter. The number of clouds with mass flux m follows a Poisson distribution: (36) $p_{N_m}\left(n\right)={\displaystyle \frac{{〈N_m〉}^n{e}^{-〈N_m〉}}{n!}\mathrm{f}\mathrm{o}\mathrm{r}n=0,1,2,\dots }$(36) where $〈N_m〉$ is the average number of clouds with mass flux m. It follows then (37) $p_{d\overline{n}\left(m\right)}\left(n\right)={\displaystyle \frac{d\overline{n}{\left(m\right)}^n{e}^{-d\overline{n}\left(m\right)}}{n!}}$(37)

For small $d\overline{n}\left(m\right)$, (38) $p_{d\overline{n}\left(m\right)}\left(1\right)\approx d\overline{n}\left(m\right)={\displaystyle \frac{〈N〉}{〈m〉}{e}^{-{\displaystyle \frac{m}{〈m〉}}}dm}$(38) and the probability of more than one cloud (i.e. $p_{d\overline{n}\left(m\right)}\left(2\right)$, $p_{d\overline{n}\left(m\right)}\left(3\right)$, …) is of higher order in $d\overline{n}\left(m\right)$ and can be neglected.

A common bias in many GCMs is that light rain occurs too frequently while the frequency of heavy rain occurrence is severely underestimated. The incorporation of stochastic convection into the ZM scheme in the NCAR CAM5 reduced a major bias in the simulated precipitation intensity in the model (Wang et al., 2016). Figure 27 shows an example of such improvements.

Fig. 27 (a) Frequency distributions of precipitation rate and (b) cumulative contribution from each binned precipitation rate based on daily mean precipitation data. The results are for the global belt of 20°S–20°N from the Tropical Rainfall Measurement Mission (TRMM) observation (black line), CTL (blue line) and EXP (red line). From Wang et al. (2016).

6 Parameterization of shallow moist convection

Shallow moist convection, which includes stratus, stratocumulus, and shallow cumulus clouds, exerts a significant impact on the radiative budget of the Earth and strongly affects the vertical transport of enthalpy, moisture, and momentum in the low troposphere. Many problems associated with the simulated variabilities of tropical convection from diurnal to intraseasonal timescales, including a poor simulation of the Madden-Julian Oscillation (MJO) and the diurnal cycle of precipitation over land, are related to the inadequate treatment of sub-grid scale (SGS) processes of shallow moist convection in numerical models (e.g. Grabowski et al., 2006; Hagos et al., 2010). Much of the uncertainty in climate change projections by general circulation models (GCMs) can also be attributed to the parametric treatment of shallow moist convection in the tropics and subtropics (Bony & Dufresne, 2005; Medeiros et al., 2008). To date, a realistic representation of the SGS processes of shallow convection in large-scale models continues to be a challenge. This section reviews how shallow moist convection is parameterized in the state-of-the-art full physics numerical models.

In large-scale models, the shallow cumulus and stratus/stratocumulus are often treated separately by different schemes or modules partially due to the different cloud processes associated with the two types of clouds and partially due to historical reasons. The development of cumulus convection parameterization can be traced back to the pioneer work by Arakawa and Schubert (1974) who constructed the dynamic view of interaction between cumulus and the environment in which the convection is embedded: Down-scaling, the environment controls the convective activities through stabilization/destabilization processes; Up-scaling, cumulus convection alters the thermodynamic structure of the environment through the subsidence induced by the convective updrafts and the detrainment of heat, water vapour, and cloud hydrometers from updrafts. This allowed Arakawa and Schubert (1974) to develop a mass-flux based cumulus scheme by assuming a quasi-equilibrium between the generation of moist convective instability by environmental processes and the stabilizing effect of cumulus convection on the environment. Since then, many convection schemes based on the quasi-equilibrium assumption have been developed. Although these schemes differ in details, they all have three components: (a) a convective cloud model that determines the vertical redistribution of heat, vapour, hydrometeor, and precipitation; (b) an entrainment and detrainment parameterization that describes how the convection interacts with the environment; and (c) a closure at the convection release height that determines when and at what height the convection initiates. Although the quasi-equilibrium assumption and three-component structure were originally proposed for deep convection parameterization, this strategy has been adopted by some sophisticated shallow convection schemes used in GCMs, such as NCAR CAM (Bretherton et al., 2004a), ECMWF (Kohler et al., 2011; Soares et al., 2004) and Météo France (Pergaud et al., 2009).

The most commonly used cloud model in convection schemes is developed under the mass-flux framework in which the kinematic vertical flux of a generic variable $\psi$ is decomposed into four parts following a top-hat distribution, (39) $\overline{w^{\mathrm{\prime }}{\psi }^{\mathrm{\prime }}}=\sigma {\overline{w^{\mathrm{\prime }}{\psi }^{\mathrm{\prime }}}}^u+\left(1-\sigma \right){\overline{w^{\mathrm{\prime }}{\psi }^{\mathrm{\prime }}}}^e+\sigma \left({\overline{w}}^u-\overline{w}\right)\left({\overline{\psi }}^u-\overline{\psi }\right)+\left(1-\sigma \right)\left({\overline{w}}^e-\overline{w}\right)\left({\overline{\psi }}^e-\overline{\psi }\right)$(39) where $w$ is the vertical velocity, overbar and overbar indexed with superscripts u and e indicate the averages over an area, positively buoyant updraft of convection in the area, and environment in which the convective updraft is embedded, respectively, and $\sigma$ is the fraction of convective updraft of the area. The first and second terms on the right-hand side (RHS) of Eq. (39) represent the fluxes induced by the turbulence within the updraft and environment, respectively. The sum of the third and fourth terms on the RHS of Eq. (39) represents the fluxes induced by the coherent updrafts and the associated environment. Assuming that a generic variable also follows the same top-hat distribution, the mean of the variable may be decomposed as $\overline{\psi }=\sigma {\overline{\psi }}^u+\left(1-\sigma \right){\overline{\psi }}^e$, then, the sum of the third and fourth terms on the RHS of Eq. (39) may be rewritten as, (40) ${\overline{w^{\mathrm{\prime }}{\psi }^{\mathrm{\prime }}}}^u=\sigma \left(1-\sigma \right)\left({\overline{w}}^u-{\overline{w}}^e\right)\left({\overline{\psi }}^u-{\overline{\psi }}^e\right)=M_u\left({\overline{\psi }}^u-{\overline{\psi }}^e\right),$(40) where ${\overline{w^{\mathrm{\prime }}{\psi }^{\mathrm{\prime }}}}^u$ is the kinematic vertical flux induced by the coherent features associated with the updraft and environment, and $M_u=\sigma \left(1-\sigma \right)\left({\overline{w}}^u-{\overline{w}}^e\right)$ is known as the kinematic updraft mass flux. The transport induced by the turbulence within the updraft and environment may be appropriately described by local mixing mechanism. Thus, in numerical simulations the vertical fluxes represented by the first and second terms on the RHS of (39) are often handled by a turbulent mixing scheme that also parameterizes the turbulence in the planetary boundary layer (PBL). In contrast, the mixing induced by convection is fundamentally non-local. Thus, a convection scheme is specifically designed to account for the non-local transport represented by ${\overline{w^{\mathrm{\prime }}{\psi }^{\mathrm{\prime }}}}^u$. Note that not all cumulus clouds are positively buoyant. For shallow convection, the negatively buoyant cumulus can be up to 40% according to observations (Taylor & Baker, 1991) and large eddy simulations (LESs, Siebesma & Cuijpers, 1995; Zhu & Albrecht, 2003). Nonetheless, Siebesma and Cuijpers (1995), Wang and Stevens (2000), and Zhu (2015) showed that the coherent updrafts account for 80–90% of the total fluxes in the shallow cumulus layer, indicating that the decomposition of total flux in an area in terms of the positively buoyant updraft is appropriate for parameterizing the vertical transport induced by shallow cumulus convection. In real practice, since cloud fraction $\sigma$ is usually smaller than 20% for shallow cumulus (Siebesma & Cuijpers, 1995), Eq. (40) is often simplified as, (41) ${\overline{w^{\mathrm{\prime }}{\psi }^{\mathrm{\prime }}}}^u\approx \sigma {\overline{w}}^u\left({\overline{\psi }}^u-\overline{\psi }\right)=M_u\left({\overline{\psi }}^u-\overline{\psi }\right),$(41) where it assumes ${\overline{\psi }}^e\approx \overline{\psi }$, ${\overline{w}}^e\approx \overline{w}\approx 0$, and $\sigma \ll 1$. From Eq. (3), it is clear to calculate the cumulus induced vertical fluxes, updraft mass flux $M_u$ and mean updraft property ${\overline{\psi }}^u$ need to be determined. From the conservation law,$M_u$ and ${\overline{\psi }}^u$ should obey the mass continuity equation and dilution equation, respectively (Betts, 1973), (42) ${\displaystyle \frac{\mathrm{\partial }{M}_u}{\mathrm{\partial }z}=E-D,}$(42) (43) ${\displaystyle \frac{\mathrm{\partial }\left(M_u{\overline{\psi }}^u\right)}{\mathrm{\partial }z}=E\overline{\psi }-D{\overline{\psi }}^u,}$(43) where $E$ and $D$ are the lateral entrainment and detrainment mass flux per unit height. Defining $ϵ$ and $\delta$ as the fractional entrainment and detrainment rate, (44) $E=ϵM_u,D=\delta M_u.$(44)

Inserting Eq. (44) into (42) and (43), it yields, (45) ${\displaystyle \frac{1}{M_u}{\displaystyle \frac{\mathrm{\partial }{M}_u}{\mathrm{\partial }z}=ϵ-\delta ,}}$(45) (46) ${\displaystyle \frac{\mathrm{\partial }{\overline{\psi }}^u}{\mathrm{\partial }z}=ϵ\left(\overline{\psi }-{\overline{\psi }}^u\right).}$(46)

Eq. (46) indicates that the change of updraft property is solely caused by entraining environmental air into the updraft.

The cloud model described by (45) and (46) is not closed. To solve the equations, $ϵ$ and $\delta$ need to be determined. For shallow cumulus convection, both observations from the Barbados Oceanographic and Meteorological Experiment (BOMEX, Esbensen, 1978) and LES (Siebesma & Cuijpers, 1995) show that the entrainment and detrainment rate decreases with height from cloud base to cloud top. One of the methods to parametrically determine $ϵ$ and $\delta$ is the buoyancy sorting approach proposed by Kain and Fritsch (1990) and Bretherton et al. (2004a). It assumes that over a small increment of depth $\delta z$ in the cloud layer, equal parts of updraft and environmental air mass flux ${\eta }_0{M}_u\delta z$ are mixed up laterally, where ${\eta }_0$ refers to a fraction of updraft and environmental air close to the cloud edge that is involved in the lateral mixing. The total mixing mass flux is, therefore, $2{\eta }_0{M}_u\delta z$. But the mixing creates a spectrum of mixtures. Some of the mixtures continue to be positively buoyant and keep rising, but the rest becomes negatively buoyant and detrains from the updraft. Thus, to describe the status of a mixture, one may define a fraction index $\chi$ ranging from 0 for the undiluted updraft to 1 for pure environmental air. The probability distribution function (PDF) of a certain mixture with $\chi$ may be described by $\xi \left(\chi \right)$. For all mixtures, there should exist a critical ${\chi }_c$. Mixtures with $\chi <{\chi }_c$ will remain positively buoyant and those with $\chi \ge {\chi }_c$ will be neutral or negatively buoyant. For a certain ${\chi }_c$, the net change of updraft mass flux $\delta M_u$ over a depth $\delta z$will be, (47) $\begin{array}{rl}\delta M_u=& {\int }_0^{{\chi }_c}2{\eta }_0{M}_u\delta z\chi \xi \left(\chi \right)d\chi -{\int }_{{\chi }_c}^{1}2{\eta }_0{M}_u\delta z\left(1-\chi \right)\xi \left(\chi \right)d\chi \\ =& {\eta }_0{M}_u\delta z{\chi }_c^2-{\eta }_0{M}_u\delta z(1-{\chi }_c{)}^2.\end{array}$(47)

The first and second term in Eq. (47) represent the environmental air entrained in the updraft and updraft air detrained in the environment, respectively. Note that to calculate the integrations in Eq. (47), the PDF $\xi \left(\chi \right)$ has been assumed to follow a uniform distribution, i.e. $\xi \left(\chi \right)=1$, (Bretherton et al. 2004). As $\delta z\to 0$, Eq. (47) becomes, (48) ${\displaystyle \frac{1}{M_u}{\displaystyle \frac{\mathrm{\partial }{M}_u}{\mathrm{\partial }z}={\eta }_0{\chi }_c^2-{\eta }_0(1-{\chi }_c{)}^2=ϵ-\delta ={\eta }_0(2{\chi }_c-1).}}$(48)

For a different PDF function $\xi \left(\chi \right)$, the integration of Eq. (47) will be different and will involve complicated calculations. But as pointed out by Bretherton et al. (2004), no observational evidence is available to show the preference of a specific PDF function to others. Thus, for simplicity here we choose the uniform distribution. From (48), it is clear that the entrainment and detrainment rate in a buoyancy sorting model is determined by the fraction of updraft and environmental air involved with the lateral mixing ${\eta }_0$ and the critical mixture index ${\chi }_c$ that makes the mixture buoyant with respect to the environment. Kain and Fritsch (1990) related ${\eta }_0$ to the inverse of the updraft diameter $R^{-1}$. Bretherton et al. (2004) further scaled $R$ to the height of cumulus convection to obtain ${\eta }_0=\frac{c_0}{H}$, where $c_0$ is am empirical coefficient taken as 15 according to the best estimate from LESs, $H$ is the height of cumulus from the previous model time step.

${\chi }_c$ may be estimated from a linear combination of the updraft and environment following Randall (1980) and Deardorff (1980)’s method. Defining liquid water potential temperature as ${\theta }_l=\theta -{\displaystyle \frac{L\theta }{C_{p}T}{q}_l}$ and total water mixing ratio as $q_t=q_v+q_l$, where $T$ is the temperature; $\theta$ is the potential temperature; $q_v$ and $q_l$ are the mixing ratio of water vapour and liquid water; $L$ is the specific latent heat of vaporization; and $C_p$ is the specific heat of dry air at constant pressure. Let $m_u$, ${\theta }_{lu}$, $q_{tu}$, $m_e$, ${\theta }_{le}$, and $q_{te}$ be mass, liquid water potential temperature, and total water mixing ratio of the updraft and environment, respectively. Not that the environment air is unsaturated, thus, ${\theta }_{le}={\theta }_e$ and $q_{te}=q_{ve}$. The liquid water potential temperature and total water mixing ratio of a mixture will be, (49) ${\theta }_{lf}=\left(1-\chi \right){\theta }_{lu}+\chi {\theta }_{le},$(49) (50) $q_{tf}=\left(1-\chi \right)q_{tu}+\chi q_{te},$(50) where $\chi =\frac{m_e}{m_u+m_e}$ is the mixing fraction index with 0 being undiluted updraft and 1 being pure environmental air, and subscript “$f$” indicates the mixture. Now we want to find ${\chi }_s$ so that the mixture is just saturated mixture, i.e. (51) $q_{tf}=q_{sf}\left({\theta }_{lf}\right)=q_{sf}\left[{\theta }_{lu}+{\chi }_s\left({\theta }_{le}-{\theta }_{lu}\right)\right],$(51) where $q_{sf}\left({\theta }_{lf}\right)$ is the saturated mixing ratio of the mixture at ${\theta }_{lf}$. Since $\left({\theta }_{le}-{\theta }_{lu}\right)\ll {\theta }_{lu}$, Eq. (51) may be expanded as the Taylor series. Taking terms of O(${\theta }_{le}-{\theta }_{lu}$), it yields, (52) $q_{tf}=q_{sf}=q_s\left({\theta }_{lu}\right)+{\chi }_s{\displaystyle \frac{\mathrm{\partial }{q}_s}{\mathrm{\partial }{\theta }_l}\left({\theta }_{le}-{\theta }_{lu}\right).}$(52)

Within the range of $[{\theta }_{le},{\theta }_{lu}]$, ${\displaystyle \frac{\mathrm{\partial }{q}_s}{\mathrm{\partial }{\theta }_l}}$ is nearly a constant, thus, Eq. (52) may be rewritten as, (53) $q_{tf}=q_{sf}=q_s\left({\theta }_{lu}\right)+{\chi }_s\left[q_s\left({\theta }_{le}\right)-q_s\left({\theta }_{lu}\right)\right].$(53)

Combining Eq. (15) with Eq. (50), one obtains, (54) ${\chi }_s={\displaystyle \frac{q_{tu}-q_s\left({\theta }_{lu}\right)}{\left[q_{tu}-q_{te}\right]-\left[q_s\left({\theta }_{lu}\right)-q_s\left({\theta }_{le}\right)\right]}.}$(54)

For the mixture that is just buoyant with respect to the environment, one needs to find, (55) ${\theta }_{vf}={\theta }_{ve},$(55) where ${\theta }_{vf}={\theta }_f\left[1+ϵq_{sf}-\left(q_{tf}-q_{sf}\right)\right]$ is the virtual potential temperature of the mixture, $ϵ=R_d/R_v$ is the ratio of gas constant of dry air to water vapour, ${\theta }_{ve}$ is the virtual potential temperature of the environment. ${\theta }_f$ may be expressed as ${\theta }_f={\theta }_{lf}+\pi {\displaystyle \frac{L}{C_p}\left[q_{tf}-q_{sf}\right]}$, where $\pi ={\displaystyle \frac{\theta }{T}={\left({\displaystyle \frac{p_0}{p}}\right)}^{\left(R_d/C_p\right)}}$, $p_0$ is a reference pressure, then, one can theoretically find a ${\chi }_c$ that makes the mixture just buoyant with respect to the environment by solving Eq. (55), but it is not a trivial mathematical problem. In practice, one may utilize the approximate linear relationship among ${\theta }_{vf}$, ${\theta }_{vu}$, and ${\theta }_{ve}$ to provide a good estimate of ${\chi }_c$ as illustrated by Fig. 28, (56) ${\chi }_c={\chi }_s\left(1-{\displaystyle \frac{{\theta }_{ve}-{\theta }_{vfs}}{{\theta }_{vu}-{\theta }_{vfs}}}\right).$(56)

Fig. 28 Virtual potential temperatures for fractional mixtures of updraft and environmental air, where ${\chi }_s$ and ${\chi }_c$ denote the fraction where the mixture is just saturated and just buoyant with respect to the environment. ${\theta }_{vfs}$ is the virtual potential temperature of the mixture that is just saturated.

With an appropriate parameterization of entrainment and detrainment rate, the convection system described by (45) and (46) now can be closed by providing the low boundary conditions at the cumulus base. It describes how convection is initiated. The cloud base may be naturally defined at the lifting condensation level (LCL). However, not all plumes that reach LCL can overcome the convection inhibition (CIN) so that they can develop freely due to their buoyancy. Only those plumes at the cloud base (or LCL) that have vertical velocities greater than a critical value, $w>w_c$, can penetrate the CIN and reach the level of free convection (LFC), from there, they may develop into positively buoyant updrafts. From energy balance perspective, $w_c$ is determined by CIN as, $w_c=\sqrt{2aCIN}$, where $a$ is the virtual mass coefficient, usually taken as $a=1$ (Simpson & Wiggert, 1969), and CIN, following its definition, may be calculate as the enclosed area by the vertical profiles of virtual potential temperature of the environment and an air parcel that undergoes an adiabatic process.

The vertical velocity at the cloud base may be estimated by assuming that the distribution of $w$ follows a Gaussian function, (57) $f\left(w\right)={\displaystyle \frac{1}{\sqrt{2\pi \overline{w^{{}^{\mathrm{\prime }}2}}}}\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \frac{w^2}{2\overline{w^{{}^{\mathrm{\prime }}2}}}}\right),}$(57) where $\overline{w^{{}^{\mathrm{\prime }}2}}$ is the vertical velocity variance. The fraction of updraft and kinetic updraft mass flux at the cloud base are the integration of$f\left(w\right)$ and $wf\left(w\right)$ from $w_c$ to infinity, respectively, (58) ${\sigma }_{cbs}={\int }_{w_c}^{\mathrm{\infty }}f\left(w\right)dw={\displaystyle \frac{1}{2}\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{c}\left({\displaystyle \frac{w_c}{\sqrt{2\overline{w^{{}^{\mathrm{\prime }}2}}}}}\right),}$(58) (59) $M_{u,cbs}={\int }_{w_c}^{\mathrm{\infty }}wf\left(w\right)dw=\sqrt{{\displaystyle \frac{\overline{w^{{}^{\mathrm{\prime }}2}}}{2\pi }}}\exp \left(-{\displaystyle \frac{w_c^2}{2\overline{w^{{}^{\mathrm{\prime }}2}}}}\right).$(59) $\overline{w^{{}^{\mathrm{\prime }}2}}$ may be scaled to either boundary layer mean turbulent kinetic energy (TKE), $\overline{w^{{}^{\mathrm{\prime }}2}}={\text{β}}_e{\overline{e}}_{pbl},$ or Deardorff convective velocity scale, $\overline{w^{{}^{\mathrm{\prime }}2}}={\text{β}}_w{w}_{\ast }$ (Stull, 1988), depending on the specific turbulent mixing scheme used in a model, where ${\text{β}}_e$ and ${\text{β}}_w$ are the empirical coefficients. The cloud base mass flux parameterization described here is based on the maritime trade-wind cumulus convection in which the sub-cloud layer buoyancy production is negligible and buoyancy production in the cloud layer is the main mechanism for maintaining the shallow cumulus layer. However, for continental shallow cumulus, the strong surface forcing causes the sub-cloud layer buoyancy production to be much larger than that in a maritime condition. Zhu and Albrecht (2003) showed that the buoyant production in the sub-cloud layer alone can be sufficient for maintaining a forced cumulus layer. In this case, the CIN and $w_c$ calculation may be adjusted to allow for the forced shallow cumulus in overland conditions (Zhu & Bretherton, 2004).

Compared to the updraft mass flux, the thermodynamic properties of updraft at the cloud base are relatively easy to determine. Observations show that the sub-cloud layer below the cloud base is commonly well mixed, thus, variables conserved for moist adiabatic process, such as, ${\theta }_l$ and $q_t$, at the cloud base may be taken their values in the surface layer with a high accuracy. Once the low boundary condition is specified, (45) and (46) can be solved upward till the height where buoyant updraft disappears, which is considered as the top of cumulus convection. For shallow cumulus convection, precipitation is commonly considered to be less important than that of deep convection. Thus, shallow convection schemes either omit the precipitation process or simply remove the amount of liquid water in the buoyant updraft in excess of an arbitrary valve value as precipitation. Lastly, it should be noted that while nearly all mass-flux cumulus schemes use the same cloud model expressed by (45) and (46), schemes (e.g. Bretherton et al., 2004a; Deng et al., 2003a, 2003b; Pergaud et al., 2009; Soares et al., 2004) do differ in model closure and the ways to determine the entrainment and detrainment rate. Interested readers may refer to specific papers for details.

Unlike a cumulus convection scheme, which works as a standalone physics module in large-scale models, there is usually no separate parameterization scheme to account for the effect of stratus and stratocumulus. In GCMs, such as NCAR CAM the parameterization of vertical transport induced by stratus and stratocumulus is incorporated in the turbulent mixing scheme (e.g. Grenier & Bretherton, 2001; Lock et al., 2000) while in ECMWF stratocumulus are mainly generated by the shallow convection scheme and maintained by radiatively driven turbulent mixing (Kohler et al., 2011). There are two scientific reasons behind this unified parameterization of stratus/stratocumulus and turbulent mixing. First, stratus and stratocumulus are intimately involved with boundary layer turbulent mixing processes and cannot be artificially separated. As illustrated by Fig. 29, in addition to the wind shear and surface sensible and latent heating often known as the bottom-up forcing, the cloud-top longwave radiative cooling and evaporative cooling due to the entrainment of the unsaturated free-atmosphere in the clouds provide an important top-down forcing for driving and maintaining the turbulence in the boundary layer topped by stratocumulus clouds. Thus, a unified parameterization allows for a tight coupling between clouds and turbulence. Second, it is relatively easy to represent the internal stratification or decoupling of stratocumulus topped boundary layer and drizzling process if there is any in a unified framework. In order for a conventional turbulent mixing scheme to appropriately represent the transport processes in the boundary layer topped by stratus/stratocumulus clouds, three modifications or improvements need to be done. First, since the dominant source of TKE in the cloud-topped boundary layer is the buoyancy production in the cloud layer, to appropriately represent the energetics, dynamics, and internal stratification, the turbulent mixing scheme needs to be formulated using the moist thermodynamics based on the variables conserved for moist reversible adiabatic processes, such as ${\theta }_l$ and $q_t$. These schemes are often known as the “moist” turbulent mixing schemes in contrast to the “dry” schemes developed based on the non-conservative variables.

Fig. 29 Schematic vertical structure and physical processes of the stratocumulus topped boundary layer.

Second, a key process that governs the evolution of the stratus/stratocumulus topped boundary layer is the cloud-top entrainment. However, conventional turbulent mixing schemes cannot appropriately account for cloud-top entrainment in large-scale simulations due to the coarse model resolution and lack of consideration of cloud-top radiative cooling and evaporative cooling. A common method to solve this problem is to implement an explicit entrainment parameterization at the diagnosed cloud top to overwrite the turbulent mixing parameterization there. For example, in the NCAR CAM “moist” TKE turbulent mixing scheme (Bretherton & Park, 2009), the eddy exchange coefficients at the cloud top determined by the turbulent mixing scheme are overwritten by the one determined in terms of the cloud-top entrainment rate (60) $K_{h,m}^E=w_E{\mathrm{\Delta }}_{E}Z,$(60) where $K_{h,m}^E$ is the eddy exchange coefficient for heat and momentum at the diagnosed cloud top, respectively; ${\mathrm{\Delta }}_{E}Z$ is the depth of the cloud top entrainment zone, which is usually taken as the vertical model layer in which the diagnosed cloud top resides. $w_E$ is the cloud-top entrainment rate, which may be parameterized by considering a combined effect of boundary layer TKE (or Deardorff convective velocity scale $w_{\ast }$ if TKE is not available), buoyancy discontinuity across the entrainment zone, cloud-top radiative cooling, and evaporative cooling. For example, in the Met Office United Kingdom (METO), and ECMWF models, the entrainment rate is parameterized as (e.g. Lock et al., 2000), (61) $w_E=A_1{\displaystyle \frac{V^3}{\mathrm{\Delta }B}+A_2{\displaystyle \frac{\mathrm{\Delta }F}{\mathrm{\Delta }B},}}$(61) where $V$ is a diagnosed turbulence velocity scaled either to TKE or $w_{\ast }$, $\mathrm{\Delta }B$ is the buoyancy jump across the entrainment zone, $\mathrm{\Delta }F$ is the cloud-top radiative flux divergence, and $A_1$ and $A_2$ are empirical coefficients. NCAR CAM incorporates $\mathrm{\Delta }F$ into the calculation of $V$, and thus, takes $A_2=0$. In this way, the processes that govern the cloud-top entrainment, such as turbulent mixing, radiative cooling, and evaporative cooling, are explicitly included in the cloud top entrainment rate parameterization.

Higher-order closure has been proposed for parametrization of the PBL in climate models (e.g. Bogenschutz et al., 2013; Bougeault, 1981; Golaz et al., 2002; Guo et al., 2015; Lappen et al., 2010; Lappen & Randall, 2001a, 2001b, 2001c; Mellor & Yamada, 1974; Miyakoda & Sirutis, 1977). Higher-order closure provides more detailed representation of PBL processes, and can predict the internal vertical structure of the PBL or a shallow cumulus layer. However, higher-order closure model needs finer vertical grids and correspondingly shorter time steps, which has hindered its application in climate models. Nevertheless, application of higher-order closure has been shown to improve many aspects of the basic state climate, such as the transition of stratocumulus to trade wind cumulus in the subtropical oceans (Bogenschutz et al., 2013).

Finally, the “moist” turbulent mixing scheme and the explicit parameterization of cloud-top entrainment discussed previously were designed to obtain the best result if the entrainment zone can be realistically resolved. However, large-scale models do not have a sufficient vertical resolution to resolve the sharp inversion that caps the stratus/stratocumulus topped boundary layer as shown in Fig. 29. To relax the problem, Grenier and Bretherton (2001) proposed numerical techniques to reconstruct the sharp inversion from the model vertical grids. Interested readers may refer to their paper for details.

7 Summary and discussions

Convective parameterization is the long-lasting bottleneck of global climate modelling and counts for half of the variance of climate sensitivity among the IPCC models. In the past 200 years, numerous observational, theoretical and modelling studies have been conducted on atmospheric convection, which provide the foundation for convective parameterization.

Convective parameterization is not simply a model-resolution problem. The non-hydrostatic high-resolution global cloud-resolving models (Kodama et al., 2015; Fig. 2c) and climate model with super-parameterization (Randall et al., 2016; Fig. 2b) still have large biases in climatological mean precipitation, which are comparable to those in climate models using traditional convective parameterization (Huang et al., 2018; Fig. 2a). Therefore, convective parameterization will be needed for a long time.

In general, the observed convective systems have four components: the always-entraining convective updrafts, unsaturated convective downdrafts, mesoscale updrafts and mesoscale downdrafts (Fig. 6d). The observed convection has self-suppression mechanisms caused by entrainment in convective updrafts, surface cold pool generated by convective downdrafts, and a warm and dry lower troposphere created by mesoscale downdrafts. The post-convection environment is often characterized by “diamond sounding”, which suggests an over-stabilization of the atmosphere rather than barely returning to the neutral state. The pre-convection environment for the future events is then characterized by slow moistening of the lower troposphere forced by surface moisture convergence and other mechanisms. The over-stabilization and slow moistening make the convection events episodic, and decouple the upper troposphere from the boundary layer, making the state-type convective quasi-equilibrium invalid.

The vertical structure of tropical deep convection has three regimes: continental, warm pool, and tropical convergence zone, with the height of heating maximum decreasing from the first one to the third one. Oceanic convection has larger stratiform precipitation fraction (40-50%) than continental convection (30-40%). Convective momentum transport is a leading term in the zonal momentum budget of the Walker circulation in both the deep convection region and shallow convection region.

Right now, the most prominent problem is that the convection schemes are too easy to be turned on, which is linked to double-ITCZ, overly weak MJO variance, precocious diurnal precipitation maximum, and missing stratocumulus clouds. We proposed the following directions for future research:

1. Remove undiluted convective updrafts and improve treatment of entrainment.

2. Add surface convergence trigger and other triggers (e.g. density current, gravity wave).

3. Change saturated convective downdrafts to unsaturated convective downdrafts, which will generate stronger decrease of boundary layer moist static energy and enhance stabilization.

4. Add mesoscale effects. The mesoscale updraft will warm up the upper troposphere and compensate any cold biases caused by suppressing convective updraft height and intensity. The mesoscale downdraft will help to generate the post-convection diamond sounding and make convection more episodic.

5. Add convective momentum transport.

6. Improve the coupling with the dynamics by possibly including some form of convective memory.

7. Develop unified parameterization that can represent PBL, stratocumulus clouds, trade wind cumulus clouds, shallow convection and deep convection in one scheme.

Traditionally, single-column models have been widely used in convection scheme development (Bechtold et al., 1996, 2000; Christensen et al., 2018; Dal Gesso et al., 2015; Davies et al., 2013; Ghan et al., 2000; Grabowski et al., 2006; Guichard et al., 2004; Lenderink et al., 2004; Randall et al., 1996; Svensson et al., 2011; Wing et al., 2018; Woolnough et al., 2010; Wyant et al., 2007; Xie et al., 2002, 2005; Zhang, 2013; Zhu et al., 2005). They are often driven by the constrained variational analysis of sounding array budgets (Zhang & Lin, 1997). However, such datasets are available only from a limited number of field experiments in selected locations and time periods. Here, we propose a new strategy for convection scheme development (Fig. 30):

1. Use reanalysis-driven GCM experiments such as the assimilation runs in weather prediction centres and the decadal prediction runs in climate modelling centres (Boer et al., 2016; Meehl et al., 2014; Merryfield et al., 2020; Smith et al., 2019). This will provide the best possible large-scale environment, and global interactive test of the scheme.

2. Use satellite simulators to directly evaluate the key characteristics such as the convective cloud-top distribution and stratiform precipitation fraction (Bodas-Salcedo et al., 2011; Hashino et al., 2013; Klein & Jakob, 1999; Masunaga et al., 2010; Matsui et al., 2009; Pincus et al., 2008; Roh et al., 2017; Roh & Satoh, 2018; Voors et al., 2007; Webb et al., 2001). In particular, the TRMM and GPM simulators will allow comparison with cloud-top distribution from precipitation radars. Other new diagnostics packages will also reduce the workload of the modellers (Eyring et al., 2016, 2019, 2020).

3. It is important to evaluate the lifecycle of convection and its feedback on the large circulation such as convectively coupled tropical waves and the MJO.

4. The convection scheme development will interact with satellite observations, field experiments and CRM/LES simulations.

5. Convection schemes should be tested and developed for the competitive 10–1 km resolution range where there is benefit from both grid-scale and parametrized representation of convective transport.

Fig. 30 Proposed strategy for convection scheme development. Schematic of convective system is from Zipser (1977). Photos courtesy of NASA and NOAA.

Convective parameterization is one of the most difficult problems in atmospheric sciences and requires long-term tedious work. Right now, most of the modelling centres are short of manpower working on convection scheme development, and the modern “bean counting” academic environment often prevent careful, deep work. It would be great if some fundings are secured to support a limited number of highly qualified researchers to focus on convection scheme development, which will reward global climate modelling in the long run.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

Jialin Lin and Taotao Qian were supported by NSF [grant number AGS-1347132]. Part of this work was done by TQ when working in the Department of Geography, The Ohio State University. Guang Zhang was supported by NSF [grant number AGS-2054697]. Ping Zhu is supported by NOAA [grant numbers NA16NWS4680029 and NA18NWS4680057] and NSF [grant number AGS-1822238].