The Impact of Microphysics Parameterization in the Simulation of Two Convective Rainfall Events over the Central Andes of Peru Using WRF-ARW

Abstract

The present study explores the cloud microphysics (MPs) impact on the simulation of two convective rainfall events (CREs) over the complex topography of Andes mountains, using the Weather Research and Forecasting- Advanced Research (WRF-ARW) model. The events occurred on December 29 2015 (CRE1) and January 7 2016 (CRE2). Six microphysical parameterizations (MPPs) (Thompson, WSM6, Morrison, Goddard, Milbrandt and Lin) were tested, which had been previously applied in complex orography areas. The one-way nesting technique was applied to four domains, with horizontal resolutions of 18, 6, and 3 km for the outer ones, in which cumulus and MP parameterizations were applied, while for the innermost domain, with a resolution of 0.75 km, only MP parameterization was used. It was integrated for 36 h with National Centers for Environmental Prediction (NCEP Final Operational Global Analysis (NFL) initial conditions at 00:00 UTC (Coordinated Universal Time). The simulations were verified using Geostationary Operational Environmental Satellites (GOES) brightness temperature, Ka band cloud radar, and surface meteorology variables observed at the Huancayo Observatory. All the MPPs detected the surface temperature signature of the CREs, but for CRE2, it was underestimated during its lifetime in its vicinity, matching well after the simulated event. For CRE1, all the schemes gave good estimations of 24 h precipitation, but for CRE2, Goddard and Milbrandt underestimated the 24 h precipitation in the inner domain. The Morrison and Lin configurations reproduced the general dynamics of the development of cloud systems for the two case studies. The vertical profiles of the hydrometeors simulated by different schemes showed significant differences. The best performance of the Morrison scheme for both case studies may be related to its ability to simulate the role of graupel in precipitation formation. The analysis of the maximum reflectivity field, cloud top distribution, and vertical structure of the simulated cloud field also shows that the Morrison parameterization reproduced the convective systems consistently with observations.

1. Introduction

The study of storms producing occasional rainfall and hail in the valleys of the central Andes is difficult, as it results from the interaction of flows from the Amazon and the Pacific at different levels, modulated by local complex orography conditions [1,2]. Particularly, in the area of the Mantaro Valley, the tropical conditions influenced by the easterly warm and humid wind from the Amazon, added to a strong heating of the soil during the afternoon hours, propitiate the convective lift, strongly interacting with local circulation, as a factor of an increase in convective development [3,4].

The Mantaro basin is surrounded by the central Andes of Peru, with altitudes of up to more than 5300 m. It has great economic and social importance for having a rich agriculture, as well as water resources, of considerable weight in the generation of electricity and the supply of drinking water. This favors the presence of numerous urban centers in the region, highlighting the city of Huancayo, located in the valley [5,6]. Sulca et al. [4] studied the meteorological conditions causing the occurrence of continuous periods with or without rainfall for the Mantaro basin, noticing that the anomalies of the easterly winds at higher levels are not always enough to produce periods of rain in the basin. In addition to this, the low level circulation from the southeast of the mountain range must provide an additional factor to guarantee the necessary orographic forcing to maintain significant periods of rain. From the conclusions of this work, it is inferred that an investigation of the factors influencing the development of precipitation in the Mantaro basin requires the application of tools that allow a deepening of the complex interactions involved. The Weather Research and Forecasting (WRF) model, in its advanced research variant (ARW) [7], is one of the main tools for the study of the mechanisms of the development of convection and the internal structure of storms in the valley, as well as for precipitation forecasting. The right configuration of the model, according to the conditions of the domain, is fundamental to achieve successful simulations, as has been shown for different parts of the world with diverse meteorological conditions [8,9,10,11,12,13,14,15,16].

Rajevan et al. [9] conducted a sensitivity study on the ability of WRF to reproduce the characteristics of a convective storm in southern India. They applied four microphysical parameterizations, among them the original version of Thompson’s scheme [17], with single predicted moments, which turned out to be the most effective. More recently, Orr et al. [15] studied the impact of four microphysics parameterization schemes on rain forecast and concentrations of different types of hydrometeors, for the Himalayan region. This work showed the good performance of Morrison’s scheme [18,19], and to a lesser extent, from the second version of the Thompson scheme [20] in complex orography conditions. Morrison also showed a good performance in reproducing rainfall in a simulation of a pre-monsoon convective storm, over the foothills of the Himalaya, Nepal, in the study by Shrestha et al. [13], who applied a configuration of WRF with an inner nested domain of a 3 km resolution. Using a similar resolution (4 km) to investigate the impact of microphysical parameterization in the simulation on an extreme precipitation event in the Himalaya, with a unique domain configuration, Karki et al. [16] found that the Thompson scheme [20] gave the best rainfall estimation. The schemes of Milbrandt-Yau [21,22] and Thompson were among the best performing in the simulations of the extreme rain event that caused flooding in California in 1997 [12]. Milbrandt and Morrison also provided good results in [11].

Moya et al. [23] developed a sensitivity study on different microphysics, cumulus, and boundary layer parameterizations for short- and medium-term precipitation forecast in the Central Andes of Peru, focusing on the Mantaro basin, for a set of three nested domains, (18, 6, and 3 km). This domain configuration has been further discussed and augmented in [24] and applied to extreme precipitation forecasting in the Central Andes in [25]. The tested configurations were Thompson [20], Lin [26], and Morrison [18,19] As a result, a configuration of the model was recommended, based on the minimization of statistical estimators of agreement with precipitation observations.

A sensitivity test using seven different MPPs for two nested domains at a 9 km and 3 km resolution was developed by Comin et al. [27] using WRF. They showed that WSM6 performs better over the highly complex topography of the Southern Andes of Chile for an extreme snowfall event, if boundary and initial conditions are taken from National Centers for Environmental Prediction (NCEP) Final Operational Global Analysis (NFL), When using ERA Interim [European Centre for Medium-Range Weather Forecasts (ECMWF) Re-Analysis] as driving data, the Goddard scheme [28,29] performed best. In another recent investigation [30], also for the Chilean Andes, three WRF MP parameterizations (MPPs) were tested for precipitation forecasting for three nested domains, with horizontal resolutions of 81, 54, and 9 km, finding the best results for the single-moment MPP WSM6 [31].

Flores-Rojas et al. (2019) [32] analyzed the dynamical mechanisms of the formation and development of three rainfall events in the Mantaro basin using three nested WRF domains, with horizontal resolutions 18, 6, and 2 km, with NCEP Global Final Analysis (FNL) data as boundary conditions and Lin as the microphysics scheme. The Advanced Regional Prediction System (ARPS) was used in an inner domain with a resolution of 0.5 km to investigate the microphysical composition of the systems, using the Schultz microphysical parameterization [33]. The events developed under the influence of the interaction of moisture fluxes from the Pacific Ocean and the thermally induced moisture fluxes coming from the South American low level jet at the east side of the Andes. The possibility of generation of inertial gravity waves in the Amazon basin was shown, creating conditions for intensification of convective cloud systems in the Mantaro valley.

The present paper is focused in the impact of the choice of microphysical scheme in the WRF reproduction of the structure of two cloud systems producing precipitation in the region of the Mantaro Basin. The conditions and dynamical mechanisms of development of these systems were discussed by Flores-Rojas et al. (2019) [32], so that in this paper, we aimed to analyze the detailed evolution of the systems and its microstructure by applying six different microphysical parameterization schemes. The WRF performance in reproducing and forecasting precipitation characteristics in the region was previously investigated by Moya et al. [23,24,25], by applying a similar domain configuration, but in this paper, an inner domain was introduced, covering the Mantaro Valley, with a 750 m resolution. The convective rainfall events (CREs) occurred in the afternoons of December 29 2015 (CRE1) and January 7 2016 (CRE2), and can be considered as typical of the Central Andes region. The CREs were documented with the observations of the Huancayo Observatory (HYO) and satellite databases.

2. Model and Data

The simulations were carried out using the WRF_ARW V. 3.7. model. The initial and boundary conditions were taken from the “Global Operational Analysis” of the National Center of Environmental Prediction (NCEP), final analysis FNL (https://rda.ucar.edu/datasets/ds083.2/), every 6 h, with a horizontal resolution of 1°. The analysis contains surface data and 26 mandatory levels. The included variables are surface pressure, sea level pressure, geopotential height, temperature, sea surface temperature, soil parameters, ice layer, relative humidity, zonal and meridional wind components, and vertical wind velocity. The SRTM digital elevation model NASA/NGA (Shuttle Radar Topography Mission [34] (https://dds.cr.usgs.gov/srtm/version2_1/)) at a 90 m resolution was used for topography.

Four domains were nested, whose characteristics are specified in

Table 1. Main parameters of the four nested domains.

Parameter	Domain 1	Domain 2	Domain 3	Domain 4
Central point	Lat: 10° S Lon: 75° W	Lat: 12.26° S Lon: 74.83° W	Lat: 12.37° S Lon: 75.03° W	Lat: −75.36° S Lon: −11.96° W
Horizontal resolution	18 km	6 km	3 km	0.75 km
Dimensions (XYZ)	115 × 140 × 28	115 × 142 × 28	127 × 163 × 28	113 × 121 × 28
Time step	50 s	16.6 s	8.3 s	2.1 s
Initial and boundary conditions	FNL 1.0 × 1.0° analysis	WRF output from Domain 1	WRF output from Domain 2	WRF output from Domain 3

, by applying the one-way nesting technique. The set of domains was designed following our previous paper [23], wherein three nested domains were used (18, 6, and 3 km), but including in this case, a fourth inner nested domain, with a resolution of 750 m, covering the Mantaro valley and its bordering mountain slopes (Figure 1a). The analysis of the evolution of the structure of the cloud systems was focused on this domain. The MPPs tested in this paper were those applied in [25], complemented by three more schemes that were recently used in similar orographic conditions.

Contour lines of the inner domain are shown in more detail in Figure 1b. The simulation period was 36 h for each experiment, extending from 00:00 UTC of January 7 2016 (19 LST of January 06) to 12:00 UTC the next day (07 LST of January 8), leaving approximately 12 h as spin-up time before the sunrise on the day of interest, so that the compatibility of the atmospheric variables of each of the interior domains with the boundary conditions of the immediate outer domain was guaranteed well before the start time of convection.

The Grell and Freitas [35] scheme was applied as convective parametrization, which consists of a modification of the “ensemble” scheme of Grell and Devenyi [36], explicitly representing the updraft and downdraft, and including the entrainment and detrainment of water and cloud ice. In domains 1 to 3, cumulus parameterization was activated, while in high resolution domain 4, only the microphysical scheme was applied [37,38]. In the case of domain 3, the resolution was higher than the 5 km limit proposed by Gilliand and Rowe [39], but it was decided to include cumulus parameterization, following the results of the experiment developed in [25]. For the boundary layer, the YSU (Yonsei University Scheme) scheme was used [40], and for the surface layer, the MM5 Similarity Scheme [41,42,43,44,45,46] was applied.

The RRTMG (rapid radiative transfer model) [45] was used as the radiation scheme. The soil model was the Unified Noah Land Surface Model, which is a result of the collaboration between NCEP, Air Force Weather Agency (AFWA), the National Center for Atmospheric Research (NCAR), and the Oregon State University (OSU), coupled to WRF with reasonable performance in its ability to capture heterogeneity in surface heat fluxes [46].

To test the impact of MPPs in the model and the ability to reproduce the cloud and precipitation evolution in the case studies, six schemes were applied, three of which had been previously tested in [23] (Lin [26] and Thompson [20], which are single-moment and Morrison [18,19], which is a two moment scheme). Considering their success in simulations of precipitation events or precipitation forecasting in orographic or tropical regions, three more schemes were tested: The two-moment Milbrandt and Yau scheme [21,22] and the single-moment schemes, Goddard [28,29] and WSM6 [31], which are associated with good performance in previous work in complex orography situations [11,12]

The single-moment schemes, used in most models up to the first decade of the present century, consist of equations for transportation and budget balance of each of the types of hydrometeors present as prognostic variables in the model. Generally, the third order moment is predicted for the distribution of each hydrometeor, as it is related to its total mass. In this type of parameterization, the concentrations are fixed a priori or are diagnosed from theoretical or empirical considerations, starting from the predicted moment. The latter is a limitation, since in many cases, it is not true that there is a one-to-one correspondence between the concentration of a type of hydrometeor and its mixing ratio, because the microphysical processes involve both the union of hydrometeors of the same type, and the diffusion growth of drops or ice particles from water vapor. These processes can imply the simultaneous growth of a moment and decrease of another. However, sensitivity studies show the best results in some cases for two-moment schemes [11,12,15,23], and in other cases for single-moment ones, as in [9,16,27,29].

In the other three tested schemes, the particle distribution of the hydrometeors is described in general by means of three-parameter gamma distribution functions:

$$N_i\left(D\right)=N_{0i}{D}^{{\mu }_i} e^{-{\lambda }_{i }D},$$

where the index, i, represents the different hydrometeors, and N₀, μ, and λ are the parameters of the distribution. The μ parameter is associated with the form of the distribution. If μ = 0, the distribution reduces to an exponential.

The different microphysical schemes deal with the determination of distribution parameters with different approaches and under different assumptions for each of the hydrometeors. The Thompson scheme [20] applied in this work is an improved version of the one-moment Thompson [17] parameterization, including improvements in both the physical processes and the calculation scheme. The assumed distribution of snow size depends on both the water and ice contents and temperature. It is represented as a sum of the exponential and gamma distributions and does not assume sphericity for snow particles. It also includes the prediction of cloud ice concentration, so it can be considered as a two-moment scheme only with respect to ice particles.

In the Milbrandt-Yau scheme, the mixing ratios and the concentrations of six types of hydrometeors are predicted separately: Cloud drops, cloud ice, raindrops, snow, graupel, and hail. For each type of hydrometeor, No and λ are predicted and the shape parameter is diagnosed, based on physical considerations, from the volumetric mean diameter, which is in turn obtained as a function of the other two parameters of the distribution.

Morrison’s scheme [18,19] is also double-moment, predicting the mixing ratios and concentrations of five types of hydrometeors: Cloud drops, cloud ice, snow, rain, and graupel. The shape parameter of the droplet gamma distribution is diagnosed, being bounded to dimensions that guarantee that the dispersion of the distribution remains within realistic limits. The scheme uses expressions for the calculation of homogeneous and heterogeneous nucleation processes as a function of temperature and supersaturation, based on Kohler’s theory [37]. For the calculation of the vertical velocity, an algorithm is used for which the sub-grid component is parameterized in terms of the turbulent diffusivity and mixing length. The droplet concentration is calculated from a quasi-stationary saturation adjustment algorithm, while the excess saturation feeds a non-stationary vapor deposition algorithm, allowing the growth of the ice crystals. The form parameter is fixed to zero for the rest of hydrometeors, including raindrops.

For the verification of the simulations, the measurements of the instruments of the Laboratory of Atmospheric Microphysics and Radiation (LAMAR) of the Geophysical Institute of Peru, and the Huayao meteorological station of the National Meteorological and Hydrological Service (SENAMHI), both located in the Huancayo Observatory (HYO) [12°02′18″ S (12.04° S), 75°19′22″ W (75.32° W), 3350 mamsl; UTC-5] were used. Radar MIRA 35-C, with a frequency of 34.85 GHz (wavelength of 8.6 mm), range resolution of 31 m, and temporal resolution of 5.6 s, is a polarimetric ka band, static Doppler radar directed vertically upwards, allowing registration of the reflectivity and vertical velocity profiles of the cloud and precipitation particles, as they move over the radar. It has a minimum valid height of measurement of 200 m a.g.l. and a beam width of 0.6° and high sensitivity, which allows detection of the reflectivity in the clouds before and during the development of precipitation, even if its relatively high attenuation for precipitation sized particles limits the accuracy of the reflectivity profile at medium and high levels in the case of high rainfall rates.

The general distribution of cloudiness was estimated from the brightness temperature of the infrared channel of the NOAA GOES satellite (National Oceanic and Atmospheric Administration, USA. Geostationary Operational Environmental Satellite; https://goes.gsfc.nasa.gov/), with a temporal resolution of 30 minutes and a spatial resolution of 4 km. The space-temporal rainfall distribution was estimated from the TRMM 3 hourly rainfall product, 3B42 [47] and the daily rainfall of the Climate Hazards Group Infrared Precipitation with Station gridded data set [48], also with 4 km resolution.

3. Results and Discussion

3.1. Description of the Rainfall Events and Their Conditions of Development

3.1.1. Evolution of Cloud Systems

To describe the general behavior of cloudiness on the dates of interest, the 3 h surface rainfall provided by the Tropical Rainfall Measuring Mission program (TRMM-3B42 product) [47] from 12:00 UTC (07 LST) to 06:00 UTC (01 LST) on the next day was used. Figure 2a shows the evolution of 3 h rainfall accumulated for a window over the Mantaro Basin on December 29 2015 (upper panels) and January 7 2016 (lower panels). The 12 h accumulated rainfall distribution over the basin for the two case studies in the period of convection of each of the study days is shown in Figure 2b. The brightness temperature observed from the infrared channel of the NOAA GOES satellite (National Oceanic and Atmospheric Administration, USA) of a space resolution of 4 km and a time resolution of 30 min was used to follow the evolution of the cloud systems developing in the basin and its neighborhood (Figure 3).

The occurrence of extended light rainfall surrounding convective cells can be noticed for December 29 2015 in the GOES data, since the morning hours at the west of the basin. Dissipation occurred between 15:00 and 18:00 UTC, and further intensification in the afternoon, as shown in the 21:00 and 00:00 UTC fields, where a strong convective cell can be seen in the basin, reaching the west border of the valley, and dissipating before 03:00 UTC. Convective cells were observed in the valley since the middle of the afternoon, immersed in the upper levels of south-eastern flow. They enter the Mantaro basin from the leeward slopes of the western mountain range, including associated stratiform precipitation. The presence of high cloud tops over the Mantaro valley can be noticed in the GOES imagery from 21:45 to 23:45, corresponding to CRE1. A convective cell with a brightness temperature of nearly −60 °C and cloud top of almost 11 km can be noticed over the valley in the 23:15 UTC image, surrounded by apparently stratiform clouds.

In the afternoon of January 7 2016, two main rainfall events developed in the Mantaro Basin region, conditioned by the interaction of the orographic flow with the general wind circulation and solar heating. At 18:00 UTC (13:00 UTC LST), there are rainfall areas at the northern extreme of the basin, while by the southwestern border, a moderate rainfall band can be observed. The maximum rainfall rate occurs at 21:00 UTC, with a peak value of ~6 mm/h, which is higher than climatology (~2 mm/day) and at 00:00 UTC, the rainfall rate has decreased, and the band has moved to the east. The 12 h accumulated rainfall in most of the basin reached 10 to 15 mm. A convective cloud system developed since the morning in the eastern border of the basin, fed by moist air from a storm area in the Amazon. This system can be observed in its maximum development in Figure 3, in the panels of 18:45 UTC and 19:45 UTC, and the 20:45 UTC shows that its remains. However, almost no rainfall is observed in the basin from this system (Figure 2). On the other hand, also since the middle of the morning, along the western border of the basin, a line of separate convective cells develops, gradually covering the boundary of the basin. The 18:45 image shows that, in the afternoon, most of the line is over the basin, developing and gradually moving to its center, so that at 20:45, the maximum top height covers the Mantaro valley, corresponding to CRE2 (Figure 3; bottom left panel), with a minimum brightness temperature reaching −60 °C, with cloud top altitudes of more than 11 km. This is consistent with the extended rainfall maximum observed at 21:00 UTC by TRMM. Afterwards, the rain band continues moving towards the east, dissipating in the early evening.

3.1.2. Large Scale Situation

The south-east Pacific subtropical anticyclone (SPSA) is a dominant force in the South Pacific basin. During the austral summer monsoon (December to March), SPSA intensifies at the coast of central-southern Chile (35° S). Because of subsidence generated by SPSA, arid and stable conditions prevail along the coasts of northern Chile and southern Peru and extend inland over the western slope of the Andes [49], which are sensitive to upper level large scale circulation [4]. The Bolivian high (BH) and the north-east low (NL) system are the main circulations over South America during the austral summer [50]. Strengthening of BH accelerates the upslope easterly flow at mid and upper levels and produces 80% of the annual precipitation (350–400 mm) over the central Andes during the austral summer monsoon. The weakening and strengthening of BH-NL systems is responsible for the creation of dry and wet spells over the Mantaro basin [4]. The NCEP-NCAR reanalysis data of the geopotential height show the meteorological situation of both case studies (Figure 4). At low levels (850 hPa), there is a weak trough extending from the southeast, warming the lower troposphere and facilitating low-level convergence, while in 200 hPa, Bolivia’s high is well represented, with a large anticyclone circulation, with divergence at high levels, stimulating convective development in the region.

The synoptic situation at low levels is similar for the two cases. On December 29 2015, low levels are characterized by very low pressure gradients, which, like on January 7 2016, warm this part of the troposphere and trigger the circulation of local winds. However, the situation at high levels is not the same, since Bolivia’s high center is displaced to the north. Here, the fundamental trigger of convection lies in the intrusion of cold air from the southwest, provided by the trough that extends from the south, with its axis almost parallel to the coasts of Chile and Peru. Together with the warming of low levels, this cold air intrusion produces high instability in the region. In both cases, the synoptic situation is favorable to the occurrence of rainfall on the Mantaro basin. A more detailed discussion of the triggering mechanisms and conditions of development of these convective events in the Central Andes can be found in [32]. Satellite-based climatology of intense rainfall events has shown that in a favorable large scale situation, such as the Andes topography condition, the geographical locations of the formation and development of rainfall events over the continent [51].

3.1.3. Radar Records of the Convective Rainfall Events

On December 29 2015, in the 23:15 UTC GOES image (Figure 3), there is a convective core over the valley, which persists in the next 30 minutes, appearing weakened after an hour in the 00:15 UTC image. This image corresponds to CRE1. This convective cell appears at the time surrounded by apparently stratiform clouds. The radar record (Figure 5) showing high reflectivity at low levels (upper panel) between 23:00 and 23:50 corroborates this. Subsequently, brief periods of intense rainfall and possibly hail persist in the context of stratiform rain, with the appearance of a bright band at almost 2 km above ground level (a.g.l.), corresponding to nearly 5 km above sea level (a.s.l.). Both phenomena are also detected in the rain gauge record, where the smaller short duration peaks are consistent with low rainfall rate accumulation in the bucket and discrete recording, characteristic of stratiform rainfall (Figure 7). To interpret the ka band radar information, it must be considered that it represents an average of its measuring volume, enclosed by its beam width of 0.6° and its range resolution of 31 m, and that its precipitation estimation corresponds to its lower valid measuring volume, beginning at 200 m a.g.l.

In the radar record of CRE2, it can be noticed that its convective core did not reach HYO, as the low Ze values correspond to very low rainfall rates. However, the height of the cloud echo reaches 12 km a.g.l. (nearly14 km a.s.l.) A melting layer at almost the same height as in CRE1 (2 km a.g.l.) is present during the whole period, with intermittent light precipitation, and a brief interval of high reflectivity near to 21:00 UTC. Despite the cloud depth, most of the record is not characteristic of a convective cloud. The radar Doppler vertical velocity is presented in Figure 5 (lower panels), and only downward (in negative) motion is presented near to the surface, with a sudden increase at the melting layer. Turbulence, which manifests by a combination of up and downwards velocities, is noticed above, up to 10 km, between 19:30 and 21:30 UTC. A brief interval of maximum values of both, Ze and vertical velocity, appears at nearly 21 h with very high velocities (~15 m/s), which are usually related to hail that could have reached the radar beam carried by a downdraft. This picture is consistent with the passage over HYO of a peripheral partially stratified sector of a convective system, which extended approximately from 19:30 to 23:30 hours.

3.3. Simulation of Temperature and Wind Speed

Figure 6 shows the comparison between the time series of the ground-based observations (HYO) and WRF simulation for surface temperature (in °C) (a,c) and average wind speed (m/s) (b,d) for six MPPs for domain 4, within a 10 km × 10 km box centered in HYO (12.0S, 75.3W) for the two cases (a,b: CRE1; c,d: CRE2). The observations are at a 1-minute time step, whereas WRF output is every 10 minutes. So, the observations were smoothed by moving averaging every 10-min to match the simulations. For CRE1, the maximum of the observed surface air temperature (T2m) recorded at 21:10 UTC shows an abrupt change of slope, due to the cooling influence of downdrafts from the neighboring storm (Figure 6a). The WRF simulations approximately reproduce this behavior with time shifts of nearly 20 minutes in advance, different across schemes, wherein Morrison, Milbrandt, and Goddard reproduce the form of the observation time series quite well. The cooling rate in the observations increases again after the detection of the rainfall, which is reproduced by Morrison, after the occurrence of the simulated precipitation at 21:50. In the case of wind speed (Figure 6b), the simulations lay within the range of the variability of the observations, and no specific trend was found. For CRE2, the simulated temperature (Figure 6c) shows peaks of less than 18 °C for the six schemes, whereas the observations’ maximum reaches almost 22 °C. The observed maximum temperature occurs at 19:30 UTC, coincident with the beginning of the influence of the storm, as can be seen from Figure 5. Observed and simulated temperature profiles show a decrease in surface temperature during the rainfall event because of evaporative cooling affecting the neighborhood of the storm. All the simulations underestimate the surface temperature during the rainfall event (between 19:30 UTC to 23:30 UTC), matching better after the event, especially in the cases of Morrison and WSM6, and with some overestimation in the rest of the configurations. As for CRE1, the simulated and observed wind speed (Figure 6d) show similar trends in the observed and simulated wind profiles.

3.4. Rainfall Characteristics: Simulations and Data

To investigate whether the model can reproduce the precipitation consistently with the available data, eight points around HYO were interpolated from the model output to obtain the rainfall simulated by the model not for a fixed point, but for a more diffused region (Figure 1b). Figure 7 shows the time series of observed and simulated rainfall. The dotted lines indicate the rain gauge (red) and ground-based radar (blue) observations, whereas solid lines show the model output for different MPPs. CRE1 (Figure 7a) was reproduced in the neighborhood of HYO, roughly an hour in advance of its detection by the radar and the rain gauge. Morrison, Goddard, WSM6, and Milbrandt reproduced the rainfall over HYO to some extent, while Thompson and Lin practically failed to reproduce it in the selected box. Morrison, Milbrandt, and Goddard simulated the precipitation event approximately from 21:50 to 22:50, while WSM6 reproduced it one hour before. Morrison shows the most significant precipitation maximum, at 21:10 UTC, 1 h and 30 min earlier than the peak in the observations. A previous precipitation event was reproduced by the model earlier in the afternoon, which was not recorded by the instruments, but could have occurred near HYO.

For CRE2 (Figure 7b), the six MPPs are able to reproduce the occurrence of low intensity rainfall at and near the station, but either in advance or with a time lag relative to the rainfall event. As an example, Thompson and WSM6 reproduce two weak rainfall events, one before and one after the observations (between 20:25 UTC to 22:25 UTC). Morrison reproduces only one rainfall maximum, nearly 2 h after the peak rainfall in the observations, whereas Goddard reproduces the maximum rainfall 2 h earlier than the observations. Lin simulated the maximum rainfall nearly 1 h earlier than the observations. The MPPs either overestimate or underestimate in nearly the same amount the surface rainfall, compared to observations, except for Milbrandt, which simulates the least rainfall.

Root mean square errors (RMSEs, (a) and differences (b) between the simulated 24n h total precipitation results and the reference gridded dataset CHIRPS were calculated. The CHIRPS data (4 × 4 km) were interpolated into the model’s grid. Figure 8 shows the spatial distribution of RMSE (a) and the difference (b) between the simulation and CHIRPS data for the six MPPs. For CRE1, all the schemes show good performance within the valley and the surrounding basin, according to these criteria, even if WSM6 and Lin show a high RMSE in some limited regions of the basin. The difference is very low for the central and southern parts of the valley, while in the north, the model underestimated precipitation in some extent. For CRE2, some of the schemes show high local RMSE values over mountain areas, and, in general, the error is limited in the valley. The best performing schemes regarding RMSE are Morrison and Goddard, for which the RMSE is lower than 12 mm for most of the valley. The spatial field of the difference in 24 h rainfall shows that the model estimations depend strongly on the MP scheme, as Goddard and Milbrandt show underestimation in the valley while the rest of the schemes show different degrees of overestimation. The highest errors correspond to the mountainous areas, as observed by Moya et al. [25]. Morrison and Thompson show the lowest differences, near to zero in most of the valley. In general, looking at both parameters, and for the two cases, Morrison seems to be the best performing scheme under the criterion of spatial reproducibility of the precipitation field.

3.5. Simulation of Averaged Vertical Velocity

The temporal evolution of the zonally averaged vertical velocity, simulated for six MPPs in the neighborhood of HYO, is shown in Figure 9 at the 5.1 km WRF output height (a.g.l) within the longitudinal belt from −72.4 to −72.5, including HYO. For CRE1, Morrison, Milbrandt, Goddard, and Lin show intense vertical velocities near the latitude of HYO and the time of occurrence of the storm shifting southwards from approximately 21:00 to 22:00 UTC, consistent with the time of occurrence of the observed storm. For CRE2, all the MPPs show the northward propagation of the updraft zones related with convective activity, but only in the case of Morrison and Lin are the updraft zones simulated near HYO at the time of the observed event, consistent with the GOES information (Figure 3).

3.6. Simulation of Hydrometeors in the Mantaro Valley and its Surroundings

Figure 10 and Figure 11 show the time–height evolution of hydrometeors in WRF simulations for a 10 km × 10 km area, centered in HYO, for the 5.1 km WRF output level. The evolution of the graupel mixing ratio (Qg (g/kg)) is shown in Figure 10, superposed with the contour lines of vertical velocity in m/s (updraft in red and downdraft in green), for different MPPs. Surface rainfall is also plotted as a thick red line, related with the right Y axis scale. Both for CRE1 and CRE2, Figure 10 reveals the importance of graupel in the formation of the surface rainfall, characteristic of convective events, which is apparent from the correlation of the maximum graupel mixing ratio with maximum rainfall, consistent with the results of Karki et al. [16] over the central Himalaya. For CRE1, which is a more intense convective system, this holds only for WSM6, Morrison, and Lin, while for CRE2, Milbrandt and, to some extent, Goddard show qualitative correlation too. For CRE1, the time of maximum graupel production and surface rainfall coincides better with the GOES information and with the radar record for Morrison and WSM6, while Lin reproduced the maximum in advance. In the case of CRE2, the averaged maximum graupel was reproduced roughly at the same interval of time by the different schemes.

Figure 11 focuses on the evolution of the mixing ratios of ice (Q_i(g/kg): red contours) and snow (Q_s(g/kg): green contours), together with rain water (Q_r(g/kg): shaded). Compared to graupel, the snow maxima show less time correlation with the regions of more intense precipitation, but in the cases of WSM6, Morrison, and Lin, a significant snow mixing ratio is present in the intervals of high precipitation, even if the snow content persists after precipitation. This seems to indicate that the aggregation processes related with snow formation plays an important role in cloud formation, particularly in the trailing stratiform clouds persisting after convection, along with the riming processes related with graupel. The ice mixing ratio is generally an order of magnitude less, and its role on precipitation development in this case is mostly through its participation in the formation of snow.

Figure 12 shows the time and area averaged vertical profiles for six MPPs for different hydrometeors in domain 4 near the HYO observatory. There are large differences in the simulated vertical profiles of the mixing ratios of different water species across MPPs. To evaluate the ability of each of the schemes to produce different types of hydrometeors, the average profiles of the five predicted mixing ratios in each of the configurations were obtained. The space averaging was performed for a 10 × 10 km box, and the time averaging was performed for the time interval of the development of the simulated storms. The diverse relationship between the mixing ratio profiles of the different hydrometeors across configurations in the middle troposphere is noteworthy. For CRE1, the height of the absolute maximum of the profiles of graupel (Q_g) and snow (Q_s) for the six schemes is not very different, at a height somewhat lower than 6 km, except for Goddard, for which the snow peak is located a little higher. The profiles of the rest of the hydrometeors is very dependent on the microphysics scheme. In the case of Thompson, the mixing ratio of hydrometeors, except snow, is very low at all heights, resulting in low surface precipitation, which may be related to the limitations found in the single moment schemes [43]. On the other hand, in the two-moment schemes, Morrison and Milbrandt, the important contribution of graupel and cloud water is observed, which is consistent with the activation of graupel growth from the riming of supercooled cloud droplets in the presence of updrafts. The good performance of the Morrison scheme seems to be a consequence of its ability to produce rainwater (Q_r) from cloud water, simulating the coalescence process properly, since the maxima of both mixing ratios are coincident near 4.5 km, which favors the efficient formation of graupel by riming. For all schemes, but more markedly for Morrison, a secondary maximum of cloud water and rainwater is observed at a height of approximately 1 km above the surface, coincident with the zero degree Celsius isotherm for a large part of domain 4. An apparent deficiency of the Milbrandt scheme for the region is its inability to reproduce significant concentrations of cloud ice (Q_i) in the area near the top of the clouds, unlike Morrison, which has a clear maximum of this magnitude near 8 km. Probably, the efficient conversion of cloud ice to snow is the source that allows the reproduction of the snow profile, which is also valid for Lin. It is remarkable that recent studies in the Himalaya have found that the Thompson scheme reproduces best the precipitation features in that region, where snow plays a very important role. This is only an apparent contradiction with our results, as in both cases, Thompson produces less graupel than Milbrandt and Morrison, and in the Central Andes, the higher temperatures propitiate the predominant precipitation formation by the coalescence-riming-graupel mechanism. Milbrandt does not reproduce the maximum of rainwater at 4.5 km, which implies that the coalescence process is not represented correctly, preferring the conversion of ice and snow to graupel by the ice-riming-graupel mechanism, based in the Bergeron-Findeisen process [52]. This may be related with the lower production of graupel.

Morrison produces the greatest graupel mixing ratio, and its peak attains the greatest height, near to 4.5 km, while the cloud water and rain water profiles reach their maxima approximately 1 km below. This suggests an efficient coalescence and riming processes, wherein rain water grows on the expenses of cloud water and graupel grows from rain water, also feeding the rain water by melting in the downdrafts. The two-peak structure of the rain water profiles show the relative importance of the liquid phase and mixed phase mechanisms in rain formation, and in the case of Morrison, both peaks are well defined. For CRE2, Morrison, Milbrandt, and WSM6 show similar graupel profiles., but their rainwater profiles show a maximum near to the ground, and the higher peak is very small. In the case of WSM6, ice and snow peaks are close in height, which suggests an important contribution of crystal aggregation processes in rain formation, which is consistent with the formation of less graupel at high altitudes. On the other hand, the Lin ice and snow profiles show absolute maxima at higher altitudes, coincident with Morrison, but much less intense, which may be related to their sedimentation mechanisms. In CRE2, in the case of Milbrandt, the graupel formation is efficient but cloud and rain water is constrained to the lower layers, while the formation of snow at relatively low levels seems to limit the ice water content, and as a result, the precipitation formation was limited, conducting to rainfall underestimation. In the case of Goddard, the ice and snow mixing ratios have high values, but the rain water is restricted to low heights and the graupel content is relatively low. The snow process is relatively efficient, but the precipitation formation mechanisms did not grant the correct simulation of precipitation in the case study.

A detailed verification of the hydrometeor characteristics from WRF simulations is not possible due to the unavailability of this kind of data from ground-based or satellite measurements. Thus, the evolution of the simulated hydrometeor fields will be tested for consistency with the GOES estimated cloud fields. Considering that the coalescence and riming processes conducting the growth of raindrops and graupel particles is inherent to the development of convective updraft [53], the mixing ratio of the sum of both hydrometeors will be used as an indicator of the presence of convective cells [54]. Figure 13 and Figure 14 show the spatial evolution of Q_g+r for most of the period corresponding with the detection of rainfall by the radar for both cases (20:00 UTC to 23:00 UTC) using the different MPPs at the WRF output height of 5.1 km above ground.

For CRE1 (Figure 13), it was found that the model reproduced convective activity in the neighborhood of the valley with a duration comparable with the observations, even if with a certain time shift, for the Morrison, Milbrandt, and Lin schemes. Convective systems were simulated within the time interval from 20:00 to 23:00 UTC, nearly two hours in advance of the observed system, though its representation differed for the three microphysical schemes. Important differences in cloud microphysics and dynamics across schemes come out. In the case of Morrison, it reproduces two convective systems with intense and extended convective cores, approaching the valley from the west and carrying high Q_g+r. The cells at the northwest develop in the western mountains, reaching the valley in a mature stage, corresponding to the first precipitation peak in Figure 7a, while the cells to the south east develop mainly on the valley, reaching a maximum Q_g+r of more than 6 g/kg, corresponding to the second group of precipitation peaks in Figure 7a, which can be matched with the observed rainfall event, with a time advance of nearly an hour. Considering the total water content (not shown), the system practically covers all the valley, and particularly HYO. Milbrandt also produces an extended system in the western mountains, but it is less developed than in Morrison and does not approach the valley, but since 20:20 UTC, a convective cell develops in the valley, which produces some early rainfall, expands and eventually merges with the system in the mountains. Also, in this case, the inclusion of snow to the mixing rate (not shown) produces the almost total overcast of the valley in the afternoon and evening. In the case of Lin, the system is also reproduced, but with less Q_g+r and less extension_, so that it does not cover HYO. The advance in the onset of WRF simulated convection and precipitation relative to the observations with explicit convection was also observed in [9], with coarser resolution, but was not observed in [16]. Lean et al. [55] developed a series of experiments, including different cases and resolutions of the Met Office Unified Model, concluding that this problem can be related with resolution, so that in the present study, it originated in the outer coarser domains, and passed to the inner domain as a consequence of the nesting process.

Although the cross sections of the other vertical levels have been omitted for a lack of space, the three-dimensional analysis of the model output evidences the existence of convective cores with a deep vertical structure, developing in the presence of local convergence and carrying significant cloud water, snow, and graupel content. The simulated system rose over the valley at 21:20 UTC and persisted until 23:00 UTC, staying on the valley until 22:40 UTC. This is consistent with the observed GOES imagery. In particular, over the Mantaro Valley, the Morrison scheme produced a clearly convective cloud, carrying an abundant graupel mixing ratio of more than 4 g kg⁻¹ and updrafts of up to 12 m s⁻¹, as can be seen from the vertical cross sections shown below.

For CRE2 (Figure 14), all the MPPs produce convective activity over the domain too, but only WSM6, Morrison, and Lin reproduce the evolution of the cloud field shown by Figure 3 to some extent, simulating a convective cloud system entering the valley from the south western border and moving through the north east, covering the valley near or over HYO. In the case of Thompson, it simulates a system entering the valley from the north east, corresponding to the other system observed in the GOES imagery, which did not enter the basin, while Milbrandt and Goddard practically did not produce convection in the valley.

As can be inferred from Figure 7, Figure 10, Figure 11, and Figure 14, for CRE2, the different configurations reproduce the convective events in different time periods. In the case of WSM6 and Morrison, with a lag relative to the observations, and in the case of Lin, in advance.

From the above analysis, it follows that the WRF configurations that reproduced the general horizontal field of convection consistently with the observations for the two cases were Morrison and Lin, even if in the case of Lin, the system formed was displaced to the north relative to HYO. The simulated local wind field, conditioning the onset of cloud development for each of these two WRF configurations, is shown in Figure 15.

In the case of CRE1, the Morrison configuration developed an extended convective system in the east of the basin in the afternoon, propagating into the valley and producing graupel and rainwater at 21:20, more than an hour before the beginning of the radar rainfall event record, as a result of the convergence triggered by the associated downdraft of its wind system, which eventually developed into the simulated system over the valley, consistent with the general pattern of convective development from east to west apparent from the GOES imagery (Figure 3). On the other hand, the onset of Q_g+r occurred in Lin as a result of a quasi-isolated cell, fed by the moisture flow from a system developing in the right slopes of the basin.

For CRE2, even if the wind field evolution in the early afternoon was similar for the two configurations, the onset of Q_g+r near HYO initiating the development of the convective cloud system was different. The analysis of the model output of CRE2 showed that Lin produced significant Q_g+r greater than the 0.1 g/kg threshold, implemented in the figure color bar since 19:40 UTC, conditioned by the convergent flow related with a convective system in the eastern slopes of the surrounding basin. This is consistent with the beginning of the CRE2 radar trace, but the simulated system did not match the whole duration of the radar record. The initiation time of Q_g+r near HYO for Morrison was only half an hour before the radar record, but the simulated rainfall on the ground only matched the last part of the record. The situation is complicated by the existence of two simultaneous systems affecting the local wind field. Morrison simulated the intense convective system approaching the basin from the northeast, mainly from 20:00 to 22:30 UTC, and after its decline, began to reproduce the system in the center of the valley. In the case of Lin, it failed to reproduce the early system in the north, but in the 20:20 cross section, several isolated cells appeared between the northern border and the center of the valley, while a narrow band was formed in the eastern border, as a result of the convergent orographic flow, later developing the simulated system reaching HYO. After 21:30, it produced an extension of the convective zone, merging with a widespread convective system from the bordering mountains.

To describe the evolution of the simulated systems as a whole during the period of maximum development of both case studies, the horizontal projection field of simulated maximum reflectivity using the Morrison parameterization is shown in Figure 16, for the time period in which this parameter attained its maximum in both cases in the neighborhood of HYO. For CRE1, a core of more than 50 dBZ formed at 21:40 UTC at the east of HYO, reaching its maximum extension at 22:00 UTC, in the neighborhood of the observatory, and slightly weakening in the next 20 min. In the case of CRE 2, different simulated systems developed in the valley, at different places and times, as the one penetrating from the western border at nearly 19:00 UTC, reaching 40 dBZ and crossing the valley at the north of HYO, but the one attaining maximum reflectivity entered the valley at approximately 22:00 UTC, and attained more than 50 dBZ at the southeast of HYO, crossing the valley in less than an hour while losing intensity. Both case study simulations are consistent with the observed ka band radar record, detecting a convective system over the radar for CRE1 and the outskirts of a system for CRE2. Unfortunately, no scanning meteorological radar data field is available for comparison with the three-dimensional simulated reflectivity fields, but the maximum reflectivity distributions are consistent with the positions of the maximum tops in the GOES output, even if with a time shift (Figure 3).

To compare the evolution of the cloud tops simulated by the Morrison configuration with the GOES satellite pictures shown in Figure 3, horizontal cross sections of the total water mixing ratio (Q_T) (including the sum of the five hydrometeors) are plotted in Figure 17 for the WRF output levels of approximately 6.9, 8.7, and 10.5 km a.g.l. Contour lines within a minimum threshold of 0.2 g/kg are plotted for the lower levels (red for the cross section at 6.9 km. and black for 8.7 km) and shaded contours indicate the Q_T distribution at the upper level of 10.5 km, within a threshold of 0.1 g/kg. The upper panels indicate the evolution of the highest tops of CRE1 in the HYO neighborhood, and a second coexisting storm over the north of the valley. At 21:00 UTC, CRE1 shows tops higher than 10.5 km with Q_T > 0.7 g/kg and at the same time, the system to the north, and the 21:20 and 21:40 cross sections show a process of slow descent of the CRE1 tops while the storm in the north develops extended high tops with inner cells of relatively high Q_T, consistent with convective cells surrounded by anvil clouds. Both storms persist in the next two cross sections at 22:00 and 22:40 UTC, with tops higher than 9.7 km and a significant total water mixing ratio. These two high top systems correspond to an extended red spot over the valley in the 21:45 UTC GOES image, with a brightness temperature lower than −55 °C and corresponding to an approximate cloud top height of 10.5 km, persisting with a smaller area only at the center of the valley in the 23:15 image. In the case of CRE2 (lower panels of Figure 16), high cloud tops are simulated by WRF from 18:30 UTC crossing the valley, with cloud tops higher than 9.7 km at 19:30, when the higher tops have displaced to the eastern boundary of the valley. High cloud tops are shown at the 20:50 and 21:30 simulated cross sections. The two high cloud top systems were observed in the GOES imagery (Figure 3), but with a time shift and reversed order, as the highest cloud tops are observed first in the north of the valley in the 19:45 UTC image as a white spot, with brightness temperatures lower than −60 °C and tops higher than 11 km, and a similar, but more extended spot was observed over the center of the valley in the 20:45 UTC image.

To investigate the characteristics of the simulated updraft–downdraft structure and its relationship with the vertical distribution of hydrometeors using the Morrison scheme, the model output latitudinal cross sections at approximately 12° S latitude, crossing HYO between 75.65 and 75.2° W, were plotted every 5 min. The plots were reviewed to detect the period of maximum vertical development of the simulated system in the neighborhood of HYO. Figure 17 and Figure 18 show the cross sections for the times of the maximum vertical velocity and maximum mixing ratio of graupel plus rain water, as indicators of the maximum development period of the clouds.

For CRE1 (Figure 18), a vigorous convective system developed, including a rapidly growing cell with an updraft of more than 12 m/s at 21:45 UTC (not shown), producing high mixing ratios of cloud water and graupel and very fast graupel and rainfall development, producing a rapidly growing cell that reached 10.5 km of cloud top height and the maximum graupel mixing ratio of more than 2 g/kg at 21:55 UTC, when the vertical velocity decreased below 10 m/s because of the hydrometeor load but strengthened again with the liberation of latent heat from graupel formation so that a second vertical velocity maximum was attained at 22:10 UTC with the formation of a new cell. It can be seen that the cloud develops towards the east, with the growth of cloud water in the updraft at the east of the graupel core, while snow is spread all over the horizontal extent of the cloud, and extending itself to the east. Figure 19 shows that the vertical cross section of CRE2 reached a maximum vertical velocity of 5 m/s, coincident with the development of cloud water and maximum graupel mixing ratio of 1 g/kg. It can be noticed that in CRE1, cloud water grows at the same time of graupel, which is a sign of development, while in CRE2, the growth of graupel partially depleted cloud water, resulting in smaller duration of this system relative to CRE1.

4. Conclusions

In the present study, cloud-resolving simulations using WRF-ARW were used to reproduce two convective rainfall events, which developed in the Mantaro valley. Its main objective was to examine the impact of the microphysical parameterization scheme to simulate the core properties of convective rainfall events in the complex topography of the Andes. The study was conducted using six microphysical parametrization schemes to understand the structure and microphysics of convective rainfall events over the Mantaro valley and its surrounding mountainous area. The main characteristics and simulation results are listed below:

• During the development of the two convective rainfall events used as case studies, the Mantaro valley was under the influence of low and medium level flow from the north and northeast. At the same time, a high-level wind from the Pacific Ocean was flowing over the mountain range of the central Andes, interacting with the humid flow from the east Amazon and causing instability, which was higher in the case of CRE1.
• Ground-based observations show that the most rainfall occurred in the afternoon (after 19:00 UTC) and mainly after 20:00 UTC for CRE2 and after 23:00 UTC for CRE1. The ka band cloud radar captured the CREs from its early stage, in its central part in the case of CRE1, but only in its periphery in the case of CRE2. However, both events left their trace in the radar record, showing significant Ze and vertical velocity profiles, consistent with deep convective clouds lasting until nearly 02:00 UTC (21 LST) for CRE1 and 23:00 UTC (17 LST) for CRE2. In both cases, the convective systems coexist in at least part of their time over the radar with stratified rainfall from the outskirts of the storm, forming an almost continuous bright band.
• In the case of CRE1, wherein the core of the system passed over the radar, the temperature record was well reproduced by most of the configurations, particularly by Morrison and Milbrandt. For CRE2, for which the periphery of the system passed over the radar, all the microphysical schemes produced underestimation of the surface temperature during the rainfall event, and matched after the rainfall, showing the temperature drop related with the evaporative cooling effect of the rainfall.
• For CRE1, all the schemes give good estimations of 24 h precipitation, with relatively low RMSE, but for CRE2, Goddard and Milbrandt underestimated the 24 h precipitation in the inner domain, representative of the valley and its surrounding mountains, while the rest of the schemes overestimated precipitation. Morrison and Thompson showed the lowest precipitation difference fields while Morrison and Goddard showed the least root mean square errors, particularly over the valley.
• For CRE1, the Morrison, Milbrandt, and Lin configurations reproduced the general dynamics of the development of cloud systems, but in the case of Lin, it did not reproduce the system over HYO, so that its comparison with radar and rain gauge data was not possible. Regarding CRE2, only Morrison, WSM6, and Lin reproduced it. Lin simulated the event in advance and with shorter duration, while Morrison and WSM6 reproduced it with a one hour lag, but the duration of the simulated storm was comparable with the radar record and GOES imagery.
• The vertical profiles of the hydrometeors simulated by different schemes show significant differences, showing that the best performance of the Morrison scheme for both case studies may be related to their ability to simulate the role of graupel in precipitation formation.
• The analysis of the vertical structure of the simulated cloud field shows that the Morrison parameterization reproduced the convective systems in a way consistent with the observations.