2.1 WRF-Chem model. Simulation domains

The WRF-Chem model was developed by the WRF Atmospheric Chemistry Working Group and simulates turbulent mixing, transport, transformation, and destiny of gases and aerosols simultaneously (^{Grell et al., 2005}).

The chemical module of WRF was used in this study with the volcanic emissions option enabled, assuming the ground-based generators as “ash-emitting small volcanoes”, and that AgI aerosol dispersion properties were similar to those of ash.

To avoid introducing large interpolation errors with a drastic decrease in resolution, a total of four domains were defined: A coarse external domain of 27 km and three nested domains of 9, 3 and 1 km, each of them identified as d01, d02, d03, and d04, respectively (Fig. 1).

Fig. 1 Simulation domains. 

The finest grid (1 km) covers the region of interest (eastern Cuba) centered on 20.334 ºN and -75.824 ºW. It has a horizontal extension of 442 × 220 grid points and 50 vertical levels distributed on an eta (ɳ) coordinate axis that decreases with an inverse dependence of the difference between the pressure of the level and the pressure at the top of the troposphere, this is between the level ɳ = 1 at the surface until ɳ = 0 at the top of the model. The pressure at the top of the troposphere in the model is 50 hPa.

Topographic data were downloaded from the website https://www.usgs.gov, with a resolution of 30 seconds of arc (approximately 900 m), which is adequate for the accuracy required in this case study, since the higher resolution domain is 1 km.

2.2 Emission parameters of the point sources

The emission parameters of the point sources are shown in Table I. The basic function of the AgI generator is the combustion of this chemical compound. The AgI solution burned by the generator is composed mostly of acetone, which acts as the primary solvent for the chemical and as fuel.

The size distribution of emitted particles (where the size is assumed as particle diameter) has a large influence on aerosol dispersion, with smaller sizes resulting in greater spatial dispersion. The distribution proposed by ^{Demott (1982)} was chosen since both the solution and the generator were similar to those proposed for installation in Cuba. A 2% AgI concentration in the solution is common for both samples. So far, no measurements have been made in Cuba for the particle size distribution of the combustion result of the AgI solution in acetone. Since measurements of different 2% AgI made by ^{Demott (1982)} solutions were similar, the distribution was selected arbitrarily, since a simulation study was made for them and results were similar. Measurements made by ^{Demott (1982)} of a 2% AgI-NH₄-acetone-water with 20 mole % NH₄ClO₄ solution (see his figure 8b) were chosen as a reasonable estimate of the proposed generators particle size distribution. From this distribution it was possible to calculate the total mass of AgI, which was split in 10 intervals of different sizes corresponding to ranges of aerosol particle diameters from 0.0075 to 0.244 μm (Fig. 2).

Fig. 2 Aerosol particle size distribution assumed from the combustion of 2% AgI-NH₄-acetone-water with 20 mole % NH₄ClO₄ distributed in 10 intervals (bins). 

2.3 Case studies

The selected days for simulations correspond to the rainy period (May to October) during 2015 and 2016. The selection was made considering the presence of convective clouds in eastern Cuba related to synoptic scale or mesoscale phenomena. This ensured the presence of atmospheric conditions suitable for convection.

The influence of the Azores High pressure system is common during the rainy season in Cuba, with the rain primarily produced by local conditions, mesoscale flux convergence and diurnal instability due to surface heating. These conditions are typical in all the cases selected for this study, as shown in Figure 3. In some cases, convection was stimulated by the presence of a cut-off low pressure system or associated trough at 200 hPa.

Fig. 3 Synoptic conditions of the selected days for simulation. Geopotential contours at the 200 hPa level are presented in panels a, c, e, g, i, k, while the corresponding patterns at 850 hPa are shown in panels b, d, f, h, j, l. 

The atmospheric variables for the selected days were taken from the output of the Global Forecast System (GFS) model, available at http://nomads.ncep.noaa.gov. These data are in 3-h intervals with a spatial resolution of 50 km. Simulations were initialized at 00:00 UTC of each selected day until 03:00 UTC of the following day, for a total of 27 h for each study day. Model outputs were recorded every 30 min.

The ground-based generators were turned on in the simulation at 13:00 UTC and remained working for 11 h. Two hypothetical spatial distributions of 13 generators were analyzed. Point sources were located on the hillside of the mountains, taking advantage of the orographic forcing, and favoring the introduction of the reagent in the cloud.

The first distribution (I) corresponds to a location of the generators in nearby positions (10 to 15 km apart, according to the criteria used in the installation plans for the generators in Cuba) at an average height of 445 masl, all located in the province of Santiago de Cuba, in the Sierra Maestra mountain range and the Cordillera de la Gran Piedra (Fig. 4a). In the second experiment (II), the generators were hypothetically distributed in the main elevations of the eastern region and the heights above sea level were higher than in the first case, with an average height of 585 masl. In this case, the generators were scattered in the Sierra Maestra, Sierra del Cristal and in the Cordillera de la Gran Piedra (Fig. 4b). In no case the heights exceed 940 masl.

Fig. 4 Spatial distribution of the ground-based generators I (a) and II (b). 

To establish if an increase in the number of generators would contribute positively to the seeding process, two case studies were simulated considering the same emission parameters but different emission rates. The first case assumed an emission rate of 10 g h^-1 and the second of 30 g h^-1, which is equivalent to assuming the presence of three generators emitting at a rate of 10 g h^-1, located close enough to the original location so that they contribute to the same node of the model grid (1 km).

2.4 Seeding threshold

According to the technical characteristics of the pyrotechnic flares PV-26-01 (^{Martínez et al., 2011a}) used in the second program of artificial increase of rain in Cuba in 2005, at -5 ºC, 1 g of AgI contributes to 10¹³ ice nuclei. Hence, for the static seeding to occur about 0.1 ng m^-3, it would be necessary to obtain 1 nucleus per liter. On the other hand, for dynamic seeding, approximately 10 ng m^-3 of AgI are necessary to obtain concentrations of the order of at least 100 nuclei per liter. Since there are no published measurements of the mixture of AgI burned in the generators, the above values are used as concentration thresholds necessary for the occurrence of glaciogenic seeding.

3.1 Analysis of the maximum concentration and area occupied by the aerosol

The maximum concentration of AgI between the 1000 and 1500 m levels above the surface was analyzed for each selected case study. The maximum daily concentrations for each of the six days simulated under experiments I and II are shown in Figure 5. For emissions of 10 g h^-1, most maximum concentrations do not exceed 1 ng m^-3, but most of the reports are above 0.1 ng m^-3 (Fig. 5a, b). The maximum average concentration over the six selected days is reached at 17 h, 4 h after the generators are turned on. Approximate values of 0.53 and 0.66 ng m^-3 were obtained for spatial distributions I and II, respectively. When the emission rates were increased to 30 g h^-1, the maximum concentrations grew proportionally to the rate of the emission increase (Fig. 5c, d).

Fig. 5 Maximum concentration of AgI between 1000 and 1500 m above the surface for an emission rate of 10 g h^-1 for the spatial distributions I (a) and II (b); and for emissions of 30 g h^-1 for the spatial distributions I (c) and II (d) for the set of days simulated. 

The area occupied by the aerosol in the studied layer was also obtained for concentrations greater than 0.1 ng m^-3 for each case study. For emissions of 10 g h^-1, the reported area for the spatial distribution I is greater than for II, with an average maximum surface of 67 and 39 km², respectively, both at 16 h (Fig. 6a, b). In the case of 30 g h^-1 emissions, the areas increase considerably, reaching more than 400 and 600 km² for distributions II and I, respectively (Fig. 6c, d). After 16:00 UTC, the AgI covered areas, with concentrations above 0.1 ng m^-3, as well as the maximum concentration value, experienced a decrease until the generators were turned off.

Fig. 6 Average area with the presence of AgI in concentrations ≥ 0.1 ng m^-3, between 1000 and 1500 m above the surface for an emission rate of 10 g h^-1 for the spatial distributions I (a) and II (b); and for emissions of 30 g h^-1 for the spatial distributions I (c) and II (d) for the set of days simulated. 

The output of the model was analyzed for greater heights between the isotherms of -5 and -10 ºC, where the aerosol particles may potentially be activated as ice nuclei. Since for 10 g h^-1 emissions, AgI aerosol concentrations were not higher than 0.1 ng m^-3 at those altitudes, the analysis was carried out only for emissions of 30 g h^-1. Maximum concentrations slightly exceed 0.1 ng m^-3 in distribution II, whereas for distribution I they are higher; also, there are more hourly reports per day above this threshold value (Fig. 7). The area occupied by the aerosol dispersed in concentrations greater than 0.1 ng m^-3 is very small, with values that do not exceed 15 km² for distribution II and 30 km² for distribution I (not shown). Therefore, for the studied days, static seeding could occur with an emission rate of 30 g h^-1. However, the areal fraction occupied by the AgI concentrations above 0.1 ng m^-3 is very small. On the other hand, dynamic seeding would not occur in any case, because concentrations above 10 ng m^-3 were not obtained for any of the studied days between the seeding levels.

Fig. 7 Maximum concentration of AgI between -5 and -10 ºC for an emission rate of 30 g h^-1, for the spatial distributions I (a) and II (b) for the set of days simulated. 

To study the cloud seeding effectiveness, the convolution of the cloud water content and aerosol concentration ≥ 0.1 ng m^-3 was calculated. The presence of cloud water above the 0 ºC isotherm is considered as supercooled water. Figure 8 shows the occupied area by the coincidence of these two variables for an emission rate of 30 g h^-1. The existence of coincidences indicates that the aerosol could serve as ice nucleus, allowing the static seeding to occur mainly for distribution I (Fig. 8a), where the area occupied is higher.

Fig. 8 Average area with the simultaneous presence of cloud water and aerosol concentration ≥ 0.1 ng m^-3, between -5 and -10 ºC for an emission rate of 30 g h^-1, for the spatial distributions I (a) and II (b) for the set of days simulated. 

3.2 Influence of the horizontal wind

Figure 9 shows the influence of wind on aerosol dispersion for levels between 1000 and 1500 m above the surface. On July 15, 2016 winds were predominantly easterly due to the influence of the North Atlantic High. Under these conditions, for spatial distribution I aerosol is advected by the prevailing wind towards the sea (Fig. 9a, c, e, g). Similarly, for spatial distribution II, aerosol is driven away from land for generators located in the Cordillera de la Gran Piedra and some of the generators in the Sierra Maestra (Fig. 9b, d, f, h). However, for this day the wind keeps within the target region the emission plumes of generators located in the Sierra del Cristal and in the westernmost part of the Sierra Maestra. In all cases, aerosol is advected towards the sea after reaching its maximum concentration over land.

Fig. 9 Average hourly concentrations of AgI ≥ 0.1 ng m^-3, between 1000 and 1500 m above the surface, on July 15, 2016 for distributions I (a, c, e, g) and II (b , d, f, h), for emissions of 10 g h^-1 at various times. Wind barbs indicate the direction and speed of the average wind at the levels between these heights. 

3.3 Vertical cross section

Figure 10 shows a longitudinal vertical cross section at latitude 20.11 °N, corresponding to the location of one of the generators at 422.9 masl and 75.95 ºW. The plume of another generator located at a different latitude can be also detected in Figure 10b, c, d. Note that the plume does not reach seeding levels in this case, as it remains under the 0 °C isotherm. The main characteristics for August 16, 2016 was the presence of a cut-off law pressure system over eastern Bahamas at high levels, inducing unstable conditions over eastern Cuba and favoring a plume ascent. Under spatial distribution I the plume reached higher altitudes than for distribution II (not shown), indicating that the local prevailing conditions on a given day vary, mainly due to local orography, which contributes to the ascent of the plume by mechanical forcing. The higher rise of the plume occurs in the first 4 h of ignition of the generators for all simulated cases.

Fig. 10 Longitudinal vertical cross section at latitude 20.11 °N, corresponding to the location of one of the generators for distribution I, on August 16, 2016, showing the concentration of AgI in ng m^-3. The continuous blue isolines represent the turbulence kinetic energy (TKE) (J kg^-1). The isotherms of 15 and 0 ºC are represented by dashed lines. Wind barbs indicate the direction and speed of the wind. The black triangle indicates the generator position. 

Aerosol dispersion is closely related to the TKE. The higher the values of TKE, the higher the probability of vertical rise of the plume. TKE is the result of the sum of several terms that, in addition to turbulent transport and wind shear terms, also consists of the buoyant production term, the TKE advection term, the dissipation term and the pressure adjustment of correlation term. TKE is an indicator of the height of the atmospheric boundary layer, therefore it can be inferred that, for some cases, the plume exceeds the height of this layer. Therefore, it is likely for the aerosol to penetrate the base of the cloud under cloudy conditions.

Figure 11 shows a comparison of the mean plume height for concentrations greater than 0.1 ng m^-3 and the mean TKE height for values greater than its mean below the PBL height (0.46 J kg^-1) between the hours where the convection is active. This analysis was made for the cross section shown in Figure 10.

Fig. 11 Mean vertical distribution of the aerosol plume height and turbulence kinetic energy (TKE) height for the cross section shown in Figure 10. 

Vertical velocity was taken as a criterion for convection initiation if values remained greater than zero in the morning and early afternoon (Fig. 12). The convection time prevailed from 15:30 to 18:30 UTC, where the TKE and the plume height show their maximum values (Fig. 11). The dependence of the aerosol plume height to the TKE height is clear in Figure 10. Maximum values for both variables were reached at 18:00 UTC.

Fig. 12 Mean vertical velocity for the cross section shown in Fig. 10. 

To visualize the AgI plume concentration values that could reach the height for effective seeding, Figure 13 illustrates a longitudinal vertical cross section at 20.08º N with the representation of the isotherms of 0, -5 and -10 ºC. The ground-based generator is located at this latitude and 75.78º W, at 383.2 masl. Note that concentrations of the order of 0.02 ng m^-3 reach and surpass the height of the -5 ºC isotherm and coincide with the presence of supercooled water (Fig. 13c, d). However, these concentrations are insufficient for static seeding to occur. Thus, a simulated emission ratio of 10 g h^-1 does not ensure that the reagent with sufficient concentration could reach the height required for seeding. Consequently, the use of ground generators with this emission rate would not accomplish the desired objective and would result in the total loss of seeding material for the selected case studies.

Fig. 13 As Figure 10, but for latitude 20.08 ºN, with the 0, -5 and -10 ºC isotherms represented by dashed lines. The continuous red isolines represent the cloud water mixing ratio (g kg^-1).