An Analysis of Long-Term Rainfall Trends and Variability in the Uttarakhand Himalaya Using Google Earth Engine

: This paper analyses the spatio-temporal trends and variability in annual, seasonal, and monthly rainfall with corresponding rainy days in Bhilangana river basin, Uttarakhand Himalaya, based on stations and two gridded products. Station-based monthly rainfall and rainy days data were obtained from the India Meteorological Department (IMD) for the period from 1983 to 2008 and applied, along with two daily rainfall gridded products to establish temporal changes and spatial associations in the study area. Due to the lack of more recent ground station rainfall measurements for the basin, gridded data were then used to establish monthly rainfall spatio-temporal trends for the period 2009 to 2018. The study shows all surface observatories in the catchment experienced an annual decreasing trend in rainfall over the 1983 to 2008 period, averaging 15.75 mm per decade. Analysis of at the monthly and seasonal trend showed reduced rainfall for August and during monsoon season as a whole (10.13 and 11.38 mm per decade, respectively); maximum changes were observed in both monsoon and winter months. Gridded rainfall data were obtained from the Climate Hazard Infrared Group Precipitation Station (CHIRPS) and Precipitation Estimation from Remotely Sensed Information Using Artiﬁcial Neural Networks-Climate Data Record (PERSIANN-CDR). By combining the big data analytical potential of Google Earth Engine (GEE), we compare spatial patterns and temporal trends in observational and modelled precipitation and demonstrate that remote sensing products can reliably be used in inaccessible areas where observational data are scarce and / or temporally incomplete. CHIRPS reanalysis data indicate that there are in fact three signiﬁcantly distinct annual rainfall periods in the basin, viz. phase 1: 1983 to 1997 (relatively high annual rainfall); phase 2: 1998 to 2008 (drought); phase 3: 2009 to 2018 (return to relatively high annual rainfall again). By comparison, PERSIANN-CDR data show reduced annual and winter precipitation, but no signiﬁcant changes during the monsoon and pre-monsoon seasons from 1983 to 2008. The major conclusions of this study are that rainfall modelled using CHIRPS corresponds well with the observational record in conﬁrming the decreased annual and seasonal rainfall, averaging 10.9 and 7.9 mm per decade respectively between 1983 and 2008, although there is a trend (albeit not statistically signiﬁcant) to higher rainfall after the marked dry period between 1998 and 2008. Long-term variability in rainfall in the Bhilangana river basin has had critical impacts on the environment arising from water scarcity in this mountainous region.


a.
To collect, tabulate and analyse rainfall data collected from i) available surface observatories managed by IMD ii) reanalysis rainfall data from CHIRPS and iii) satellite rainfall data PERSIANN-CDR b.
To map annual and seasonal time series of rainfall dynamics and juxtapose the performance of reanalysis and satellite-based data products compared to ground station data c.
To evaluate annual, seasonal, and monthly distributions and trends in rainfall and rainy days and compare the results from the three sources through the application of non-parametric tests d.
To briefly discuss the regional socio-economic implications of observed rainfall variability

Environmental Characteristics of the Bhilangana River Basin
The Bhilangana river basin, a sub-basin of the upper Bhagirathi, is located in Uttarakhand State, western Himalaya, India. The basin, which occupies a total area of 1482 km 2 , lies between 30 • 32 71"N to 30 • 87 79"N latitude and 78 • 48 22"E to 79 • 03 53"E longitude ( Figure 1) at elevations ranging between 616 to 6635 m amsl. The headwaters originate at the terminus of the Khatling glacier (3718 m asl) and flow~80 km southwest to confluence with another river and feeds into the Tehri dam (India's highest elevation river dam), after which it is known as the Bhagirathi river. The dam supports a 2400-megawatt hydropower system which provides electricity to adjacent states (Punjab, Haryana, Uttar Pradesh, Jammu and Kashmir) including the national capital (New Delhi). Indeed, the Bhilangana river is the lifeline for approximately 1.3 million people (around 13 percent of Uttarakhand's population) and 295 villages [24] in terms of fresh water, irrigation, household activities and social wellbeing. Geologically, this area is primarily underlain by the Garhwal formation [25]. The major litho-tectonic units in the area are a) Central Crystalline and b) Garhwal groups [26]. Tectonically, this area is located in an earthquake zone "V", (i.e., seismically more active) and is frequently impacted by mass movements including landslides and debris flow following monsoonal precipitation in summer months (June-September) [27,28]. The area is prone to extreme rainfall events, such as cloudbursts and torrential downpour, some of which together with other extreme events such as glacier lake outburst floods have resulted in loss of life, property, and infrastructure in the region during recent past [18,29].
Remote Sens. 2020, 12, x FOR PEER REVIEW 3 of 25 a. To collect, tabulate and analyse rainfall data collected from i) available surface observatories managed by IMD ii) reanalysis rainfall data from CHIRPS and iii) satellite rainfall data PERSIANN-CDR b. To map annual and seasonal time series of rainfall dynamics and juxtapose the performance of reanalysis and satellite-based data products compared to ground station data c. To evaluate annual, seasonal, and monthly distributions and trends in rainfall and rainy days and compare the results from the three sources through the application of non-parametric tests d. To briefly discuss the regional socio-economic implications of observed rainfall variability

Environmental Characteristics of the Bhilangana River Basin
The Bhilangana river basin, a sub-basin of the upper Bhagirathi, is located in Uttarakhand State, western Himalaya, India. The basin, which occupies a total area of 1482 km 2 , lies between 30°32′71″N to 30°87′79″N latitude and 78°48′22″E to 79°03′53″E longitude ( Figure 1) at elevations ranging between 616 to 6635 m amsl. The headwaters originate at the terminus of the Khatling glacier (3718 m asl) and flow ~ 80 km southwest to confluence with another river and feeds into the Tehri dam (India's highest elevation river dam), after which it is known as the Bhagirathi river. The dam supports a 2400-megawatt hydropower system which provides electricity to adjacent states (Punjab, Haryana, Uttar Pradesh, Jammu and Kashmir) including the national capital (New Delhi). Indeed, the Bhilangana river is the lifeline for approximately 1.3 million people (around 13 percent of Uttarakhand's population) and 295 villages [24] in terms of fresh water, irrigation, household activities and social wellbeing. Geologically, this area is primarily underlain by the Garhwal formation [25]. The major litho-tectonic units in the area are a) Central Crystalline and b) Garhwal groups [26]. Tectonically, this area is located in an earthquake zone "V", (i.e., seismically more active) and is frequently impacted by mass movements including landslides and debris flow following monsoonal precipitation in summer months (June-September) [27,28]. The area is prone to extreme rainfall events, such as cloudbursts and torrential downpour, some of which together with other extreme events such as glacier lake outburst floods have resulted in loss of life, property, and infrastructure in the region during recent past [18,29].

Observed Data
Observed monthly rainfall and rainy days data for the period from 1983 to 2008 were obtained from the India Meteorological Department (IMD, Pune) ( Table 1). Since both rainfall and rainy-day records of IMD have data gaps, we carefully selected suitable stations for this study based on a) more than 30 years of data, b) up to 80% of data availability, and c) the most recent records. In some cases, we interpolated the missing data by applying long-term monthly averages. Eventually, three surface observatories were selected for both rainfall and rainy-day parameters. As per the IMD, a rainy-day is defined as a day with a minimum of 2.5 mm rainfall. Reliable and reasonably continuous monthly data for these ground stations were only available for the period from 1983 to 2008.

Gridded Rainfall Data
Along with observed data, one reanalysis (CHIRPS) and one satellite-based (PERSIANN-CDR) gridded rainfall product was derived using GEE. Conspicuously, the CHIRPS gridded daily rainfall product is a combination of observed and satellite rainfall and has low systematic bias [30][31][32][33][34][35] with an extended period of record. CHIRPS data is accessed from the website of the Climate Hazard Group (https://www.chc.ucsb.edu) (USAID, NASA, and NOAA).PERSIANN-CDR data is satellite-derived daily rainfall estimate using artificial neural networks (ANN), with high spatio-temporal resolution and accuracy [36][37][38][39], and retrieved from the website of NCAR and UCAR (https://climatedataguide. ucar.edu). The ANN estimates precipitation from cold cloud pixels and adjacent features using gridded satellite brightness temperature details. In addition to this, two intensive field visits were conducted in 2016 to assess the environmental and socio-economic impacts of rainfall changes.
To estimate long-term trends, observed monthly rainfall and rainy-day data were collected from 1983 to 2008, as per the existing observation records of IMD. This was further compared with gridded data available between 1983 and 2008, based on high spatial resolution and a similar period of record. Thereafter, the more recent spatial and temporal trends were carried out using reanalysis data for the period 2009 to 2018.

Data Handling and Statistical Applications in GEE
Annual, seasonal and monthly gridded rainfall products were derived from the GEE platform [40] using the ee.ImageCollection algorithm and applying a filter command (ee.Filter.calendarRange) to cover the study period . After processing the images (.sum), the clip function was used to confine the boundaries to the Bhilangana study region (Supplementary code 1). Both high-resolution modelled datasets were downloaded and analysed in the GEE code editor. To ensure spatial similarity in both gridded datasets, bilinear interpolation resampling technique was used [41]. To match the similar spatial information with station data, nearest point-pixel rainfall values of both the gridded datasets Remote Sens. 2020, 12, 709 5 of 24 were extracted in GEE platform for station localities (Supplementary code 2) [33,36]. Subsequently, pixel-based long-term spatio-temporal trend for both gridded data sets was performed in GEE using ee.Reducer.senSlope algorithm (Supplementary code 3). The Sen's slope (magnitude of trend) was tested through Mann-Kendall (alpha ≤ 0.05 and Z statistics = 1.96, obtained from the standard normal table) expression.
Additionally, gridded daily rainfall records were then aggregated to monthly and seasonal totals viz., pre-monsoon (MAM), monsoon (JJAS), post-monsoon (ON), and winter (DJF), following the scheme of Basistha et al. [13] and Singh and Mal [42]. Thereafter, Lag-1 autocorrelation technique (alpha ≤ 0.05) was employed to check the presence of positive and negative autocorrelation in rainfall data, since this may affect the observed magnitude of trends in a dataset [13,[43][44][45]. As observed rainfall records in a high mountainous region like IHR are poorly available and remain a challenge for the scientific community, several descriptive statistics were then applied to check data distribution and relative performance of the gridded datasets to ensure high-quality satellite rainfall data with respect to station observed data be used for further investigation. This statistical metrics includes determination of (a) mean-variance between observed and satellite/model rainfall values (Bias), (b) multiplicative bias (MBias) to measure the mean magnitude of bias, (c) relative bias (RBias) to reveal systematic standard error (d) mean absolute error (MAE) to obtain the average slope of error (e) root mean square error (RMSE) to assess the average magnitude of error, and (f) correlation coefficient (CC) to evaluate the degree of relationship between observed and satellite data, (g) coefficient of variation used as a estimation of variability when the average is proportional to standard deviation, (h) uncertainty to measure the amount of error in an average of large scale estimation, together with an (i) empirical cumulative distribution frequency (CDF) and (j) the Taylor diagram [22,33,34,46] (Appendix A). A normalised standard deviation method was applied to the Taylor diagram to exhibit relative performance of both gridded and station datasets [47]. Furthermore, Sen's slope estimator was used to calculate the long-term magnitude of change in station observed data, the Mann-Kendall (MK) test was applied using MATLAB programming language to evaluate their statistical significance (Appendix B).

Spatio-Temporal Distribution of Rainfall and Rainy Days in the Bhilangana River Basin (1983-2008)
Mean annual rainfall  in the Bhilangana river basin is observed as 1196 mm across an average of 66 rainy days per year, although with marked spatial variability and seasonality. Rainfall is highest in the mountainous northeast (1473.6 mm), and somewhat less (710.8 mm) in the lower elevation areas in the southwest ( Figure 2). Maximum rainfall is recorded during the monsoon season (78.6%), along with an average of 44 rainy days, followed by the winter (10.1%), pre-monsoon (8.4%), and post-monsoon (3.0%) along with around 8, 11, and 3 rainy days correspondingly. Within the Indian Summer Monsoon (ISM) monsoon season, August exhibits maximum rainfall (28.5%) followed by July (24.7%), September (16.8%) and June (10.5%) with 15, 13, 8, and 7 rainy days respectively. Station observed data on the other hand revealed mean annual rainfall of the basin is 1313.4. Monsoon is the highest rainfall season followed by pre-monsoon, winter, and post-monsoon ( Table 2). Out of the two gridded datasets, CHIRPS reveals similar spatial information with station observed data in compare to PERSIANN-CDR. Mean annual rainfall revealed by the CHIRPS was 1119.7 mm, where PERSIANN-CDR underestimate rainfall counts by 998.1 mm. Seasonal and monthly distribution highlights similar spatial integrity with observed and CHIRPS dataset, where PERSIANN-CDR overestimate rainfall counts during pre-monsoon, post-monsoon and winter season (Figure 3a,b, Table 3). In both station and gridded datasets, there is greater annual and seasonal spatial variability at lower elevation stations such as Tehri and Narendra Nagar, than the higher elevation observatory (Mukhim) (Figure 3c,d). Both gridded datasets indicate similar spatial variability, although CHIRPS has better match with all the observed stations at seasonal and monthly scale (Figure 3c,d).
PERSIANN-CDR overestimate rainfall counts during pre-monsoon, post-monsoon and winter season (Figure 3a,b, Table 3). In both station and gridded datasets, there is greater annual and seasonal spatial variability at lower elevation stations such as Tehri and Narendra Nagar, than the higher elevation observatory (Mukhim) (Figure 3c,d). Both gridded datasets indicate similar spatial variability, although CHIRPS has better match with all the observed stations at seasonal and monthly scale (Figure 3c,d).   PERSIANN-CDR overestimate rainfall counts during pre-monsoon, post-monsoon and winter season (Figure 3a,b, Table 3). In both station and gridded datasets, there is greater annual and seasonal spatial variability at lower elevation stations such as Tehri and Narendra Nagar, than the higher elevation observatory (Mukhim) (Figure 3c,d). Both gridded datasets indicate similar spatial variability, although CHIRPS has better match with all the observed stations at seasonal and monthly scale (Figure 3c,d).

Detailed Outcomes of Descriptive Statistics
The autocorrelation assessment indicates that no significant homogeneity between both observed and gridded data sources (alpha ≤ 0.05). According to the results of descriptive statistics, the PERSIANN-CDR, for example, overestimates rainfall in the higher elevation region (Mukhim) with a positive MBias, but underestimates it in the low elevation part of catchment (higher spatial bias and lower CC (Table 3). Figure 4a-c illustrates linear regression and CC for both gridded and individual observational data, revealing strong spatial agreement with CHIRPS and somewhat less correspondence in the case of PERSIANN-CDR. Both the datasets were well matched with lower elevation station data, although CHIRPS has better agreement (r = 0.85) than the PERSIANN-CDR (r = 0.74) with an average lower uncertainty value (±6.2 mm). Higher elevation surface observatories have a similar spatial relationship with gridded datasets. Statistical indices indicated better performance of CHIRPS data in corresponding with direct observation, showing less bias at Tehri (4.59 mm) and Nagendra Nagar (6.26 mm), and a little higher, albeit not beyond the acceptance level, in Mukhim (bias acceptance level = 10 mm) (Figure 4d,e). The positive bias for each station indicates that CHIRPS sometimes underestimates rainfall amount. However, MBias and RBias values are low suggesting that the existing error is not systematic. PERSIANN-CDR, on the other hand, overestimates rainfall with both greater bias and lower spatial agreement (CC), therefore signifying systematic errors (Table 3), exception for Tehri. The maximum RMSE was detected by the PERSIANN-CDR in the high elevation region, i.e., Mukhim (111.30 mm), followed by Nagendra Nagar (96.2 mm) and Tehri (81.24 mm) whereas CHIRPS shows better performance with all the observed stations. Moreover, empirical CDF reveals spatial distribution for mean monthly rainfall ( Figure 4f); notably, more than 80% of the data are evenly distributed among CHIRPS and station records, where PERSIANN-CDR show a relatively higher distribution in respect to observed data due to greater RMSE and systematic error. Data (1983Data ( -2008 This study reveals annual and seasonal rainfall distribution and trend for the period from 1983 to 2008. Maximum annual rainfall occurred in highest elevation surface observatory (Mukhim) followed by Narendra Nagar and Tehri. Seasonally, Narendra Nagar received maximum rainfall than the Mukhim and Tehri ( Table 4). The analysis also reveals that the three surface observatories experienced decreasing rainfall trends, some of which are statistically significant. Narendra Nagar showed the maximum decrease (−32.82 mm/decade) followed by Tehri (−12.26 mm/decade) and Mukhim (−2.04 mm/decade) between 1983 and 2008 (Table 4). Stations at a lower elevation (Tehri) witnessed rainfall changes of lower magnitude (statistically insignificant). All surface observatories experienced a decreasing trend in rainfall throughout the year, except Mukhim during the winter season. The most marked rainfall changes occurred during the monsoon season (−11.23 mm/decade) followed by winter (−1.1 mm/decade), pre-monsoon season (0.63 mm) and post-monsoon (0.53 mm) seasons. Narendra Nagar recorded the maximum decrease in monsoon season rainfall (−22.37 mm/decade) and also Remote Sens. 2020, 12, 709 8 of 24 for the winter (−2.56 mm/decade) followed by post-monsoon (−0.82 mm/decade) and pre-monsoon (−0.676) rainfall. The station at Tehri also witnessed significant rainfall decline during monsoon (−7.45 mm/decade), although the reductions during winter (−0.71 mm/decade) and pre-monsoon (−0.31 mm/decade) were not statistically significant. This was also the case, for Mukhim in the monsoon (3.79 mm/decade) and winter (−0.009 mm/decade); the negative trend (0.91 mm/decade) in the post-monsoon season was, however, statistically significant at this locality.

Annual and Seasonal Trends in Observatory Rainfall
Remote Sens. 2020, 12, x FOR PEER REVIEW 8 of 25 This study reveals annual and seasonal rainfall distribution and trend for the period from 1983 to 2008. Maximum annual rainfall occurred in highest elevation surface observatory (Mukhim) followed by Narendra Nagar and Tehri. Seasonally, Narendra Nagar received maximum rainfall than the Mukhim and Tehri ( Table 4). The analysis also reveals that the three surface observatories experienced decreasing rainfall trends, some of which are statistically significant. Narendra Nagar showed the maximum decrease (−32.82 mm/decade) followed by Tehri (−12.26 mm/decade) and Mukhim (−2.04 mm/decade) between 1983 and 2008 (Table 4). Stations at a lower elevation (Tehri) witnessed rainfall changes of lower magnitude (statistically insignificant). All surface observatories experienced a decreasing trend in rainfall throughout the year, except Mukhim during the winter season. The most marked rainfall changes occurred during the monsoon season (−11.23 mm/decade) followed by winter (−1.1 mm/decade), pre-monsoon season (0.63 mm) and post-monsoon (0.53 mm) seasons. Narendra Nagar recorded the maximum decrease in monsoon season rainfall (−22.37 mm/decade) and also for the winter (−2.56 mm/decade) followed by post-monsoon (−0.82 mm/decade) and pre-monsoon (−0.676) rainfall. The station at Tehri also witnessed significant rainfall decline during monsoon (−7.45 mm/decade), although the reductions during winter (−0.71 mm/decade) and pre-monsoon (−0.31 mm/decade) were not statistically significant. This was also the case, for Mukhim in the monsoon (3.79 mm/decade) and winter (−0.009 mm/decade); the negative trend (0.91 mm/decade) in the post-monsoon season was, however, statistically significant at this locality.

Monthly Trends in Observatory Data (1983-2008)
Analysis of monthly distribution and trend of rainfall offers a more detailed picture of rainfall change in the region over the study period. Maximum rainfall was observed in July, August, and September in all the surface observatories. However, Narendra Nagar experienced maximum rainfall in the last 26 years followed by Mukhim and Tehri (Table 5). Furthermore, statistical trend analysis shows marked decreases in Narendra Nagar for July (−8.4 mm/decade), August (−12.4 mm/decade) and October (−0.5 mm/decade) ( Table 5). Tehri also recorded significant decreases but mainly in July (0.8 mm/decade) and September (0.1 mm/decade). Monthly values for Mukhim show no statistically significant trends, although rainfall for June, July, October, and November did decrease over the study period. Annually, all surface observatories witnessed a statistically significant reduction in the number of rainy days ( Table 6). Narendra Nagar recorded the most marked frequency reduction during the monsoon season, followed by Mukhim and Tehri (Table 6). Mukhim and Tehri stations experienced fewer rainy days in post-monsoon and winter seasons respectively. The negative trend in rainy days during monsoon season indicates a declining rainfall pattern for June through September. Significant changes in annual and seasonal rainy days have direct influence on long-term rainfall patterns in the study area. Statistically significant declining annual rainfall trends were recorded in both gridded data products, averaging 10.9 mm/decade (Table 7). CHIRPS exhibited the most prominent changes (−14.8 mm/decade), followed by PERSIANN-CDR (−6.9 mm/decade). Modelled rainfall declined in all seasons across the study period, with CHIRPS revealing a greater magnitude of change during monsoon, pre-monsoon, winter, and post-monsoon (significant) seasons compared with PERSIANN-CDR (not significant). The magnitude of rainfall reduction revealed by PERSIANN-CDR is lower, due to higher spatial bias and greater mean absolute error. Due to the underestimation of rainfall counts (beyond the acceptance level), PERSIANN-CDR showed very fewer changes in seasonal and monthly values in comparison to observed station datasets, whereas CHIRPS has better agreement over the study period exhibited by the higher correlation coefficient (r 2 = 0.74). Spatially, the statistically significant reduced rainfall trend is confined to the south-western parts of the study area ( Figure 5). The decreased rainfall trend in the north-western windward parts of the area over the area is not statistically significant. Low elevation regions exhibit less rainfall variability throughout the decades, because both models (CHIRPS and PERSIANN-CDR) underestimate rainfall amounts (bias = 4.6 and 16.88 mm respectively). counts (beyond the acceptance level), PERSIANN-CDR showed very fewer changes in seasonal and monthly values in comparison to observed station datasets, whereas CHIRPS has better agreement over the study period exhibited by the higher correlation coefficient (r 2 = 0.74). Spatially, the statistically significant reduced rainfall trend is confined to the south-western parts of the study area ( Figure 5). The decreased rainfall trend in the north-western windward parts of the area over the area is not statistically significant. Low elevation regions exhibit less rainfall variability throughout the decades, because both models (CHIRPS and PERSIANN-CDR) underestimate rainfall amounts (bias = 4.6 and 16.88 mm respectively).   Lower elevation parts of the basin record greater magnitude changes both annually and seasonally supported by lower systematic errors and bias with high spatial integrity. CHIRPS data also indicated more change in low-lying areas with higher spatial agreement between station observed data, whereas PERSIANN-CDR show no significant trend due to its underestimation of rainfall counts (Bias = 16.88). The eastern side of the basin showed maximum changes in both the datasets, which suggests lower monsoonal rainfall ( Figure 6). The north and north-western areas exhibit lower winter rainfall changes from 1983 to 2008.

Monthly Trends: CHIRPS and PERSIANN-CDR (1983-2008)
Markedly reduced rainfall is exhibited by CHIRPS for the months May through November (Table 8). In contrast, the PERSIANN-CDR output showed a (statistically insignificant) slightly increasing trend for May, June, and October. Maximum changes are recorded in July and August (CHIRPS), whereas PERSIANN-CDR data reflect a lower magnitude of change over all due to higher spatial bias and systematic errors compared to observe rainfall. Lower elevation parts of the basin record greater magnitude changes both annually and seasonally supported by lower systematic errors and bias with high spatial integrity. CHIRPS data also indicated more change in low-lying areas with higher spatial agreement between station observed data, whereas PERSIANN-CDR show no significant trend due to its underestimation of rainfall counts (Bias = 16.88). The eastern side of the basin showed maximum changes in both the datasets, which suggests lower monsoonal rainfall ( Figure 6). The north and north-western areas exhibit lower winter rainfall changes from 1983 to 2008.

Monthly Trends: CHIRPS and PERSIANN-CDR (1983-2008)
Markedly reduced rainfall is exhibited by CHIRPS for the months May through November (Table 8). In contrast, the PERSIANN-CDR output showed a (statistically insignificant) slightly increasing trend for May, June, and October. Maximum changes are recorded in July and August (CHIRPS), whereas PERSIANN-CDR data reflect a lower magnitude of change over all due to higher spatial bias and systematic errors compared to observe rainfall. Bold numbers are statistically significant (alpha ≤ 0.05).

Inter-Annual Rainfall Variability and Trends (1983 to 2018) from Station Observed and CHIRPS Data
Inter-annual rainfall variability was tested using station observed and gridded data sets for the period from 1983 to 2008. Due to the unavailability or unreliability of ground station records for the basin in the post-2008 period, annual, seasonal, and monthly rainfall trends and variability need to be assessed using modelled data. The study suggests an overall increasing tendency for the period from1983 to 2018 of approximately 8.76 mm/decade, although there is substantial inter-annual variability across the entire study period. CHIRPS is shown here to produces values that are more consistent with observed data, and we selected this to estimate long-term inter-annual rainfall variability and trends as illustrated in Figure 7a indicate three distinct phases are evident in the annual, seasonal and monthly records, such that 1983 to 1997 (phase I) exhibits mean annual

Inter-Annual Rainfall Variability and Trends (1983 to 2018) from Station Observed and CHIRPS Data
Inter-annual rainfall variability was tested using station observed and gridded data sets for the period from 1983 to 2008. Due to the unavailability or unreliability of ground station records for the basin in the post-2008 period, annual, seasonal, and monthly rainfall trends and variability need to be assessed using modelled data. The study suggests an overall increasing tendency for the period from1983 to 2018 of approximately 8.76 mm/decade, although there is substantial inter-annual variability across the entire study period. CHIRPS is shown here to produces values that are more consistent with observed data, and we selected this to estimate long-term inter-annual rainfall variability and trends as illustrated in Figure 7a indicate three distinct phases are evident in the annual, seasonal and monthly records, such that 1983 to 1997 (phase I) exhibits mean annual precipitation fluctuating around 1270 mm with lower uncertainty (±18 mm); phase II (1998-2008) is associated with significantly lower annual totals (less than 620 mm), strongly suggested by station observed data (Figure 7d,e); phase III indicates sharply increased annual rainfall of around 1586mm with high uncertainty values. More recently, there is a suggestion of an increasing trend (4.11 mm/decade), although this is not statistically significant. This investigation further revealed that monthly and seasonal rainfall has increased in recent years (2009-2018) in comparison to previous decades, whereas phase II experienced severely less rainfall and the significant reduction rate was −4.24 mm/decade (Figure 7a). Bold numbers are statistically significant (alpha ≤ 0.05).  This recent (phase III) increase is primarily due to higher rainfall during the pre-monsoon season in particular (Table 9; Figure 7b). However, there are significant decreasing changes in rainfall during monsoon, pre-monsoon, and winter season in the north-western part of the basin in particular. The eastern and south-eastern parts on the other hand experienced increasing rainfall during these same seasons (not statistically significant) (Figure 8). Analysis of monthly trends indicate statistically insignificant decreasing rainfall with lower uncertainty values in September (−8.2 mm/decade) followed by February (−2.1mm/decade) and October (−1.9 mm/decade) with higher uncertainty. However, January exhibits a significant increasing trend (0.79 mm/decade) with relatively greater uncertainty (±11.4 mm), while March to July recorded an increasing but not significant trend in rainfall with an average of ± 2.9 mm uncertainty (Table 9). This study also highlights that the highest increases are for June (7.5 mm/decade) and July (6.1 mm/decade) (Figure 7c, Table 9). This increasing tendency in monthly rainfall is most likely to be a result of extreme rainfall events in the study area.

Discussion
This study presents an analysis of annual, seasonal, and monthly trends and variability in rainfall and rainy days using three surface observatories and two high-resolution gridded datasets at different elevation ranges (from 6635 to 616 m amsl) over the last four decades. These three datasets were calibrated and compared using several non-parametric statistical techniques which reveals significance levels (alpha ≤ 0.05) for a wide range of data characteristics. The present study reveals that maximum annual rainfall is observed during the monsoon season (79%) followed by the winter, pre-monsoon, and post-monsoon, which is corroborated by the studies in the adjacent regions of Uttarakhand and Nepal Himalaya, during the presence of ISM (monsoon) and Westerly Disturbances (winter) [13,14,17].

Temporal Trends in Observed Rainfall (1983-2008)
Significance changes in annual, seasonal, and monthly rainfall are observed using station-based records. Spatial variability in rainfall trends and orographic effects on the distribution of rainfall are apparent in the study area. Narendra Nagar witnessed the significant reductions in annual and monsoonal rainfall (−32.82 and −22.37 mm/decade) followed by Tehri and Mukhim. The trends align well with the previous studies in the western Himalaya and in the abject Nepal Himalaya [13,42]. In such circumstances, topography, relief, and aspects are dominant factors influencing rainfall distribution, for instance, Tehri is located on the western side of the Gangotri range and receives more rainfall from frontal systems in the winter compared to Narendra Nagar [48,49].
Significant inter-seasonal variations are also apparent, whereby average monsoon rainfall changes (−11.23 mm-decade) are most considerable across the basin. Similar outcomes have been reported by others, for example Ramesh and Goswami [15] described the increased occurrence of "monsoon-breaks" across India between 1951 and 2008, which is attributed largely to less atmospheric water vapour reaching the subcontinent. Shekhar et al. [5] and Kumar and Jain [16] highlighted decline in frequency and intensity of westerly disturbances in winters along with shrinkage of monsoonal cloud cover in the Himalayan region from 1971 to 2000. A growing body of work concerning rainfall distribution and trends in the IHR describes declining thermal contrast between the Indian Ocean and Tibetan Plateau, resulting in less moisture reaching the Indian landmass, and contributes to decreased rainfall in higher elevation regions [16,17,50]. Besides, anthropogenic factors such as large scale land use transformation from forest to agricultural land and settlement in low elevation regions under the 'Tehri Dam Rehabilitation Program' [51,52] and the 'Terai Colonization Scheme' also play a dominant role in altering local rainfall patterns [53].
Significant reduction in the annual frequency of rainy days was identified from observational data in the Bhilangana river basin. Several other studies on rainfall in India and Himalayan region have reached similar conclusions [13,16,17,54,55]. This temporal pattern is perhaps due to a combination of factors, including substantial decrease in the frequency and intensity of westerly disturbances from 1977 to 2007 [16], population migration from high to low elevation regions with associated land use change and weakening of the land-ocean thermal contrast between the Tibetan plateau and Indian Ocean [50]. In general, large scale deforestation [56], global climate shift [57], dynamics of the El Niño Southern Oscillation [58] and North Atlantic Oscillation [10] are also the possible drivers of the rainfall trend.

Appropriateness of Gridded Datasets in Compared to Observed Data (1983-2008)
There remains controversy, however, around long-term climate changes in the Himalaya. This is largely a consequence of the complexity of monsoon system itself, coupled with heterogeneity of surface topography [59]. Long-term observation stations are few and far between due to the inhospitable climate and inaccessibility, and remotely sensed and reanalysis products must be relied upon to estimate uncertainties that cannot be addressed with station rainfall data. Using the facility of GEE cloud computing, reanalysis (CHIRPS) and satellite-derived (PERSIANN-CDR) rainfall products may be compared and statistically analysed in terms of relative performance against the limited observational data.
Satellite products do have limitation in regard to estimate rainfall in mountainous regions due to the presence of cold cloud, ice, and snow [32,33,59,60]. Moreover, the role of local topography in orographic rainfall is widely described [61][62][63][64] and this certainly impacts on the resolution and accuracy of satellite-derived rainfall data. Prakash [34] and Beck et al. [35] demonstrated the suitability of CHIRPS rainfall data through several statistical methods and concluded that the product adequately estimates annual and seasonal long-term rainfall trends with respect to observed data. Which seems in general to underestimate changes relative to the observational data. Ullash et al. [33] noted the proficiency of CHIRPS data to estimate seasonal and monthly rainfall counts in different climate zones of Pakistan, with low bias.
However, Ashouri et al. [38] established the proficiency of PERSIANN-CDR in terms of long-term rainfall-runoff modelling with low bias. In recent years, there has been growing literature in the performance of PERSIANN-CDR daily rainfall product to acquire rainfall information in Asian mountainous regions and assess its efficiency in terms of long-term performance for rainfall and hydrological modelling [32,37,[65][66][67][68][69].
Our study further highlights the performance of these two gridded data products whereby CHIRPS a greater degree of correspondence with observed rainfall compared to PERSIANN-CDR. The monthly, annual, and seasonal distribution from 1983 to 2008 depicts similar temporal trends and spatial patterns to the three surface observatories (Figure 9). Statistical analysis supports this interpretation of rainfall, in that CHIRPS shows lower systematic bias and spatial distributional errors (RMSE). Annual, seasonal, and monthly trends from 1983 to 2008 reflect the ground station data well, whereas the trends in PERSIANN-CDR data are variable.
Remote Sens. 2020, 12, x FOR PEER REVIEW 15 of 25 (RMSE). Annual, seasonal, and monthly trends from 1983 to 2008 reflect the ground station data well, whereas the trends in PERSIANN-CDR data are variable. Figure 9. Taylor diagram compares observed and gridded datasets (Notations, "Jan to Dec" are used for defining months, "Ann to Win" for annual and seasons, "obs" for the observation stations and "C" and "P" denote CHIRPS and PERSIANN-CDR respectively).

Estimation of 1983 to 2008 Rainfall with Gridded Data
Observed rainfall variability based on the gridded data presented here exhibit similar temporal trends to station data, with slight differences in annual, seasonal, and monthly rainfall, with a mean annual reduction of around 11 mm/decade across the basin. CHIRPS with a greater degree of agreement with station data, shows reduction in annual and seasonal rainfall by comparison with PERSIANN-CDR., because this product appears to systematically underestimate rainfall counts in the region. The study shows that most of the low-lying regions in the Bhilangana river basin have experienced a significantly decreasing trend in both annual and monsoon season rainfall, whereas PERSIANN-CDR showed much less change. Similar results are noted by Ullah et al. [33] in the mountainous regions of Karakoram and found CHIRPS data has better performance than PERSIANN-CDR to estimate rainfall counts. Taylor diagram compares observed and gridded datasets (Notations, "Jan to Dec" are used for defining months, "Ann to Win" for annual and seasons, "obs" for the observation stations and "C" and "P" denote CHIRPS and PERSIANN-CDR respectively).

Estimation of 1983 to 2008 Rainfall with Gridded Data
Observed rainfall variability based on the gridded data presented here exhibit similar temporal trends to station data, with slight differences in annual, seasonal, and monthly rainfall, with a mean annual reduction of around 11 mm/decade across the basin. CHIRPS with a greater degree of agreement with station data, shows reduction in annual and seasonal rainfall by comparison with PERSIANN-CDR., because this product appears to systematically underestimate rainfall counts in the region. The study shows that most of the low-lying regions in the Bhilangana river basin have experienced a significantly decreasing trend in both annual and monsoon season rainfall, whereas PERSIANN-CDR showed much less change. Similar results are noted by Ullah et al. [33] in the mountainous regions of Karakoram and found CHIRPS data has better performance than PERSIANN-CDR to estimate rainfall counts.

Estimation of 2009-2018 Rainfall with CHIRPS
The CHIRPS analysis of rainfall pattern from 2009 to 2018 indicates a tendency towards increasing annual rainfall, although this is not statistically significant, however, other studies have suggested increased rainfall over the western Himalaya after 2008 due to fluctuation in monsoonal rainfall, the formation of deep and wide convective cores in the Arabian seas and release at the foothills of western Himalaya, might be the cause of substantial rainfall at the south-eastern area due to orographic uplifting [70,71].
An increased frequency of extreme events in the recent past, including cloudbursts, lake avulsion floods and avalanches have been recorded for the Bhagirathi and Bhilangana basins, more especially in the monsoon season [18,29]. Although Dimri et al. [29] indicate that gridded data products may not capture such short term events in high mountain regions, CHIRPS data in this study do successfully record extreme rainfall events during the monsoon seasons of July and August 2012 and 2013.

Variability of Rainfall (1983-2018) from CHIRPS Data
Temporal variability in annual and seasonal rainfall was analysed using the CHIRPS dataset for the period from 1983 to 2018. Based on our analysis, we identified three distinctive rainfall phases (1983 to 1997, 1998 to 2008, and 2009 to 2018). This pattern is notably due to a particularly dry period, although this has not previously been shown for this or other parts of the Himalaya. However, Bhutiyani et al. [11] conducted a spatial-temporal analysis of teleconnections between two meteorological phenomena (NAO, ENSO) and monsoon rainfall and found that these systems brought less moisture to the upland areas between 1983 and 2005. They further suggested that reduced rainfall during 2004 and 2005 was due to rapidly increasing temperatures in the mountainous areas of the Himalaya. Archer and Fowler [48] highlighted the importance of local topography in influencing the amount and distribution of moisture during the onset of the ISM.

Regional Environmental Impacts of Rainfall Changes
Rainfall change is certain to have significant environmental and socio-economic consequences, especially in vulnerable high mountain regions. Aside from the historical changes highlighted in this paper, has been predicted that future surface temperatures will increase [72] and that surface runoff of Himalayan major river basins will increase due to melting glaciers and early melting of snow [73]. Higher temperatures impact rainfall mechanisms, which will, in turn, affect spatial and temporal patterns of rainfall in the study region. Such scenarios are potentially very important for people living in mountainous areas. Appollo, [74] stated that in the last few decades (1961-2011), the Himalayan population has increased by 2.5 times from 19.9 to 52.8 million with a growth rate of 3.3 percent, which was three times the global average. The temporal patterns of rainfall reported here suggests that the population has already faced considerable challenges in this regard and will be increasingly exposed. Almost 70 percent of the population of the basin is dependent on agriculture and small-scale industries for their livelihoods and is therefore entirely reliant on the river water supply of the Bhilangana river [24].

Effects of Decreasing Rainfall: 1983 to 2008
The long-term trend analysis of the Bhilangana river basin has revealed a decreasing rainfall pattern, averaging roughly 15 mm/decade using both stations observed and gridded datasets. Of the total area of the Bhilangana community development block (21350.67 hectares), more than 40 percent of the land is cultivable and, within this, more than one-quarter is irrigated [24]. This changing rainfall pattern presents a major challenge to the people of the study area, since less water availability impacts agriculture and other related activities directly. Several other studies have noted that rainfall in the Himalayan mountains has become more strongly seasonal, whereby the pre-and post-monsoon seasons exhibit prominent drying trends [13,42]. The negative trends in the pre-monsoon rainfall have significant implications for the economy and local livelihoods. Since the pre-monsoon season marks the beginning of the growing season, when the demands for water for irrigation is high. Decreasing rainfall has led to abandonment of agricultural lands (Figure 10a). Moreover, small-scale household industries, such as small hydropower stations have been discontinued due to lack of water in the streams and rivers (Figure 10b). movement events that affected the area between 800 to 1800 m amsl and resulted in losses of livestock and destruction of houses during the monsoon months. Banerjee et al. [28] reported that almost 10 percent of the total basin area falls within the high landslide risk zone and increased rainfall further heightens the risk of such disasters. Dramatic evidence of the impact of such events changes can be seen in the damage to buildings and infrastructure (Figure 10d). During a flood event in 2012 flood, a medium-scale (24 mega-watts) hydropower station was partially destroyed (Figure 10c).

Conclusions
The prime aim of this study was to evaluate the spatio-temporal dynamics of rainfall and rainy days for the Bhilangana river basin during the period from 1983 to 2018 using both in situ and gridded rainfall data. Non-parametric statistical indices were applied coupled with advance web-based cloud computing platform (GEE). This basin was studied because of its diverse topography, large extent of population dynamics, heterogeneous land use and limited availability of ground data. The present study adequately uses big data applications using GEE platform, alongside MATLAB programming language to delineate the monotonic trend of gridded and observed data. The main conclusions of the study can be summarised as follows. Recent trends in inter-annual rainfall patterns suggest the long-term reduction in rainfall has been reversed and, while this is not statistically significant at the annual scale, an increase in the frequency and magnitude of extreme events has become an additional threat for local people as well as infrastructure. Dimri et al. [29] describes several recent major cloudburst and associated mass movement events that affected the area between 800 to 1800 m amsl and resulted in losses of livestock and destruction of houses during the monsoon months. Banerjee et al. [28] reported that almost 10 percent of the total basin area falls within the high landslide risk zone and increased rainfall further heightens the risk of such disasters. Dramatic evidence of the impact of such events changes can be seen in the damage to buildings and infrastructure (Figure 10d). During a flood event in 2012 flood, a medium-scale (24 mega-watts) hydropower station was partially destroyed (Figure 10c).

Conclusions
The prime aim of this study was to evaluate the spatio-temporal dynamics of rainfall and rainy days for the Bhilangana river basin during the period from 1983 to 2018 using both in situ and gridded rainfall data. Non-parametric statistical indices were applied coupled with advance web-based cloud computing platform (GEE). This basin was studied because of its diverse topography, large extent of population dynamics, heterogeneous land use and limited availability of ground data. The present study adequately uses big data applications using GEE platform, alongside MATLAB programming language to delineate the monotonic trend of gridded and observed data. The main conclusions of the study can be summarised as follows.
1) This study provided a comprehensive understanding of the long-term historical trend of observed rainfall and rainy days, together with satellite and reanalysis data sets. The statistical indices indicated CHIRPS data show relatively better performance and accuracy than PERSIANN-CDR, when compared to observational data. However, associating the two results, it can be seen that both station and gridded data products record statistically significant decreasing trends in rainfall from 1983 to 2018. CHIRPS showed similar tendency corresponding with station data, averaging 11 mm rainfall/decade over the period of record.
Notably, The PERSIANN-CDR showed less alteration in annual rainfall (6.9 mm/decade) due to high RMSE and systematic bias, resulted in miscounting rainfall totals. From the observational data exposed, Narendra Nagar has exclusively witnessed very substantial reduced annual rainfall (32.82 mm/decade). Seasonally there is a statistically significant decreasing trend in monsoon, winter, and post-monsoon rainfall of all the data products. Areas at lower elevation were subject to a greater degree of rainfall decrease while the magnitude of change was highest in the monsoon season.
2) The temporal pattern of rainfall change can bring about substantial impacts on the regional economy, including reduced water supply for crop cultivation, hydropower generation, tourism, and livelihoods in general. In addition, the increasing recent rainfall pattern from 2009 to 2018, has caused loss of life and damage of property and infrastructure.
3) Understanding spatial and temporal patterns, trends, and variability of rainfall in the Himalaya remain a challenge to assess due to a combination of complex topography and low spatial coverage of observation stations and inhospitable climatic conditions. Eventually, the systematic bias of gridded rainfall in mountainous regions specify that precise estimation by satellite-based rainfall products remains a big task for the scientific community. Sometimes low spatial resolution of satellites is also another cause to detect short term changes at micro level regions. Nevertheless, high-resolution reanalysis data products (CHIRPS) can be analysed through powerful cloud computing applications such as GEE and provide valuable information about the dynamics of rainfall in vulnerable high mountain regions.

Appendix A.
Appendix A.1. Auto-Correlation Occurrence of positive and negative autocorrelation may affect the analysis of trend in dataset [13,43]. The autocorrelation coefficient was tasted using the interval of Lag-1 (at 95% significance level) and computed as follows −1 − 1.645 * (N − 2) 0.5

. Bias
Bias denotes mean variance between in situ and satellite or model data. The positive bias indicates underestimation and negative bias specifics overestimation of rainfall.
whereas, R S i and R O i depicts satellite and observed rainfall values respectively. N is showing number of observations.

Appendix A.3. Multiplicative Bias
Multiplicative Bias is the average magnitude of gridded precipitation over observed counts, where underestimation seems to count less than unity and overestimation is greater than unity [22]. This is actually estimate probability mass function to distinct kernels of a set of observations including variance and mean integrated standard error.
where R S i and R O i depict satellite and observed rainfall values respectively. N is the number of observation days over a particular grid.

Appendix A.4. Relative Bias
Relative bias illustrates the systematic bias of satellite precipitation and acts as parallel to bias. This actually estimates the standard error between the observed and modelled datasets.
where R S i and R O i depict satellite and observed rainfall values respectively. N is the number of observations.

Appendix A.5. Mean Absolute Error
Mean absolute error is employed to demonstrate average magnitude of error in satellite and observed rainfall values.
where R S i and R O i depict satellite and observed rainfall values respectively. N is the number of observations days over a particular grid.
Appendix A.6. Root Mean Square Error Root mean square error is similar to MAE. However, it gives higher priority to calculate the mean magnitude of error in satellite and in situ measurement [33]. where R S i and R O i depict satellite and observed rainfall values respectively. N is the number of observations.

Appendix A.7. Correlation Coefficient
The correlation coefficient is widely used to portrait spatial agreement between observed and gridded data sets. Perfect relationship specifies towards +1.

. Coefficient of Variation
The coefficient of variance is used to determine the degree of variability from one data series to another by simply the ratio of mean and standard deviation.
where R S i and R O i depict satellite and observed rainfall values respectively. N is the number of observations.