Probabilistic Precipitation Estimation with a Satellite Product

: Satellite-based precipitation products have been shown to represent precipitation well over Nepal at monthly resolution, compared to ground-based stations. Here, we extend our analysis to the daily and subdaily timescales, which are relevant for mapping the hazards caused by storms as well as drought. We compared the Tropical Rainfall Measuring Mission (TRMM) Multi-satellite Precipitation Analysis (TMPA) 3B42RT product with individual stations and with the gridded APHRODITE product to evaluate its ability to retrieve different precipitation intensities. We ﬁnd that 3B42RT, which is freely available in near real time, has reasonable correspondence with ground-based precipitation products on a daily timescale; rank correlation coefﬁcients approach 0.6, almost as high as the retrospectively calibrated TMPA 3B42 product. We also ﬁnd that higher-quality ground and satellite precipitation observations improve the correspondence between the two on the daily timescale, climatology and show that the quantitative precipitation estimates produced by this model are well calibrated compared to APHRODITE.


Introduction
Precipitation products based on remote sensing offer the potential for improving hazard response and water resource management in mountainous areas with inadequate near-real-time ground-based data [1].We focus on Nepal and its surroundings (26 • -31 • N, 79 • -89 • E), a region that encompasses the Himalaya range and its foothills in the north of the Indian subcontinent, has wide geographic and seasonal ranges of precipitation frequency and intensity [2,3], and whose population is largely agrarian and highly vulnerable to climate-related hazards, including flooding and drought [4,5].We have previously compared the performance of several remote sensing based precipitation products over Nepal relative to station observations on the monthly timescale, finding that the Tropical Rainfall Measuring Mission (TRMM) Multi-satellite Precipitation Analysis (TMPA) 3B43 precipitation product, which combines data from TRMM and other satellites with calibration to surface measurements, performed best and showed little bias [6].Here, we extend this work by (a) considering the reliability of satellite precipitation at daily and subdaily temporal resolutions, which are better suited than monthly resolution for most hydrological hazard assessment work; (b) comparing the performance of the near-real-time TMPA 3B42RT product to the research TMPA 3B42 product (the 3B43 product previously assessed is the monthly-resolution version of 3B42); (c) estimating rainfall rate probabilities conditional on the satellite data [7], which enables, for example, the mapping of areas that experienced heavy rainfall with high probability.We validate and calibrate the TMPA products against the Asian Precipitation-Highly Resolved Observational Data Integration Towards Evaluation of Water Resources (APHRODITE) gridded daily precipitation product, which is based on station precipitation observations, and against automatic weather stations that provide precipitation time series at high (sub-hourly) time resolutions.

Satellite-Based Precipitation
TMPA Product 3B42 is released by the Goddard Space Flight Center of the National Aeronautics and Space Administration (NASA) [8].Product 3B42 is a 3-hour precipitation product with a spatial resolution of 0.25 • , available for the period since 1998, which incorporates microwave and infrared observations from multiple satellites including TRMM [9].3B42 is a retrospective research product that is generally updated with a lag of several months because it uses monthly rain gauge accumulations in order to correct the estimated precipitation field.Product 3B42RT [10] has the same spatial and temporal resolution as 3B42 but is processed in near real time, available since March 2000, with updated fields typically posted about 3 hours after the end of each 3-hour window.Both products are available globally over 50 • S-50 • N.While 3B42RT currently does not come with quantitative error estimates, it does include source information that distinguishes between pixels where the precipitation estimate is based on microwave sounders (more accurate) or only on infrared imagery (less accurate).Here, we consider the 3B42RT product because its near real time availability makes it more suitable than 3B42 for most hydrological applications, but also compare the two products to assess the extent to which after-the-fact calibration improves the correspondence with gauge-based precipitation.The product versions used were the most recent available as of early 2015, i.e., 3B42 Version 7, which has generally been found to improve on the older Version 6 [11][12][13].Several studies have previously evaluated 3B42RT against ground measurements and other satellite-based products in different regions [14][15][16][17][18][19][20][21].

Gauge-Based Precipitation
Our primary "ground-truth" dataset in this study is from the Asian Precipitation-Highly Resolved Observational Data Integration Towards Evaluation of Water Resources (APHRODITE) project, which has produced a daily precipitation dataset over Asia based on an extensive compilation of rain gauge measurements [22], primarily station data obtained from national meteorological agencies, including Nepal's Department of Hydrology and Meteorology (DHM).Station data were interpolated using an algorithm that takes topography into account and uses climatology to estimate missing values.We used the APHRODITE V1101 monsoon Asia product at the best available resolution of 0.25 • , which is available for 1951-2007.APHRODITE includes information on the percent of 0.05 • cells in each pixel containing gauge data.Unless otherwise specified, we employed only 0.25 • pixels for which this quantity was positive, i.e., which are directly based on at least some gauge measurements, in calibrating and comparing with the satellite products.The satellite products were aggregated to daily resolution for comparison with APHRODITE.
Note that one possible contributor for inconsistency between these data sets is that the days in the aggregated satellite products follow Universal Time (Greenwich Mean Time) while the reading times for the daily values recorded from gauges and used as the basis for the APHRODITE interpolation may vary [22].However, in preliminary testing we found that shifting the satellite product aggregation window to allow for time zone and gauge reading time differences did not greatly change the correlations with APHRODITE daily values.
To characterize the performance of the TMPA satellite products in reproducing the precipitation diurnal cycle, we also compared them with precipitation measured at high temporal resolution at 3 automatic weather stations (AWSs) in Nepal installed by our team.These stations are located in farm fields representing different agro-ecological zones at (a) Baireni, Valtar in Dhading district, 27 83"E (115 m).These stations record accumulated precipitation at 5 minute intervals and temperature and relative humidity at 15 minute intervals.The rain gauges use tipping buckets to measure precipitation at 0.2 mm resolution; the gauge models are TR-525S (manufactured by Texas Electronics, Inc., Dallas, Texas, USA) at the Dhading site and TB3 (Hydrological Services, Ltd., Sydney, NSW, Australia) at the other two sites.The gauges are not heated, since the site elevations are low enough that they do not experience snow.Data were available beginning June 2013 and ending September 2014 (Syangja and Kapilvastu) or December 2014 (Dhading).

General Approach
Probabilistic precipitation forecasts are common in the short-term synoptic context [23,24] and also for long-term seasonal prediction [25][26][27][28].An ensemble of precipitation fields consistent with observations has been used to represent the uncertainty of radar-based precipitation estimates [29,30].The uncertainty represented by such probability distributions or ensembles can be propagated to streamflow and water supply quantities using hydrological models [31].Bellerby and Sun [32] is one study that adopted a probabilistic approach to precipitation retrieval from remote sensing, fitting gamma distributions to represent precipitation intensity probabilities at different measured cloud-top infrared brightness temperatures.Authors in [33] used a nonparametric (kernel density) approach to estimate 3-hourly precipitation occurrence and intensity probabilities over part of the United States conditional on retrievals from the satellite precipitation product CMORPH, but noted that sampling the resulting nonparametric probability distributions is not straightforward.Here, we model the probability distribution p(x) of APHRODITE precipitation amount x by bringing together separate models for precipitation occurrence (p(x > 0)) and precipitation intensity conditional on occurrence (p(x|x > 0)).(An alternative approach would be to subsume both in one distribution that has a point mass at zero precipitation amount [34,35].)We use flexible parametric functional forms to facilitate fitting the conditional probability distributions to the observed precipitation amount data.

Precipitation Amount Modeling
For precipitation intensity, we assume a generalized linear model [36] where N denotes the normal distribution, Z is an 1 × l vector of predictor values, which can include climatological (observed spatial and seasonal variability), and 3B42RT terms, and β is a k × 1 vector of coefficients.Values for β and the model standard deviation σ are determined by linear regression using the data from 2001-2007.Given this training data, standard statistical results show that g(p(x|x > 0)) for unobserved data with a known set of predictor values is given by a t distribution [37].Note that, for simplicity and to make use of all the available valid data, β, σ were fitted over the entire study domain and not for each grid point separately.We compared 4 models, which differ in the predictor values considered: (a) no predictors (baseline probability distribution, which is the same for all times and places); (b) predictors based on the 3B42RT precipitation amount; (c) predictors based on spatial location and Fourier components of the seasonal cycle (climatology); (d) both 3B42RT and climatology predictors.For each of models (b)-(d), which predictors to retain out of a candidate set was decided using stepwise linear regression [38] as implemented in the stepwisefit function of the Statistics package in GNU Octave [39].
The transformation g(p), which serves as the link function in the generalized linear model, is based on approximating the APHRODITE precipitation intensity probability density as a sum of exponential distributions, also known as a hyperexponential distribution: with the component distributions arranged in order of scale, so that 0 < b 1 < b 2 < . . .< b N , and additional constraints a i > 0, N i=1 a i b i = 1.The coefficient values are chosen to maximize the likelihood of observed APHRODITE precipitation over 1981-2000, and the number of components N is chosen by applying the corrected Akaike information criterion [40,41] with the same data.(See Appendix for more details.) A variety of probability distributions have been used to model daily precipitation intensity series, including the gamma distribution [42][43][44], a stretched exponential distribution based on the idea that precipitation intensity can be thought of as the product of three independent Gaussian variables representing mass flux, specific humidity, and precipitation efficiency [45], or a hybrid exponential and generalized Pareto distribution [46].However, these distributions have few (2-3) adjustable parameters and cannot be made to fit daily precipitation series well over the entire range of precipitation intensity.By contrast, the hyperexponential distribution provides an analytically tractable approximation to the observed precipitation intensity distribution that can approximate well even long-tailed positive-valued distributions with monotone decreasing densities, such as those for queuing times [47].Wilks [48] previously considered a hyperexponential distribution with N = 2 for precipitation intensity, and [49,50] use an exponential distribution (i.e., N = 1) to model hourly precipitation intensity.
Given the fitted hyperexponential distribution, the link function is a normal quantile transform [51]: where P H (x) = x 0 p H (u)du denotes the cumulative distribution function obtained by integrating the hyperexponential density function p H above, and P inv N is the inverse of the standard normal cumulative distribution function.
The skill of the fitted probability distribution was evaluated using leave-one-month-out cross-validation with root mean square error (RMSE) and mean negative log likelihood (NLL) metrics.NLL is a particularly useful metric for probabilistic models because it evaluates not only the quality of the model's central estimate, as RMSE does, but also whether the model standard deviation is consistent with the model-observation spread [27].The difference in NLL between an augmented model and a baseline can be interpreted as the information gain in bits (assuming base-2 logarithms are taken) from the data (such as 3B42RT) included in the augmented model.RMSE and NLL averages were computed over all APHRODITE rainy pixels (x > 0) for the 2001-2007 period.

Precipitation Occurrence Modeling
For precipitation occurrence, we assumed a logistic regression model, one form of a generalized linear model: where L(p) denotes the logistic function log p 1−p , Z is a 1 × k vector of predictor values identical with that used in the precipitation amount model, and γ is a k × 1 vector of coefficients.Given the predictors and APHRODITE data over 2001-2007, the coefficients γ were determined via maximum likelihood, using the trust region Newton method implemented in LIBLINEAR [52,53].
Skill measures for the precipitation occurrence model included both NLL and a probabilistic form of RMSE closely related to the popular Brier skill score [54,55], namely the square root of the average value of (p(x > 0) − p o ) 2 , where the observation p o is equal to 1 if the day was rainy and 0 if not.

Results
We begin with a comparison of precipitation patterns between the satellite products (TMPA 3B42RT and 3B42) and ground-based products (APHRODITE and the AWSs).First, both satellite products have broadly similar regional spatial patterns of mean precipitation as APHRODITE when evaluated for the overlap period of 2001-2007 (Figure 1).3B42 shows the benefits of retrospective calibration against rain gauges in a sharper precipitation gradient across the Himalayas that compares better with APHRODITE, whereas, as previously shown [20], 3B42RT overestimates precipitation over the Tibetan Plateau.On the other hand, 3B42RT better shows the double heavy-precipitation band just south of the Himalayas also seen in APHRODITE and missed by 3B42 because its calibration at 1 • resolution imposes smoothing [6].
We next consider how the correspondence of satellite-derived with ground-based precipitation changes with timescale (Figure 2).We use rank correlation as an appropriate overall measure of correspondence, as this is less sensitive to the non-normal distribution of precipitation and to outlying values than Pearson product-moment correlation.The correlation coefficients of 3B42RT and 3B42 with APHRODITE are comparable on the daily timescale (0.58 and 0.61 respectively), and increase with aggregation time to 0.90 and 0.96 respectively at 30 days.Correlations with the AWSs behave similarly on the daily to monthly scales (coefficients of 0.53 and 0.62 at 1 day and 0.70 and 0.90 at 30 days-the somewhat lower correlations are expected given that we are comparing an individual station and not an interpolated product with the gridded satellite precipitation field), but show a sharp drop-off at shorter time scales than daily, down to 0.26 and 0.30 at 3 hours.On the other hand, the satellite products do roughly reproduce the pronounced precipitation diurnal cycle seen in the AWS observations, which features a maximum in early morning (local time) and a minimum in the afternoon (Figure 3).Previous studies suggest that the precipitation diurnal cycle in Nepal varies by season [56] and differs for precipitation occurrence versus amount [57], which we may consider in a future study.
At the daily timescale, 3B42RT was better correlated with APHRODITE when the satellite data quality was higher (as measured by more 3-hour periods with microwave precipitation rate determinations) and when the ground-based data quality is higher (as measured by more 0.05 • squares with reporting gauges within the 0.25 • pixel), with the satellite data quality having a more pronounced impact (Figure 4).This highlights the potential role of better observations in improving satellite-based precipitation products.The regional APHRODITE daily precipitation amount distribution turned out to be fit very well by the sum of N = 5 exponentials (Figure 5).The exponential function was a good basis for this probability distribution because the empirical probability density decreased for increasing precipitation amounts.(APHRODITE's spatial interpolation procedure favors showing light precipitation amounts rather than exactly zero precipitation, with about 60% nonzero daily precipitation amounts in the region.) This baseline hyperexponential distribution was modified based on location, season, and 3B42RT precipitation to give spatiotemporally varying predictive distributions.An example of how this affects the base distribution is shown in Figure 6, illustrating a situation when both the monsoon season and the detection of heavy precipitation by 3B42RT implies a higher probability of heavy precipitation compared to the baseline.Probabilities for any precipitation level of interest being exceeded can be mapped over the entire region based on our combined climatology and satellite based model, and can be seen to reflect both climatological precipitation patterns and precipitation reported in 3B42RT (Figure 7).Skill measures for the different probabilistic models show that both the climatology and 3B42RT measures yield improvements over the baseline probability distribution, and their combination yields the best performing models for both precipitation occurrence and precipitation amount (Table 1).Note that while the satellite product yields some 0.19 bits of information gain over baseline for precipitation amount (S versus B NLL in Table 1), satellite plus climatology only outperforms climatology by 0.04 bits (Comb versus C).Results for the other skill measures show the same trend, presumably because much of the correlation of the satellite product with APHRODITE reflects spatial and mean seasonal patterns of variability that can be captured in climatology, whereas it does not offer as much useful information on interannual variability [6].
Well calibrated probabilities are an important feature for a probabilistic monitoring model, and the models developed here indeed generally give probabilities that match the observed (APHRODITE) frequency of occurrence for precipitation over a wide range of threshold values and model probabilities (Figure 8), despite their linearity in transformed precipitation amount.

Discussion and Conclusions
Likely the most useful extension of the model offered here would be adding other sources of information available in near real time.These include numerical weather prediction model output fields, regional circulation pattern variables and global modes of variability that affect precipitation expectation [58][59][60], and any ground-based radar or weather stations available in near real time.Improvements in the availability and quality of past ground-based precipitation datasets for model calibration (such as extending APHRODITE past 2007) would also be expected to improve precipitation retrievals.While most precipitation in the region occurs in summertime as rain, special treatment of snow may be explored to improve high-altitude winter precipitation retrievals [61] considering that the current algorithms for converting satellite sensor values to precipitation rates have been developed primarily for liquid precipitation and have been shown to underestimate snowfall water equivalent [62].
Including spatial and temporal correlations is another possible direction for extension.How to efficiently specify mutivariate probability distributions in highly non-normal quantities such as daily precipitation amounts is an active area of research [29,46,63,64].
Our work here presents a comprehensive approach to constructing a probabilistic model of precipitation given an imperfect, deterministic satellite precipitation product without explicit error information and ground-based observations from past years.We anticipate that this method could be adapted to different regions where near-real-time precipitation observations are scarce.

Figure 1 .
Figure 1.Mean precipitation (mm•d −1 ) over 2001-2007 from the gauge-based APHRODITE product and from the satellite-based 3B42RT and 3B42 products for Nepal and vicinity.

Figure 2 .
Figure 2. Rank correlation of precipitation of the satellite-based 3B42RT and 3B42 products as a function of aggregation timescales: (a) Correlation with the gauge-based APHRODITE product for 2001-2007 over aggregation timescales from 1 to 90 days; (b) Correlation with automated weather stations for 2013-2014 over aggregation timescales from 3 hours to 90 days.

Figure 3 .
Figure 3. Mean precipitation diurnal cycle, 2013-2014, averaged across 3 stations in Nepal (5-minute values and a Fourier series smoothed fit are shown) and for satellite products with 3-hour resolution subsampled at the same grid points.Nepal Standard Time is 5:45 hours ahead of Universal Time.

Figure 4 .
Figure 4. Rank correlation of daily precipitation between the satellite-based 3B42RT product and APHRODITE as a function of both station data coverage (expressed as the number of 0.05 • subcells in a 0.25 • grid cell with precipitation gauges) and satellite data coverage (availability of higher-quality microwave (MW) sounder precipitation estimates for at least 4 of the daily 3-hour windows versus only lower-quality estimates based on infrared (IR) imagery).Error bars are 95% confidence intervals for the correlation coefficients.

Figure 5 .
Figure 5. Empirical probability density for APHRODITE daily precipitation amount at the 0.25 • grid scale (conditional on occurrence) for the region 26 • -31 • N, 79 • -89 • E and for two time periods, compared with a hyperexponential distribution fit only to the 1981-2000 data.

Figure 6 .
Figure 6.Example probability density functions for daily precipitation amount at the grid point containing Kathmandu, Nepal (27.7 • N, 85.3 • E) for a heavy-rainfall monsoon day (2 September 2013): background PDF, climatology PDF taking into account location and season, and probabilities that incorporate data from the satellite precipitation product 3B42RT.

Table 1 .Figure 8 .
Figure 8. Frequency of precipitation above different threshold values ((a) 1; (b) 10; (c) 50 mm/day) as a function of modeled probability (combined satellite and climatology model).1-1 lines indicating an ideally calibrated model are also shown.