Assessment of the Weather Research and Forecasting (WRF) Model for Extreme Rainfall Event Simulations in the Upper Ganga Basin

Reliable estimates of extreme rainfall events are necessary for an accurate prediction of floods. Most of the global rainfall products are available at a coarse resolution, rendering them less desirable for extreme rainfall analysis. Therefore, regional mesoscale models such as the Advanced Research version of the Weather Research and Forecasting (WRF-ARW) model, are often used to provide rainfall estimates at fine grid spacing. Modelling heavy rainfall events is an enduring challenge, as such events depend on multiscale interactions, and the model 15 configurations such as grid spacing, physical parameterization and initialization. With this background, the WRFARW model is implemented in this study to investigate the impact of different processes on extreme rainfall simulation, by considering a representative event that occurred during 15 – 18 June 2013 over the Ganges basin in India, which is located at the foothills of the Himalayas. This event is simulated with ensembles involving four different microphysics (MP), two cumulus (CU) parameterizations, two planetary boundary layer (PBL), and two 20 land surface physics options; and different resolutions (grid spacing) within the WRF model. The simulated rainfall is evaluated against the observations from 18 rain gauges and the Tropical Rainfall Measuring Mission MultiSatellite Precipitation Analysis (TMPA) 3B42RT version 7 data. From the analysis, it is noted that the selection of MP scheme influences the spatial pattern of rainfall, while the choice of PBL and CU parameterizations influence the magnitude of rainfall in the model simulations. Further, WRF run with Goddard MP, Mellor–Yamada–Janjic 25 PBL and Betts–Miller–Janjic ́ CU scheme is found to perform ‘best’ in simulating this heavy rain event. The model performance improved through incorporation of detailed land surface processes involving prognostic soil moisture evolution in Noah scheme as compared to the simple Slab model. To analyze the effect of model grid spacing, two sets of downscaling ratios – (i) 1:3, Global to Regional (G2R) scale; and (ii) 1:9, Global to Convection-permitting scale (G2C) are employed. Results indicate that higher downscaling ratio (G2C) causes higher variability and 30 consequently, large errors in the simulations. Therefore, G2R is opted as a suitable choice for simulating heavy rainfall event in the present case study. Further, the WRF simulated rainfall is found to exhibit least bias when compared with that of the Coordinated Regional Climate Downscaling Experiment (CORDEX) data and the NCEP FiNaL (FNL) reanalysis data. Hydrol. Earth Syst. Sci. Discuss., https://doi.org/10.5194/hess-2017-533 Manuscript under review for journal Hydrol. Earth Syst. Sci. Discussion started: 30 August 2017 c © Author(s) 2017. CC BY 4.0 License.


Introduction 35
Indian Summer Monsson Rainfall (ISMR) is often associated with very heavy (124.5 to 244.4 mm/day) to extremely heavy (more than 244.5 mm/day) rainfall (Ray et al., 2014), particularly during June to September (Srinivas et al., 2013).The extremely heavy rainfall events often occur due to the presence of organized meso-convective systems (MCSs) embedded in large scale monsoonal features such as offshore troughs and vortices, depressions over the Bay of Bengal /Arabian Sea, and mid tropospheric cyclones (Sikka and Gadgil, 1980;Benson Jr and Rao, 1987;Zipser et 40 al., 2006).
Extremely heavy rainfall at shorter time scales are particularly difficult to predict in mountainous terrains, and continue to be a challenge to operational and research community (Das et al., 2008;Li et al., 2017).Because of the multiscale features associated with these events, there has been an ongoing effort to implement mesoscale Numerical Weather Prediction (NWP) models for the ISMR simulations.The operational and research community 45 has widely adopted the Advanced Research version of the Weather Research and Forecasting (WRF-ARW) model (hereafter referred as the WRF model), to simulate a variety of high impact meteorological events, such as rainfall (Vaidya and Kulkarni, 2007;Deb et al., 2008;Kumar et al., 2008;Chang et al., 2009;Routray et al., 2010;Mohanty et al., 2012), tropical cyclones (Raju et al., 2011;Routray et al., 2016;Osuri et al., 2017b) and thunderstorms (Madala et al., 2014;Osuri et al., 2017a).However, setting up the WRF model, that simulates extremely heavy 50 rainfall over the ISMR region is still considered as a challenging task, which involves consideration of several aspects such as forcing data, model grid spacing/resolution, land surface parameterization and choice of an appropriate physics scheme.
Earlier studies such as by Krishnamurthy et al., (2009); Misenis and Zhang (2010); Rauscher et al., (2010);Mohanty et al., (2012); Chevuturi et al., (2015) indicated that heavy rainfall predictions can be improved through ensemble 55 model techniques and fine grid resolution.However, the influence of the model parameterization schemes on mesoscale rainfall simulations over India is still an understudied issue.In particular, heavy rain simulations studies have reviewed the impact of individual parameterization options such as the Microphysics (MP) scheme (Rajeevan et al., 2010;Raju et al., 2011;Kumar et al., 2012), Cumulus (CU) parameterization scheme (Deb et al., 2008;Mukhopadhyay et al., 2010;Srinivas et al., 2013;Madala et al., 2014), Planetary Boundary Layer (PBL) scheme (Li 60 and Pu, 2008;Hu et al., 2010;Hariprasad et al., 2014), and Land Surface Model (LSMs) options (Chang et al., 2009).However, the ensemble analysis that reviews the relative impact of different configurations is lacking.It is important to study the impact of different parameterizations in an ensemble mode because it is often likely that the performance of one scheme depends on other model configurations considered.For example, the conclusions regarding which CU scheme performs best would be intimately tied to the choice of the MP or land surface options 65 considered in conducting the numerical experiments.With this perspective, an attempt is made in this paper, to assess sensitivity of the WRF model to predict an extremely heavy rainfall episode, that occurred from 15 June through 18 June 2013, over the Ganges basin in India, which is located at the foothills of the Himalayas.Additional details regarding the event are presented ahead.The performance of the WRF model is evaluated against the global reanalysis and downscaled data.70

Description of the heavy rainfall event
The 2013 summer monsoon had a normal onset but the trough advanced rapidly, covering whole of India by mid-June, instead of mid-July (Ray et al., 2014).This large-scale setting is thought to have created a platform for 80 interaction of two synoptic scale events -northwest moving depression from Bay of Bengal and preexisting westerly trough in mid-troposphere.Meteorological studies conducted over the region (Kotal et al., 2014;Ray et al., 2014;Rajesh et al., 2016) established that there was a monsoon low pressure system during this period.The longitudinal time section for 850 hPa geopotential height along with anomaly averaged over 20° -26° N showed high negative anomaly on 14 June, which migrated to west, over 75° E by 17 June.The meridional wind anomaly 85 within the belt of 35° -45° N showed a westerly wave, moving from 10° E on 12 June to 70° E on 17 June.These two anomalies are found to be in phase, consequently causing interaction between the eastward moving trough in mid-upper troposphere and westward moving monsoon low in the lower troposphere.The monsoon low provided the moisture feed and the upper level westerly trough provided the divergence to lift the moisture.This whole system eventually led to an unanticipated heavy rainfall during 15 -18 June 2013 in the Kedarnath valley and 90 adjoining areas in the state of Uttarakhand, India (Kotal et al., 2014;Ray et al., 2014;Rajesh et al., 2016).The region received rainfall greater than 370 mm in one day (17 -18 June 2013), which is 375% above the daily normal rainfall (65.9 mm) during the monsoon season (Ray et al., 2014).Consequently, heavy floods occurred in the region, causing unprecedented damage to life and property.
The synthesis of the synoptic setting of the event has been carried out in a number of studies (Dube et al., 2014;95 Kotal et al., 2014;Ray et al., 2014;Shekhar et al., 2015;Rajesh et al., 2016), but the mesoscale assessment pertaining to the simulation of this rainfall event is still lacking.Therefore, the present study emphasizes on quantitatively evaluating and conducting sensitivity analysis on the prediction capability of the WRF model on rainfall simulations.
Since the epicenter of the heavy rainfall was Kedarnath, study region comprising of the upstream part of the Ganga 100 Basin in India, referred as Upper Ganga Basin (UGB) hereafter, is selected in this paper.Figure 1 presents the topography of the UGB as described for the three Domains of the WRF (Domain 1, Domain 2a and Domain 2b with 27, 9 and 3 km grid resolution respectively) along with the 18 rain gauge stations located within the region.The region is of social, cultural and economic importance to India, further making this study necessary.It is noticed that the northwest part of the region received higher rainfall compared to the northeast, with stations such as Tehri and Dehradun showing 327% and 210% (respectively) more rainfall than their historic means.A few stations like Chamoli in the northeast region, received 250 mm of cumulative rainfall over the 4 days' period, which is 144% higher than the historic mean.In general, most of the stations in the southern part of the basin, which are 120 located at lower elevation, recorded relatively less rainfall with a cumulative range of 445 mm, in comparison to the northern part (at higher altitude) having a rainfall range of 515 mm.Additionally, three stations in the southeast region, i.e., Mukteshwar, Haldwani and Nainital received extremely heavy rainfall with a cumulative average of 498 mm.From the above analysis, it is evident that the system moved from east to west direction with two distinct regions in the UGB -southeast and northwest, receiving extremely heavy rainfall.125

Observed Data 110
The region has complex topography and a limited number of rain gauges because of the difficulty in operating a network in this region.To further capture the spatial variability in rainfall, Tropical Rainfall Measuring Mission Multi-Satellite Precipitation Analysis (TMPA) 3B42RT (version 7) product, which is available at 0.25° resolution at daily scale is analyzed (Fig. 3).It is to be noted that, since the focus area for the analyses is the upstream region of the UGB (Fig. 1 (iii) and Fig. 2), results are presented in this paper with respect to the geographical extent of 130 Domain 2b throughout this paper.From Figure 3, it can be noted that the TMPA data is able to capture the spatial variability in the rainfall -with distinct clusters corresponding to heavy rainfall in the northwest and southeast regions of the study area.However, the rainfall amount is significantly underestimated by the TMPA product, with the maximum value of 265 mm against the recorded 650 mm.This under reporting for gridded satellite product versus rain gauge in the ISMR region is a well-known feature (Rahman et al., 2009;Kneis et al., 2014).The TMPA 135 estimates are verified against the IMD station observations for baseline quality check.Mean absolute error (MAE), root mean square error (RMSE) and bias (β) are computed using the nearest neighborhood mapping approach and are presented in Table 1.TMPA data is observed to behave differently for different ranges of rainfall values over the study domain.TMPA overestimated the rainfall at stations with cumulative rainfall less than 200 mm (e.g., Pauri and Almora).In contrast, 145 rainfall at stations receiving more than 250 mm of cumulative rainfall is underestimated.Stations that received rainfall of 200 -250 mm are well represented in the TMPA data (e.g., Mussorie, Pithogarh and Champawat).From the analysis, it could be inferred that in the TMPA data rainfall values are clustered towards the mean value.Errors noticed in the TMPA data could be attributed to two factors: first, large spatial coverage and coarse resolution of the TMPA data, and second, for comparison with the observed data a simple approach of selecting nearest grid point is 150 implemented.

CORDEX Data
The Coordinated Regional Climate Downscaling Experiment (CORDEX; http://www.cordex.org/)aims at providing high-resolution climate projections for historic and future time periods.To achieve this, climate scenarios from the Atmosphere-Ocean coupled General Circulation Models (AOGCMs) under Coupled Model Intercomparison Project 155 Phase 5 (CMIP5) are dynamically downscaled.For downscaling limited domain, relatively finer resolution Regional Climate Models (RCMs) are used with lateral boundary conditions from the coarse resolution AOGCMs.RCMs resolve the topographic details and land surface heterogeneity in order to obtain climate variables at finer spatial scales, in contrast to the driving AOGCMs.CORDEX provides data for historic simulations from 1971 to 2005 and future projections from 2006 to 2099/2100.For the present paper, CORDEX data corresponding to dynamically 160 downscaled projections from six CMIP5 AOGCMs at 0.5° (~ 50 km) resolution (Table 2), for two Representative Concentration Pathways (RCPs), RCP 4.5 and RCP 8.5 are procured from Centre for Climate Change Research (CCCR), India (http://cccr.tropmet.res.in/home/cordexsa.jsp).The CORDEX projections significantly underestimate the heavy rainfall event.Further, negligible variability is observed across the models for all the days except 18 June, indicating that none of the models is able to capture the realistic magnitude of the event.A few models (NorESM1-M with RCP 4.5 and ACCESS1.0 with RCP 8.5) 175 simulated rainfall close to the observed value of 18 mm for 15 June.Cumulative rainfall across all the models has relatively higher variability, primarily due to the variability in heavy rainfall on 18 June.From the analysis, it can be concluded that the CORDEX data is capable of estimating the qualitative features of the rainfall and has significant under-prediction, indicating that care must be exercised while using the data for applications involving heavy rainfall events, such as flood modelling.180

Model Configuration and Experimental Setup
The simulation experiments in this paper are conducted using the Advanced Research Weather Research and Forecasting (ARW-WRF, or simply WRF) model, version 3.8.WRF is a widely used, NWP non-hydrostatic, mesoscale model, available with several advanced physics and numerical schemes, designed for better prediction of atmospheric processes.The model description and updates can be found from (Skamarock et al., 2005) and the WRF 185 user webpage (http://www2.mmm.ucar.edu/wrf/users/).
The parent domain provides lateral boundary conditions to the inner domains, resulting in the downscaling ratios for simulations as 1:3 and 1:9.The three domains use 30 vertical pressure levels, with the top fixed at 50 hPa.Model 205 time steps were function of grid spacing: 135 s, 45 s and 15 s respectively for the three domains.
The model configuration used default parameterization options following (Osuri et al., 2012).For example, shortwave radiation is based on Dudhia (Dudhia, 1989), and long wave based on Rapid Radiative Transfer Model (RRTM; (Mlawer et al., 1997)) scheme.Other physical parameterization options such as Microphysics (MP), Cumulus (CU) parameterization schemes, Planetary Boundary Layer (PBL) and Land Surface Models (LSMs) were 210 selected as outlined ahead.There is currently no known unique configuration that can best simulate an extremely heavy rainfall event.Therefore, based on literature (e.g.Kumar et al., 2008;Hong and Lee, 2009;Misenis and Zhang, 2010;Mukhopadhyay et al., 2010;Argüeso et al., 2011;Cardoso et al., 2013;Efstathiou et al., 2013), four MP schemes, two CU schemes, two PBL schemes and two LSMs are considered representative to obtain an ensemble of rainfall simulations.The two PBL schemes considered are the Yonsei University (YSU) scheme (Hong 215 et al., 2006) and the Mellor-Yamada-Janjic (MYJ) scheme (Janić, 2001).YSU is a non-local scheme, wherein fluxes are calculated at certain height in the PBL considering the profile of the entire domain.MYJ scheme, on the other hand, is a local scheme in which fluxes are calculated at various heights within the PBL and are related to vertical gradient in the atmospheric variables at the same height.Further details regarding the difference between the YSU and the BMJ schemes can be obtained from (Misenis and Zhang, 2010;Efstathiou et al., 2013).The two CU 220 schemes considered are the Kain -Fritsch (KF) scheme (Kain, 2004) and the Betts-Miller-Janjic (BMJ) scheme (Janjić, 1994(Janjić, , 2000)).KF is a shallow convection scheme, based on entrainment and detrainment plume model with updrafts and downdrafts of mass flux.Potential energy is removed in the convective time scale within this scheme.
Furthermore, it includes cloud, rain, snow and ice detrainment at cloud top.BMJ, on the other hand, considers convection at both shallow and deep levels.However, there is no updraft and downdraft of mass flux and no cloud 225 detrainment.Domain 2b is configured without any CU scheme, assuming MP to explicitly solve the convection at finer resolution (Sikder and Hossain, 2016).Four MP schemes considered are, the Purdue Lin (PLin) scheme (Lin et al., 1983;Chen and Sun, 2002), the Eta Ferrier (Eta) scheme (NOAA, 2001), the WRF Single-Moment 6-class (WSM6) scheme (Hong and Lim, 2006) and the Goddard scheme (Tao et al., 1989).Although both, the PLin scheme and the WSM 6 scheme are based on the parameterization from Rutledge and Hobbs (1984) Hong et al., (2009).The Eta scheme was designed primarily for computational efficiency in NWP models, wherein the total condensate and the water vapors are directly advected into the model.The Goddard scheme is a slight modification from the PLin scheme for ice-water saturation.In general, all the MP schemes are known to influence the rainfall simulations at fine grid resolution by influencing the water phase component (Li et 235 al., 2017).In addition to the physics options mentioned above, two LSMs considered in the present work are, Simple 5-layer Soil Model (Slab; (Dudhia, 1996)) and the Noah LSM (Chen and Dudhia, 2001;Tewari et al., 2004).Slab is based on simple thermal diffusion in the soil layers that has constant soil moisture availability but a prognostic soil temperature term.The Noah LSM is a modestly detailed model, which includes explicit land surface parameterization with prognostic soil moisture and soil temperature evolution and snow cover prediction.240 (Skamarock et al., 2005) Since each scheme is associated with a distinct feature, it is important to examine the effect of their interactions on the rainfall simulations.Table 3 provides the summary of the WRF physics schemes considered to simulate the extremely heavy rainfall event.scheme option as compared to other MP schemes.Further, most of the model runs are able to reproduce the spatial gradient in the rainfall amount, which is perhaps primarily due to the topographical variation in the region.For 265 locales below 1000 m, observations show distinctly lower rainfall as compared to the high elevation regions (> 1000 m).Further, distinct clusters corresponding to heavy rainfall event are observed in the northeast and northwest areas of the study region.These clusters are found to be consistent with the TMPA data; however, due to lack of surface rain gauge observations, amount of rainfall in these regions could not be verified at this stage.Incidentally, the observed heavy rainfall event in southeast part of the region is seen in a few WRF configurations, such as 270 configuration (b) and (c).In general, WRF simulated rainfall fields show a similar spatial pattern as that of TMPA rainfall product.However, the magnitude of WRF rainfall is significantly high as compared to TMPA and is attributed to the negative bias in TMPA for heavy rains.Figure 6 summarizes the comparison of WRF rainfall with rain gauge observations for the three domains.Figure 6 indicates that Domain 1 captures rainfall within the range of 150 to 400 mm for most of the WRF configurations.For Domain 2a and Domain 2b, increase in the predicted rainfall amount is noted, particularly; for 280 small rainfall thresholds.Further, the WRF runs still under predict extremely heavy rainfall and each of the configuration considered (across all the three domains) underestimated the rainfall amount more than 400 mm.However, the underestimation of rainfall is less in Domain 2b (G2C scale) compared to others, indicating the necessity of finer grid spacing as the first-order requirement for simulating the magnitudes of the extremely heavy rainfall events.The bias in the WRF simulations is typically due to number of interactive factors: (i) scale feedback 285 between mesoscale convection and large-scale processes within the model (Bohra et al., 2006); (ii) lack of local observations that can add mesoscale features (Osuri et al., 2012;Osuri et al., 2015); and (iii) lack of proper land surface processes (Niyogi et al., 2006;Chang et al., 2009;Osuri et al., 2017a).To assess performance of the WRF simulations, quantitative scores (MAE and RMSE) with respect to the observed data are computed for cumulative rainfall over the 4 days' period.Results pertaining to these are presented in Figure 7. 290  Figure 7 indicates that there is more error at the stations Dehradun and Haldwani, which received higher rainfall.295 The highest rainfall obtained in different WRF configurations for these stations was less than 500 mm and this underestimation is highlighted in the error statistics.The model results show higher error and variability in the simulations for the northern part of the domain as compared to southern.This is likely due to the complex terrain in the northern part of the domain.
To identify the 'best' and the 'worst' performing configurations, temporal errors across all the stations are summed 300 up and reviewed (Appendix B, Table B.1).From this, the configuration (b) that is, YSU PBL, KF CU and Eta MP, produces maximum error, whereas configuration (p) with MYJ PBL, BMJ CU and Goddard MP, gives minimum error.This was also the 'best' performing configuration across all the locations for the three domains.In addition to this, spatial analysis is also conducted, wherein MAE across the 4 days' is computed and averaged across the station locations.Corresponding results are presented in Appendix B (Fig. B.1).These results are consistent to the temporal 305 analysis and again configurations (b) and (p) give maximum and minimum error respectively.Therefore, through this analysis, it can be inferred that the WRF model with MYJ PBL, BMJ CU and Goddard MP schemes is the 'best' in simulating the spatial and temporal variability of the extremely heavy rainfall over the upstream region of the UGB.Why this combination emerged as the best performing is an intriguing but difficult question to address, and needs to be studied through more cases and observational analysis as it becomes available.Note that the rainfall 310 prediction is the combination of many nonlinear, interactive factors including behavior of each configuration and cannot be realistically studied with the sparse rainfall data and absence of vertical sounding observations.Some possible factors that could contribute would be that local boundary formulation in MYJ may be more appropriately capturing the vertical environment in the complex terrain as compared to the nonlocal YSU scheme which seeks to simulate vertical mixing and boundary layer evolution using averaged and grid representative fields (Alapaty et al., 315 1997).As regards to the BMJ CU emerging in the top configuration, there are a number of studies for the ISMR where it has emerged as performing "overall best" (Ratnam and Kumar, 2005;Vaidya, 2006;Rao et al., 2007;Kumar et al., 2010;Mukhopadhyay et al., 2010;Srinivas et al., 2013).As for the MP scheme, there are limited studies in comparison to those that have studied the CU configuration for the ISMR.Further, the MP scheme performance has been evaluated for tropical cyclone cases because of the warm versus cold pool processes that are 320 critical in the simulation of the cyclone intensity.Of those available in literature, studies such as Sing and Mandal (2014) found that the Goddard scheme has a "slightly better" performance than other schemes.This conclusion is also supported by studies such as Choudhury and Das (2017) and has been used in hailstorm studies such as Chevuturi et al., (2014).Note that in citing these studies there is no claim being made about proving or even explain why the configuration that has emerged as the 'best' is indeed such.What these studies do provide is a reasonable 325 basis to support the notion that the 'best' configuration that has emerged is realistic and plausible to be considered as such.
The impact of downscaling ratio on the rainfall simulations is addressed next.On comparing the simulations from G2R and G2C domains with the rain gauge data, it is noted that the former gives less error for most of the locations (Appendix C).G2C scale has large resolution (grid spacing) gap from outer to inner domain in comparison to G2R, 330 which could result in less accurate initial and lateral boundary conditions, and consequently, more simulation errors in G2C.Another possibility is that the metric being used, which is the rainfall observation from in-situ data, itself is more conservative with regards to the grid in which rainfall occurs in the coarser domain and may slightly favor the G2R.However, on reviewing the overall structure of rainfall fields and the amounts across the domain, results suggest that the G2R scale with moderate downscaling ratio may be better suited for simulation of the extreme 335 rainfall event as in the present case study.The results are found to be consistent with other studies, such as by Liu et al., (2012), wherein moderate ratio of 1:3 is found to perform best.However, it is to be noted that errors corresponding to the grid point nearest to the rain gauge are considered here for comparison.Result may vary upon selection of another grid point.

Impact of Different Parameterization Schemes 340
Although configuration (p), with MYJ PBL, BMJ CU and Goddard MP, appears to be the 'best' physics configuration for the study region, significant variability exists among the simulations pertaining to different configurations of the WRF model.This variability causes significant uncertainty across different runs, which is quantified through computation of SE and CV (Fig. 8 and 9).Deviation in model simulations with respect to observed data provides the SE, however, CV gives variation within different model simulations.345  Most of the model configurations have SE value clustered around 1 (Fig. 8), indicating that the variability in simulated rainfall is similar to the observed rainfall.However, variability in northeastern part of the domain is 355 observed to be high compared to others.Same is reflected in the CV plot (Fig. 9), wherein grid points around the Chamoli station (on the northeastern side) have CV between 41-51 %, whereas stations closer to Uttarakashi and Tehri have values ranging between 11 -30 %.Further, grid points closer to Dehradun have low CV value, which could be due to the models consistently underestimating the rainfall in this subdomain.Southern part of the region, which received low rainfall, also exhibited high variability.In general, it can be inferred that uncertainty in rainfall 360 is more in the northeastern part compared to northwest.The regions that received very high or very low rainfall during this period also displayed higher uncertainty.Uncertainty in rainfall simulations varies between the domains, with Domain 2b having maximum uncertainty.This could be attributed to high variability in the simulated values at higher spatial resolution.
Since consideration of different parameterization schemes is the reason for variability in rainfall simulations, it is of 365 interest to understand how former influences the model output.For this, the average cumulative rainfall over the region, across different configurations is considered.The differences between various configurations is evaluated to assess the influence of PBL, CU and MP parameterization schemes on the rainfall simulations.Results for the same are presented in Figure 10.It is noticed that, in general, the WRF configurations with KF convective scheme produce rainfall of higher magnitudes.This result is consistent with previously conducted studies (Gallus Jr, 1999;Fonseca et al., 2015;Pieri 375 et al., 2015).For PLin MP, it is noted that considering YSU PBL along with KF CU scheme has a synergistic effect, leading to maximum amount of rainfall over the region.This additive effect could be attributed to the YSU being a non-local scheme making it suitable to convective, unstable PBL conditions (Bright and Mullen, 2002).Upon changing the PBL scheme (from YSU to MYJ), and maintaining the convective scheme as KF, notable difference in the fields is simulated (as shown by red circles in Fig. 10).This difference obtained for changing the PBL (with 380 PLin MP and KF CU), is found to be either equivalent to or more compared to the case when only CU is changed (from KF to BMJ) under either YSU (shown by blue plus sign in Fig. 10) or MYJ PBL (shown by yellow star in Fig. 10) conditions.BMJ CU, irrespective of the PBL scheme, results in less simulated rainfall across the region.
Therefore, it can be inferred that with PLin MP, CU plays a dominant role in determining the amount of rainfall over the region.For WSM6 MP, within YSU PBL scheme, changing the CU (from KF to BMJ) parameterization 385 produces significant variability (displayed by blue plus sign) in rainfall than changing the PBL scheme itself (from YSU to MYJ PBL with KF CU).However, with MYJ PBL, the effect of changing CU scheme is insignificant (yellow star in Fig. 10).Furthermore, with BMJ the difference in rainfall produced due to changing the PBL is minimal.With Eta and Goddard MP, CU schemes have significant influence on rainfall irrespective of the PBL condition.It can be concluded from this section that, for Eta and Goddard MPs, the choice of PBL and CU schemes 390 dominates the rainfall simulation.The relationship between PBL and CU for PLin and WSM6 MPs is interlinked and the choice of CU appears to have dominant impact on the simulation of rainfall over the region.

Impact of Land Surface Boundary Condition
The main differences in the two LSMs -Slab and Noah, considered here are related to (i) the soil depths along with the inclusion of land surface processes and (ii) the temporal evolution of soil moisture.Slab is a relatively simple 395 LSM with 5 soil layers (at 1, 2, 4, 8 and 16 cm depths) and uses a thermal diffusion equation to compute surface fluxes based on a surface temperature and drag coefficient formulations.Noah LSM is modestly detailed (compared to Slab) LSM with 4 soil layers (at 10, 30, 60 and 100 cm depths) and explicit representation of land surface parameters, which includes the effect of soil moisture changes, snow cover, evapotranspiration and hydrologic  In a number studies, the differences in the surface energy fluxes simulated by the choice of different LSMs i.e.Slab versus Noah has been discussed (see (Niyogi et al., 2016) for a review).The main reason being that the surface processes affect the boundary layer feedbacks which in turn create zones of mesoscale convergence that can affect the location and intensity of convection.These convective systems then contribute to the simulated rainfall.The 420 results obtained in this study emphasize this feature with differences in the rain amounts and locations through the domain in response to the change in LSM.Better performance of using Noah model could be attributed to temporal evolution of soil moisture fields and thus, the importance of soil moisture initialization over India for extreme weather conditions is highlighted through this work (Osuri et al., 2017a).

Comparison between Rainfall from the WRF, CORDEX and FNL datasets 425
Simulated rainfall from the WRF model runs is assessed with respect to the CORDEX downscaled data and the NCEP FNL reanalysis dataset.To achieve this, it is necessary to bring all the datasets to a common spatial resolution.Therefore, WRF simulated rainfall is upscaled, through averaging of the grids, to match the resolution of NCEP FNL (1°×1°) and CORDEX (0.5°×0.5°) data.For the analysis, simulations pertaining to the 'best;' performing configuration (p), are only considered.Bias (β) in rainfall simulations from the three datasets 430 corresponding to 18 rain gauge locations is obtained, results for which are presented in Figure 12.From Figure 12 (a), it can be observed that the NCEP FNL data overestimates the rainfall for most of the locations.435 Upon dynamic downscaling of the FNL data through the WRF model, rainfall simulations improved over the UGB region.Locations such as Mussorie, Pauri and Roorkee, which have shown β between 2.5 to 3.5 in the FNL data, reduced to 0 to 0.25 in the WRF simulations.Uttarakashi and Pithoragarh locations having small bias in the FNL data, show similar small bias in the WRF simulations.Dehradun along with three stations from south-eastern region, (Mukteshwar, Haldwani and Nainital), which recorded heavy rainfall (Section 2.1), are observed to have small bias 440 in the FNL data, and the rainfall at these location is underestimated by the WRF model.Overall, rainfall simulations from the WRF model (for all the three domains) have less β compared to the FNL data even after upscaling to the resolution of 1°×1°.As expected, upon upscaling, the spatial variability between the domains is reduced due to averaging across several grid points.Nainital and Haldwani, which are not simulated well by the WRF model (Section 3.1.1),have shown maximum underestimation in the range of -0.6 to -0.8.It is noticed that despite upscaling, spatial variability is preserved in the WRF simulations which is not seen in the CORDEX data, wherein rainfall is consistently underestimated, within the same range, across all the locations.In addition to this, slight variability in rainfall across the three domains is noticed unlike the earlier case (a), wherein the WRF simulations are compared with the NCEP FNL data.This is 455 attributed to the fact that in this case the simulated data is upscaled to 0.5°, whereas in the previous case it was upscaled to 1°.
From the above analysis, it is evident that the WRF model can simulate extreme precipitation better than the CORDEX data and the reanalysis data.This can be attributed to increase in spatial resolution, and better representation of surface and meteorological features, with respect to the lateral boundary conditions as suggested in some of the previous works such as by Argüeso et al., (2011); Mishra et al., (2014); Giorgi and Gutowski Jr (2015); Singh et al., (2017).

Summary and Conclusions
The main aim of this study is to provide a general guideline for setting up the WRF model configuration to simulate heavy rainfall events.For the analysis, an extremely heavy rainfall event, which occurred from 15 to 18 June 2013, 465 over the Ganges basin, in the foothills of Himalayas in the Uttarakhand State of northern India is considered.
Ensemble experiments are conducted with the WRF model with different grid spacing, four microphysics schemes, two cumulus parameterization schemes, two planetary boundary layer schemes and two land surface model conditions.The rainfall simulations are evaluated against the observed rain gauge data and the TMPA precipitation data.The WRF configuration with Goddard microphysics, Mellor-Yamada-Janjic planetary boundary layer 470 condition and Betts-Miller-Janjic´ cumulus parameterization scheme is found to perform 'best' in simulating this extremely heavy rain event.Although complex interactions are observed between different physics options, microphysics schemes are noticed to influence the spatial pattern of the rainfall, while the choice of cumulus scheme is found to modulate the magnitude of the simulated rainfall.Upon analyzing the impact of downscaling ratios on rainfall simulations, it is concluded that downscaling from global to regional scale with moderate downscaling ratio 475 may give least model errors and thus, be considered as suitable for reproducing the extreme rainfall event.In addition to this, effect of land surface models (LSMs) on rainfall simulations is also assessed in this paper.The Slab LSM significantly underestimates the rainfall values, and incorporating Noah helped improve the performance.
In addition to the sensitivity experiments, the WRF simulated rainfall is also compared with the CORDEX downscaled data and the NCEP FNL reanalysis data.The NCEP FNL data is found to overestimate the rainfall 480 whereas, the CORDEX downscaled data significantly underestimated this event, indicating that care must be taken while employing global datasets for regional analysis.The WRF simulated rainfall, on the other hand, has least bias.Through this, it can be established that the rainfall values obtained from the high-resolution mesoscale model can be effectively used in hydrologic models for realistic streamflow estimates.
Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2017-533Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 30 August 2017 c Author(s) 2017.CC BY 4.0 License.Thus, the tasks undertaken in this work are: (i) quantitative verification of the WRF model to simulate an extremely heavy rainfall event; (ii) assessment of sensitivity of the model simulated rainfall to different parameterization options, downscaling ratios and land surface models; and (iii) comparison of the WRF simulated rainfall with the global reanalysis data and the Coordinated Regional Climate Downscaling Experiment (CORDEX) downscaled data to investigate the impact of local versus global factors on rainfall simulations.A related objective is to use the model 75 results and provide recommendations on a possible optimal choice for model configuration to simulate such events in the region.

Figure 1 .
Figure 1.Topography of the study region (shown with black outline) as represented in the WRF model for (i) Domain 1 -27 km grid spacing; (ii) Domain 2a -9 km grid spacing (downscaling ratio -1:3); and (iii) Domain 2b -3 km grid spacing (downscaling ratio -1:9).Locations of the rain gauge stations within the UGB are presented as black dots in Figure 1 (iii).

Figure 2
Figure 2 presents daily and cumulative rainfall data from 15 to 18 June 2013 (obtained from the Indian Meteorological department (IMD) and literature (Ray et al., 2014) for the 18 official rain gauges located within the UGB.

Figure 3 .
Figure 3. Cumulative rainfall during 15 -18 June 2013 in the upstream region of the UGB obtained from the Tropical Rainfall Rainfall values from 15 to 18 June 2013 (consistent with the observed data) extracted from 12 climate scenarios (6 models × 2 RCPs), are presented as boxplots in Figure 4. CORDEX grids falling within the geographic extent of Domain 2b, along with the rain gauge data are spatially averaged to obtain single rainfall value for each day.Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2017-533Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 30 August 2017 c Author(s) 2017.CC BY 4.0 License.

Figure 4 .
Figure 4. Boxplot presenting the variability in daily and cumulative rainfall obtained from the CORDEX downscaled data for six These large-scale conditions are regridded by the model domain considering the grid spacing, and local topographical as well as other terrain conditions.As is common for most WRF studies over the Indian region, National Centers for Environmental Prediction (NCEP) global FiNaL (FNL) analysis dataset, based on GlobalData  190    Assimilation System (GDAS) with Global Forecast System (GFS) is considered.The FNL data is available at a coarse resolution of 1° × 1°, at every six hours' interval -00, 06, 12 and 18 UTC (Coordinated Universal Time) and is used to provide initial and boundary conditions to the model.The WRF model is initialized using the FNL dataset Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2017-533Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 30 August 2017 c Author(s) 2017.CC BY 4.0 License.from 14 June 2013:00 UTC through 19 June 2013:00 UTC, for 121 hours of forecast that was output at 1-hour interval.The lateral boundary conditions in the WRF model are updated at 6-h intervals.Considering the short 195 duration of the run, the model was forced with fixed Sea Surface Temperature (SST) throughout the integration, and no regional data assimilation is carried out.The land surface boundary conditions are taken from the Moderate Resolution Imaging Spectroradiometer (MODIS) International Geosphere-Biosphere Programme (IGBP) 21category landuse/cover fields that are available with a horizontal grid spacing of 10 min.Three telescopically-nested domains are used in this study -the parent domain (Domain 1) is fixed between 60°E and 100°E with grid-spacing 200 of 27 km; the first nested domain (Domain 2a) covers 70-85°E, 22-37°N with 9 km grid spacing and is indicative of scheme and the WSM 6 scheme are based on the parameterization fromRutledge and Hobbs (1984), former has 5-230

Figure 6 .
Figure 6.Scatter plots between the rainfall data from the rain gauges and the WRF simulations for (i) Domain 1; (ii) Domain 2a; and (iii) Domain 2b for (a) to (p)* WRF configurations.*Refer to Appendix A (TableA.1) for the list of the WRF configurations.

Figure 7 .
Figure 7. Root mean square error (top panel) and mean absolute error (bottom panel) computed temporally for (i) Domain 1; (ii) Domain 2a; and (iii) Domain 2b for (a) to (p)* WRF configurations.*Refer to Appendix A (TableA.1) for the list of the WRF configurations.

Figure 10 .
Figure 10.Difference in simulated rainfall due to PBL, CU and MP parameterization schemes corresponding to (i) Domain 1; (ii) Domain 2a; and (iii) Domain 2b over the UGB region.
Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2017-533Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 30 August 2017 c Author(s) 2017.CC BY 4.0 License.processes such as runoff and drainage (to sub surface layers).Further, in the Noah LSM soil moisture and 400 temperature is prognostically computed for each of the 4 soil layers, whereas in the Slab LSM only soil temperature is prognostic and moisture is considered as a constant value based on the land-use.To understand the influence of each LSM on the rainfall estimation, simulations using the Slab LSM are conducted for the 'best' and the 'worst' performing configurations from the Noah LSM case ((p) and (b), respectively).Results comparing the two land model based runs are presented in Figure 11.405

Figure 11 .
Figure 11.Scatter plot (top panel) and mean absolute error (bottom panel) for the observed rainfall data and the WRF simulations (for (b) and (p) configurations) pertaining to Noah and Slab LSMs corresponding to (i) Domain 1; (ii) Domain 2a; and (iii) Domain 2b.The Slab LSM based run significantly underestimates the rainfall in comparison to the Noah LSM.For example, the 410 locations which recorded rainfall greater than 400 mm have the Slab LSM based simulated values in the range of 100 -150 mm.As stated earlier, although Noah LSM also underestimated the rainfall for such stations, the bias with the Noah LSM is significantly less than the Slab LSM (-26 % with the Noah LSM, in contrast to -64 % with the Slab LSM for domain 2a and (p) configuration).Further, MAE in rainfall is also found to be higher with the Slab LSM in comparison to the Noah LSM.This is essentially due to significant underestimation of rainfall during 16 and 17 June 415 2013 by the Slab LSM.

Figure 12
Figure 12 (b) presents comparison between the CORDEX data and the WRF simulations upscaled to the resolution 445 of 0.5°.Rainfall is underestimated across all the locations by the CORDEX downscaled data, with most of the models having β in the range of -0.9 to -1.The WRF simulations are observed to have less β compared to the CORDEX data.Locations in the northeast along with Roorkee, Ranikhet, Almora and Pithoragarh from the southern part of the region are noticed to have negligible β in the WRF simulations.Northwestern locations (except Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2017-533Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 30 August 2017 c Author(s) 2017.CC BY 4.0 License.The analyses presented in this paper are subjected to a few limitations: first, results are limited to the physics 485 parameterization schemes considered in this paper, and may vary upon inclusion of other schemes; second, the best configuration obtained needs to be evaluated for reproducing another extremely heavy rainfall event; third, only two sets of downscaling ratios, i.e., 1:9 and 1:3 are tested in the current work.The sensitivity of simulations pertaining to other downscaling ratios should be tested in future; and fourth, only G2R and G2C sensitivity is assessed in

Table 2 .
List of six CMIP5 AOGCMs considered to obtained CORDEX dynamically downscaled RCM simulations

Table 3 .
Configuration of the WRF model considered for simulation of rainfall