Rainfall events (RE) are among the phenomena that cause significant negative impacts on the economy and society. This kind of events is quite frequent in Panama, where national and regional important activities such as the operations at both the Panama Canal and the international flight hub, have often been affected. Therefore, the implementation of a numerical weather forecast system becomes a high priority contribution to the improvement of early warning mechanisms and the development of new climate services oriented to energy, agriculture and transport sectors.
Panama, a relatively small country located in the extreme eastern part of Central America is frequently affected by RE of different levels of magnitude. In a general sense, high amounts of precipitation occur throughout the whole year but they are more frequent in summer. Large-scale weather systems such as disturbances associated with the Intertropical Convergence Zone (ITCZ) and cold fronts penetrating from northern high latitude appear as two of the major circulation systems that trigger convective activity  . The narrow and elongated position of the country between the Pacific Ocean and the Caribbean Sea, as well as the moisture that is transported from both bodies of water, is also a key factor . Finally, the influence of orography over the low-level convergence tends to enhance the convection activity and favors the occurrence of heavy rainfall episodes. As has been pointed by   a significant convective activity takes place over the Panama Bight with strong afternoon rainfall over the narrow isthmus of Panama.
Precipitation episodes and particularly RE are very difficult to predict over short timescales in a complex region like Panama, where several different mechanisms interact. Despite Global Numerical Weather Prediction Models have continuously improved during the last years, their low resolution prevents them from producing a reasonably good forecast in such complex regions where convection is a critical factor producing RE. The Regional Models, which represent explicitly many of these complexities like the convection and the interactions with large-scale weather phenomena appear as a good option for implementing a forecast precipitation system   that can produce useful and relevant information for a wider user community and to create the base for future development of end-user services   .
This work presents an assessment of the ability of the advanced research core of the Weather Research and Forecast model (hereafter referred to as WRF) to forecast precipitation over Panama by exploring different model configurations based on various combinations of cumulus and microphysics parameterizations. The model outputs are thoroughly evaluated using observations from different available sources to find physically-based answers to the produced forecasts.
WRF is a very popular regional model widely used for research    and real-time operational forecasting     of a variety of meteorological events, such as extreme rainfall. Despite the large amount of published work that can be found on the use of this model around the world, there are just a few studies in Central America   and in particular over Panama .
The paper has been organized as follows. Section 2 offers the details of the WRF version employed in this study as well as the different model configurations considered. It introduces the different sources of rainfall records and other data used for synoptic and processes analysis. Methods applied for verification and sensitivity assessments are also explained. Section 3 is dedicated to present and to discuss the results obtained.
2. Model Configurations, Data and Methods
WRF-ARW V3.8.1 was used to carry out the simulation of the precipitation events. The initialization and boundary conditions every 3 hours, were provided by NOAA’s Global Forecast System (GFS) with 0.5˚ resolution. Figure 1 shows the three simulation domains defined, with resolutions of 27, 9, and 3 km. The verification was focused in the innermost domain.
The main details of the physical configuration used for the model are presented in Table 1. Static data selected were provided by the USGS, with 33 categories of land use and a topography resolution of 30 s. The vertical grid was composed of 29 sigma levels from the surface to ~10 km using a smaller spacing in the low levels. This configuration was the base to conduct 15 different experiments resulting from the combination of five bulk microphysics parameterizations (BMP) and three cumulus parameterizations (CU). The experiments were developed following the results from    during the implementation of the operational forecast system established in Cuba. Simulations were carried out for a period of 48 hours starting on January 16 2018 0000 UTC.
Figure 1. Simulation domains over Panama. The red, blue and green squares correspond with resolutions of 27, 9 and 3 km respectively. The area covered by the green square signals the region identified as the core verification area (CVA).
Table 1. Physical configuration details for the WRF-ARW model.
The BMP’s include three single-moment schemes (Lin , WSM6 , and Thompson  ) and two double-moment schemes (WDM6  and Morrison ). In the case of CU, the schemes selected were the convective adjustment scheme of Betts-Miller-Janjic (BMJ)  and the two mass flux schemes Grell-Freitas (Gr)  and Kain-Frisch (KF) .
Single-moment BMP’s prognoses the third-order moment of the size distributions (mass mixing-ratio variables) of different water categories. Lin scheme simulates five classes of hydrometeors (cloud water, cloud ice, rain, snow, and hail). The exponential size distribution and terminal speed formulation are used to represent only the precipitation particles: rain, snow, and hail, including ice sedimentation. WSM6 is based on WSM5 , adding processes related to graupel. The terminal velocity for graupel is calculated using a mass-weighted expression following  . Thompson has implemented a generalized gamma distribution shape for each hydrometeor, a variable gamma distribution shape parameter for cloud water droplets based on observations and the calculation of the intercept variable of rain and graupel depending on the respective mixing ratios. Double-moment BMPs were developed to improve the explicit representation of clouds and precipitation in mesoscale atmospheric models. The scheme predicts the evolution of mass, as well as number densities of the five hydrometeor types, cloud droplets, raindrops, cloud ice, snow and graupel  . Since the number concentrations of all these hydrometeors are calculated explicitly, the relevant homogeneous and heterogeneous nucleation processes have been parameterized including the activation of cloud condensation nuclei, which is not predicted in most state-of-the-art cloud-resolving models.
The BMJ scheme of CU introduces the parameter “cloud efficiency” to calculate the relaxation time and the deep convection using specific reference states in the Betts-Miller parameterization  . The Kain-Fritsch parameterization scheme uses a one-dimensional entraining/detraining plume model assuming that any mixture that becomes negatively/positively buoyant detrains/entrains from/into the cloud. On the other hand, the Grell scheme implemented in the WRF model follows the method developed by , to represent an automatic smooth transition as the resolution increases.
2.1. Available Data
Table 2 summarizes the relevant characteristics of the different sources of available data that were used in this work. ERA-Interim reanalysis  was employed for the description of the synoptic conditions during the selected day. Precipitation data within the CVA comes from a group of 38 weather stations (WS) operated by the Meteorological Service of Panama and the Panama Canal Authority (PCA). Also precipitation data from the Global Precipitation Mission (GPM)  were used, considering their high spatial and temporal resolution as well as their availability for the time period of this study, this allowed these data to be very conveniently compared to the model’s precipitation output fields. A set of other
Table 2. Meteorological observations and data used for verification.
meteorological data such as upper air soundings provided by the PCA was also used to assess the WRF skill in relation to precipitation characteristics and other processes associated with its formation. Figure 2 shows the location of the WS.
The methodological approach used involved the following steps. First, a comparison of 24 hours accumulated precipitation from WS and their corresponding underlying grid-point cells of GPM was made. The aim of the comparison was to assess the accuracy of GPM to represent the rainfall behavior over Panama. WRF experiments were then bilinearly regridded to the GPM grid and statistically compared, adopting different widely used metrics. These metrics included Pearson’s correlation coefficient (r) and a set of categorical indexes  , that were calculated to explore the rainfall detection performance of WRF experiments. The indexes were: the probability of detection (POD), false alarm ratio (FAR), and critical success index (CSI). POD is used to describe the extent of rainfall events that are correctly forecasted; FAR reflects the extent of rainfall events forecasted that are false alarms, and CSI illustrates the overall proportion of rainfall events that are correctly detected in the forecast. Table 3 presents the mathematical expressions for computing the categorical indexes.
3. Results Discussion
3.1. Synoptic Situation Analysis
The surface chart (Figure 3) from Jan 17 2018 0000 UTC, shows a stationary
Figure 2. Geographical location of the different weather stations (black points) and upper air sounding (red star) used as sources in the verification.
Figure 3. Synoptic surface analysis and contours of mean sea level pressure from ERA-Interim reanalysis.
Table 3. Categorical indexes used in the assessment of the experiments . Capital letters H, M, FA corresponds with hits, misses and false alarms respectively.
front over the Caribbean Sea in dissipation stage, spanning from the eastern part of Cuba to the northern border area between Costa Rica and Panama. A trough axis is projected on the Pacific coast of Panama, which is coupled with the ITCZ. An increase and development of cloudiness due to daytime heating is observed through the day, more relevant in the center and east of the country. This was influenced by a stationary frontal system and two low pressure systems over the southwest of Panama and the north of Colombia.
The analysis of the vertical profile (Figure 4), given by an upper air sounding, shows a high water content in the low levels mainly from the surface up to 3 kilometers high and a wind of 10 to 15 knots on average from the northeast. In
Figure 4. Upper air sounding launched at 8.98 N, 79.57 W Jan 17 2018 0000 UTC.
medium levels a layer with a relative humidity of 50% and 10 to 15 knots east winds is observed, while high levels show southwestern and western winds with a relative humidity of 50% to 60%. The wind vertical profile shows little shear and the trigger temperature is 29.6 degrees.
3.2. Correspondence between GPM Estimated Precipitation and Surface Stations Reports
As it was mentioned in the previous section, in order to assess the performance of the experiments with a better spatial and temporal resolution, the precipitation estimated by the GPM product was used. However, as it is not a direct measurement, a characterization of the product is needed. Figure 5 presents the 24 hour accumulated rainfall registered by the WS (Figure 5(a)) and the same variable estimated by GPM (Figure 5(b)). As can be seen, despite showing a low density of stations, accumulated rainfall appears well represented for this particular case study. The spatial distribution of precipitation events that occurred in the area of Panama City and the Panama Canal, shows that there were two nuclei of these events, one over the Caribbean Sea area where the highest rainfall
Figure 5. Daily rainfall (mm/24h) for Jan 16 2018: (a) Registered at weather stations; (b) Estimated by the GPM product.
was recorded of about 60 mm/24h and the other area on the Pacific slope of the canal with accumulates very close to 30 mm/24h; both measured at the stations. On GPM these two maxima are observed, showing agreement with the observations. Other places where this correspondence is observed are the western part near the coast and the border with Costa Rica, and south of the Central-American Sierra Madre; where small nuclei of precipitation are observed, also reflected in GPM. An important feature is that GPM data did not reflect precipitation at places where the stations did not report it, which is a good indication of agreement. Despite the good spatial representation of the precipitation areas, an overall underestimation of the amount of precipitation is observed on GPM product. Precipitation values are approximately between 20 and 30 mm below those registered at the weather stations.
Although there are quantitative differences in the precipitation obtained from GPM for the case studied in this investigation, this data was useful in the spatial sense, for the verification of the experiments developed. Also it was important for covering the regions where no weather stations are present, for example near the border with Colombia. Hence, these data were used to illustrate rainfall at sub-daily values in the absence of both stations and hourly stations data.
3.3. Precipitation Verification
Diurnal cycle graphics of precipitation for the 15 combinations of BMP and CU on a nine boxes mosaic that divides the WRF inner domain over Panama is presented in Figure 6. Each graph is located over its corresponding geographic region and a shaded map of the 24 hours GPM accumulated precipitation is shown in the background. Wherever significant total rainfall occurred, most of the combinations properly represented the corresponding peaks, although in general with slight displacements in time (boxes 1, 2, 5, and 7). For these boxes there stand out, with better results the combinations Mor-BMJ, Thom-Gr, WSM6-BMJ and WSM6-KF.
The rainfall estimated by GPM for regions 6 and 9, that cover the Caribbean Sea was not significant and it was over these areas where more than half of the combinations approached best the accumulates from GPM. A similar behavior occurred at box 1 which is also located over the Pacific Ocean. This indicates that most of the combinations presented greater ability to forecast precipitation
Figure 6. Hourly diurnal cycle of precipitation on Jan 16 2018, starting at midnight. The CVA was divided into nine boxes where the total area precipitation every one hour was computed.
over sea for the case studied. This is an expected performance, considering that over sea, the surface friction is lower leading to a more uniform wind field and less pronounced temperature changes, allowing a better representation of convective processes. Lin-KF, Mor-KF and Thom-KF presented the mayor overestimation in these areas.
At box 4 however, despite having a large part of the area occupied by sea, the model had difficulties representing the rain in the earliest and latest hours of the day for all configurations, with a huge overestimation. Something similar occurred with box 3 but less pronounced and with greater dispersion in forecasts. Over zone 8, which presents irregular coastal areas, the model presented difficulties, overestimating precipitation. It also occurred at box 7 where overestimation occurred in the afternoon. Once more combinations with the KF cumulus parameterization produced the higher values of accumulated rainfall been the ones with mayor differences from GPM.
Regions 2 and 5, also enclosed water-land interface zones and the mountainous system (Central-American Sierra Madre). In spite of it, most of the configurations were able to forecast the precipitation with low errors, representing well, though earlier, the peaks of precipitation. Mor-BMJ achieved better representation of the diurnal cycle over the box number 5, where weather stations recorded the most significant accumulated precipitation.
In general, the ability of each combination of BMP-CU varies significantly from one region to another, so it can’t be easily identified, which ones achieved better accuracy. As a common feature, all experiments overestimate the hourly total area of precipitation in each box, presenting higher precision at regions where sea is predominant. It was found that experiments that use KF as CU strongly overestimate the rainfall values while the ones with BMJ fit better the daily behavior of the precipitation.
From the point of view of categorical verification indexes, in which the skill of the experiments is evaluated for the prognosis of the occurrence or not of precipitation, Figure 7 shows the probability of detection for each BMP-CU combination. The POD is a measure of whether precipitation events were correctly identified in terms of occurrence, not amount. Notice that the 15 numerical experiments presented a similar spatial pattern of the index, having higher values in the areas where GPM estimated the most significant accumulated rainfall. This is over the Pacific Ocean region near the Panama-Colombia border and over the southwestern part of the CVA. In these regions POD reached values of 0.8 while at zones where the greatest episodes of precipitation were recorded (values greater than 50 mm/24h), the ability of detection drops to 0.5. On the other hand, the performance of event detection decreases, having values below 0.4 at locations where weather stations registered values between 30 and 60 mm/24h, such as south-west and north-west of the Central Panama province.
The best skill is presented by the group of combinations that use BMJ as CU parameterization. For the area comprising the Panama Canal, BMJ presented an improvement that could be quantifiable up to 10% in the POD index when compared to the other two cumulus parameterizations. The Lin-BMJ experiment stood out as it presented values close to 1 in the northern and southern parts of Central Panama, decreasing to 0.5 in zones where the greatest episodes of precipitation were recorded by GPM.
The CSI index, which gives the proportion of events being predicted correctly is shown in Figure 8. In general, the values are very low in comparison with POD, this is due to the consideration of false alarms in the calculus of the index. The above suggests that all the experiments give a considerable amount of areas that are indeed FA. The FAR index (not shown) rises to values very close to 1 in most of the core verification area which explains the low values of CSI. Particular interest deserves the regions covered by the Caribbean Sea, where CSI barely achieves values of 0.1. Over this region no precipitation event occurs according to Figure 5(a). The poor behavior of the combinations of BMP-CU indicates that, although the ability of the model is good when rainfall is registered, it fails when precipitation does not occur. These results are highly related with the
Figure 7. POD index for all experiment conducted. Rows are ordered by cumulus parameterizations and columns are organized by microphysics parameterizations.
Figure 8. CSI index for all experiment conducted. Rows are ordered by cumulus parameterizations and columns are organized by microphysics parameterizations.
overestimation encountered over the CVA in the diurnal cycle analysis. For this particular case, just like with POD, the best combination of parameterizations was BMJ with Lin, while the experiments that have KF as CU showed the worst skill.
Pearson correlation was computed for both GPM (Figure 9(a)) and WS (Figure 9(b)) 24 hours accumulated precipitation. This correlation is spatial in the sense that was calculated for each combination, for the entire CVA and all weather stations. All experiments here showed a poor skill, especially when compared with GPM. Values lower than 0.4 indicate that the combinations of BMP-CU failed, possibly because areas of precipitation are shifted in space, misshapen or not forecasted. WDM6-BMJ, stands out slightly in the verification with GPM with 0.3 as correlation and WSM6-BMJ highlights as the best with values around 0.8 in the comparison with WS. Focusing on the experiments assessments with weather stations, another two combinations that used BMJ as CU, reached high correlations values. These are Mor-BMJ and Thom-BMJ, both over 0.5. Again BMJ stands out as the better CU selection in three of the combinations studied. On the contrary, those that use Gr cumulus parameterization showed no correlations or negative correlations with real values.
As a verification of the displacement and/or deformation of the areas of precipitation, Figure 10(a) presents the one hour accumulated rainfall from GPM for Jan 17 2018 0000 UTC. Figure 10(b) and Figure 10(c) show the accumulated precipitation for the same period from Mor-BMJ and Mor-KF respectively. Figure 10(d) and Figure 10(e) show the categorical maps resulting from the comparisons with GPM. It can be observed that the forecasted areas for both combinations are quite smaller than those shown on the GPM chart. It can be seen as well that there are zones of precipitation registered by GPM that the experiments did forecast, but spatially shifted (see for instance the precipitation area over Veraguas). Other precipitation events shown in GPM were missing in the experiments, for example the one southwest of Panama City. All of this could explain the presence of large areas of missing and false alarm categories, and therefore the low areal correlations between GPM and BMP-CU combinations. KF overestimation of precipitation is also made evident in the numerous areas of false alarm that appear in Figure 10(e).
Figure 9. Pearson correlation: (a) Computed with GPM; (b) Computed with weather surface stations.
(a) (b) (c) (d) (e)
Figure 10. Spatial categories verification for Jan 17 2018 0000 UTC; H, CN, M, FA mean hits, correct negatives, misses and false alarms respectively: (a) GPM estimated precipitation; (b) Mor-BMJ accumulated rainfall; (c) Mor-KF accumulated rainfall; (d) Mor-BMJ categories; (e) Mor-KF categories.
3.4. Verification with an Upper Air Sounding
In order to achieve a better understanding of the results, an evaluation with the sounding at Jan 17 2018 0000 UTC, was carried out. The numerical vertical profiles for wind velocity, relative humidity and temperature were extracted from each experiment following the real sounding trajectory. Figure 11 shows these vertical profiles. The discontinuities in the forecasted profiles are due to the fact that the sounding device track left the model domain area at about 400 hPa and reentered at about 200 hPa.
Regarding the wind profiles (Figure 11(a)), whereas for BMJ, the prognosticated profile underestimated the measured wind at lower levels, the rest of the
Figure 11. Profiles from all experiments and real sounding on Jan 17 2018 0000 UTC. (a) Wind speed (m/s); (b) Relative humidity.
combinations and WSM6-BMJ exceeded the intensity values. The rapid but light changes in wind velocity in medium levels were better represented by BMJ, except WSM6-BMJ. In general, most of the combinations successfully followed the sharp peak above 200 hPa level reported by the sounding, as well as the significant wind intensity decrease experienced at higher levels.
All experiments were able to represent the vertical humidity profile close to the sounding measurement (Figure 11(b)). However, between 700 and 600 hPa, they showed difficulties in the estimation of the variable, this could be related to the amount of vertical levels defined in the model grid. It is appreciable that BMJ and KF schemes combined with all BMP schemes are the ones that represented best the vertical profile of relative humidity. Even though BMJ does not account for any changes that occur below the cloud base , it showed the most appropriate behavior of humidity at lower levels (except when combined with WSM6). In the BMJ scheme, moisture profiles are a determinant factor in the identification of the different states of convective equilibrium. The scheme also succeeded to represent well the profiles of humidity and temperature (not shown), which are used to calculate the entropy and “cloud efficiency”, that measures the ability of the convective column to transport the enthalpy in the vertical producing the least possible precipitation. This could explain the ability presented by BMJ in well representing the total amount of precipitation as well as the behavior of the diurnal cycle, for this case studied. The rain in the sounding area was mainly caused by diurnal heating. The temperature and humidity profiles reflected the adjustments made by the scheme that forced them to climatological reference profiles. Those schemes managed to adequately represent the vertical profile of both variables. The BMJ parameterization, which is supported by the concept that convection is a purely thermodynamic process, was therefore the most capable of reproducing precipitation in the area.
On the other hand, Gr presented less skill predicting the variables, overestimating the precipitation in most combinations. A study conducted in the tropical zone describing the African monsoon presented similar difficulties with this parameterization . All KF combinations also overestimated the amount of rainfall, which may be related to the fact that this scheme tends to leave an extremely deep saturated layer after convection  in order to better represent the sounding profiles. This is a deficiency, since as these conditions remain in the environment, they activate the microphysics scheme that tends to produce stratiform precipitation.
Figure 12 represents the Taylor’s diagram for the wind velocity and relative humidity. Note that once again the experiments that use BMJ as cumulus parameterization had the higher correlations while the ones configured with Gr and KF had the worst performance as shown in the wind profile diagram. The contrary occurred when comparing moisture profiles, where combinations that use KF correlated slightly higher.
The main purpose of this work was to assess the ability of the WRF model to forecast rainfall events in Panama. In this regard, the performance of WRF was tested through verifying the precipitation simulated by the model as well as analyzing the sensitivity of the forecasts to a group of different parameterization schemes. The investigation includes the evaluation of model capacity to reproduce features of processes associated with rainfall occurred in the study region. The model simulations were carried out for a period of 24 hours over three nested domains with a resolution of 27, 9 and 3 km, respectively. Experiments were developed including combinations of Lin, WSM6, Thompson, WDM6, and Morrison microphysics parameterizations and Betts-Miller-Janjic, Grell-Freitas and Kain-Fritsch cumulus schemes.
The ability of the selected microphysics-cumulus combinations to forecast the rainfall varies significantly across different areas, making difficult the clear identification of those which better simulate the diurnal cycle of rainfall. In areas where high daily precipitation amounts occurred, most of the experiments produced an acceptable representation of the peaks, although with slight displacements in time. The Mor-BMJ emerged as the best combination representing the diurnal cycle in various areas of the CVA, followed by Thom-Gr, WSM6-BMJ, and WSM6-KF combinations. Many of the combinations present a good skill to forecast precipitation over sea areas. The largest overestimation for these zones is produced by combinations including the KF cumulus scheme.
Figure 12. Taylor’s diagram for wind profiles and relative humidity profiles.
The categorical evaluation of the occurrence of precipitation demonstrates the best ability of those schemes that combine BMJ as cumulus parameterization. This happens mainly over the area comprising the Panama Canal where Lin-BMJ presented a clear improvement over all other cumulus schemes. The overall WRF performance of rainfall/no-rainfall events forecast was fairly poor and it explains the overestimation found in the diurnal cycle within the validation area. The worst combinations were those that include the KF parameterization. For all combinations, there were low spatial correlation values of predicted rainfall with respect to GPM estimates, suggesting a model deficiency to represent the spatial distribution of precipitation for the selected study case. When spatial correlations were computed versus the station data, higher values emerged. Here also BMJ combinations were better. The experiments that use the Gr scheme stood out as the worst, showing no correlations or negative correlations with observations from stations.
All experiments represented the vertical humidity profile close to the sounding. It is appreciable that BMJ and KF schemes combined with all BMP schemes yielded the best humidity profiles. BMJ showed the most appropriate behavior of humidity at lower levels. While most BMJ experiments underestimated the wind profiles at lower levels, the rest of the BMP-CU combinations and also WSM6-BMJ exceeded the intensity values. For medium levels BMJ represented better the changes in wind velocity with the exception also of WSM6-BMJ. In general, most of the combinations followed the profile peak above 200 hPa level reported by the sounding.
A general feature that stands out from all the tests performed is that changes in cumulus parameterizations have a higher impact on the model’s performance statistics than those in microphysics, with BMJ scheme yielding the best results.
This research was funded thank to 2018-4-IOMA17-011 project “Analysis of the WRF-ARW Numerical Model for the prediction of rainfall at the basin scale in Panama” financed by the National Secretariat of Science Technology and Innovation of Panama (SENACYT), coordinated by the Water Center for the Humid Tropics of Latin America and the Caribbean (CATHALAC) and with the participation of the Center for Atmospheric Physics, Institute of Meteorology of the Republic of Cuba (INSMET). Authors want to thank Johisy Yasury Bethancourt from Dirección de Hidrometeorología, Empresa de Transmisión Eléctrica S.A.
 Nakaegawa, T., Pinzon, R., Fabrega, J., Cuevas, J.A., De Lima, H.A., Cordoba, E., Nakayama, K., Batista Lao, J.I., Lau Melo, A. and Gonzalez, D.A. (2019) Seasonal Changes of the Diurnal Variation of Precipitation in the Upper Río Chagres Basin, Panamá. PLoS ONE, 14, e0224662.
 Durán-Quesada, A.M., Reboita, M. and Gimeno, L. (2012) Precipitation in Tropical America and the Associated Sources of Moisture: A Short Review. Hydrological Sciences Journal, 57, 612-624.
 Mapes, B.E., Warner, T.T., Xu, M. and Negri, A.J. (2003) Diurnal Patterns of Rainfall in Northwestern South America. Part I: Observations and Context. Monthly Weather Review, 131, 799-812.
 Warner, T.T., Mapes, B.E. and Xu, M. (2003) Diurnal Patterns of Rainfall in Northwestern South America. Part II: Model simulations. Monthly Weather Review, 131, 813-829.
 Moya-álvarez, A.S., Gálvez, J., Holguín, A., Estevan, R., Kumar, S., Villalobos, E., Martínez-Castro, D. and Silva, Y. (2018) Extreme Rainfall Forecast with the WRF-ARW Model in the Central Andes of Peru. Atmosphere, 9, 362.
 Patel, P., Ghosh, S., Kaginalkar, A., Islam, S. and Karmakar, S. (2019) Performance Evaluation of WRF for Extreme Flood Forecasts in a Coastal Urban Environment. Atmospheric Research, 223, 39-48.
 Baklanov, A., Schlünzen, K., Suppan, P., Baldasano, J., Brunner, D., Aksoyoglu, S., Ca michael, G., Douros, J., Flemming, J. and Forkel, R. (2013) Online Coupled Regional Meteorology-Chemistry Models in Europe: Current Status and Prospects. Atmospheric Chemistry and Physics, Discussions, 13, 12541-12724.
 Givati, A., Gochis, D., Rummler, T. and Kunstmann, H. (2016) Comparing One-Way and Two-Way Coupled Hydrometeorological Forecasting Systems for Flood Forecasting in the Mediterranean Region. Hydrology, 3, 19.
 Simpson, C.C., Sharples, J.J., Evans, J.P. and McCabe, M.F. (2013) Large Eddy Simulation of Atypical Wildland Fire Spread on Leeward Slopes. International Journal of Wildland Fire, 22, 599-614.
 Bhomia, S., Kumar, P. and Kishtawal, C.M. (2019) Evaluation of the Weather Research and Forecasting Model Forecasts for Indian Summer Monsoon Rainfall of 2014 Using Ground Based Observations. Asia-Pacific Journal of Atmospheric Sciences, 55, 617-628.
 Sierra-Lorenzo, M., Ferrer-Hernández, A.L., Valdés-Hernández R., González-Mayor, Y., Cruz-Rodríguez, R.C., Borrajero-Montejo, I., Rodríguez-Genó, C.F., Quintana-Rodríguez, N. and Roque-Carrasco, A. (2015) Sistema Automático de Predicción a Mesoescala de Cuatro Ciclos Diarios.
 Sierra-Lorenzo, M., Borrajero-Montejo, I., Ferrer-Hernández, A.L., Morfa-ávalos, Y., Morejón-Loyola, Y. and Hinojosa-Fernández, M. (2017) Estudios de Sensibilidad del SisPI a Cambios de la PBL, la Cantidad de Niveles Verticales y, las Parametrizaciones de Microfisica y Cúmulos, a muy Alta Resolución.
 Mitrani-Arenal, I., Perez-Bello, A., Cabrales-Infante, J., Povea-Perez, Y., Hernandez-Gonzalez, M. and Diaz-Rodriguez, O.O. (2019) Coastal Flood Forecast in Cuba, Due to Hurricanes, Using a Combination of Numerical Models. Revista Cubana de Meteorología, 25, 121-138.
 Hsiao, L.-F., Chen, D.-S., Kuo, Y.-H., Guo, Y.-R., Yeh, T.-C., Hong, J.S., Fong, C.-T. and Lee, C.-S. (2012) Application of WRF 3DVAR to Operational Typhoon Prediction in Taiwan: Impact of Outer Loop and Partial Cycling Approaches. Weather Forecast, 27, 1249-1263.
 Poleo, D., Vindas, C. and Stoltz, W. (2017) Comparación y Evaluación de Diferentes Esquemas de Parametrización de Cúmulos con WRF EMS Aplicadas al Caso del Huracán Otto. Tópicos Meteorológicos y Oceanográficos, 16, 28-40.
 Lin, Y.-L., Farley, R.D. and Orville, H.D. (1983) Bulk Parameterization of the Snow Field in a Cloud Model. Journal of Applied Meteorology and Climatology, 22, 1065-1092. https://doi.org/10.1175/1520-0450(1983)022<1065:BPOTSF>2.0.CO;2
 Thompson, G., Rasmussen, R.M. and Manning, K. (2004) Explicit Forecasts of Winter Precipitation Using an Improved Bulk Microphysics Scheme. Part I: Description and Sensitivity Analysis. Monthly Weather Review, 132, 519-542.
 Thompson, G., Field, P.R., Rasmussen, R.M. and Hall, W.D. (2008) Explicit Forecasts of Winter Precipitation Using an Improved Bulk Microphysics Scheme. Part II: Implementation of a New Snow Parameterization. Monthly Weather Review, 136, 5095-5115.
 Lim, K.-S.S. and Hong, S.-Y. (2010) Development of an Effective Double-Moment Cloud Microphysics Scheme with Prognostic Cloud Condensation Nuclei (CCN) for Weather and Climate Models. Monthly Weather Review, 138, 1587-1612.
 Morrison, H., Thompson, G. and Tatarskii, V. (2009) Impact of Cloud Microphysics on the Development of Trailing Stratiform Precipitation in a Simulated Squall Line: Comparison of One- and Two-Moment Schemes. Monthly Weather Review, 137, 991-1007.
 Janjić, Z.I. (1994) The Step-Mountain Eta Coordinate Model: Further Developments of the Convection, Viscous Sublayer, and Turbulence Closure Schemes. Monthly Weather Review, 122, 927-945.
 Grell, G.A. and Freitas, S.R. (2014) A Scale and Aerosol Aware Stochastic Convective Parameterization for Weather and Air Quality Modeling. Atmospheric Chemistry and Physics, 14, 5233-5250.
 Lim, J.-O.J., Hong, S. and Dudhia, J. (2004) The WRF Single-Moment-Microphysics Scheme and Its Evaluation of the Simulation of Mesoscale Convective Systems. 20th Conference on Weather Analysis and Forecasting/16th Conference on Numerical Weather Prediction, Seattle, 1-15 January 2004, 1-4.
 Rutledge, S.A. and Hobbs, P.V. (1984) The Mesoscale and Microscale Structure and Organization of Clouds and Precipitation in Midlatitude Cyclones. XII: A Diagnostic Modeling Study of Precipitation Development in Narrow Cold-Frontal Rainbands. Journal of the Atmospheric Sciences, 41, 2949-2972.
 Betts, A.K. (1986) A New Convective Adjustment Scheme. Part I: Observational and Theoretical Basis. Quarterly Journal of the Royal Meteorological Society, 112, 677-691.
 Betts, A. and Miller, M. (1986) A New Convective Adjustment Scheme. Part II: Single Column Tests Using GATE Wave, BOMEX, ATEX and Arctic Air-Mass Data Sets. Quarterly Journal of the Royal Meteorological Society, 112, 693-709.
 Arakawa, A., Jung, J.-H. and Wu, C.-M. (2011) Toward Unification of the Multiscale Modeling of the Atmosphere. Atmospheric Chemistry and Physics, Discussions, 11, 3731-3742.
 Iacono, M.J., Delamere, J.S., Mlawer, E.J., Shephard, M.W., Clough, S.A. and Collins, W.D. (2008) Radiative Forcing by Long-Lived Greenhouse Gases: Calculations with the AER Radiative Transfer Models. Journal of Geophysical Research: Atmospheres, 113, D13103.
 Ruiz-Arias, J.A., Dudhia, J., Santos-Alamillos, F.J. and Pozo-Vázquez, D. (2013) Surface Clear-Sky Shortwave Radiative Closure Intercomparisons in the Weather Research and Forecasting Model. Journal of Geophysical Research: Atmospheres, 118, 9901-9913.
 Ek, M.B., Mitchell, K.E., Lin, Y., Rogers, E., Grunmann, P., Koren, V. and Tarpley, J.D. (2003) Implementation of Noah Land Surface Model Advances in the National Centers for Environmental Prediction Operational Mesoscale Eta Model. Journal of Geophysical Research: Atmospheres, 108, 1-16.
 Sukoriansky, S., Galperin, B. and Perov, V. (2005) Application of a New Spectral Theory of Stably Stratified Turbulence to the Atmospheric Boundary Layer over Sea Ice. Boundary-Layer Meteorology, 117, 231-257.
 Dee, D.P., Uppala, S.M., Simmons, A., Berrisford, P., Poli, P., Kobayashi, S., Andrae, U., Balmaseda, M., Balsamo, G. and Bauer, P. (2011) The ERA-Interim Reanalysis: Configuration and Performance of the Data Assimilation System. Quarterly Journal of the Royal Meteorological Society, 137, 553-597.
 Hou, A.Y., Kakar, R.K., Neeck, S., Azarbarzin, A.A., Kummerow, C.D., Kojima, M., Oki, R., Nakamura, K. and Iguchi, T. (2014) The Global Precipitation Measurement Mission. Bulletin of the American Meteorological Society, 95, 701-722.
 Dorninger, M., Friederichs, P., Wahl, S., Mittermaier, M.P., Marsigli, C. and Brown, B.G. (2018) Forecast Verification Methods Across Time and Space Scales—Part I. Meteorologische Zeitschrift, 27, 433-434.
 Gbode, I.E., Dudhia, J., Ogunjobi, K.O. and Ajayi, V.O. (2019) Sensitivity of Different Physics Schemes in the WRF Model during a West African Monsoon Regime. Theoretical and Applied Climatology, 136, 733-751.
 Gallai, I., Giaiotti, D.B., Gladich, I. and Stel, F. (2008) On the Onset and Evolution of Deep Moist Convection over Areas Characterized by Complex Orography: The Case of Friuli Venezia Giulia; Università di Trento. Dipartimento di ingegneria civile e ambientale.