Assessment of WRF microphysics and cumulus parameterizations in simulating heavy rainfall events over Badulu Oya catchment, Sri Lanka

The accurate prognosis at different temporal scales ranging from short and medium range (a few hours to some days ahead) can reduce the impacts of water-related hazards and provide the basis for more effective water resources management (Emmanouil et al. 2021). However, the forecasting task is more rigorous due to the complex interactions between the Earth and the atmosphere (Piciullo et al. 2018). There are well-known Global Climate Models (GCMs) that often have oversimplified physics, parameterized structures and coarser horizontal spatial resolution to simulate the interactions among various macro- and microphysical processes of a weather system. Regional Climate Models (RCMs) are applied to overcome the GCM's limitations which are dynamically forced by GCMs and have high spatial and temporal resolution and improved physics (such as convection) (Bliznak et al. 2019). With the advent of recent high-performance computing facilities, they can also resolve the regional features related to topography, local climate and the associated precipitation. All-weather centers and meteorological services employ the GCMs to simulate weather and climate over the area of interest at desired resolutions (Srinivas et al. 2013; Ji et al. 2014).

The WRF model (Weather Research and Forecasting model) is one of the most widely applied, cutting-edge RCMs designed to combine the expertise and experience of mesoscale meteorology, land surface and climate science over the last few decades (Skamarock et al. 2019; Scaff et al. 2020). The model provides a comprehensive range of options for selecting physics schemes: microphysics, cumulus, land surface, boundary layer and long/short wave radiation schemes that can be modified and fine-tuned to suit its regional environment (Tian et al. 2017; Skamarock et al. 2019). Advanced Research WRF (ARW) and Nonhydrostatic Mesoscale Model (NMM) are two dynamical cores of the model which have Eulerian mass dynamical cores with terrain-following vertical coordinates (Skamarock et al. 2019). The WRF-ARW is widely used in regional climate and seasonal time-scale research, chemistry applications, global simulations and idealized simulations at many scales (convection, baroclinic waves, extensive eddy simulations, etc.). This model performs better than the fifth-generation Penn State/NCAR (National Center for Atmospheric Research) Mesoscale Model (MM5), while rainfall is still one of the most challenging variables to simulate and predict (Tian et al. 2017; Skamarock et al. 2019). Distinct physical parameterization combinations within the WRF model must be applied in different cases because each physics scheme emphasizes a different physical process and has its unique structure and complexity to storm formation and development (Tian et al. 2017).

Previous studies have proven that the microphysics and cumulus parameterization schemes directly affect precipitation predictions of the WRF model, whereas the precipitation simulation is less susceptible to microphysics in contrast to the cumulus convection scheme (Chawla et al. 2018; Liu et al. 2018). The microphysics scheme includes particular water vapor, cloud and precipitation processes. The cumulus parameterization scheme is responsible for the subgrid-scale effects of convective or shallow clouds and it is theoretically valid only for coarser grid sizes. Activating the cumulus scheme within the high-resolution inner domain improves the precipitation forecasting skill in explicit microphysics alone (Lee et al. 2011; Sun & Barros 2014). Model horizontal resolutions between 2 and 10 km are referred to as the ‘gray zone’, and users must either avoid the grid scale or perform some basic testing to decide whether utilizing a cumulus scheme is ideal for their specific grid scale beyond certain limits (On et al. 2018).

Sri Lanka is an island nation in the Asia-Pacific region prone to severe tropical storms that cause widespread flooding and landslides in the country, destroying dwellings and submerging entire villages. In Sri Lanka, landslides threaten several districts and spread over a land area with the Badulla district being the most vulnerable (Jayasinghe et al. 2017). Rainfall is a significant cause of landslides in Sri Lanka, and sloped agricultural lands in the Badulla district might experience landslides after two days of continuous rainfall of 75–100 mm (Perera et al. 2017). As a result, the Badulu Oya catchment area in the Badulla district is a significant area where large-scale landslides are observed during intense rainfall events that have disastrous effects on infrastructure and result in significant monetary and human losses in a short period.

The numerical studies employing the WRF-ARW model are more rigorous in Sri Lanka due to its complex topographical features, notably its central mountainous area. Previous research utilizing the WRF model to forecast heavy rainfall events in Sri Lanka discovered that the model's results did not accurately estimate rainfall amounts and that rainfall patterns varied depending on the model's configurations (Nandalal et al. 2012; Darshika & Premalal 2015; De-Silva et al. 2016; Rodrigo et al. 2018; Samarasingha et al. 2021). Hence, there is no definitive conclusion regarding physics systems in Sri Lanka, necessitating additional research. Hopefully, this study will address this research vacuum with appropriate content for operational usage and future research into the WRF model's use in rainfall forecasting.

The purpose of this study is to simulate extreme rainfall events that occur during the North-East monsoon over the Badulu Oya basin with an emphasis on evaluating the performance of the physics schemes combinations to predict the spatial distribution of daily accumulated rainfall accurately. In an initial desk study, more accurate precipitation amounts are obtained with activated cumulus parameterization schemes inside the finer grid resolution consistent with previous investigations (Chawla et al. 2018; Rodrigo et al. 2018). The model outputs on model grid points were interpolated and extracted to each observation location to perform an inter-comparison. The deterministic approach was used to examine the accuracy of the amount and position of precipitation using a categorical verification method and two continuous verification methods by comparing the simulated 24-h precipitation output with the observed 24-h precipitation levels.

The structure of this paper is as follows: in Section 2, the methodology including the study area, selected precipitation events, model set-up, data and results comparison methods used for all test events are described. The results are discussed in Section 3 and the conclusions in Section 4.

In December 2019, extensive rainfall events were recorded in the Badulu Oya catchment area with 100–140 mm large cumulative rainfall per day during the North-East monsoon which instigated severe damage for residential people in that area (Ministry of Defence 2019). The dates and recorded maximum daily rainfall amounts within the study area are 1–8 December (130 mm) (First Rainfall Event) and 15–22 December (135 mm) (Second Rainfall Event) (DMSL 2020).

These two recorded extreme rainfall events which caused considerable damage within the 2019 year were selected as model testing events for this study: 1–8 December (First Rainfall Event) and 15–22 December (Second Rainfall Event). To model the validation process, another two extreme precipitation events were selected within the North-East monsoon: 21–28 December 2018 and 17–24 December 2020 (DMSL 2020).

Fifteen physics combinations (CM1–CM15) were configured to evaluate a series of mesoscale simulations with various parameterization schemes for the study. The WRF-ARW uses Yonsei University planetary boundary layer scheme (YSU) (Hong et al. 2006), RRTM longwave radiation scheme (Mlawer et al. 1997), Dudhia shortwave radiation scheme (Dudhia 1989), Unified Noah land surface scheme (Mukul Tewari et al. 2004), Revised MM5 surface layer scheme (Jimenez et al. 2012) as default schemes which were not changed within the selected combinations. The microphysics and cumulus schemes were changed among the physics combinations which are sensitive to the numerical weather prediction over tropical areas (Nandalal et al. 2012; Chawla et al. 2018; Liu et al. 2018). Eta (Ferrier) scheme (Rogers et al. 2001), WRF single–moment 3–class scheme (WSM3) (Han & Hong 2018), WRF single–moment 5–class scheme (WSM5) (Hong et al. 2004), WRF single–moment 6–class scheme (WSM6) (Hong & Lim 2006), Kessler scheme (Kessler 1969), WRF double-moment 6–class scheme (WDM6) (Lim & Hong 2010) were selected as microphysics schemes for the study. Kain–Fritsch scheme (KF) (Kain 2004), Betts–Miller–Janjic scheme (BMJ) (Janjic 1994), Multi–scale Kain–Fritsch scheme (MKF) (Zheng et al. 2016), Grell–Freitas Ensemble scheme (GF) (Grell & Freitas 2014) were used as cumulus schemes for the tested combinations. Table 1 provides a list of the tested physics scheme combinations that interchange the selected microphysics scheme and cumulus schemes with the default remaining schemes. Moreover, the same microphysics scheme and cumulus scheme were inserted into all the nested domains for a given model simulation in an attempt to minimize inconsistencies at the interface of the computation grids (Warner et al. 1997).

Selected heavy rainfall events were simulated in WRF-ARW software for an 8-day period where the first day (24 h) was considered as the spin-up period. The models were initialized at 0600 UTC and rainfall periods of 24 h ending at 0600 UTC on each day within 7 days were investigated during this study. This is the same as the 24-h rainfall criterion used in the Department of Meteorology, Sri Lanka (DMSL) while considering slight deviations of the actual recording time (0300 UTC) of field precipitation measurements. Boundary conditions for the finer domain were updated every 6 h. It should be noted that each model combination used identical configurations and initialization data, other than marginally different parameterization schemes. All the proposed combinations (CM1–CM15) were run for the selected two rainfall periods separately within December 2019. For verification and validation of the selected final five best models, these models were tested for the two North-East monsoon rainfall events: 21–28 December 2018 and 17–24 December 2020.

NCEP FNL (Final) high-resolution analyses data were used as initial and lateral boundary conditions to the atmospheric model which are on 1-degree by 1-degree grids available at 6 hourly intervals (National Center for Atmospheric Research 2019). The initial and lateral boundary conditions were forced onto the outer domain and the model was left to solve its nested domain numerically using the forcing from the mother domain. Required terrestrial datasets were applied from the WRF users’ website which is the most commonly used dataset for the WRF model simulations (National Center for Atmospheric Research 2020). Surface observations of DMSL stations were utilized for the investigated periods of the events with the intention of assessing the model performance. Five rain gauging stations within the Badulu Oya catchment including Kandaketiya (7.167 N, 81.017 E), Ledgerwatte Estate (7.033 N, 81.017 E), Badulla (6.983 N, 81.050 E), Telbedda Estate (6.978 N, 81.084 E) and Lower Spring Valley (6.917 N, 81.100 E) were selected for the study (Figure 1) considering their data availability. Among these stations, only Kandaketiya and Badulla stations are auto record type gauges (3 h period), while the rest are manually operated rain gauges. According to the Department of Meteorology, all rainfall observations are for the cumulative rainfall at 0300 UTC (08.30 a.m. in Sri Lanka time). However, the actual recording time of field precipitation measurements has few deviations. The manually operated rain gauges have an accuracy of 0.1 mm. These surface observations were used as ground-truth observations in this study.

The quality of rainfall predictions was compared in this study based on the deterministic approach considering the categorical verification method, as well as continuous verification methods after the WRF model grid point values were interpolated and extracted to each observation location.

The joint distribution of forecast and observation can be described using the dichotomous categorical verification technique (contingency tables). A 2 × 2 contingency table based on forecast categories that were established on 24-h total precipitation accumulated thresholds of above 10 mm, above 25 mm and above 50 mm was used to evaluate a total of 35 (5 × 7) data pairs that were relevant to the 5 rainfall gauging stations for a period of 7 days. The observed daily rainfall data and the simulated data were evaluated to determine the proportion correct (PC), probability of detection (POD), false alarm ratio (FAR), bias (B) and Equitable Threat Score (ETS) for each threshold category across the entire domain (Wilks 2011). The PC determines the fraction of all correct forecasts, which is sensitive to both hits and correct rejections. The value ranges between 0 and 1, and the perfect score is 1. The POD displays the percentage of stations at which WRF correctly simulated precipitation when precipitation was observed, which is also called the hit rate. The value ranges between 0 and 1 and the perfect score is 1. The FAR determines the percentage of falsely simulated events of precipitation when compared to the total number of simulated precipitation events. The value varies between 0 and 1, and the perfect score is 0. B determines whether WRF simulated precipitation at more stations than observed, or fewer than observed. The value can range between 0 and infinity, and the perfect score is 1. The ETS or Gilbert Skill Score determines overall skill when simulating precipitation. It includes a correction term that reduces the effect of a correct precipitation simulation by chance. An ETS score of 0 indicates the same accuracy as a random precipitation simulation and positive ETS scores indicate some level of accuracy. In contrast, a perfect precipitation simulation would have an ETS score equal to 1. To summarize all results related to the three rainfall thresholds, average values of the ETS and B were employed to represent the overall performance of the model rainfall predictions. Appendix A contains formulae for the above precipitation statistics indices. The study primarily focused on the higher thresholds of precipitation due to the hazardous consequences of high precipitation rates.

Precipitation forecasts can be compared by direct comparison with observed and calculation of various statistical parameters which include the measures of distributional properties (mean, variance), forecast confidence (sample size), correlation statistics (scatter plots, correlation coefficient) and accuracy statistics (Mean Absolute Error, Root Mean Square Error (RMSE)). The most common option for illustrating the phase association (correlation) and amplitude (standard deviation) errors in the field of weather research forecasting is the normalized version of Taylor diagrams (NCAR 2019). After obtaining the forecasted 24-h accumulated rainfall at the five rainfall stations from the WRF output data, a total of 35 (5 × 7) data pairs pertinent to the 5 rainfall gauging stations for a period of 7 days were used to calculate the statistics mentioned above indices which are listed in Appendix B. These indices were presented graphically using the Taylor diagram after reviewing technical details at the NCAR website (NCAR 2019). In the diagram, WRF models that agree well with observations lie nearest to the point marked ‘REF’ on the x-axis. These WRF models have relatively high correlation and low normalized RMSE. WRF models lying on the reference have the correct Normalized Standard Deviation (which indicates that the model variations are of the correct amplitude). Furthermore, lower bias values indicate that the model's outputs do not significantly deviate from the observations.

The WRF output data and observed rainfall data on the day of the highest daily rainfall already recorded were spatially distributed within 3 km × 3 km horizontal grids (a total of 44 grids) over the Badulu Oya catchment and compared using the Inverse Distance Weighting (IDW) technique of the ArcGIS software. The discrepancy between the daily rainfall prediction and the observed rainfall was illustrated using Mean Absolute Model Error percentage (MAME, Appendix C) (Samarasingha et al. 2021). The tested models’ overall skill was compared using the cumulative correctly predicted area (CCPA) percentage within the −50 to 50% range of MAME over the basin.

The computed statistical indices in each accumulated precipitation threshold for the combinations (CM1-CM15) are displayed in Table 2. Overall, the highest PC values are in the largest threshold; 50 mm. CM9 displays the highest PC values in the 10 mm (PC value is 0.74) and 25 mm (PC value is 0.83) thresholds, while CM3 and CM11 record the highest PC value in the 50 mm (PC value is 0.94) threshold. CM12 has the highest POD value in the 10 mm (POD value is 0.95) and 25 mm (POD value is 1.00) thresholds. In the largest threshold, CM8, CM10 and CM12 have the best POD value which is 1.00. The lowest FAR value in the 10 mm (FAR value is 0.08) threshold belongs to CM5, while CM7 has the lowest FAR value in the 25 mm (FAR value is 0.00) threshold. CM3 and CM11 both record the lowest FAR value in the 50 mm threshold (FAR value is 0.17). The FAR value relevant to the CM7 in the 50 mm threshold could not be calculated due to both hits and false alarms being zero in the largest threshold category. CM1 and CM2 display the lowest skill in each category with Ferrier and WSM3 microphysics, respectively. Considering the overall statistical indices, CM9 and CM12 which include the WSM6 microphysics scheme display positive forecast skills among others in the 10 and 25 mm threshold categories. The highest threshold category has the best overall results from CM3 and CM11, including the KF and MKF cumulus schemes. The overall performance of the first rainfall event analysis was identified using the average values of the ETS and B related to the three threshold categories (Table 2). CM3 displays the best performance with the largest ETS value (ETS value is 0.38) and bias value of 0.95. Next, CM4 records positive skill similar to CM3 which has an ETS value of 0.36 and a bias value of 1.04. After the CM3 and CM4 combinations, CM11, CM9 and CM8 record good skill among the remaining combinations which have ETS values of 0.31, 0.31, 0.34 and bias values of 1.03, 1.45, 1.50, respectively. CM1 and CM2 are less-performed models which also agree with the above threshold categories’ findings. These two combinations have Ferrier and WSM3 microphysics, respectively. Considering the best-performing two models, WSM5 and WSM6 microphysics schemes with KF cumulus scheme simulated the first rainfall event more accurately than the other schemes.

The statistical indices in each accumulated precipitation threshold category for the second rainfall event are indicated in Table 3. CM12 displays the highest PC value in the 10 mm (PC value is 0.80) threshold category, while CM10 records the highest PC value in the 25 mm (PC value is 0.69) and 50 mm (PC value is 0.77) threshold categories. CM12 has the highest POD value in the 10 mm (POD value is 0.91) and 25 mm (POD value is 0.63) thresholds, while CM14 has the highest POD value in the 50 mm (POD value is 0.46) threshold category. CM10 records the lowest FAR value in all three threshold categories (FAR value is 0.0), while CM8 also has zero FAR value in the 25 and 50 mm threshold categories. The FAR value relevant to the CM5, CM6 and CM9 in the 25 mm threshold and CM3, CM4, CM5, CM6 and CM9 in the 50 mm threshold could not be calculated due to hits and false alarms being zero during the calculation steps. CM1 and CM2 display the lowest skill in each category with Ferrier, WSM3 microphysics scheme and KF cumulus scheme, respectively. Considering the overall statistical indices, CM9 and CM5 display positive forecast skill among others in the 10 and 25 mm threshold categories, including WSM6 and Kessler micro physic schemes with BMJ and KF cumulus schemes, respectively.

CM3 and CM11 with KF and MKF cumulus schemes have the overall best performance in the largest threshold category. Considering the obtained ETS values, CM10 indicates the best performance with an ETS value of 0.19 and bias value of 0.30. Although CM12 has a higher ETS value (ETS value is 0.18), it exhibits a considerably higher bias value (bias value is 0.96) than CM10. Next, CM2, CM7 and CM11 have good skill among the remaining combinations which have ETS values of 0.13, 0.13, 0.15 and bias values of 0.39, 0.46, 0.54, respectively. Considering the calculated ETS values in the second rainfall event, the first rainfall event has higher ETS values than the second event. This resulted in a higher predicting skill of the models during the first rainfall event than the second rainfall event. However, the first rainfall has higher bias values which exceed 1.0 (over forecasting) than the second rainfall event (under forecasting). CM9 and CM6 are less performing combinations with Ferrier and WSM6 microphysics with BMJ cumulus scheme, respectively. Considering the two best-performing configurations, WSM6 and Kessler microphysics schemes with MKF and BMJ cumulus schemes could simulate the second rainfall event more accurately than the other schemes.

The simulation of rain events over complex orography is still challenging, as WRF is influenced by the albedo, wind speed and direction (Lorente-Plazas et al. 2016). The surface wind was calm and varied throughout the first rainfall event, while it changed to a north-easterly direction during the second rainfall event with wind speeds ranging from 5 to 20 kmph over the country (DMSL 2020). The differences in wind characteristics and tropical weather features of the two tested rainfall events may have influenced the first rainfall event in this study to have better statistical parameters than the second rainfall event. Besides, the physics sensitivity of the parameterization schemes might need further optimization and fine tuning based on many more event tests. This could be an extended work of this study.

The model training was conducted for the selected two rainfall events in December 2019 and predicting skills of the tested models were evaluated based on the deterministic approach considering a categorical verification method and continuous verification methods in the above sections. Table 4 has a summary of the findings of the conducted evaluation process. The first four best-performing model combinations identified in each verification method relevant to the tested two rainfall events are presented in Table 4. The selected combinations are listed in descending order considering their identified performance skills in each verification method.

Considering the obtained summary of model training results, CM3, CM4 and CM11 display more positive skills than the remaining tested model combinations for the first rainfall event simulation considering the outcomes of the three deterministic approaches. In the second rainfall event simulation, CM10, CM11 and CM12 performed better than other models based on the three deterministic approaches. Among these mentioned best model combinations, CM11 performs well in first rainfall and second rainfall simulations. Finally, CM3, CM4, CM10, CM11 and CM12 model combinations were filtered to achieve an improved station-wise forecast within the North-East monsoon season over the Badulu Oya catchment.

The selected best-performing five models were WSM5, WSM6, Kessler and WDM6 microphysics schemes, and KF, BMJ and MKF cumulus schemes. The WSM5 is an enhanced microphysics scheme developed from the simple WSM3 scheme. The mixed phase WSM5 includes prognostic water substance variables of water vapor, cloud, ice, rain and snow (Hong & Lim 2006). Hong & Lim (2006) stated that increasing the number of hydrometers increases the rainfall amount for the finer grid resolution. Moreover, their study specified that sedimentation of cloud ice is an essential factor for better simulations of monsoonal precipitation and large-scale features within the East Asian summer monsoon (Hong & Lim 2006). The WSM6 includes graupel as another predictive variable added to WSM5. These schemes better represent the ice cloud properties considering the cloud radiation feedback. The ice cloud affects longwave heating more significantly than shortwave heating in terms of upper-level mean temperature (Hong & Lim 2006). The Kessler microphysics scheme is the oldest microphysics scheme in the WRF model. This liquid-only scheme includes water vapor, cloud water and rain. The production, fall and evaporation of rain, the accretion, auto-conversion of cloud water and the production of cloud water from condensation are the main processes within this scheme (Mielikainen et al. 2013). The Kessler scheme uses a lower threshold to convert cloud water to rain which occurs early automatically, increases rainfall rates and shifts precipitation upstream (Mielikainen et al. 2013). The WDM6 microphysics scheme is a double-moment warm-rain microphysics scheme that predicts the mixing ratio of the hydrometeors and their number concentration of cloud and rainwater. These have been added to the WSM6 scheme. The WDM6 scheme can allow flexibility in variable raindrop size distribution by predicting the number of concentrations of clouds and rain (Lim & Hong 2010).

The KF cumulus scheme is a mass flux parameterization scheme that uses the Lagrangian parcel method including vertical momentum dynamics. It can be generally grouped into three parts: (1) convective trigger function, (2) mass flux formulation and (3) closure assumptions which can represent the subgrid scale feature of the updraft and downdraft process that lead to the development of strong convection (Kain 2004). Notably, experiments with KF cumulus parameterization often perform better with WSM5 and WSM6 schemes. This statement is in line with the results of Pennelly et al. (2014), who examined three flooding events in Alberta, Canada, using five different cumulus parameterization schemes (Pennelly et al. 2014). However, the KF scheme does not perform much better when the spatial distribution of the catchment is investigated. This concurs with Gallus Jr.'s (1999) findings which showed that more precisely simulating excessive rainfall happening at high resolutions with the ascent facilitates activation of the scheme and shifts the maximum by several hundred kilometers (Gallus 1999). The BMJ cumulus scheme is a convective adjustment scheme that includes deep and shallow convection. The deep convection uses empirically based quasi-equilibrium thermodynamic profiles as a reference state rather than a moist adiabat. The construction of the reference profiles and the specification of the relaxation time scale are two significant components of this scheme (Vaidya & Singh 2000). Earlier research revealed that the BMJ scheme is more consistent with land-based data during the spatial analysis when compared to the other cumulus schemes (Ratna et al. 2014; Remesan et al. 2015). Additionally, Vaidya & Singh (2000) found that the spurious rainfall spread in the Arabian Sea and in the Bay of Bengal is suppressed in the BMJ scheme (Vaidya & Singh 2000). Despite the BMJ scheme not explicitly defining subgrid updrafts and downdrafts in addition to the mesoscale microphysical processes, it performs better in this assessment with the Kessler microphysics scheme (CM10). The MKF cumulus scheme is an updated version of the KF scheme which includes subgrid-scale cloud radiation interactions, a dynamic adjustment time scale, the impact of subgrid-scale cloud updraft mass fluxes on grid-scale vertical velocity and an entrainment methodology based on the lifting condensation level. These changes introduce scale dependency for some of these critical parameters in the KF scheme, and are expected to improve weather forecasts at 9 and 3 km grid spacings (Zheng et al. 2016).

In the evaluated two extreme rainfall events, the WSM5-KF combination (CM3) records excellent statistical indices among other combinations. Although the WSM6 scheme generally outperformed the WSM5 scheme (Kar & Tiwari 2016), the WSM6 scheme performed less well in this evaluation even when a graupel particle was introduced to make it more complicated than the WSM5 with the KF scheme. This phenomenon might be caused by a slowdown in the direct conversion of snow and cloud particles into rainwater (Akinola & Yin 2019). The KF scheme outperformed the BMJ system overall during the tested first rainfall event. This statement is consistent with a study by Bukovsky & Karoly (2009) that examined several WRF modeling systems including two land surface models and two cumulus schemes over the United States (Bukovsky & Karoly 2009).

The location and topography greatly influence the effectiveness of short-term rainfall forecasting and the WRF can capture the intricate interactions between land and lower atmospherics in generally flat terrains but tends to over- or under-simulate in mountainous terrains (Maussion et al. 2011). In our case, a similar observation was made regarding the deviations in the rainfall magnitudes estimated by the model combinations. This also influences the selected five model combinations that outperformed others in the mountainous Badulu Oya catchment. Yet again, the model and physics uncertainties could be reasons for such performance.

In the above physics selection process, findings of the categorical verification and Taylor diagram were considered more than the analysis of spatial distribution. The reason for this decision is that the conducted analysis approach using the IDW technique of the ArcGIS software does not consider the inhomogeneity in the topography or circulation over the catchment. The filtered best five model combinations were evaluated in verification simulations.

The obtained average values of ETS and B over the three threshold categories are indicated in Table 5 relevant to the December 2018 rainfall event (from 21st to 28th of December 2018) and December 2020 rainfall event (from 17th to 24th of December 2020). During the analyzed December 2018 rainfall event, CM12 had the highest ETS value (ETS value is 0.20). Thereafter, ETS values of CM4 and CM3 are 0.17 and 0.13 in descending order, respectively. CM10 and CM11 are less-performed models which have zero ETS value. The highest bias value (bias value is 1.61) belongs to the CM3 model. Only the CM10 model has zero value for both hits and false alarms within the three threshold categories which produced zero value of average ETS and bias. Considering the tested December 2020 rainfall event, CM10 has the highest ETS value (ETS value is 0.18) and the lowest bias value (bias value is 0.26). After CM10, CM11 has a higher ETS value (ETS value is 0.11), while the remaining tested combinations have lower ETS values of less than 0.10. CM12 and CM4 have higher bias values among the tested models (bias values are 1.46 and 1.42).

The model outputs showed acceptable predictability for all the tested cases. In contrast, model predictions displayed high sensitivity for the physical parameterization and the atmospheric features during the tested seasons because each model's physical parameterization has unique methods for resolving weather features at different scales. The corresponding results of the heavy rain cases showed seasonal and physical parameterization effects. Considering obtained overall results of the three verification methods, CM12, CM4 and CM3 models can be selected as the best-performed models in the December 2018 rainfall event, while CM10, CM11 and CM3 models are the best-performed models in the December 2020 rainfall event. Hence, CM3, CM4, CM10, CM11 and CM12 models have noticeable forecasting skills relevant to the North-East monsoon season over the Badulu Oya catchment.

Overall, the sensitivities of the model physics to different spatial resolutions and areas of interest might need more scrutiny. Nevertheless, this study remains essential as this has helped to highlight and shortlist the good/poor combinations of physics options which can be studied in detail under more numerical simulations and examinations as necessary. Such an exercise would lead to better optimizing the WRF model over this tropical island of Sri Lanka.

The simulations of heavy rainfall events relevant to the North-East monsoon season performed in this research over the Badulu Oya catchment with different parameterization schemes demonstrated significant sensitivity for microphysics and cumulus schemes.

CM3, CM4, CM10, CM11 and CM12 model configurations with the WSM5, WSM6, Kessler and WDM6 microphysics schemes along with the KF, BMJ and MKF cumulus parameterization schemes displayed the overall best performances during the North-East monsoon season.

Among these models, the CM3 model configuration including the WSM5 microphysics scheme, KF cumulus scheme along with the remaining default parameterization schemes of WRF-ARW records the best statistics; ETS value of 0.38, bias value of 0.95, correlation coefficient of 0.76, NRMSE of 0.72, normalized standard deviation of 1.06 and CCPA value of 75%.

The output of the five developed model combinations displayed tremendous sensitivity with seasonal and physical parameterization effects during the model validation due to their unique methods for resolving weather features at different scales. The ETS value of 0.20, correlation coefficient value of 0.35 and CCPA value of 88.6% are the highest statistical indices among the two rainfall events assessed during the model validation.

Although the model outcomes indicated large deviations in the rainfall magnitudes estimated by the model combinations, several combinations (CM3, CM4, CM10, CM11 and CM12) could capture spatial distribution similar to the spatial distribution obtained from the ground observation network.

To conclude, the filtered five model combinations can be used to predict the future significant rainfall events over the Badulu Oya catchment within the North-East monsoon season and those data can be used to deliver early warnings for the landslides within the study area. Moreover, several studies can also use the five model combinations to predict soil moisture content, wind characteristics and radiation characteristics in the area.

The results of this study imply that more precise predictions can be produced by taking into account other parameterization schemes, including radiation, boundary layer and land surface physics. Increasing the number of simulations of extreme rainfall events within different seasons to find better combinations of physics for heavy rainfall forecasting for all seasons would be valuable. Additional sensitivity testing may be desired to determine the sensitivity of global model data such as NOAA and ECMWF due to the WRF-ARW model predictions being highly dependent on the lateral boundary and initial conditions provided for initialization.

P. G. S. G. conceptualized the whole work, developed the methodology, performed simulations using the software, validated the model results, conducted a formal analysis, investigation, data curation, wrote the original draft, visualized the article, wrote the review and edited the article. P. N. conceptualized the whole work, brought resources, wrote the review, edited the article and supervised the work. R. A. A. supervised the work. S. V. R. wrote the review, edited the article and supervised the work.

The authors acknowledge Computer Department of University of Peradeniya for giving accessibility to their high-performance computing (HPC) for conducting research simulations. Moreover, the authors are grateful to Department of Meteorology, Sri Lanka (DMSL) for providing their recorded precipitation data throughout the study.

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.