Bulk cloud microphysical properties as seen from numerical simulation and remote sensing products: case study of a hailstorm event over the La Plata Basin

2.1. Atmospheric model

The BRAMS model (Freitas et al. 2005, 2009) derives from the Regional Atmospheric Modelling System (RAMS, Pielke et al. 1992; Cotton et al. 2003), and includes new features to improve the representation of fundamental physical processes in tropical and subtropical regions (Silva et al. 2012). Recent BRAMS model versions include a set of state-of-the-art physical and chemical parameterisation schemes that extend RAMS original functionalities towards a fully integrated environmental model (Moreira et al. 2013; Freitas et al. 2017). BRAMS has been mainly applied for operational forecasts and research related to severe weather (e.g. Freire et al. 2022), urban heat islands (e.g. Souza et al. 2016), urban and remote air pollution (e.g. Bela et al. 2015), aerosol–radiation–cloud interactions, carbon and water cycles over the Amazon (e.g. Moreira et al. 2013), and volcanic ash dispersion (e.g. Pavani et al. 2016). The model source code and related preprocessors are open source and are available under the GNU General Public License at http://brams.cptec.inpe.br.

The BRAMS cloud microphysics schemes able to account for hail formation include those of RAMS Colorado State University (RAMS/CSU) (Walko et al. 1995; Meyers et al. 1997) and Thompson (Thompson and Eidhammer 2014), both two-moment bulk schemes. The cloud microphysics scheme used in this study (RAMS/CSU) predicts the mixing ratio and number concentration for seven hydrometeor categories: cloud droplets, rain, pristine ice, snow, aggregates, graupel and hail. The general gamma distribution is the basis function used for hydrometeor size in each hydrometeor category (Walko et al. 1995). Cloud droplets and rain are liquid water, but may become supercooled. Pristine ice, snow and aggregates are assumed to be pure ice, whereas graupel and hail are mixed-phase categories (ice only or a mixture of ice and liquid). Ice–rain interactions usually result in most of the collected mass and number concentrations being converted to the hail category; however, since hail is assumed to be a mixed phase hydrometeor, a look-up table for melting is applied to determine the number tendency of hail particles that completely melt into rain. This number concentration tendency is similar to the look-up table for melting of snow and aggregate hydrometeors. Although the hail category is permitted to contain liquid, excess liquid is shed to the rain category, with shedding from hail being based on the work of Rasmussen et al. (1984). Depending on atmospheric conditions, riming of graupel can also serve as a primary source of hail (Meyers et al. 1997). A detailed description of the main BRAMS physics schemes can be found in Freitas et al. (2017).

2.2. Experiment design

Three BRAMS simulations were performed for a thunderstorm event observed in the eastern LPB on 14–15 July 2016. The experimental setup consists of a single domain with 250 × 280 (latitude × longitude) horizontal grid points at a resolution of 6 km. The criteria for choosing a resolution of 6 km were to align, in terms of spatial resolution, the model outputs with satellite- and radar-based observations. For the vertical coordinate, 45 sigma-type levels were set up with grid spacing ranging from 20 to 1000 m, and the model top extending up to 23 km (~33 hPa). Initial and boundary conditions to drive the simulations were based on global data at 0.25° horizontal resolution from the National Centers for Environmental Prediction (NCEP) GFS and CFSv2 models (Saha et al. 2014), and the European Centre for Medium-Range Weather Forecasts (ECMWF) Reanalysis v5 (ERA5) model (Hersbach et al. 2020), each driving a BRAMS simulation. The GFS data are provided at 26 levels from 1000 to 10 hPa, whereas CFSv2 and ERA5 provide meteorological data at 37 levels from 1000 to 1 hPa. Meteorological fields, including temperature, relative humidity, geopotential height, zonal wind and meridional wind, are interpolated to the model grid from the original pressure-level data sets, with boundaries updated every 6 h.

The main physical parameterisations include the Rapid Radiative Transfer Model for General Circulation Model applications (RRTMG) scheme (Iacono et al. 2008) for shortwave and longwave radiation, the Mellor–Yamada scheme (Mellor and Yamada 1982) for turbulence and the RAMS/CSU scheme (Meyers et al. 1997) for cloud microphysics. Cloud droplets are assumed to nucleate from constant CCN vertical profiles fixed at 450 nuclei per kilogram of air – sensitivity tests, which involved changing the CCN concentration to evaluate its impact on ice cloud formation, were previously performed for this particular study (not shown). Overall, new ice particles are formed when CCN is enhanced, whereas hail shows a slight decrease in concentrations with increasing CCN. This result agrees well with that of Camponogara et al. (2018). However, although the BRAMS model developers recommend coupling convection schemes for horizontal resolutions usually coarser than 2 km, other studies suggest that convective schemes should be scale-dependent (Arakawa et al. 2011), with assumptions that may vary with horizontal resolution (Bryan et al. 2003; Grell and Freitas 2014). Sensitivity tests with the standard RAMS convection scheme (Tremback 1990) turned on and off were conducted for verification. However, the simulations with the convection scheme turned on yielded slightly worse outcomes; therefore, they are not presented here for the sake of simplicity. Model results presented in this study derived from the microphysics configuration that best reproduced the observed ice cloud conditions. Static geographical data (e.g. topography, land use and sea surface temperature) and masked surface fields (e.g. soil moisture) needed for the model initialisation can be found at http://brams.cptec.inpe.br/input-data/. Model simulations were integrated over time for a simulation period of 54 h, starting at 18:00 hours UTC on 13 July 2016. The main physical parameterisations and other simulation attributes are listed in Table 1.

2.3. Observational data

Observational data used for the model evaluation include surface accumulated precipitation from MERGE and Climate Hazards group Infrared Precipitation with Stations (CHIRPS) merging techniques, reflectivity fields from a Doppler weather radar (DWR) instrument, cloud optical properties derived from Moderate Resolution Imaging Spectroradiometer (MODIS) retrievals and vertical profiles from sounding data. MERGE combines 0.25° Tropical Rainfall Measuring Mission (TRMM) satellite precipitation estimates with in situ station data to create daily gridded rainfall time series over South America (Rozante et al. 2010). The TRMM values at locations that coincided with the observation locations were interpolated to the TRMM grid points using successive corrections of the Barnes scheme (Barnes 1973; Koch et al. 1983). CHIRPS is a 0.05° land-only precipitation database available at different time scales that incorporates satellite, model and in situ data sets (Funk et al. 2015). Both MERGE and CHIRPS provide 24-h accumulated precipitation data. Although MERGE covers the period from 12:00 to 12:00 hours UTC, CHIRPS does it from 00:00 to 00:00 hours UTC. Precipitation data from CHIRPS lack uncertainty information at the moment. Reflectivity data were collected by a DWR instrument operating in the 2.7–3.0-GHz frequency band with a coverage extending out to a range of 400 km. The DWR is installed on the top of a 1.8-km peak in southern Santa Catarina State (28.13°S, 49.47°W), the second-highest peak within the modelling domain. MODIS is an optical sensor onboard the NASA Terra and Aqua satellites that takes images for the entire Earth’s surface every 1 or 2 days, and acquires data in 36 spectral bands in a wavelength range of 0.4–14.4 μm. Statistical fields from Aqua MODIS L3 daily cloud products at 1.0° resolution (Platnick et al. 2019) are used to evaluate the model skill in reproducing the cloud water path for liquid and ice phases. Sounding data for the Florianópolis station, located off the coast in western Santa Catarina State (27.67°S, 48.55°W), are provided by the Department of Atmospheric Science at the University of Wyoming (https://weather.uwyo.edu/upperair/sounding.html). Fig. 1 shows the locations of sounding and radar stations.

2.4. Model evaluation

A number of statistical parameters can be used to evaluate the performance of atmospheric models, including the correlation coefficient (r), mean bias error (MBE) and root-mean-square error (RMSE); r is a measure of the strength and direction of the linear relationship between simulation and observation, MBE measures the mean difference between simulation and observation, and RMSE represents the square root of the mean squared error between simulation and observation. All three parameters are appropriate for evaluating model performance over multiple time and space scales and can be calculated as follows:

where, P and O represent the predictions and observations respectively; k represents the kth grid point; and N is the total number of grid points. Overbars signify means over site, domain or time.

To compare against the DWR reflectivity observations, the model’s equivalent radar reflectivity factor (Z_e) was estimated following the approach of Gao and Stensrud (2012). Additionally, in order to understand the strengths and weaknesses of the model in terms of the representation of atmospheric water and ice content, the Vertically Integrated Liquid (VIL) and Vertically Integrated Ice (VII) were calculated for the model and radar reflectivity factors, following the formulations of Greene and Clark (1972) and Mosier et al. (2011) respectively. The Echo Top of 20 dBZ (ET20), indicating the highest level of a storm detected by radar, was calculated and evaluated against the corresponding value from the model. Therefore, it is possible to visualise if the vertical development of the storm observed with the DWR system is well represented by the model. The fractional skill score (FSS, Roberts and Lean 2008) was applied for evaluating reflectivity, VIL and VII for each model simulation and radar observation:

where P_F(k) and P_O(k) are the fractional coverages of reflectivity for the kth grid point that exceeds a given threshold, and N is the total number of grid points in the domain. Thresholds of 20 dBZ for reflectivity and 0.5 kg m⁻² for VIL/VII were chosen, along with a radius of influence of 10 grid points. These thresholds were selected based on our goal to determine if the model can accurately reproduce any amount of water (reflectivity and VIL) or ice (VII) caused by convection. In terms of the radius of influence, conducting multiple tests showed that employing 10 grid points yielded results with less noise. Consequently, the overall outcome did not significantly differ when using fewer grid points, although the FSS was smaller in all cases. Furthermore, it is important to note that, due to characteristics of the radar system such as air refraction and Earth curvature, the radar’s beam moves away from the surface as it extends farther from the radar. As a result, valid data are not available at lower levels in areas distant from the radar. Therefore, as VIL and VII are calculated by integrating data across atmospheric levels, the error in these calculations increases in regions far from the radar, resulting in an underestimation of the true values.

3.1. Case-study description

As mentioned in the introduction section, the border region between north-eastern Argentina and southern Brazil has been identified as one of the largest hail-prone areas in the Southern Hemisphere, with high hail fall frequencies during late winter and early spring. The hailstorm event selected for numerical experiments occurred in the eastern LPB on 14–15 July 2016, a month with a number of hail days significantly less than those usually reported in October, the peak hail month (Beal et al. 2020). Fig. 2 shows the evolution of the large-scale atmospheric conditions associated with the hailstorm event. Prior to the hailstorm event (12:00 hours UTC on 12 July 2016), there was strong moisture transport from the Amazon to the south-east of Brazil and the adjacent ocean due to a cyclogenesis taking place at 40°S and 55°W. Cyclogenesis over south-eastern South America is more often observed during austral wintertime (Crespo et al. 2020a, 2020b; de Jesus et al. 2021). This extratropical cyclone moved and developed over the South Atlantic and had an associated cold front that affected the weather over the continent. Circa 00:00 hours UTC on 14 July 2016, the environment over the study region still had a high moisture content (Fig. 2d), which was sustained by the cold front and a low-pressure system that were acting over the continent (20°S, 65°W). The low-pressure system deepened over the continent, giving rise to a new cyclone centred at 30°S and 50°W at 00:00 hours UTC on 15 July 2016 (Fig. 2c). From 00:00 hours UTC on 14 July 2016 onwards the jet stream accelerated while crossing over the Andes at the same time that two intense jet streaks took place near the region of interest: one over the continent (~25°S) and the other over the South Atlantic near the continent (~35°S) (Fig. 2h, i). In this way, the study region was under a jet stream bifurcation, i.e. an upper-level divergence zone corresponding to the polar exit and the equatorial entrance of these two different jet streaks. Over time, the ascent (air motion in the upward direction) intensified, reaching values below −1 Pa s⁻¹ (Fig. 2h). The synoptic environment described here configures an SALLJ event, as previously introduced in Section 1.

Fig. 2.

(a–c) Mean sea level pressure (hPa, contour) and 500–1000-hPa thickness (decametres, dam, shaded), (d–f) horizontal wind (m s⁻¹, vectors) and specific humidity (g kg⁻¹, shaded) at 850 hPa, (g–i) horizontal wind (m s⁻¹, vectors) and vertical velocity (Pa s⁻¹, shaded) at 250 hPa at 00:00 hours (left) and 12:00 hours UTC on 14 July 2016 (centre) and at 00:00 hours UTC on 15 July 2016 (right) derived from ERA5 reanalysis. The red marker indicates the hailstorm location.

Hail fall events during the study period were reported in south-western Paraná State (e.g. Beal et al. 2020) and across various areas in southern Santa Catarina State between the end of 14 July and the beginning of 15 July 2016. South-eastern South America has a little more than half of its storms overnight, with a broad peak between 15:00 and 00:00 hours local solar time (Cecil and Blankenship 2012). Nocturnal convection over LPB is primarily triggered by orographic lifting (Romatschke and Houze 2010), highlighting the important role of the terrain in the development and timing of convection (Rasmussen and Houze 2016). According to the Civil Defence of Santa Catarina State, the hail fall event was registered in nine municipalities of the state: Lages, São Joaquim, Praia Grande, São Miguel do Oeste, Xanxerê, Peritiba, Ipira, Xaxim and Concórdia. Although the Civil Defence agency did not provide valuable information on the hailstorms, such as duration and hail size, it did report significant economic damage in São Miguel do Oeste and Praia Grande, where citizens were rendered homeless (source: https://estado.sc.gov.br/noticias/sc-registra-granizo-em-pelo-menos-nove-cidades/). Most of the affected regions are occupied by small farms and relatively poor people, for whom any hail-related damage causes several problems. The simulated spatial distribution of hail covers locations where hail fall was reported, namely São Joaquim and Praia Grande, both in southern Santa Catarina State (see Supplementary Fig. S1).

3.2. Surface accumulated precipitation

Fig. 3 shows the spatial distributions of the surface accumulated precipitation derived from MERGE, CHIRPS and BRAMS estimates. Fig. 3a, c, e and g correspond to MERGE, BRAMS/GFS, BRAMS/CFSv2 and BRAMS/ERA5, representing the 24-h surface accumulated precipitation during the period from 12:00 hours UTC on 14 July 2016 to 12:00 hours UTC on 15 July 2016. Similarly, Fig. 3b, d, f and h correspond to CHIRPS, BRAMS/GFS, BRAMS/CFSv2 and BRAMS/ERA5, representing the 48-h surface accumulated precipitation during the period from 00:00 hours UTC on 14 July 2016 to 00:00 hours UTC on 16 July 2016. CHIRPS provides accumulated precipitation data from 00:00 to 00:00 hours UTC, requiring a 48-h time period to encompass the entire hail fall event between the end of 14 July and the beginning of 15 July 2016. In contrast, MERGE covers the hail fall event with a 24-h period, from 12:00 hours UTC on 14 July 2016 to 12:00 hours UTC on 15 July 2016. In accordance with observations, rainfall is predicted in the south-eastern part of the modelling domain. However, slight to moderate differences in accumulated precipitation exist between observations and model simulations, with BRAMS/CFSv2 and BRAMS/GFS performing best against MERGE and CHIRPS, respectively. Table 2 summarises domain-wide performance statistics over land. The r between the BRAMS simulations and CHIRPS range within 0.4–0.5, whereas those between BRAMS and MERGE are 0.38–0.69. Notably, the mean bias and RMSE between BRAMS and MERGE are smaller compared to those between BRAMS and CHIRPS. Additional comparisons at the National Institute of Meteorology (INMET) site-specific locations are shown in Supplementary Fig. S2 and Supplementary Tables S1 and S2. Even though areas with maximum accumulated precipitation are well reproduced by the model in most cases, it is clear that BRAMS overestimates rainfall, with high RMSE values probably due to errors in both initial and boundary conditions as well as to model simplifications. For instance, it is initially assumed that cloud droplets nucleate from a prescribed CCN vertical profile, whereas in nature, there exists spatial variability in CCN concentrations due to prior activation in clouds (Zhao et al. 2017). The study conducted by Camponogara et al. (2018), also over the LPB, revealed that the number and intensity of downdrafts and, consequently, their total covered area increase as CCN concentration is enhanced.

Fig. 3.

Spatial distribution of surface accumulated precipitation (mm) during the periods from 12:00 hours UTC on 14 July 2016 to 12:00 hours UTC on 15 July 2016 (left panels) and from 00:00 hours UTC on 14 July 2016 to 00:00 hours UTC on 16 July 2016 (right panels). Precipitation data correspond to (a) MERGE and (b) CHIRPS merging techniques, and (c, d) BRAMS/GFS, (e, f) BRAMS/CFSv2, and (g, h) BRAMS/ERA5 simulations.

3.3. Vertical profiles and stability

The status of a thunderstorm development can be monitored using Skew-T Log-p diagrams produced from soundings and atmospheric models. Fig. 4 shows the vertical profiles of temperature, dew-point temperature and winds derived from the sounding launched at 00:00 hours UTC on 15 July 2016, from the Florianópolis station. Vertical profiles inferred from BRAMS simulations, at the sounding launching time and for the closest grid point to the sounding station, are also shown in Fig. 4. The simulated vertical profiles of temperature agree well with those derived from sounding data. However, all simulations have difficulties in capturing strong changes in dew-point temperature profiles. In terms of winds, the BRAMS/GFS simulation better captures the vertical profile of wind direction. As indicated in previous studies, systematic model errors in representing sharp inversions of temperature and wind shear are often related to the misrepresentation of land-use and stability (Mölders and Kramm 2010). However, the main contributing factor is that vertical profiles from atmospheric models, such as BRAMS, do not properly account for the impact of advection on thermodynamic parameters (Anand and Pal 2023). In addition, distances between model grid-points and stations elevations contribute to biases, as temperature and wind vary with height. The Florianópolis station is located in Santa Catarina Island, a small rugged island which makes the task of modelling the atmosphere over this location challenging for atmospheric scientists.

Fig. 4.

Skew-T diagrams showing the profiles of temperature (red lines, °C), dew-point temperature (green lines, °C) and wind (barbs, knots) derived from the (a) sounding (CAPE = 71.4 J kg⁻¹; CIN = −132.6 J kg⁻¹; 0–8-km bulk shear of 49.8 kt, ~25.6 m s^–1) and the (b) BRAMS/GFS (CAPE = 439.1 J kg⁻¹; CIN = −87.2 J kg⁻¹; 0–8-km bulk shear of 54.4 kt, ~28.0 m s^–1), (c) BRAMS/CFSv2 (CAPE = 279.5 J kg⁻¹; CIN = −87.7 J kg⁻¹; 0–8-km bulk shear of 58.3 kt, ~30.0 m s^–1) and (d) BRAMS/ERA5 (CAPE = 160.0 J kg⁻¹; CIN = −174.1 J kg⁻¹; 0–8-km bulk shear of 56.1 kt, ~28.9 m s^–1) simulations at 00:00 hours UTC on 15 July 2016. Vertical profiles inferred from BRAMS simulations are extracted at the sounding launching time and for the closest grid point to the Florianópolis station.

Because atmospheric instability indices are used primarily for diagnosing thermodynamic instability when they exceed certain critical values, there is the possibility of severe convection, and thus, the outbreak of a severe thunderstorm. The CAPE and the convective inhibition (CIN) are two important parameters commonly used for assessing the stability of the atmosphere in terms of moist convection. As illustrated in Fig. 4, both sounding data and model simulation depict the presence of CAPE (red areas in the diagrams) and CIN (blue areas in the diagrams), with the layer of CIN acting as a barrier to convection even when higher CAPE values are observed. Using the approach of Hobbs and Wallace (1977), the simulated CAPE and CIN are 439.1 and −87.2 J kg⁻¹ for BRAMS/GFS respectively, and correspondingly 279.5 and −87.7 J kg⁻¹ for BRAMS/CFSv2, and 160.0 and −174.1 J kg⁻¹ for BRAMS/ERA5, against 71.4 and −132.6 J kg⁻¹ obtained from sounding data. Although severe thunderstorms are usually associated with high values of CAPE (>2500 J kg⁻¹), low values of CAPE (<800 J kg⁻¹) can indicate that convection is possible. Based on the severe weather indices for CAPE proposed by Miller (1972) and Gordon and Albert (2000), the atmospheric condition around the sounding site falls into the weak thunderstorm group, and thus the atmospheric instability is considered low. As discussed later in Section 3.5, convection activity was dissipating over this region with clouds slowly vanishing in the process. All simulations perform similarly in terms of a 0–8-km bulk shear, i.e. the magnitude of the vector difference between the surface wind and the winds at 8 km, with respectively 54.4, 58.3 and 56.1 kt (~28.0, ~30.0 and 28.9 m s^–1) for BRAMS/GFS, BRAMS/CFSv2 and BRAMS/ERA5 simulations against 49.8 kt (~25. m s^–1) obtained from sounding data.

It is important to mention that, despite point-based vertical profiles from atmospheric models being often used to evaluate the development of storms, they do not usually reflect the track of radiosondes, which are continuously influenced by advection. As a result, even minor spatial or temporal deviations in the model’s representation of hailstorms can lead to the penalisation of its performance. To mitigate this issue, the adoption of a more flexible evaluation approach that takes into consideration both spatial and temporal aspects would be beneficial. For instance, establishing thresholds to delineate an acceptable range of deviations can be particularly crucial for variables like CAPE, where even minor spatial shifts may still indicate valid model performance. Other than that, point-based vertical profiles from atmospheric models can be useful in the absence of radiosondes. According to modelling results, vertical profiles from BRAMS/CFSv2 and BRAMS/ERA5 over locations where hail fall was reported and simulated by these simulations show relatively low values of CAPE of 700–1000 J kg⁻¹, combined with a strong 0–8-km bulk shear of 60–70 kt (~30.9–36.0 m s^–1; Supplementary Fig. S1 and S4). At the same time, vertical profiles from BRAMS/GFS – simulation that did not produce any surface hail concentrations – show values of CAPE even lower with 50–400 J kg⁻¹, although a 0–8-km bulk shear of roughly the same order of 70–75 kt (36.0–38.6 m s^–1; Supplementary Fig. S5). In terms of spatial distribution, regions in the north-western parts of Paraná and western parts of Santa Catarina states reach CAPE values exceeding 1600 J kg⁻¹; however, no surface hail is produced in these areas (Supplementary Fig. S6).

3.4. Cloud optical properties

In order to analyse simulated cloud optical properties, unavoidable constraints of satellite observations such as microphysical assumptions (Yang et al. 2013), spatial and temporal resolution (Eliasson et al. 2013), and other limitations determined by instrument type (Rybka et al. 2021), must be carefully evaluated. Unfortunately, there were no high-resolution satellite data available for the study period, and model–satellite comparison is performed using MODIS cloud water path (CWP, g m⁻²) and cloud water effective radius (CWER, μm) data at 1.0° resolution. The CWP and CWER together help in understanding the mechanisms behind the formation of raindrops or ice particles, such as collision and coalescence (Freud and Rosenfeld 2012). Additionally, they are often used to characterise the scattering and absorption properties of clouds (Slingo 1989). The MODIS cloud optical properties data derive from multispectral reflectances for pixels identified as probably cloudy or cloudy by the cloud mask during the daytime portions of each orbit. The basic physical principle behind the calculation of MODIS cloud optical properties is the bi-spectral solar reflectance method described in Nakajima and King (1990).

Fig. 5 shows the spatial distributions of CWP for liquid and ice phases derived from MODIS products and BRAMS simulations. Owing to the limited availability of MODIS data over the thunderstorm development region, simulated CWP fields for the ice phase are calculated for 14 July 2016, whereas simulated CWP fields for the liquid phase are calculated for 15 July 2016, both based on the MODIS crossing time over the modelling domain c. 16:30 hours UTC. In general, the model simulations reproduce the observed north-west to south-east CWP fields for the liquid phase across the centre of the domain, except for some regions where simulations show significant biases. These regions include the south-central area of the modelling domain, located on the border between Brazil and Uruguay, as well as the northern regions of Paraná State (Fig. 5a, c, e and g). Regarding CWP fields for the ice phase, both observed and simulated fields show intense convection over the south-east part of the modelling domain, extending east beyond LPB up to the South Atlantic (Fig. 5b, d, f and h). In addition, CWP fields derived from BRAMS/ERA5 simulation spread over a bigger area and are slightly displaced to the north compared to the fields derived from BRAMS/CFSv2 and BRAMS/GFS simulations. Still, for the ice phase, simulated maximum CWP fields coincide with the locations of hail fall according to model results. Only simulations BRAMS/CFSv2 and BRAMS/ERA5 reproduce ground-level hail concentrations. Simulated vertical profiles from these two simulations between 16:00 and 17:00 hours UTC on 14 July 2016 over the maximum CWP areas in Fig. 5f and h (not shown) indicate that the updrafts reach the equilibrium level at ~12 km (anvil cloud), with values of up to 810.0 J kg⁻¹ for CAPE and 61.4 kt (~31.6 m s^–1) for 0–8-km bulk shear. Despite this low value of CAPE, hail fall areas derived from BRAMS/CFSv2 and BRAMS/ERA5 simulations agree well with some of the reported hail fall areas such as São Joaquim and Praia Grande in southern Santa Catarina State respectively (see Supplementary Fig. S1). Previous studies suggest that CAPE and hail size do not present a clear correlation (e.g. Johnson and Sugden 2014). Additionally, there are no conceptual models able to explain why some storms produce large hailstones whereas others produce large amounts of small hailstones (Kumjian et al. 2019). Including ice microphysics in atmospheric models in a realistic way represents a major challenge because of the wide variety of ice particle sizes, shapes and types in the atmosphere (Morrison et al. 2020).

Fig. 5.

Spatial distribution of cloud water path for liquid (left panels) and ice (right panels) phases (g m⁻²) as retrieved from (a, b) Aqua MODIS instrument, and as simulated with (c, d) BRAMS/GFS, (e, f) BRAMS/CFSv2, and (g, h) BRAMS/ERA5 on 14 July 2016 for ice phase and on 15 July 2016 for liquid phase. Model–satellite comparison based on the MODIS crossing time over the modelling domain (c. 16:30 hours UTC).

The BRAMS CWER fields are calculated following the generalised power-law expression of Liu et al. (2008):

where LWC represents the cloud liquid water content (g m⁻³), determined as the cloud water path divided by the cloud depth; ρ_w is water density (g cm⁻³); N is the cloud droplet number concentration (cm⁻³); and β is a dimensionless parameter that depends on the spectral shape of the cloud droplet distribution, set based on observation as follows:

where a_β = 0.07 and b_β = 0.14 (Liu et al. 2008; Freitas et al. 2017). Fig. 6 shows the spatial distributions of CWER derived from MODIS products and BRAMS simulations. As for MODIS, BRAMS simulations are able to reproduce qualitatively the salient features of the observed CWER pattern. Both MODIS and BRAMS show a CWER range of 6–12 μm over most cloudy regions, with the exception of some regions in the south-east part of the modelling domain where large amounts of surface accumulated precipitation are observed (Fig. 3). Over these regions (yellowish areas in Fig. 6), effective radii of up to 16 μm were both observed and simulated. Although no observations of CCN activity were available for the study period, the fact that an aerosol-sensitive quantity is predicted adequately, namely effective cloud-droplet size, is consistent with the prescribed CCN concentration set up in the model configurations. Model bias in simulation of cloud optical properties depends mainly on the aerosol loading (number concentration and chemical composition) and microphysical assumptions (Twomey 1977; Albrecht 1989). Although not shown here, our sensitivity test shows that an increase in the CCN concentration can lead to an increase in the concentration of ice hydrometeors, in accordance with previous BRAMS studies (Martins et al. 2009; Camponogara et al. 2018). Owing to the coarse resolution of the MODIS data, performance statistics were not performed in this case.

Fig. 6.

Spatial distribution of cloud water effective radius (μm) as retrieved from (a) Aqua MODIS instrument, and as simulated with (b) BRAMS/GFS, (c) BRAMS/CFSv2, and (d) BRAMS/ERA5 simulations on 15 July 2016. Model–satellite comparison based on the MODIS crossing time over the modelling domain (c. 16:30 hours UTC).

3.5. DWR-derived measurements

Doppler weather radars are remote sensing instruments widely used for monitoring intense precipitation associated with severe weather systems such as thunderstorms and hurricanes. Radar reflectivity (Z, mm⁶ m⁻³) is often used as the basis to quantify precipitation rate (R, mm h⁻¹), with empirical Z–R relationships of the form Z = aR^b being determined for any particular DWR instrument. In addition, parameters a and b are strongly dependent on the drop size distribution, and therefore on the physical processes that affect it, such as coalescence, breakup, or evaporation (Verrier et al. 2013). In this work, the predicted mixing ratios are converted to reflectivities using the approaches of Smith et al. (1975) for rainwater, and of Lin et al. (1983) and Gilmore et al. (2004) for ice, as follows:

where Z(q_r) and Z(q_i) represent the reflectivities for rainwater and ice respectively; ρ is the atmospheric density (kg m⁻³); and q_r and q_i are the mixing ratios of rainwater and ice (g kg⁻¹) respectively. The model’s equivalent radar reflectivity factor (Z_e) is then obtained by summing the contributions from rainwater and ice, following the approach of Gao and Stensrud (2012):

The Z_e, VIL, VII and ET20 derived from the radar and model simulations at 00:00 hours UTC on 15 July 2016 are shown in Fig. 7. Convective cells over continental and maritime regions evolve differently from each other. Convection over the continent is influenced by the combined action of the SALLJ and regional topography (Fig. 1). The low-level convergence of water vapour flux and the consequent upward movement favour extreme convection and the development of a horizontally extensive (eastward-propagating), mature mesoscale convective system, lasting well into the night. Over time, convective cells over the South Atlantic region just off the coast of Brazil became weak and broke down into discrete cloud patches. This rainfall pattern is fairly well reproduced by the BRAMS simulations (not shown), but there are some differences worth mentioning. First, the storm evolution from the BRAMS/GFS and BRAMS/CFSv2 simulations is delayed compared to radar reflectivity, whereas the simulation initiated with ERA5 places the storm much closer to the observation (Fig. 7a). However, in the south-east of Santa Catarina State, where the storm was more intense according to the observation, the BRAMS/ERA5 simulation shows a much weaker reflectivity field, despite being the only simulation with reflectivity greater than 25 dBZ over the coast. Overall, ET20 values from all simulations are smoother and more than 2 km lower compared to that from radar, which shows several peaks greater than 13 km in height (Fig. 7b). This result shows that the model is not able to represent the vertical development of the system as high as that observed by the radar. In addition, the model shows an underestimation of rainwater production over the continent, with peaks ~12 kg m⁻² compared to over 17 kg m⁻² from observations (see Fig. 7c). By contrast, despite being placed at different locations, the ice peaks from model and radar are both comparable, of the order of 35 kg m⁻². This result shows that the simulations do not accurately represent the vertical extent of the storm, although they can produce hail, as observed by the radar and reported by the media. Supplementary Fig. S7–S14 show the spatial distributions of Z_e, ET20, VIL and VII derived from the DWR instrument and model simulations at 20:00, 21:00, 22:00 and 23:00 hours UTC on 14 July, and 01:00, 02:00, 03:00 and 04:00 hours UTC on 15 July 2016.

Fig. 7.

Spatial distribution of (a) reflectivity (dBZ), (b) Echo Top of 20 dBZ (km), (c) Vertically Integrated Liquid (VIL, kg m⁻²), and (d) Vertically Integrated Ice (VII, kg m⁻²) derived from the DWR instrument (first column) and from BRAMS/GFS, BRAMS/CFSv2 and BRAMS/ERA5 simulations (second, third and fourth columns respectively), at 00:00 hours UTC on 15 July 2016.

To better characterise the evolution of the quality of the simulations, the statistical measures FSS and RMSE were calculated for Z_e, VIL and VII for the period from 20:00 hours UTC on 14 July to 04:00 hours UTC on 15 July 2016 (Fig. 8). In the early stage of the storm, the BRAMS/CFSv2 simulation performs better, as indicated by higher FSS values for Z_e, VIL and VII (Fig. 8a, c, d respectively). Later, at 23:00 hours UTC on 14 July, the BRAMS/ERA5 simulation better locates the storm and produces more accurate predictions of rain and hail (Fig. 8c, d respectively). The BRAMS/GFS simulation has the worst performance among all simulations. The same behaviour observed in the FSS for Z_e is also found in the RMSE, with a greater error for the BRAMS/GFS simulation (Fig. 8b). The RMSE values for VIL and VII (Supplementary Fig. S15 and S16) are similar for the three simulations. Possible sources of model–radar discrepancies include errors in initial and boundary conditions, mismatches between model and radar grid points, and the inherent limitations and assumptions of the radar, as discussed in Section 2.4.

Fig. 8.

Domain-wide fractional skill score for reflectivity (a), Vertically Integrated Liquid (VIL) (c), and Vertically Integrated Ice (VII) (d), and domain-wide root-mean-square error for reflectivity (b), for the period from 20:00 hours UTC on 14 July 2016 to 04:00 hours UTC on 15 July 2016.