A comprehensive view on trends in extreme precipitation in Nepal and their spatial distribution

We investigate trends in monsoon and extreme precipitation in Nepal based on rain gauge measurements. We find that precipitation amounts in Nepal vary considerably in space and time. The number of occurring extremes and the amount of precipitation are controlled mainly by the Indian summer monsoon. Almost all extreme precipitation events, recorded by 98 considered meteorological stations, occur during the Indian summer monsoon with the maximum in July. For Nepal in general, we find that the amount of precipitation and number of extreme events per station are neither significantly increasing nor decreasing between 1971 and 2010. However, on a regional scale we identify areas with positive and negative trends. A comparison of the combined precipitation time series with the ENSO 3.4 index reveals a connection between ENSO and the variability in monsoon precipitation. The correlation with the number of monsoon extremes vanishes for increasing percentiles. We investigate trends of upper percentiles of daily precipitation which pinpoint regions of increasing and decreasing extremes. These patterns are similar to spatial patterns in mean monsoon precipitation trends, whereas the median of the precipitation distribution undergoes only minor changes. Further analysis using extreme value theory confirm the prevailing trends from quantile regression for most stations and depict strong changes in return levels. Especially for Far‐West Nepal, we find robust evidence for a systematic increase in extreme precipitation.


Introduction
The Special Report on Extreme Events of the Intergovernmental Panel on Climate Change (Field et al., 2012) and the Fifth Assessment Report (Stocker et al., 2013) state that the magnitude of precipitation extremes and their frequency of occurrence have been increasing over most of the globe. They further conclude that a shift in the distribution, for instance a shift in the mean, can affect the extremes. For south Asia a change in return periods to more frequent extreme precipitation is projected for different emission scenarios. Bookhagen (2010) finds that the mountainous Himalaya has almost twice as many extreme events as the Ganges Plain or the Tibetan Plateau. Nepal, a country with complex topography, home to the highest peaks in the world, large glacier systems, and consequently great exposure to natural hazards like flash floods and landslides needs good estimates of changes in extreme temperature and precipitation for risk management. In this study, we comprehensively assess systematic changes in monsoon and extreme precipitation in Nepal. We hope that these findings will be useful to risk management which relies on correct estimations of changes in return values. * Correspondence to: P. Bohlinger, University of Bergen, Geophysical Institute, Allégaten 70,  Natural variability has to be considered and explored to understand changes in precipitation and extremes. Several studies agree on the connection between the El Niño-Southern Oscillation (ENSO) and the magnitude of the monsoon precipitation. The strength of this connection has been subject to continuous debate (Trenberth, 1997;Kumar et al., 1999;Shrestha, 2000;Ichiyanagi et al., 2007;Sigdel and Ikeda, 2013). The general opinion condenses to a negative correlation between ENSO and the precipitation amount with a varying strength over time.
Recently, also a lagged relationship between the Pacific quasi-decadal oscillation and the monsoon precipitation in Nepal was found (Wang and Gillies, 2013).
In order to homogenize the efforts to assess extremes worldwide and implement a common practice approach, the Expert Team on Climate Change Detection, Monitoring and Indices (ETCCDI) published a list of indices of climate extremes in the report WCDMP-No. 72 in 2009. In their report, they state that extreme weather phenomena are often not directly tied to environmental disasters. The impact of extremes depends to a large degree on the vulnerability and resilience of a system. A systematic change in extremes, however, is accompanied by a systematic change in hazardous weather situations leading to disasters with high socio-economic impact. An important piece of advice in the ETCCDI report is that changes in very rare extreme 1834 P. BOHLINGER AND A. SORTEBERG  events are not well represented by these indices and should be treated using the extreme value theory.
Utilizing ETCCDI indices for station-based research over central and south Asia, Tank et al. (2006) concluded with a warming signal but non-coherent precipitation trends. Zooming in on the Himalayas, limited work has been done assessing precipitation trends and possible future changes. This is partly due to limited availability of rain gauges compared to other surrounding countries (Dhar and Nandargi, 2000;Nandargi and Dhar, 2011). However, studies have shown evidence for a systematic increase in temperature and natural hazards (Shrestha et al., 1999;Khanal, 2001, 2002;Shrestha, 2005). Studies on changes in annual, seasonal, and extreme precipitation in Nepal produced less univocal results. Both increasing and decreasing trends have been reported at different sites across the country as well as a sensitivity to elevation, while no consistent trend could be found for the monsoon season over the entire country Shrestha, 2005;Ichiyanagi et al., 2007). Nandargi and Dhar (2011) studied extremes along the Himalayas and stated that there was a sudden increase in the frequencies of extreme rainfall events in the 60s, but a decrease in the first decade from 2000. Nepal, however, was excluded since the respective dataset did not include the dates of occurrence of extreme events. Caesar et al. (2011) investigated changes in ETCCDI indices across the Indo-Pacific region. Nepal was only represented by seven stations which indicated an increase in east and west Nepal and a decrease in central Nepal. Due to data sparseness, however, caution should be taken when associating trends at single stations with larger regions of Nepal.
Few studies exist that focus on changes in extreme precipitation in Nepal. For instance, Shrestha (2005) found indications of an increasing trend in extreme precipitation in Nepal. Other studies computed least-square fits on ETC-CDI indices (Baidya et al., 2008;Shrestha et al., 2017). A quite comprehensive attempt on assessing rainfall and its extremes in Nepal was published in 2009 by the Practical Action Nepal Office in form of a report (hereafter N09) with the title 'Temporal and spatial variability of climate change over Nepal (1976Nepal ( -2005' (Marahatta et al., 2009). Spatial maps of precipitation, return levels and changes of seasonal and annual rainfall were produced using Ordinary Least Square Regression (OLSR) for temporal changes and a Gumble type 1 distribution for computing stationary return levels. However, this report only computed trends for the means and not for extreme precipitation. The studies including Nepal are summarized in Table 1.
We highlight three issues which we complement in our study. First, we extend trends estimated with OLSR from N09 to different quantiles of daily precipitation. The reason being that OLSR is very sensitive to questionable values and outliers especially at the beginning and end of time series. Secondly, we estimate changes in return values for block maxima which has not yet been done. Thirdly, instead of a Gumble type 1 distribution we use the full Generalized Extreme Value (GEV) distribution, which, includes a shape parameter. This is because marginal rainfall distributions at the stations show shape parameters which despite large uncertainty, are likely to be different from zero. This in turn can have considerable impact on the distribution and consequently the estimation of return levels.

Data and methods
We start out with 278 meteorological stations from the Department of Hydrology and Meteorology (DHM) in Nepal (Figure 1). Purchased station data consist of 40 years of 24 h precipitation amounts (9 am previous day to 9 am local time, UTC/GMT +5:45 hours) for the years 1971 to 2010. Daily rainfall data analysed in this study was obtained using US standard 8 'diameter manual rain gauges and undergoes basic quality control like removing negative values and outliers.' More information about the different measurement networks in Nepal can be found in Talchabhadel et al. (2016). Data is not continuously available and some stations are not measuring precipitation for the entire period but for some years only. Stations are well distributed over the country although the number of stations decreases with increasing elevation. Relatively densely populated areas like Kathmandu and Pokhara have a higher number of stations close-by while other areas, in particular the Mid-West and Far-West Nepal stations are more sparse. For describing regional characteristics in our results we use the names of the different regions in Nepal as shown in Figure 1.
Since there is a vast amount of missing values in the dataset, only stations offering a set of sufficiently representative data are chosen for further investigations (see Section 2.1 for details). We provide a complete table with station names and their basic meta-data as supporting information (Table S1, Supporting Information). In Section 3.4, two arbitrary stations are used to illustrate the applied concepts and findings. These two stations with the identification numbers 512 and 1111 are also marked in Figure 1.

Selection of appropriate stations
We selected stations based on a minimum amount of observations available and a desire to cover most of Nepal's climatic zones defined by elevation as sketched in Shrestha et al. (2017). Considering all available stations, strong precipitation events happen mostly during the Indian summer monsoon from June to September (Figure 2). Constraints for selecting suitable stations are therefore mainly based on these four months rather than the entire period. The challenge of choosing representative sites is approached by allowing only stations which have at least 75% data coverage for the period  in each monsoon month. Since too many stations disqualify when applying this constraint to the entire period of 40 years (100%) we ease the restriction to 30 years (75%) which is the second constraint. An optimal time period of 30 years from the total 40 years is chosen where the number of available stations fulfilling the criteria is maximized. The maximum number of stations is 112. Of these 112 stations at least 75% (84 stations) have to fulfil the quality criteria over the same time period (third constraint). The choice of 75% instead of 100% is due to the fact that with 75% we cover all regions in Nepal and all climatic zones but the High Himalayas which are higher than 2500 m.a.s.l. (Shrestha et al., 2017). Choosing 30 years as a suitable reference period is mainly to retain time series that are somewhat representative for climate and long enough to level out variations that are associated with e.g. ENSO. In our case the reference period lasts from 1979 to 2008. The availability of stations is sensitive to the constraints as illustrated in the following example. Keeping the second constraint untouched since we want to retain a 30 year data record, we can test the sensitivity to the remaining constraints. For instance, using 70% or 80% for the first criterion results in 112 and 111 stations, respectively. Changing only the third constraint to 70% or 80% would change the number of remaining stations to 122 and 100, respectively. The position of the optimal time period varies only by a few years when changing the constraints above as described. We note that our results are not very sensitive to the number of stations included. Ichiyanagi et al. (2007) demanded availability of 80% of the entire period (1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996). Choosing a certain degree of coverage for an entire period could mean less coverage in the monsoon months where most of the extremes are found and more coverage in the rest of the year. This is why we introduced an additional availability constraint for each monsoon month to avoid that potential bias on our results.
These constraints result in 112 stations which are rather uniformly spread over the country ( Figure 1) and range from 72 m.a.s.l. to 2742 m.a.s.l. Higher mountain regions have unfortunately no stations fulfilling the requirements and were excluded from further investigations. It would be of course ideal to have close to full data coverage instead of 75%. This choice has to be made specifically for every study. We aimed for a balance between spatial and temporal availability. A uniform spatial spread gives confidence in capturing the total picture of precipitation in Nepal. Using the constraints we have chosen, we end up with a set of 112 stations that is largely a subset of the 166 stations used for the N09 report.

Homogeneity of time series
We tested the remaining 112 stations for homogeneity to reduce the risk to include artificially induced trends, e.g. by stations changing location. Due to the sparseness of stations in highly complex topography we did not find it feasible to test against a reference time series or nearby stations as suggested by Alexandersson (1986). Instead we used a penalized maximal F test (PMFT) (RHtestsV4, Wang (2008a,b)). This software for testing homogeneity is freely available and described in detail on the ETCCDI homepage (http://etccdi.pacificclimate.org/software.shtml). We did not adjust time series when a significant (5% significance level) break point was found but rather rejected the respective records. Additionally, time series were visually inspected and three more stations excluded. Two of them depicted break points that were not significant but to the eye a shift in mean seemed to be probable. The third station clearly depicted a strong, sudden change in variance and was excluded due to its heteroscedasticity. In total we rejected 14 stations and proceeded with a subset of 98 both plain area and hill stations (Table S1).

Trend analysis for extreme events
We used two methods to assess trends in extreme precipitation. Our approach using extreme value theory and specifications for the Bayesian estimation of the GEV distribution parameters are described in Appendix.
In order to determine the trend of different quantiles we performed a quantile regression on daily precipitation for each station following the method from Koenker and Hallock (2001). Advantages of quantile regression versus OLSR can be summarized as follows. Quantile regression is more flexible as it is possible to explore the impact of covariates on all quantiles within a given distribution, meaning that a different behaviour in different quantiles can be observed. Additionally, quantile regression is robust to outliers and leads to reasonable estimates when the error distribution is non-normal. Based on these advantages quantile regression is well suited for this study to investigate changes in different quantiles of daily precipitation. Fan and Chen (2016) illustrate some advantages of this method when calculating trends in extreme precipitation indices across China.
Block maxima values like annual or seasonal maxima follow the GEV distribution which is defined by the three parameters for location ( ), scale ( ), and shape ( ). Given a set of values these parameters can be estimated for instance by maximum likelihood or Bayesian approaches which are both used in this study. We computed changes of the 100 year return level with 1 where Q 99 represents the 99th quantile from a GEV with the location parameter for year 2010 versus 1971 for a single station s: (1)

Results and discussion
The results presented in this chapter are computed with respect to the reference period explained in Section 2.1. For percental trends the respective climatological value is used. This allows us to directly compare stations and regions without biasing the results with records of different length. For the following analysis a day is considered a wet day when precipitation is recorded, irrespective of the amount.

Station monsoon climatology
Over a variety of high quantiles, almost all extremes were measured during the Indian summer monsoon period ( Figure 2). The higher the percentiles the more likely it is that it is recorded during the monsoon. The distribution flattens out for lower percentiles. July claims most extreme events with a substantial difference compared to June, August, and September. Most of the annual precipitation is measured during the Indian summer monsoon, ranging from 412 mm up to 4575 mm (Figure 3(a)). The contribution to the annual precipitation ranges from 54% to 88% depending on the station (Figure 3(b)). The monsoon mean precipitation indicates local maxima in the station cluster close to Pokhara in the north-west (Lumle), north-east of Kathmandu (Gumthang), and north in East Nepal (Num) (Figure 3(a)). The contribution of the monsoon precipitation to annual precipitation is highest close to the Indian border, close to Pokhara, Kathmandu and north-east of Kathmandu (Figure 3(b)). There are local minima in the contribution at the border between Far-West Nepal and Mid-West Nepal, and central in East Nepal. The spatial distribution of summer monsoon precipitation and the percentage of monsoon rainfall compared to annual rainfall is in agreement with findings from Shrestha (2000).
The spatial distribution for median values of wet days depict a similar pattern as described for mean monsoon precipitation (Figure 3(c)). The 99.5-percentiles of daily precipitation has the highest values along the Indian border, close to Pokhara, and north in East Nepal with maximum values of over 190 mm d −1 (Figure 3(d)). Using only wet days for the extreme percentile calculation does not change the pattern qualitatively. Extreme values become higher towards India and seem to be rather constant along the Indian border where the topography is less complex.
The number of wet days combined with the rainfall amount indicates whether it rains often or seldom and if the events are strong or weak. The number of wet days increases towards China and further into the Himalayas (Figure 3   and East Nepal along the northern border, it rains almost every day on average, whereas at the southern border and Far-West and Mid-West Nepal there are stations that register rain on only half of the 122 monsoon days. The standard deviation of the number of wet days (Figure 3(f)) is probably highly influenced by local topography and spatially very variable.

Observed time series and trends
A shift in mean precipitation can affect extremes which is why we first explore possible changes in means of monsoon precipitation in Nepal. Precipitation records vary strongly in time and space. Variability ranges from strong fluctuations between single days to multi-year and multidecadal oscillations.  Figure 4(a) shows trends of the 42 remaining stations. Although the spatial coverage of stations in Nepal is sparse when using 42 stations, we still get an impression of the rainfall trends in various regions. The precipitation records show a high spatial variability with positive and negative trends. There are hardly any statistical significant trends (Figure 4(b)). The significance is estimated using the student t-test for the 0.05% significance level. The occurrence of a similar trend at several stations can indicate real changes in precipitation for a larger region. As an example for this, around Kathmandu, several negative trends coincide indicating a decrease of monsoon precipitation in that area whereas trends close to the Pokhara region and in Far-West Nepal indicate an increase. We merge the seasonal time series of all 98 stations to overcome data gaps and to obtain a precipitation trend that can be associated with Nepal. This is done by creating standard z-normal scores (z-scores) which make precipitation sums at different stations comparable. z-Scores are calculated for 40 years (1971-2010) using a reference period of 30 years to calculate a mean and a standard deviation for each station. The reference period consists of the 30 years where most stations fulfil the data availability criteria at the same time as explained in Section 2.1. z-Scores consist of precipitation anomalies, calculated from the difference between the seasonal value and the 30 years mean standardized by the standard deviation of the respective stations. Monsoon precipitation z-scores from different stations are merged, resulting in a time series for all of the 40 years .
The estimated trend for the combined time series is not significantly different from zero at a 95% confidence level using a student t-test. This means that there are no indications of a decline or increase in precipitation ( Figure 5(a)). This is consistent with the high spatial variability of station-based trends and the existence of regions with a negative trend (around Kathmandu) and regions with positive trends like Far-West Nepal and close to Pokhara.
We explore the trends in the number of extremes (multiple upper percentiles) and count all extremes registered by stations from 1971 to 2010. The number of extreme events is normalized by the number of active stations at the time of the event. Normalization is necessary since the number of active stations varies substantially between 1971 and 2010 and thereby also the potential of measuring an extreme which would bias the result. Figure 5(b) shows time series of the amount of recorded extreme events with extreme defined as values higher than the displayed percentiles. We show results for the three percentiles 99.9, 99.5 and 99. The similarity between the time series of number of extremes and the combined precipitation time series increases as the percentiles decrease. We do not find a statistically significant trend. This suggests that occurrences of extremes have neither increased nor decreased. This conclusion is consistent with findings from Shrestha et al. (2000).

Variability in the observed time series
The merged time series show strong inter-annual variability. In particular ENSO has been associated with variability in rainfall over India, Bangladesh and other places around the Indian subcontinent by changing large scale circulation patterns in the south Asian region (Pant and Parthasarathy, 1981;Rasmusson and Carpenter, 1983;Parthasarathy et al., 1988;Kripalani and Kulkarni, 1997;Lau and Wu, 2001;Ashok et al., 2001). The El Niño phase is usually associated with less precipitation whereas the La Niña phase coincides with heavy precipitation and floods. However, this correlation seems to be subject to temporal variability as described by Kumar et al. (1999) indicating the complexity of this relationship. Using the Niño 3.4 region bounded by 120 ∘ W-170 ∘ W and 5 ∘ S-5 ∘ N (Trenberth, 1997), Chowdhury (2003) investigated seasonal flooding in Bangladesh and shows that for strong El Niño events this relationship is valid, whereas an increase of precipitation is also observed during moderate El Niño years underlining the spatial differences in the impact of ENSO.
In this study, we investigate the relationship between ENSO and precipitation amounts and the number of extremes in Nepal. We use the cold and warm ENSO periods calculated by the NOAA Center for Weather and Climate Prediction (Climate Prediction Center Internet Team, 2016) over the Niño 3.4 region using the extended reconstructed sea surface temperature version 4 (ERSST. v4) (Huang et al., 2015). We find a significant correlation (r ≈ −0.6) between our combined seasonal time series and the ENSO signal (mean June to August, 3.4 region) (Figure 6(a)). This is consistent with other findings (e.g. Shrestha, 2000;Shrestha et al., 2000) that related ENSO to the rainfall in Nepal using for instance the Southern Oscillation Index (SOI) which is derived from the mean sea level pressure difference measured in Tahiti and Darwin. Shrestha (2000) found a similar correlation of r ≈ 0.58 between percentage departure of Nepal monsoon rainfall (PNMR) and Monsoon Southern Oscillation Index (MSOI). For the monsoon months June to September a correlation between all-Nepal monsoon precipitation and SOI of r ≈ 0.43 to r ≈ 0.63 depending on the time period was found by Shrestha et al. (2000). However, the correlation between ENSO and the number of extreme events per season requires a closer look (Figure 6(b)). As the percentile increases the correlation vanishes whereas for events lower than 99th percentile, the correlation starts to resemble the cross-correlation function of ENSO and the combined precipitation time series (Figure 6(a)). This shows that there is in fact a connection between ENSO and moderate extremes (here: lower than 99th percentile) but no linear correlation for percentiles even further out in the tail of the distribution.

Trends for extreme events
In this chapter, we use two different methods to assess changes of extremes. Comparing results from different methods gives a notion of sensitivity of observed trends to the method which goes beyond the obtained standard errors of a single method. We extend existing work using quantile regression to investigate and illustrate changes across a variety of high percentiles of daily precipitation. Finally, we use extreme value theory with the block maximum approach to make inference for very rare extreme events (100 year event) where the ETCCDI descriptive climate indices are not considered reliable anymore.

Investigating the evolution of quantiles over time using quantile regression
Temporal trends in extreme precipitation amounts are investigated by exploring changes of high quantiles of daily precipitation over the 40 year period from 1971 to 2010. This will answer questions like, does the value of the 99.5th percentile shift over time towards stronger or weaker extremes? Are the trends similar across multiple quantiles?
Considering station based trends of quantiles, spatial variability is strong (Figure 7). It is noteworthy that only negative trends are statistically significant. In a total Nepal perspective there is no uniform increase or decrease of percentile values although some regions have several stations with similar trends indicating the presence of increasing or decreasing values for a given percentile in a larger area. The percentage of stations with positive and negative trends is summarized in Table 2. For example stations around Pokhara depict an increase of higher percentile values. A second region showing an increase of high percentile precipitation values is in most Far-West Nepal. Around Kathmandu a group of stations show negative trends. This pattern of regions of prevailing trends is consistent with the pattern in seasonal precipitation trends discussed in Section 3.2.
A tendency of stronger trends in the higher percentiles becomes visible when comparing trends of the 99.9, 99.5, 99, 95, 75 and 50th percentile. The median seems to be quasi constant for almost all stations which is consistent with a trend not significantly different from zero for the combined seasonal precipitation time series of all stations ( Figure 5). Some stations even have a shift in sign for different percentiles. The result of the quantile regression reveals that the increase and decrease in seasonal precipitation ( Figure 4) is not due to changes in the bulk of the distribution but rather caused by changes in the upper tail. Furthermore, this shows that the dispersion of precipitation amounts increases and decreases with time and that the conditional distribution of daily precipitation is positively skewed with a long upper tail. Figure 8 illustrates the behaviour of a time-dependent increase and decrease of higher percentiles as described above for stations 512 and 1111. Both stations show different trends in their quantiles. For station 512, trend curves of different quantiles converge with increasing time. Quantiles at station 1111 are diverging over time, resulting in a positive trend in seasonal precipitation which is visible in Figure 4. The fact that the pattern for the trends in the higher quantiles is similar to the trends in seasonal precipitation for July to September indicates that seasonal precipitation trends are strongly influenced by trends in higher quantiles.  . Trends from quantile regression for different stations based on quantiles obtained from daily precipitation. Stations that do not indicate a trend significantly different from zero at the 95% confidence level are marked using a circle, significant stations with a square. Confidence intervals are obtained using bootstrapping with 500 times of re-sampling. Results are little sensitive to the number of bootstrap samples. For trends in the median fewer stations are plotted because some have a median of zero which leads to division by zero when calculating the percental change. Their respective absolute trends are zero (of the order |10 −8 | mm decade −1 or smaller) and therefore neglected.

Modelling trends of extremes using extreme value theory
We extend the previous analysis using a block maximum approach to detect return level changes and their spatial distribution over Nepal for 100 year events. For the block maximum approach we use annual and monsoon maxima taken from our set of daily precipitation values. Results are obtained using Bayesian and maximum likelihood parameter estimation for the location, scale and shape parameters of the GEV distribution. Due to strong similarities between  Bayesian and maximum likelihood estimated parameters we only use the results from the Bayesian estimation for illustrating consequences of using time as a covariate for the location parameter. The systematic change over time which we described in the previous sections can also be seen by eye for the block maximum time series of the same stations 512 and 1111 ( Figure 9). For both stations, the annual maximum daily precipitation depicts a change over time. At station 512 the annual maxima decrease and at station 1111 the extremes seem to increase over time, thus matches our assumption of an underlying linear trend in the location parameter. Starting from the year 1971 at station 512 a former 12-year event is now a 100-year event meaning that the probability of occurrence of such a high value has decreased. At station 1111 a 580-year event has become a 100-year event suggesting a higher probability of occurrence. The results from the block maximum approach are consistent with the trends in the high percentiles using quantile regression.
After illustrating our approach we fit a GEV distribution (Equation (A1)) to each station to assess the spatial distribution of trends in annual maxima. Since we are dealing with 98 stations we proceed with the same simple covariate for each station. We allow for a linear change in the location parameter at each station for the temporal evolution. The scale parameter is held constant as well as the shape parameter, which is notoriously difficult to estimate and comes along with large uncertainties. A linear change over time was chosen since many station-specific time series exhibit a somewhat linear trend, whereas for the variance we did not discover such systematic changes. Ideally a more physical covariate could be chosen. However, this would require knowledge about the underlying mechanisms for the change which are not yet established. ENSO as a more physical covariate does not significantly improve the fit. Figure 10 illustrates the change of the 100 year return level in 1971 compared to the same return levels in 2010. Data from each station have been fitted separately with the approach described above. The respective change of the 100 year return level in percent was calculated using Equation (1) in Section 2.3. For the maximum likelihood estimated models we performed a likelihood-ratio test to determine whether the trend contributes to improve the model. The null hypothesis is that the model without a linear trend in the location parameter is the same as the model including the trend. Rejecting this hypothesis therefore means that the models are significantly different from each other and the trend improves the model. Additionally, we compare the magnitude of the trend estimates and the respective standard errors. For the Bayesian estimation we decide on using the 95% credible intervals that come along with the MCMC sampling of the trend parameter. Since for the maximum likelihood estimation (MLE) the standard error approach and the likelihood-ratio test result in the same significant stations, except one, we just display the likelihood ratio test results with respect to 95% confidence interval. Bayesian and MLE estimations agree on the same significant stations, except one. This adds confidence to the robustness of the results which is supported by the fact that when considering only monsoon maxima instead of annual maxima the significant stations are the same.
In both estimation methods, the distribution parameters are very similar resulting in return values that are  (c) and (d) but for annual block maxima. Squares depict stations with significant changes considering a 95% credible/confidence interval. For maximum likelihood figures there is one station indicated by a transparent circle and square for which the parameters could not be estimated and the change is set to NaN. Topography is displayed in grey shading for orientation purposes with darker grey indicating higher elevation.
much the same which in turn is reflected in the change of the return values over time. Comparing these methods reflects the sensitivity of the results to the use of different methods. The spatial pattern is consistent with the pattern obtained from quantile regression. Computed changes are substantial and typically range from −30% to +30% increase of the 100 year return level. Roughly 55% (Bayes: 53%, MLE: 56%) of the stations indicate a positive trend whereas 45% (Bayes: 47%, MLE: 44%) show a decline of annual maxima. Considering monsoon maxima instead of annual maxima, both methods agree on a ratio of ca. 61% (Bayes: 60%, MLE: 61%) of the stations indicating a positive trend. This is comparable with ratios obtained from quantile regression. The trend pattern shows an inhomogeneous trend distribution over Nepal where only smaller regions agree on the sign of a trend. For instance between Pokhara and Kathmandu there is an area with significant negative trends which is also visible in the results from quantile regression. The negative trend in the higher quantiles around Kathmandu is much less univocal in the block maxima. In Far-West Nepal there are almost exclusively positive trends. Another important finding is the robustness and significance of the positive trend in Far-West Nepal. A homogeneous region in terms of positive trends where some show statistical significance leads to strong indications for a systematic change towards higher extremes in western Nepal. Our results differ partly from the monsoon mean trends in N09 which underlines the importance of considering extremes in an own framework. The trends do not seem to depend on the geographical parameters elevation, latitude and longitude. Using the parameters from the marginal distributions, we checked dependencies on geographical or physical parameters. Scatter plots reveal dependencies for the location offset and the scale parameter on latitude and elevation for both the MLE and Bayesian estimation (not shown here). Unfortunately, the observed dependencies are not enough to obtain a reliable fit of a spatial GEV using the mentioned covariates and creating a map of return values. The spatial model over-and underestimates return values compared to the marginal distributions at the station sites (not shown here). We conclude that a well-working spatial GEV-model probably has to be more complex, with spatially variable dependencies on physical variables.

Summary and conclusions
In this study, we investigated previously not considered elements of monsoon precipitation and extremes in Nepal on basis of rain gauge stations. We used quantile regression and extreme value theory to obtain a comprehensive view on trends in extreme precipitation.
Most events in the tail of the daily rainfall distribution occur during the Indian summer monsoon between the months June and September. The highest monsoon precipitation amounts were measured in the station cluster close to Pokhara in the north-west (Lumle), north-east of Kathmandu (Gumthang) and north in East Nepal (Num). Most of the monsoon contribution to annual rainfall is recorded close to the Indian border, around Pokhara, Kathmandu and north-east of Kathmandu. The highest absolute upper tail daily precipitation values can be found along the Indian border, close to Pokhara, and north in East Nepal.
Monsoon precipitation trends are heterogeneous across the country where a combined time series of 98 stations does not exhibit a significant trend. Seasonal monsoon precipitation is anti-correlated with ENSO (r ≈ −0.6). The number of extremes throughout a variety of higher quantiles do not depict a significant trend. There is linear dependency between ENSO and the number of moderate extremes (lower than 99th percentile). Percentiles above are not correlated with ENSO.
Quantile regression reveals a spread in the trends of different upper quantiles which can be responsible for trends in the seasonal precipitation time series at the same stations. This means that the upper tail undergoes a change rather than the bulk of the distribution. This in turn reflects in the characteristics of the distribution of daily rainfall. By computing sums or means of strongly skewed distributions, a trend in seasonal monsoon precipitation can occur.
Results using a block maximum approach from extreme value theory show a pattern similar to the one obtained by quantile regression. Especially Far-West Nepal depicts homogeneously positive trends throughout a variety of higher quantiles and block maxima where some show statistical significance at the 5% significance level. This concludes with robust evidence for a systematic change towards higher extremes in Far-West Nepal. It is important to note that trends in the annual and monsoon daily maxima differ in parts from the respective mean trends in the N09 report. This underpins the importance of considering extremes in an own framework and under non-stationarity which, due to the observed magnitude (change from −30% to +30%), should also be included in a spatial model. parameter to capture the visible changes of block maxima for the observation period. The GEV including a time dependent location parameter looks as follows: As priors for Bayesian estimation we used uniform distributions for the location offset and trend as well as for the scale and shape parameter. Additionally, the scale parameter is not allowed to become negative since that would make no physical sense. The posterior is obtained through sampling with a Markov chain Monte Carlo (MCMC) method where for the MCMC method the Metropolis-Hastings (MH) algorithm is used. Since station specific distributions can be very different the MH algorithm is initialized with the sample mean for the location parameter and the sample standard deviation for the scale parameter. The trend on the location parameter and the starting point for the shape parameter are both initialized with zero. For convergence we sample 90 000 times and use a burn-in period of 30 000 values. From the remaining values we skim every 10th value which is the basis for our final posterior distribution.

Supporting information
The following supporting information is available as part of the online article: Table S1. Metadata of 98 rain gauge stations in Nepal including the 90th, 99th and 99.9th percentile of daily precipitation.