A Parameterization for Cloud Organization and Propagation by Evaporation-Driven Cold Pool Edges

Abstract

When the negatively buoyant air in the cloud downdrafts reaches the surface, it spreads out horizontally, producing cold pools. A cold pool can trigger new convective cells. However, when combined with the ambient vertical wind shear, it can also connect and upscale them into large mesoscale convective systems (MCS). Given the broad spectrum of scales of the atmospheric phenomenon involving the interaction between cold pools and the MCS, a parameterization was designed here. Then, it is coupled with a classical convection parameterization to be applied in an atmospheric model with an insufficient spatial resolution to explicitly resolve convection and the sub-cloud layer. A new scalar quantity related to the deficit of moist static energy detrained by the downdrafts mass flux is proposed. This quantity is subject to grid-scale advection, mixing, and a sink term representing dissipation processes. The model is then applied to simulate moist convection development over a large portion of tropical land in the Amazon Basin in a wet and dry-to-wet 10-days period. Our results show that the cold pool edge parameterization improves the organization, longevity, propagation, and severity of simulated MCS over the Amazon and other different continental areas.

Introduction

Cold pools are formed by negatively buoyant air arriving at the surface in the cloud downdrafts, spreading out horizontally as a gravity currents. These drafts usually originate in the lower to mid-level parts of the convective cells and are driven by rainfall evaporation. The denser and energetic gravity currents can lift moist and warm conditionally unstable environmental air producing new convective cells downstream. Additionally, cold pools are not only able to trigger new convective cells but also upscale them into larger mesoscale convective systems (MCS), if certain environmental conditions are present, such as strong vertical wind shear (Doswell, 2001). Cold pools and wind shear have fundamental roles in the organization, longevity, propagation, and severity of MCS in any part of the Earth (Houze, 2018; Rotunno et al., 1988). The MCSs, which cover hundreds to thousands of km², significantly impact the global scale circulation, energy budget, hydrological cycle, and population safety.

In the last decade, parameterizations that account for the sub-grid scale effects of the cold pools have been developed. A prominent example is the density current parameterization described in Grandpeix and Lafore (2010), developed to be coupled with any convection parameterization. The parameterization includes prognostic equations for the temperature and moisture in the sub-grid scale cold pools and has its own formulations for the trigger function and the closure. One of the great of benefit from this scheme was the realistic shift produced in the diurnal cycle of precipitation in GCM LMDz, as shown in Rio et al. (2009). Del Genio et al. (2015) describe a cold pool parameterization based on the combination of the elements from previous works. The authors then applied the parameterization to test the sensitivity of the Goddard Institute for Space Studies GCM in simulating the interaction of the large-scale circulation and the Madden–Julian oscillation (MJO) phenomena. More recently, new formulations have been proposed for GCMs with coarse spatial resolutions unable to resolve convection (Rooney et al., 2022) and for convection permitting ones (Hirt & Craig, 2021). We propose here to focus on designing as simple as possible parameterization, looking at essential aspects, from observational and modeling perspectives, of the function of cold pools or, more precisely, the associated gust front properties in triggering and organizing new convection.

Cloud-resolving model simulations, as described by Feng et al. (2015), can be also another source of information for the design. The authors performed simulations of cold pools over the warm tropical ocean. An important aspect of this work was the separated analysis of the thermodynamical settings in two environments: the center and the edge of the cold pools. Robust statistical evaluation of the model simulations confirmed that the environment at the edge of cold pools is more favorable for developing new convection. Along the edges of the cold pools, in contrast to the centers, the temperature depressions are typically less than −1 K, and the moist static energy (MSE) is, on average, 1.6 kJ kg⁻¹ higher. The edges also show clearly higher convective available potential energy (CAPE), lower convective inhibition (CIN), and ∼25% lower LCL height than the center.

A relevant aspect that emerges from CRM/LE simulations is the manifestation of rings of higher MSE around the edge of the cold pools, which are primarily associated with a positive anomaly of water vapor (Feng et al., 2015; Fuglestvedt & Haerter, 2020; Langhans & Romps, 2015; Tompkins, 2001; Torri et al., 2015). The origin of the water vapor anomalies in the rings of moist enclosing the cold pools is still unclear. For example, Tompkins (2001) described the formation of water vapor rings by the rainfall evaporation in the sub-cloud layer during the initiation phase of the deep convection. Then, the occurrence of downdrafts pushes the moist band at the edges of the cold pool. The author concludes by stating that the temperature recovery is faster than the moisture at the edges, so the MSE positive anomaly is primarily due to the excess water vapor. Other authors (e.g., Langhans & Romps, 2015) found that water vapor comes majority from the latent heat fluxes from the surface ocean. With an alternative focus, Torri and Kuang (2016) showed that despite the surface fluxes being the primary source of the water vapor, the associated positive anomalies are equally produced by the rain evaporation and the surface fluxes.

In terms of the capability of triggering new convection by density currents, Torri et al. (2015) advocated that besides the mass convergence in the gust front, which provides an efficient mechanical lift for the environmental air parcels to reach their level of free convection (LFC), the existence of the water vapor rings also plays an essential role. This role is called thermodynamical lift and is associated with the larger buoyancy of the air parcels convecting from the rings.

Some effects of cloud organization by cold pools are described as follows. Kurowski et al. (2018) used large-eddy simulations (LESs) to study the transition of shallow-to-deep convection over land, focusing on the role of cold pools and the interplays with the surface fluxes. A simulation with interactive surface fluxes compared to the control one where the cold pools are eliminated shows a substantial increase in the near-surface temperature and moisture variability. A more profound vertical development of convection was also revealed, associated with a reduction of lateral entrainment. They found that with cold pools, a systematic reduction of about 30%–40% of the lateral entrainment rate (see Figure 16 therein) is observed. The authors explain the less diluted updrafts and taller plumes around the cold pools because of enhanced convection organization. Also, Feng et al. (2015) showed that intersecting cold pools are much more efficient in triggering convection with enhanced updraft velocities than isolated ones. However, the main reason for these convection characteristics from intersecting cold pools was the closer spacing between the convective elements (aggregation) and the consequent much moister local environment. Similarly, Becker et al. (2018) found that aggregated updrafts have 40% higher entrainment because of stronger wind flows. However, the larger entrainment does not translate to a more substantial dilution of the updrafts because they are embedded in a moist shell. Colin et al. (2019) analyzed the role of sub-grid scale boundary layer structures chasing for temporal convective memory and its role in spatial cloud organization. The authors showed that the most critical sub-grid boundary layer structures in that regard are the MSE structures, which are dominated by moisture structures.

On the other side, observations show that strong transients in thermodynamic quantities usually accompany the passage of cold pools associated with precipitating lines. For example, Silva Dias et al. (2002) reported a drop of about 10 K and 3 g kg⁻¹ in near-surface equivalent potential temperature (θe) and water vapor mixing ratio, respectively, during the passage of a precipitating line in the southwest of Amazon. On the other hand, a gust front wind speed of up 20 m s⁻¹ was observed. The authors also reported a decline of about 40 hPa in the height of the lifted condensation level (LCL). A comprehensive report of those events, dividing mesoscale systems versus unorganized convective cells, was presented by Schiro and Neelin (2018) for a site near the city of Manaus in the Amazon. On average, both isolated and MCSs events were characterized by similar decreases of ∼1.5 g kg⁻¹ in specific humidity, ∼4.2 K in temperature, and ∼9.6 K in θe, with an increase of 4.2 m s⁻¹ in wind speed at the surface. The estimated lifetime of the cold pools was about 1 hr for isolated events and larger than 3 hr for MCSs.

Tobin et al. (2013) point out the need for the process of convective aggregation, and the consequent mesoscale organization, to be included in convection parameterizations. The authors justify it to account for couplings between moisture, large-scale circulation, and radiative balance not represented by the traditional parameterizations used by low resolution global circulation models (GCM). However, although we are now (10 years after) in an era of the so-called Global Storm-Resolving Models running with a few kilometers of nominal horizontal resolution (Freitas et al., 2020; Stevens et al., 2019), the need for improving the convection parameterization framework is still valid, as we discuss below.

Given the broad spectrum of scales of the atmospheric phenomenon involving the gravity currents and the MCS life cycle, computational limitations will require parameterizations for processes not explicitly resolved at the model computational grid. This is especially true for Earth System or climate models where the need for low spatial resolution due to the computational constraints will still exist for a long time, and even more for applications involving probabilistic forecast via ensemble of global model runs. Herein lies the purpose of this work: how to account for the sub-grid scale interplay between gravity currents, mainly their associated gust fronts, wind shear, and MCS formation and propagation in the context of coarse-resolution atmospheric models?

At this point, we will set our framework to design the parameterization:

1. The parameterization will focus only on the impacts of the gust front and the associated moist rings. It will not explicitly consider the sub-grid scale region populated by the cold pools.

2. The scheme should prognostic the MSE anomalies, including their transport and environmental mixing, at the outer region enclosing the gust front.

3. The gust front propagation speed should be diagnostic.

4. The MSE anomalies should be applied to encourage the triggering of new convection and make it stronger.

5. The entrainment rate of the convection parameterization should be sensitive to the presence of cold pools in a way that it further helps the cloud organization. However, this aspect will be left for future work and will not be discussed here.

Section 2 starts with cloud-resolving simulations to help the parameterization development. Then we introduce the parameterization formulation for roles of the edges of cold pools driven by convective scale downdrafts. Then, in Section 3, we will discuss its application in simulating Amazon Basin's squall lines with an atmospheric model lacking enough spatial resolution to explicitly resolve their formation and development. Section 4 expands the parameterization evaluation by applying it to different land regions and environments. We will discuss our results and the performance of the cold pool edge parameterization in Section 5.

Developing the Parameterization Using Cloud-Resolving Simulations to Help the Design of the Parameterization

Convection resolving modeling (CRM) was performed for a large portion of the Amazon Basin to explicitly simulate the developing of cold pools over this tropical land. The simulations were produced with the BRAMS regional atmospheric modeling system (Freitas et al., 2017) configured with 500 m horizontal grid spacing. On the vertical, model grid spacing started at 70 m, increasing to 500 m at the upper levels, with the model top of about 20 km in height. The simulations used the non-hydrostatic, Boussinesq compressible (Cotton et al., 2003) dynamic core, integrated with the Runge-Kutta third order time scheme with the third (vertical) and fifth order (horizontal) advection operators (Wicker & Skamarock, 2002). The PBL parameterization was the Mellor–Yamada level 2.5 scheme (Mellor & Yamada, 1982). The cloud microphysics applied the hybrid single and double moment scheme from Thompson and Eidhammer (2014), with all the water species (vapor, cloud and ice, and hydrometeors) transported using a monotonic advection formulation (Freitas et al., 2012). Solar e terrestrial radiative transfers used the RRTMG (Iacono et al., 2008) radiation scheme. The JULES surface scheme (Moreira et al., 2013) handled the atmosphere-surface fluxes exchange.

The horizontal domain comprised a 1,000 × 1,000 km² region enclosed by latitudes approximately 11–2 S and longitudes of 62.5–53.5 W. The ECMWF Global Reanalysis ERA-5 (Hersbach et al., 2019), 0.25-degree horizontal resolution, provided the initial and boundary conditions for the zonal and meridional wind, air temperature, water vapor mixing ratio, and geopotential height. The initial condition for soil moisture and temperature originated from ERA-5 as well. The model time integration was 48 hr starting from 00 UTC 21 MAY 2017. For the analysis of the model results, the first four hours are discarded for spin-up, as well as any data distant up to 100 model grid boxes from the domain borders. Figure 1 shows true-color images from the Aqua/MODIS sensor for the afternoon (∼14 UTC) of 21 MAY 2017. The image covers a large portion of the Amazon Basin, representing the model domain. Note the presence of an easterly predominant ambient wind in the lower-mid troposphere.

Typical model results for a time in the late afternoon are shown in Figure 2. This figure introduces the near surface (∼34 m above local terrain) moist static energy (panel A), water vapor mixing ratio (panel B), and air temperature (panel C), all as instantaneous values simulated for 19:40 UTC 21 MAY 2017. Panel D shows the accumulated precipitation in the next 30 min (from 19:40 to 20:10 UTC). The simulation shows a large population of cold pools spread over this model sub-domain, denoted by relatively colder and drier regions (panels C and B). Those regions are also typically enclosed by rings of comparatively moist air (panel B) and higher moist static energy (panel A). In panel D, we call attention to a convective system organized as a straight precipitating line traveling to the southwest of the domain. The following 30 min of accumulated precipitation (panel D) seems to be developed over an area of a strong gradient of temperature (panel C), just over the positive anomalies of water vapor and MSE (panels B and A).

We now look for average characteristics of the MSE deviation from a large-scale value using the results from this CRM. The main objective is to determine the differences between the MSE anomalies in the region just outer the cold pool and the anomalies in the cold pool itself. To do that, we examine with some detail the dynamics and thermodynamics of a cold pool, as shown in Figure 3. Panel (A) contains three fields: in shaded is the instantaneous (21:40 UTC 21 MAY 2017) near-surface MSE anomaly relative to a large-scale value defined as the instantaneous MSE horizontal model domain average. In contour, are shown the previous (green) 30 mn accumulated precipitation (i.e., accumulated from 21:10 to 21:40 UTC) and the subsequent 30 mn accumulated precipitation (from 21:40 to 22:10 UTC, black). Observe that the previous precipitation falls mainly in the region with negative MSE anomalies, while the subsequent one falls mostly over the area with positive MSE anomalies. Also, the MSE anomalies typically are in the range of ∼ ±10 kJ kg⁻¹, as we more firmly demonstrate below.

Panels (B) and (C) show the instantaneous near-surface vertical velocity and the pressure gradient anomaly, respectively. The black dashed curves denote the position of an arc of upward motion (panel B) which coincides with the position a strong horizontal pressure gradient (panel C). In this way, it broadly characterizes the position of a gust front.

Based on these typical characteristics, the following filter for automatic discrimination between areas of the positive (prone to forge new convection) and negative (prone to inhibit it) MSE anomalies is proposed. The filter follows the criteria:

1. if the near-surface vertical velocity is >0.1 m s⁻¹ and the subsequent 30 mn accumulated precipitation is >4 mm => may indicate a region of positive MSE deviation or the region of the moist rings and just outer region of the cold pools

2. if the vertical velocity is <−0.1 m s⁻¹ and the previous 30 mn precipitation is >4 mm => may indicate a region of negative MSE deviation or the cold pool inner region.

Figure 4 shows the probability density function (PDF) for the MSE anomalies obtained from applying the filter to the CRM data. Two curves are shown: the PDF MSE anomaly associated with the moist rings (outer of the cold pool, in red color) defined by criterium 1 and with the inner of the cold pools (blue) defined by criterium 2. The MSE for the entire model domain (not shown) oscillates between approximately ±10 kJ kg⁻¹ around the mean 344.7 kJ kg⁻¹ and with a standard deviation of 3.3 kJ kg⁻¹. Perhaps remarkably, are the similarities between the PDFs associated with the blue and red regions, in which the main disparate difference is the displacement of its centers by about 4–5 kJ kg⁻¹. Feng et al. (2015) found similar results from their simulations over tropical warm ocean (see their Figure 7b) with similar PDFs associated with the center and edge of the cold pools and displaced by about 1.6 kJ kg⁻¹. Here we found a displacement about twice which can be conceivable as the simulation covered a tropical land during the transition from wet to dry season of Amazonia.

Table 1 shows some statistical indicators of the two PDFs of Figure 4. Observe that the two MSE anomaly PDFs associated to the just outer and in the inner of cold pools are 2.0 and −2.1 kJ kg⁻¹, respectively, with similar standard deviation of about 2.8–3.0 kJ kg⁻¹.

Table 1 Statistical Indicators of the MSE Anomalies PDF Shown in Figure 4

Region	Mean (kJ kg−1)	STD (kJ kg−1)	Number of data	Percentage of the total data
Red (outer)	2.0	2.8	1.1 106	0.1%
Blue (inner)	−2.1	3.0	1.8 106	0.2%

The Cold Pool Edge Parameterization

We designed the cold pool edge parameterization to work with the Grell-Freitas convection parameterization (Grell & Freitas, 2014, hereafter GF) using this scheme's readily available physical quantities. A short description of the GF is given here and elsewhere for the reader's benefit. The scheme follows the mass flux approach and includes the convective scale downdraft, besides a formulation for the cloud updraft. Initial entrainment and detrainment rates are prescribed but are re-scaled after the normalized mass fluxes profiles are determined. The scheme is scale-aware and has a set of different closures to determine the updraft mass flux at the cloud base. The mass flux of the downdraft at its initiation level is proportional to the updraft mass flux after applying the scale-awareness factor. The proportionality constant and the process of the downdraft development are more described at the end of this section.

We start with the definition of a new three-dimensional positive-definite scalar field named buoyancy-excess ($${\beta }_{x}$$) and given by: [Image Omitted. See PDF]where $${H}_{d}$$ and $$\tilde{H}$$ are the downdraft and environment moist static energy (MSE), respectively. We will further assume that the detrained β_x by downdrafts scales the typical local spatial variability of the near-surface environment MSE in a domain populated by cold pools. In other words, based on the PDF's similitudes of the MSE anomalies shown in Figure 4 and Table 1, we assume that the negative value of the MSE deficit ($${\beta }_{x}$$) in the inner of the cold pool approximates the MSE excess in the moist rings.

Several equations were developed for the estimation of gust front horizontal speed (V_gf) relative to the ambient flow. In most cases, the equations are expressed in terms of the gradient of some quantity across the gust front, where the quantity can be the air density (Charba, 1974), density potential temperature (Tompkins, 2001), virtual potential temperature (Grandpeix & Lafore, 2010; Hirt & Craig, 2021) or the potential temperature (Reif et al., 2020; Rotunno et al., 1988). To express in terms of the prognostic buoyancy-excess quantity, we follow the Arakawa and Schubert (1974, see also Plant, 2010) formulation for the buoyancy, assuming the cold pool is nearly air saturated. In this case, the gust front horizontal speed in estimated by [Image Omitted. See PDF]where κ is the internal Froude number (adopted as 0.9), $$\tilde{T}$$ is the environmental air temperature, g is the acceleration due to gravity, c_p the specific heat of dry air, and z is the cartesian height. The parameter $$\gamma $$ is given by the Equation 57 in Arakawa and Schubert (1974).

Following Reif et al. (2020), we adopt the Equation 3 to estimate the vertical velocity at the leading edge of the cold pool ($${\mathrm{W}}_{\text{gf}})$$: [Image Omitted. See PDF]

The quantity α is the slope angle of the cold pool head and the value of 45⁰ is set following the approximated mean value of this parameter as shown in Reif et al. (2020, Table 1 therein). The last quantity D is the depth of the cold pool and is defined by 2 km or the vertical height where $${\beta }_{x}$$ < 0.1 kJ/kg, whichever occurs first.

To determine the horizontal propagation of the gust front, we first calculate the mean cloud layer horizontal velocity by the following equation: [Image Omitted. See PDF]where p₁ = 900 hPa, p₂ = 600 hPa, $${\left(u,v\right)}_{env}$$ is the horizontal environment wind and p is the atmospheric pressure. Then, the two-dimension propagation velocity of the gust front is estimated by: [Image Omitted. See PDF]where the last term on the RHS accounts for propagation through environments with an ambient wind (Parker, 1996).

Herein cold pool edge parameterization proposes that the following prognostic equation governs the time evolution of β_x: [Image Omitted. See PDF]where:

1. The first term on the RHS is the two-dimension advection by the horizontal propagation wind. The advection follows the forward upstream of the second-order scheme (Tremback et al., 1987) and is applied only in the horizontal by setting the Cartesian vertical velocity to zero. We assume that the cold pool spreading is essentially horizontal following the terrain.

2. The second term is the three-dimension diffusion. In BRAMS, as in most the atmospheric moldels, unresolved transport is parameterized using K-theory, which involves evaluating covariances as the product of an eddy mixing coefficient (EMC) and the gradient of the transported quantity.

   • 2a.

     For the proposed horizontal grid spacings for the application of the cold pool edge parameterization, the horizontal EMC is determined by selecting the larger value between a scheme that connects diffusivity coefficients to the deformation rate of the horizontal wind (Smagorinsky, 1963), and the minimum value necessary to maintain model numerical stability (which is proportional to grid_lenght^4/3). Moreover, for the β_x prognostic, we reduce the EMC to 10% of its actual value, leaving a minimum of lateral spread of the cold pool edge.

   • 2b.

     In vertical, any of the planetary boundary layer (PBL) schemes available in the Brazilian developments on the Regional Atmospheric Modeling System (BRAMS, Freitas et al., 2017) may be used.

The third term represents the source of β_x with δ_d, the downdraft detrainment mass flux, which is readily available in the GF convection parameterization. The saturated downdrafts parameterization in GF is described in Grell (1993). Like the updraft, a one-dimensional, steady-state cloud model with entrainment and detrainment is formulated for downdrafts. The evaporation of condensate in downdraft is limited to keep it saturated and its mass flux at the level of free sinking is proportional to the updraft mass flux at the cloud base, as stated before. The proportionality factor is 1–P_ef, with the last term being the precipitation efficiency given in terms of the vertical shear of the horizontal wind as formulated by Fristch and Chappel (1980). In this case, higher wind shear environments will result in lower P_ef and, consequently, stronger cold pools.

The fourth term represents the sink associated with processes other than mixing with the environment, which act to destroy the cold pools (e.g., surface fluxes). The tunable parameter τ_cp is the cold pool lifetime and is 1–3 hr long. We currently apply 2 hr based on Schiro and Neelin (2018) for Amazon Basin. However, it can be made dependent on the underneath surface type (land or ocean) to better represent the role of surface fluxes.

Figure 5 anticipates the results of the model simulations described in the upcoming Section 2.3 to give a glimpse to the readers of the spatial distribution of $${\beta }_{x}$$. It shows a typical output of the near-surface buoyancy-excess simulated for 21 UCT 16 FEB 2014, which corresponds to the summer season's late afternoon (∼5 p.m.) in the Amazon. The units are kJ kg⁻¹, and the maximum simulated values range from 10–20 kJ kg⁻¹. The model domain horizontal average is ∼2.5 kJ kg⁻¹, if one masks out values below 1 kJ kg⁻¹. Perhaps not surprisingly, the spatial distribution of the $${\beta }_{x}$$ does not resemble an individual closed-cell gust front of a cold pool. Instead, it resembles a combination of several gust fronts resulting in a large mesoscale one covering a contiguous area on a much ampler spatial scale, as described in Houze (2004, 2018). The sharp distribution and the larger values of the $${\beta }_{x}$$ are also noticeable in the leading edge and in the direction of propagation of the mesoscale cold pool, as denoted by the yellow to red colors. This is a desirable property as it may represent the water vapor anomalies rings in the leading edge of a cold pool gust front. We will further discuss this aspect in next section.

How to Connect the Buoyancy-Excess With the Nearby and Near-Future Convection?

In this section, we design how the cold pool edge parameterization disturbs the environment and convection properties emulating the roles of the mass convergence and the moist rings of the gust front of a cold pool. The three formulations are presented in the next three subsections.

Modifying the Cloud Work Function of the Downstream Convecting Air Parcels

Although the GF CP has a set of closures to determine the updraft mass flux at the cloud base, in this work, we limit the application of the gust front parameterization to only one. We choose the closure that combines the stability removal formulation with the non-equilibrium diurnal cycle originally proposed by Bechtold et al. (2014). Freitas et al. (2021, see Section 2.4 therein) describe the implementation of this combined closure in GF. The closure is based on the ratio of the cloud work function (A) by the convective adjustment time scale (Γ) to determine the updraft mass flux at the cloud base, being it proportional to $$\sfrac{A}{{\Gamma }}$$. The cloud work function is defined by [Image Omitted. See PDF]where z_b and z_t are the height of the cloud base and cloud top, respectively, $$g$$ is the gravity, c_p the specific heat of dry air, $${Z}_{u}$$ is the normalized updraft mass flux, $$\tilde{T}$$ is the environment air temperature, and $${H}_{u}$$, $${\tilde{H}}^{\ast }$$ are the updraft and environment-saturated moist static energy, respectively. The parameter $$\gamma $$ was already introduced in Equation 2. The presence of the moist rings enclosing the cold pool is taken into account in the boundary condition for the MSE of the updraft in the propagation direction, serving as an additional source of buoyancy for the convecting air parcels: [Image Omitted. See PDF]So, $${\beta }_{x}$$ in Equation 8 represents the sub-grid scale rings of positive anomalies of MSE. In this way, the parameterization further assumes that the MSE excess is entirely in the form of water vapor, and so, the boundary condition for the water vapor mixing ratio in the updraft at the cloud base is given by [Image Omitted. See PDF]Where $${q}_{u}$$ and $$\tilde{q}$$ are the updraft and environment water vapor mixing ratio, respectively, and L_v the latent heat of vapourization.

The approach of adding a surplus with respect to the environment to the convecting air parcel in a convection parameterization to account for sub-grid scale variability of thermodynamical quantities is not new. For example, based on Jakob and Siebesma (2003), the ECMWF IFS (IFS, 2021) convection parameterization adds a temperature ∆T_u and moisture ∆q_u excess to the test air parcel depending on the surface sensible and latent turbulent heat fluxes.

Modifying the Convective Adjustment Time Scale

The convective adjustment time scale ($${\Gamma }$$) in the selected GF CP closure is given by [Image Omitted. See PDF]where $${w}_{u}$$ is the sub-grid scale updraft vertical velocity. The gust front parameterization impacts this time scale by setting the vertical velocity at the leading edge of the cold pool (W_gf, see Equation 2) as the boundary condition for the diagnostic equation of sub-grid scale updraft vertical velocity in GF: [Image Omitted. See PDF]

Equation 11 implies on stronger convection (larger updraft mass flux at the cloud base) by a resulting shorter time scale $${\Gamma }$$ in Equation 10.

Adding a New Trigger Function

The last implementation in GF related to the cold pool edge parameterization is an optional trigger function based on the kinetic energy ($${E}_{k}=\frac{1}{2}{W}_{\text{gf}}^{2}$$) of the air parcels at the leading edge of the gust front. In this case, deep convection is allowed in a model column if [Image Omitted. See PDF]where A_cin is the cloud work function given by Equation 7 but with the integral being computed from the air parcels source level to the cloud base. A_cin scales the energy necessary for the air parcels to overcome a potential convective inhibition layer capping the PBL. In practical, the kinetic energy in Equation 12 is given by $${E}_{k}=\frac{1}{2}{\left[\max \left({W}_{\text{gf}},{w}^{\ast }\right)\right]}^{2}$$, with $${w}^{\ast }$$ the boundary layer convective velocity scale.

Grandpeix and Lafore (2010) apply this kinetic energy in their closure for the mass flux at the cloud base. In our approach, $${E}_{k}$$ is used only as a condition for allowing deep convection to be developed in each model column. We assume that the initial kinetic energy of the convecting air parcels at the edge of the gust front is somewhat accounted for in the shorter convective adjustment time scale used in our closure, as described in Section 2.3.2.

Model Results Over the Amazon Basin Methodology

We performed limited-area model simulations over a large portion of the Amazon basin with the BRAMS modeling system using the cold pool parameterization coupled with the GF convection parameterization. The model was configured with a nominal horizontal grid spacing of 8 km × 8 km covering an area of approximately 2,800 km × 2,200 km. The chosen grid spacing is at the upper range of the gray zone for deep convection and is close to the current horizontal resolution of the GCMs currently used for the deterministic, mid-range numerical weather forecast in the major operational centers, for example, the ECMWF. On the vertical, model grid spacing started at 90 m, increasing to 750 m at the upper levels, with the model top of about 20 km in height. The physics configuration follows the description of Section 2.1 with two differences as follows. We used the cloud microphysics WSM 5-class single moment (Hong et al., 2004) scheme for the grid-scale clouds, and for the sub-grid scale convection, the GF CP was applied, being configured with deep and shallow plumes (Freitas et al., 2018, 2021).

The ECMWF Global Reanalysis ERA-5 (Hersbach et al., 2019), 0.25-degree horizontal resolution, provided the initial and boundary conditions for the zonal and meridional wind, air temperature, water vapor mixing ratio, and geopotential height. The boundary conditions were applied during the BRAMS time integration using the nudging (Newtonian relaxation) technique on the lateral and top of the model grid. The initial condition for soil moisture and temperature originated from ERA-5 as well. The experiments involved two periods of 2014, spanning 10 days each, representing the wet season (15–24 FEB) and the transition from dry to wet (01–10 OCT) of the Amazon basin. For each day, the model time integration was 48 hr, starting at 00 UTC, with the first 24 hr being discarded as the time for the spin-up. The results we will discuss here correspond to the last 24 hr (the second forecast day).

Simulated Properties of the Cold Pool Edge Parameterization

Figure 6 introduces the model simulation of a bow-shaped convective system, one of the modes of mesoscale convective organization (Weisman, 2001) called bow echo, from its initiation to its declining phase. As described in Weisman (2001), a bow echo typically begins as a large single convective cell (upper panels of the second column). Because of the higher propagation speed at the north part of the convective cell (see wind vectors in the last column), the cell acquires a bow-shaped configuration. A few hours later, the system transforms into the so-called comma-shaped echo. Observe that during its mature and declining phases, the vertical velocity at the leading edge of the gust front is well above 5 m s⁻¹, reaching up to 10 m s⁻¹, while the maximum magnitude of the horizontal propagation speed is above 22 m s⁻¹. The second column also shows the bow echo moved about ∼330 km in 6 hr, giving a mean propagation speed of ∼15 m s⁻¹. However, the first column of Figure 6 indicates that the development of a single cell into a bow echo and, finally, a comma-shaped echo is not simulated in the control experiment.

Figure 7 further illustrates how the cold pool parameterization improves the GF CP in producing or enhancing cloud organization. At least three MCSs archetypes are reasonably well simulated by the model experiment with the cold pool parameterization (Panel A). Following Schumacher and Rasmussen (2020, see Figure 1 and references therein), the simulated convective systems in Figure 7a resemble the MCSs types identified as training line/adjoining stratiform, bow echo, and squall line. Each of one MCSs is composed of convective and stratiform precipitation parts, and the respective propagation velocity vectors are also shown. In contrast, Figure 7b shows that the associated events simulated in the control experiment, at the same time and places, are not clearly organized and recognizable as the ones simulated with the cold pool parameterization. Also, with the cold pool edge parameterization, large precipitation rates (>10 mm hr⁻¹) are more frequent, mainly of a sub-grid scale origin. At the same time, the precipitation rates are slighter without it, with a substantial part of the rainfall coming from the grid-scale cloud microphysics (not shown).

The cold pool edge parameterization is statistically characterized with the information displayed in Table 2 in terms of its main properties as follows: the buoyancy-excess (Equation 1), the vertical velocity at the leading edge of the cold pool (Equation 2), the gust front speed (Equation 3) and the cold pool propagation speed (Equation 5). The first two columns show the experiment and its corresponding time period. Subsequent columns inform the domain area average, standard deviation and the upper 5% range of values, considering a filter of β_x > 1 kJ kg⁻¹. The mean β_x ranges from 2.2 to 3.2 kJ kg⁻¹ with 5% of the model domain having values above 5.5–8.8 kJ kg⁻¹. Similarly, the mean W_gf varies between 3.3 and 3.9 m s⁻¹, with 5% with values greater than 6.4–9.3 m s⁻¹. Observe that, as expected, the simulated cold pools properties tend to be stronger during the dry-wet period of the year. In general, the max simulated values for $$\left\vert {\left(u,v\right)}_{prop}\right\vert $$ were ∼23 and 30 m s⁻¹ for wet and dry-wet periods, respectively.

Table 2 Statistical Values for the Main Simulated Cold Pools Physical Quantities

Exp	Days of year 2014	βx (kJ kg−1)			Wgf (m s−1)			Vgf (m s−1)			|(u,v)prop| (m s−1)
		Mean	Std Dev	Upper 5%	Mean	Std Dev	Upper 5%	Mean	Std Dev	Upper 5%	Mean	Std Dev	Upper 5%
Wet	15–24 FEB	2.2	1.5	5.5	3.3	1.7	6.4	4.7	2.4	9.1	11.4	4.0	18.0
Dry-Wet	1–10 OCT	3.2	2.3	8.8	3.9	2.7	9.3	5.7	4.3	14.2	13.5	5.1	23.1

The mean simulated β_x (∼2–3 kJ kg⁻¹) is compatible with the Feng et al. (2015) finding that the moist static energy is, on average, 1.6 kJ kg⁻¹ higher along the edges of their simulated cold pools, as we discussed in the introduction section. Also, our simulated mean vertical velocity around 3–4 m s⁻¹ is compatible with the observational and modeling data shown in Reif et al. (2020, see their Table 1 and Figure 17). For the propagation speed, Garstang et al. (1994) reported Amazon coastal squall lines moving across the central Amazon Basin at speeds of 50–60 km hr⁻¹, which are consistent with the mean (41–49 km hr⁻¹) and the maximum (65–83 km hr⁻¹) speeds informed in Table 2.

Effects on the Mean Precipitation

We shall discuss how the cold pool edge parameterization impacts the spatial mean model precipitation over the Amazon basin. Figure 8 introduces the spatial distribution of the 10-day mean rainfall for the two analyzed periods: wet and dry-wet cases. The two upper rows show the rainfall observation-based estimation (GPM-IMERG, Huffman et al., 2019) and the modeling-based ERA-5 ECMWF reanalysis. The spatial resolution of these two data sets is 0.1 and 0.25°, respectively. The two lower rows depict the correspondent mean precipitation as simulated in the control experiment and with the cold pool edge parameterization. The two BRAMS model simulations (rows 3 and 4) follow the seasonal cycle of the observed precipitation, varying from ∼9 to ∼5 mm day⁻¹, consistent with GPM-IMERG and ERA-5. They also show a realistic spatial distribution of the rainfall compared to the GPM-IMERG. It is also noticed that the introduction of cold pool edge parameterization did not degrade the model results of the control experiment but offered some improvements regarding producing more events with higher precipitation rates. Also, the precipitation occurrence with the CPool ON is more discrete and less continuously spatially distributed, better resembling the GPM-IMERG.

Effects on the Intensity of the Convective Systems

The frequency distribution of the precipitation rate in mm day⁻¹ is shown in Figure 9 for two cases: the 15–16 FEB 2014 (wet, left) and 09–10 OCT 2014 (dry-wet, right) periods. The 1-hourly GPM-IMERG precipitation rate (black lines) is applied as the reference. For both cases, the cold pool edge parameterization intensifies the upper-range distribution of the convective precipitation (not shown). As a result, there is an increase of the number of severe convective systems and impacting the overall frequency of distribution of the total precipitation at higher amounts, as demonstrated in Figure 9.

Impacts on Storm's Propagation

Aspects of storm propagation are described in Figure 10, which corresponds to a Hovmöller diagram with the vertical axis corresponding to the time in UTC and the horizontal axis the longitude of the model domain. In the first row (Figures 10A1 and 10B1) is the reference data based on the GPM-IMERG estimation for the dry and dry-wet cases, 10 days average each. The quantity displayed is the mean diurnal cycle of the total precipitation using the 10 days data for each case and then averaged again over the latitude range from 9 S to 6 N. The storms' propagation across the Amazon basin is usually associated with instability waves that form on the northern coast of Brazil, near the state of Para (Cohen et al., 1995). Most of the propagation is toward to west or southwest of the Amazon basin.

The model results for the control experiment (CPool OFF) and the CPool ON are shown in the subsequent rows. The mean horizontal cold pool edge wind speed is also pictured in the lowest row (white arrows). During the wet season, the wind direction is primarily toward the southwest, while in the dry-wet season, it is mainly toward the west, following the trade wind.

For the wet season (Figure 10A1), GPM-IMERG shows storms predominantly propagating to the southwest, following the wind vector shown in Figure 10A3. Both model simulations for this case (Figures 10A2 and 10A3) capture the main features displayed by the GPM-IMERG product adequately. And so, clearly distinguishing improvements from the cold pool edge parameterization is challenging in this case.

During the dry-wet case (Figure 10B1), storms start to pop-up around 16–18 UTC and propagate for a few hours, predominantly toward to west due to the easterly wind flow (see Figure 10B3 for the wind vector propagation). They are predominantly concentrated in the late afternoon and early night, with a few initiating and propagating late night and morning. The model simulations (Figures 10B2 and 10B3) also are consistent with the timing of the convection and propagation. Both simulations seem to consistently reproduce the main characteristics of the storm's propagation tracks.

Mesoscale Convective Systems Cases

We are now exploring into the model results for four cases of mesoscale convective systems across various continental regions and environments. This is aimed at thoroughly assessing the effectiveness of the cold pool edge parameterization. In Table 3, we present the geographical locations where the studied MCSs occurred (column 1), the initial starting time of the model simulations and the length of its time integration (column 2), the horizontal and vertical spacings of the model grid (column 3), and, finally, the extension of the model domain in the last column.

Table 3 The Mesoscale Convective Systems Cases Studies and the Correspondent Model Configurations

Region	Time of start -integration length	Model grid spacing	Domain extension (# of grid points)
CONUS	12UTC 02 OCT 2014	H: 12 km × 12 km	H: 350 × 300
	48 hr	V: 70–750m	V: 45 (top @ 20.3 km)
Amazon Basin	12UTC 04 OCT 2020	H: 8 km × 8 km	H: 360 × 280
	48 hr	V: 90–750m	V: 43 (top @ 20.6 km)
Rio de la Plata Basin	12UTC 15 FEB 2023	H: 10 km × 10 km	H: 550 × 320
	48 hr	V: 90–600m	V: 49 (top @ 18.6 km)
Equatorial Africa	00UTC 06 AUG 2016	H: 12 km × 12 km	H: 555 × 370
	48 hr	V: 90–750 m	V: 45 (top @ 21.1 km)

Figure 11 introduces the MCS case study that occurred over the continental United States (CONUS) between October 2nd and third, 2014. This event involved the development of an intense line of thunderstorms that originated in the central part of the USA and propagated toward the East Coast. This squall line spanned more than 1200 miles, covering the entire latitudinal extent of the CONUS. The squall line formation was a result of the convergence of cold air from the north with warmer and more humid air from the south.

In Column 1 of Figure 11, the progression of the squall line is displayed using a 2-hr averaged rainfall representation from the GPM-IMERG product. The series starts from approximately the initiation time at 15 UTC on October 2nd (Panel A1) and continues until it approaches the East Coast around 13 UTC on October 3rd (Panel F1).

The corresponding simulations conducted using the BRAMS model are presented in columns 2 and 3. In column 2, the simulations were performed with the control setup where the cold pool edge parameterization was turned off (CPOOL OFF). In column 3, the simulations were run with the cold pool edge parameterization activated (CPOOL ON). Upon visual examination of the model simulations, it becomes evident that the inclusion of the CPOOL ON parameterization enhances the model's performance in terms of both organization and intensity. This improvement is notable from the time of squall line initiation and persists throughout its entire development phase.

Figure 12 provides visual snapshots of the four mesoscale convective system (MCS) case studies outlined in Table 3. These snapshots focus on a specific time when the convective systems are fully mature. As before, the first column of Figure 12 exhibits the 2-hr average rainfall observed from the GPM-IMERG product in mm h⁻¹. The corresponding simulations using the BRAMS model are presented in the subsequent columns. In column 2, the simulations were conducted with the control setup (CPOOL OFF), while in column 3, the simulations utilized the cold pool edge parameterization (CPOOL ON).

The left side of the first column provides information regarding each case's region and the specific time period, indicated by the hour interval and date. This arrangement allows for a direct comparison between the observed rainfall and the model simulations with and without the cold pool edge parameterization.

The first row (panels A) shows a fully developed squall line over the Amazon Basin on 19–21 UTC 04 OCT 2020. This squall line extends by approximately 2,000 km. The maximum GPM-IMERG 2-hr average precipitation rate ranged from around 10 to 15 mm hr⁻¹. In the control experiment (Figure 12A2), the simulation shows a weaker and less organized convection, with maximum precipitation rates ∼5 mm hr⁻¹.

However, the implementation of the cold pool edge parameterization (Figure 12A3) notably enhances both the strength and organization of the squall line. Precipitation rates reaching up to 10–15 mm hr⁻¹ are now evident, and the squall line's definition is much improved, although its length isn't as far-reaching as observed GPM-IMERG.

The second row illustrates an MCS moving westward across the western portion of Equatorial Africa. This is the same case previously investigated by Becker et al. (2021), as seen in Figure 8 of their work. The GPM-IMERG retrieval shows a prominent MCS, spanning approximately 1000 km in length. It extends from around latitude 5°N to 15°N near longitudes 6°W to 7°W. The maximum GPM-IMERG 2-hr averaged precipitation rates ranges between 10 and 20 mm hr⁻¹.

The control experiment (panel B2) attempts to replicate the observed precipitation pattern in terms of location and general form. However, the intensity is notably underestimated, with maximum values reaching only around 2–3 mm hr⁻¹.

Conversely, the experiment with CPOOL ON significantly enhances the representation of that MCS. The spatial organization of the MCS takes the form of a bow-echo, indicating a better simulation of its structure. Additionally, the precipitation rates are improved are closer to the observations, although they still underestimate the peak values.

Interestingly, a similar enhancement in model performance is highlighted in the work by Becker et al. (2021). In their study, they expand the mass flux closure in the IFS convection parameterization to explicitly include the advective moisture tendency (as shown in Equation 2 of their paper). This enhancement is visually evident in their Figure 8, particularly when comparing the panels labeled by “Qadv off” and “Qadv on.” This suggests that incorporating the cold pool edge parameterization produces a similar performance gain achieved by explicitly accounting for advective moisture tendencies in their study.

The third one is the extremely long squall line over the CONUS that took place on 02 and 03 October 2014 discussed the begin of this section. This case is revisited here for the purpose of comparing it with the other MCS cases for an overall performance evaluation and includes a time interval that is not depicted in Figure 11.

Lastly, we present the case of a large MCS that took place over the Rio de la Plata Basin in the central south portion of South America. On February 15th and 16th, 2023, a complex interplay of meteorological systems unfolded, resulting in the development of a large MCS that impacted the southern region of the South American continent. Multiple factors converged to bring about the formation of this system. First, there was a notable transport of warmer and more humid air at the surface, originating from the Amazon region. This phenomenon was linked to the presence of a low-level jet on the east side of the Andes that was active in the area. Besides, at higher levels and over the southern part of the continent, an intense and well-structured trough was associated with the coupling between the Polar Jet and the Subtropical Jet. This resulted in the advection of cyclonic vorticity to lower levels, intensifying a low-pressure system already present over the central north portion of Argentina. This intensification of the low-pressure system, along with the abundance of moisture and high instability, gave rise to the MCS that later moved northward following the upper-level trough, affecting areas of central and northern Argentina and, subsequently, Paraguay.

Figure 12 D1 shows that MCS covering nearly the entire Paraguay region and extending southeastward, crossing the Brazilian state Rio Grande do Sul and reaching the northern portion of Uruguay. The GPM-IMERG 2-hr averaged precipitation rates surpass the 30 mm hr⁻¹ scale, indicating its severity. The subsequent two columns (D2 and D3) display simulations without (D2) and with the cold pool edge parameterization (D3). Even in the CPOOL OFF (D2) simulation, the model captures the general location and structure of the MCS reasonably. However, precipitation rates are significantly underestimated, and the southeastern branch of the MCS is largely absent.

The inclusion of the cold pool edge parameterization (D3) enhances the MCS representation. While precipitation rates remain somewhat underestimated, they are larger, and the organization of the MCS improves, with the southeastern branch better depicted.

Conclusion

We designed and tested a parameterization to account for critical interplays between evaporation-driven cold pools, wind shear, and mesoscale convective systems in the context of atmospheric models with an insufficient spatial resolution to resolve their interactions explicitly. Motivated by previous works (e.g., Mapes & Neale, 2011), a new 3-day scalar field named buoyancy-excess (β_x) was created to assist in parameterizing certain effects of cold pools edges. The source for β_x comes from the downdrafts' detrained deficit of the MSE with respect to the environment. The scalar β_x is horizontally, terrain following transported between grid columns, three-dimensionally mixed with the environmental air, and has a generic sink term to account for cold pool dissipating processes other than the turbulent mixing. The definition of β_x readily allows the estimation of two essential properties of the cold pool gust fronts: the propagation speed and the vertical velocity of the environmental air at its edge. Inspired by CRM simulations and observations, we assumed that β_x scales the typical local spatial variability of the near-surface environment MSE in the presence of cold pools. Our model results show a suitable consistency between β_x and the typical MSE anomalies in the leading edge of the gust fronts (see Section 2.1 and, e.g., Feng et al., 2015), where the environmental air parcels are lifted to their level of free convection forging new updrafts. This consistency holds well in spatial distribution and magnitudes, as shown in Figures 5 and 6 and discussed in Section 3.1.

Assuming this framework, the cold pool edge parameterization clearly improves the organization, longevity, propagation, and severity simulation of Amazon Basin MCSs. At the same time, there is no degradation in the mean spatial and diurnal cycle of modeled precipitation, with some cases even offering improvements. This work also shows in Section 4 (see Figure 12) that the parameterization improves the MCS modeling representation over a set of contrasting continental regions and environments.

We are considering the possibility of exploring and adding more characteristics of the cold pool edge parameterization. Investigating how the environmental entrainment rate responds to cloud organization holds potential as an area of model further improvement. Also, we anticipate that this parameterization could establish a valuable connection between the trimodal plume formulation in our convection parameterization. To illustrate, the gust front could play a role in shaping the development of shallow and congestus plumes near the area enclosing the cold pool. This phenomenon is often observed in remote sensing imagery but is currently not accounted for in our formulation here described. Additional ideas involve incorporating a direct interaction between the surface and the gust front. This approach could potentially lead to a scenario where the duration and strength of the cold pool edge are influenced by the actual surface moisture and temperature. Worth trying is also a suggestion of an anonymous reviewer about the exact definition of the buoyancy excess quantity. The reviewer suggested directly using the buoyancy excess defined as the MSE of downdraft minus the saturation MSE of the environment.

One potential limitation of our current design is that it does not consider the detailed descriptions of individual cold pools. As a result, it cannot capture the effects of the collisions between cold pools, which are known to play a significant role in triggering stronger convective cells, as demonstrated in studies such as Feng et al. (2015). Our future work also will board this aspect.

Lastly, the assessment of the performance of this parameterization will extend to global-scale simulations, incorporating a variety of different types of mesoscale convective systems in additional environments, such as maritime and continental settings, and using different spatial resolutions. This evaluation will also encompass mid-latitude cold fronts, the Madden-Julian Oscillation, the development of haboobs, tropical cyclones, and severe weather events in general. These comprehensive evaluations are part of our planned future research efforts.

Acknowledgments

The first author acknowledges the support from the Brazilian National Council for Scientific and Technological Development (CNPq) process # 309565/2021-9 and from the State of Goiás Research Support Foundation (FAPEG) process # 202110267001027.

Data Availability Statement

The BRAMS data simulations herein shown are available in Freitas (2023a). The code of the Brazilian developments on the Regional Atmospheric Modeling System (BRAMS) applied for the here published model results, including the cold pool edge and Grell-Freitas convection parameterizations, is available in Freitas (2023b). The ERA5 output used in this study (Hersbach et al., 2018) has been provided by the Climate Data Store. The GPM-IMERG data used here has been provided by Goddard Earth Sciences Data and Information Services Center (Huffman et al., 2023).