A Spatio-Temporal Analysis of Rainfall and Drought Monitoring in the Tharparkar Region of Pakistan

: The Tharpakar desert region of Pakistan supports a population approaching two million, dependent on rain-fed agriculture as the main livelihood. The almost doubling of population in the last two decades, coupled with low and variable rainfall, makes this one of the world’s most food-insecure regions. This paper examines satellite-based rainfall estimates and biomass data as a means to supplement sparsely distributed rainfall stations and to provide timely estimates of seasonal growth indicators in farmlands. Satellite dekadal and monthly rainfall estimates gave good correlations with ground station data, ranging from R = 0.75 to R = 0.97 over a 19-year period, with tendency for overestimation from the Tropical Rainfall Monitoring Mission (TRMM) and underestimation from Climate Hazards Group Infrared Precipitation with Stations (CHIRPS) datasets. CHIRPS was selected for further modeling, as overestimation from TRMM implies the risk of under-predicting drought. The use of satellite rainfall products from CHIRPS was also essential for derivation of spatial estimates of phenological variables and rainfall criteria for comparison with normalized di ﬀ erence vegetation index (NDVI)-based biomass productivity. This is because, in this arid region where drought is common and rainfall unpredictable, determination of phenological thresholds based on vegetation indices proved unreliable. Mapped rainfall distributions across Tharparkar were found to di ﬀ er substantially from those of maximum biomass (NDVI max ), often showing low NDVI max in zones of higher annual rainfall, and vice versa. This mismatch occurs in both wet and dry years. Maps of rainfall intensity suggest that low yields often occur in areas with intense rain causing damage to ripening crops, and that total rainfall in a season is less important than sustained water supply. Correlations between rainfall variables and NDVI max indicate the di ﬃ culty of predicting drought early in the growing season in this region of extreme climatic variability. Mapped rainfall and biomass distributions can be used to recommend settlement in areas of more consistent rainfall.


Introduction
The Thar Desert, located to the northwest of the Indian subcontinent, is one of the largest subtropical deserts [1]. It is regarded as the only fertile desert in the world, with a population dependent on rain-fed agriculture and livestock rearing, and it was declared by the World Food Program as the most food-insecure region of Pakistan [2]. Pakistan itself is the fifth most affected country in the long-term Climate Risk Index [3]. The climate of the Thar is characterized by low and erratic rainfall, high temperature, and long spells of dry weather, which threaten agricultural livelihoods. Tharparkar is the largest district in Sindh province, but climatic stations in Sindh are sparsely distributed between 100 and 150 km apart. Mithi, the only continuous rain gauge station within Tharparkar (Figure 1), has 277 mm of rainfall annually, but this varies greatly from year year. Pressure on water resources is compounded by population growth, as the population of Tharparkar increased from 0.4 million in the 1960s to 0.9 million in 1998 [4], and 1.649 million in 2017 [5], almost doubling over the last two decades (Table 1).  Table 1. Human population growth in Tharparkar district [5]. The increase in frequency and severity of droughts in recent decades reduced agricultural yield, leading to increased poverty and food insecurity. In fact, in February 2014, drought due to extremely low rainfall in Tharparkar resulted in 167 human fatalities and great loss of animal life [2]. A recent study [1] challenged prevailing notions about aridity and changing rainfall pattern in Tharparkar. They found that lower agricultural yield, leading to food scarcity in the last four decades, was actually accompanied by increased annual rainfall at the rate of 6.35 mm/year. Thus, they attribute recent hardships of food insecurity and water scarcity in Tharparkar to erratic rainfall rather than overall rainfall deficit. With an average winter temperature of approximately 20 °C, temperature is never a limiting factor to plant productivity.
The scenario of increasing hunger and starvation in recent years demands a better response from authorities in the future, for which identifying the spatial and temporal dimensions of drought in any season is required. A number of studies used satellite datasets to identify and characterize drought by comparing rainfall products [6,7] with vegetation indices over time [8][9][10]. However, the ability to detect patterns of past drought using time-series satellite data does not automatically translate into drought prediction ability. Indeed, Aziz et al. [11] working in Pakistan demonstrated satellites' ability to detect drought over large areas in Pakistan using vegetation indices, but did not show how drought or its impacts could be predicted in advance.
More recent developments using satellite data for phenological analysis show more promise if spatial relationships between total crop yields and seasonal phenological variables can be established. These variables include the start, length, and end of the growing season, drought periods, and the nature of rainfall, whether early, late, frequent, intense, or well-distributed. Such variables can potentially be mapped using satellite images and their trends identified [12][13][14]. However, different methodologies have been shown to give significantly different seasonal transition dates [14] and  Table 1. Human population growth in Tharparkar district [5]. The increase in frequency and severity of droughts in recent decades reduced agricultural yield, leading to increased poverty and food insecurity. In fact, in February 2014, drought due to extremely low rainfall in Tharparkar resulted in 167 human fatalities and great loss of animal life [2]. A recent study [1] challenged prevailing notions about aridity and changing rainfall pattern in Tharparkar. They found that lower agricultural yield, leading to food scarcity in the last four decades, was actually accompanied by increased annual rainfall at the rate of 6.35 mm/year. Thus, they attribute recent hardships of food insecurity and water scarcity in Tharparkar to erratic rainfall rather than overall rainfall deficit. With an average winter temperature of approximately 20 • C, temperature is never a limiting factor to plant productivity.
The scenario of increasing hunger and starvation in recent years demands a better response from authorities in the future, for which identifying the spatial and temporal dimensions of drought in any season is required. A number of studies used satellite datasets to identify and characterize drought by comparing rainfall products [6,7] with vegetation indices over time [8][9][10]. However, the ability to detect patterns of past drought using time-series satellite data does not automatically translate into drought prediction ability. Indeed, Aziz et al. [11] working in Pakistan demonstrated satellites' ability to detect drought over large areas in Pakistan using vegetation indices, but did not show how drought or its impacts could be predicted in advance.
More recent developments using satellite data for phenological analysis show more promise if spatial relationships between total crop yields and seasonal phenological variables can be established. These variables include the start, length, and end of the growing season, drought periods, and the nature of rainfall, whether early, late, frequent, intense, or well-distributed. Such variables can potentially be Remote Sens. 2020, 12, 580 3 of 16 mapped using satellite images and their trends identified [12][13][14]. However, different methodologies have been shown to give significantly different seasonal transition dates [14] and algorithms effective in one study area may not apply in others. White et al. [15] compared 10 different techniques for defining the start of the growing season (SOS), finding average differences of ±60 days and standard deviation of ±20 days among the different methods. Since temperature is the most common limiting factor to growth, many studies used temperature thresholds to define growing season days such as daily average above 5 • C [16]. Landscapes in more arid regions and/or those with low seasonal amplitude of vegetation indices were shown by several studies to give the lowest agreement between satellite normalized difference vegetation index (NDVI) and ground-measured phenology [12][13][14]. Kang et al. [12] also demonstrated uncertainty in phenology measurement arising from the spatial resolutions of different satellite sensors.
Thus, satellite-based precipitation estimates may be able to provide a viable alternative to sparsely distributed climate station data, as well as provide regular indication of vegetation productivity over the whole region. Satellite rainfall products have high spatial and temporal resolution. Thus, they provide timely, repetitive, and cost-effective information of rainfall at different time scales from daily to annual [6]. However, large differences were observed in algorithm performance [17,18], and, for drought monitoring, assessment of low rainfall situations is essential [17,19]. This would enable spatially detailed assessment of the relationship between vegetation productivity and rainfall for advanced drought warnings and relief measures.
The specific objectives of this study are as follows: (1) To evaluate satellite-based rainfall estimates against ground station data in their ability to provide realistic and useful estimates of on-ground rainfall conditions and adequacy throughout the growing season; (2) To determine if satellite NDVI data can be used along with satellite rainfall data to provide timely crop forecasts, thus enabling rapid response to crop failure in drought situations.

Study Area
Tharparkar is the largest district in the arid Sindh province of Pakistan within latitude 24 • -26 • north (N) and longitude 69 • -71 • east (E), with Mithi as its district headquarters ( Figure 1). Tharparkar has a population of approximately 1.65 million and an area of 22,000 km 2 , 97% of which is desert. The natural vegetation of grasses and scattered, drought-resistant shrubs and trees was replaced in many areas by croplands, and cultivation occupies 42% of the land area, 98.4% of which is rainfed and 1.6% of which is irrigated from canals. The desert landscape of Tharparkar is dominated by Quaternary parabolic dunes, in the form of sand ridges, deposited by pre-monsoon winds from the southwest [20]. Although the landscape rises gradually from below 50 m in the west to over 150 m in the east, local topography is controlled by the ridge and hollow pattern of the dune-inter-dunes with cultivation taking place mainly in inter-dune depressions. Dug wells are the only source of drinking water in the area, with depths ranging from 18 m in the west to 80 m in the east, but the supply is increasingly brackish, and the water table depth is falling [21]. Local people often construct small ponds lined with clay to store rainwater for domestic and animal use. The subsistence nature of farming in this arid region where rainfall is highly variable introduces considerable risk for livelihoods. The main crops are millet, mung beans, and guar (cluster beans). These crops are planted in July/August soon after the first rains and harvested in October/November following the last rains. In recent years, below average rainfall adversely affected livestock and agriculture, as well as human health, resulting in the death of newborn children and nutrition issues in pregnant and lactating women [5]. Malnutrition is cited as a cause of stunted growth in approximately 50% of children under five years in the southern Sindh Province of Pakistan [22]. Livestock rearing in Thar compounds the water deficiency situation, as the animal population is increasing at almost twice the rate of population growth, from 0.4 million in 1976 to 1.26 million in 1985 to four million today (Table 1). The average water requirement for a cow in Thar is 57 liters per day and 7.5 for a sheep or goat, and, in the 1970s, the ratio of cows to goats was 1:1, whereas it was 1:6 at the 1998 census [4]. Following a severe drought in 2018, with only 58 mm of rainfall received at Mithi, the Sindh government declared a drought calamity in six districts and decided to provide 50 kg of wheat per month to a total of 323,435 families in Tharparkar [23]. During drought years, migration away for alternate livelihood is common. Although drought vulnerability may have reduced due to improved road infrastructure and other economic opportunities, the frequency and intensity of drought appears to have increased since 1990 due to climate change. Thus, there is a trend of increasing exposure to drought [20].

Rainfall in Tharparkar
Annual rainfall varies from 50 to 300 mm, and 87% of Thar's total rainfall is in the monsoon months of June-September [4]. Mean annual rainfall at Mithi station is 277 mm, with a maximum in August, and the only other months with significant rainfall are July and September ( Figure 2). However, rainfall is highly variable, and the rainy season may commence in early June to mid-July ( Figure 3). Any rainfall deficit during the rainy monsoon months severely impacts the socio-economy of the region [24], and Tharparkar was officially declared a drought calamity area 15 times since 1968, with the most recent times being the three years 2013, 2014, and 2018.
Remote Sens. 2020, 12, x FOR PEER REVIEW 4 of 16 alternate livelihood is common. Although drought vulnerability may have reduced due to improved road infrastructure and other economic opportunities, the frequency and intensity of drought appears to have increased since 1990 due to climate change. Thus, there is a trend of increasing exposure to drought [20].

Rainfall in Tharparkar
Annual rainfall varies from 50 to 300 mm, and 87% of Thar's total rainfall is in the monsoon months of June-September [4]. Mean annual rainfall at Mithi station is 277 mm, with a maximum in August, and the only other months with significant rainfall are July and September ( Figure 2). However, rainfall is highly variable, and the rainy season may commence in early June to mid-July ( Figure 3). Any rainfall deficit during the rainy monsoon months severely impacts the socio-economy of the region [24], and Tharparkar was officially declared a drought calamity area 15 times since 1968, with the most recent times being the three years 2013, 2014, and 2018.  alternate livelihood is common. Although drought vulnerability may have reduced due to improved road infrastructure and other economic opportunities, the frequency and intensity of drought appears to have increased since 1990 due to climate change. Thus, there is a trend of increasing exposure to drought [20].

Rainfall in Tharparkar
Annual rainfall varies from 50 to 300 mm, and 87% of Thar's total rainfall is in the monsoon months of June-September [4]. Mean annual rainfall at Mithi station is 277 mm, with a maximum in August, and the only other months with significant rainfall are July and September ( Figure 2). However, rainfall is highly variable, and the rainy season may commence in early June to mid-July ( Figure 3). Any rainfall deficit during the rainy monsoon months severely impacts the socio-economy of the region [24], and Tharparkar was officially declared a drought calamity area 15 times since 1968, with the most recent times being the three years 2013, 2014, and 2018.

Weather Station Rainfall data
Daily rainfall data for four stations ( Figure 1) were obtained from the Pakistan Meteorological Department (PMD) ( Table 2). Daily rainfall data were accumulated to form dekadal (10 days) and monthly rainfall for comparison with satellite-based rainfall estimates.

Satellite-Based Rainfall Products
Satellite-based rainfall products typically exploit a combination of data from thermal infrared (TIR), passive microwave (PMW), and ground-based gauge observations, and these datatypes are often combined to create an optimal product. Various rainfall datasets were produced using convective cloud top temperature by applying the cold cloud duration (CCD) technique [25]. In this study, two satellite-based rainfall products Climate Hazards Group Infrared Precipitation with Stations (CHIRPS) and Tropical Rainfall Monitoring Mission (TRMM) were selected for evaluation against rainfall gauge data, because of their long time series, near-real-time data availability, and free access.

Tropical Rainfall Measuring Mission (TRMM)
The Tropical Rainfall Measuring Mission (TRMM) is a joint mission between the National Aeronautics and Space Administration (NASA) and the Japan Aerospace Exploration Agency (JAXA) aimed at improving observations of precipitation across the globe between 45 • N and 45 • south (S). The most widely used outputs are the TRMM Multi-Satellite Precipitation Analysis (TMPA) three-hourly (3B42) product accumulated to daily and monthly (3B43) products, which are available from 1998 to 2014 at spatial resolution of 25 km [25,26]. The TMPA product depends on input from a combination of optical, thermal, and microwave sensors, as well as gauge data [18]. Daily TRMM 3B42 V7 and monthly 3B43 V7 products were used in this study.

Climate Hazards Group Infrared Precipitation with Stations (CHIRPS)
The CHIRPS Version 2.0 rainfall dataset was developed by the United States (US) Geological Survey (USGS) and Climate Hazard Group at the University of California, Santa Barbara. It is available from 1981 onward at a spatial resolution of 5 km. The CHIRPS algorithm (i) incorporates satellite thermal infrared (IR) data to represent sparsely gauged locations, (ii) blends station data to produce a preliminary information product with a latency of about two days and a final product with an average latency of about three weeks, and (iii) uses a novel blending procedure incorporating the spatial correlation structure of CCD estimates to assign interpolation weights. CHIRPS also uses the TRMM's TMPA product, which includes several microwave sources, to calibrate global cold cloud duration (CCD) rainfall estimates [27]. images at 250-m resolution [28], based on compositing daily images to select the highest-quality, highest-value pixels over 16-day periods.

Statistical Measures
Satellite rainfall estimates were compared with gauge rainfall using pairwise comparison statistical measures, Pearson product moment coefficient of correlation (R), bias, mean error (ME), and root-mean-square error (RMSE).
The Pearson correlation coefficient (R) measures the strength of the linear relationship between satellite and gauge rainfall. Values of R close to 1 indicate a perfect relationship between satellite and gauge rainfall estimates. The statistical significance of correlation (R) is represented by asterisks (** p < 0.01 and * p < 0.05).
where G is the gauge rainfall amount, − G is the average gauge rainfall amount, S is the satellite rainfall estimate, − S is the average satellite rainfall estimate, n is the total number of data. Bias indicates how well the average of variable 1 corresponds with the average of variable 2. Thus, values close to 1 show that cumulative satellite rainfall estimates are close to cumulative gauge rainfall measures, while values greater than 1 indicate that the satellite overestimates the measure, and values less than 1 indicate that the satellite underestimates the measure.
Mean error (ME) is the measure of average difference. Thus, a positive value reflects an overestimation of satellite rainfall, whereas a negative value indicates an underestimation of satellite rainfall.
Root-mean-square error (RMSE) is the standard deviation of the difference. Thus, a higher value of RMSE indicates a large difference between the satellite and gauge rainfall measures.

Smoothing of NDVI Time Series Using Savitzky-Golay filter
Due to cloud cover, varying atmospheric conditions, and bi-directional effects, there are some disturbances in the NDVI time series. During the rainy season, clouds are a major problem and often prevail for more than two weeks. Although NDVI datasets are an Maximum Value Composite product (selection of maximum NDVI value over 16 days), there is still noise, represented by negatively biased NDVI values. To overcome this noise, TIMESAT, which uses a Savitsky-Golay filter [29,30], was applied to smooth the NDVI time series pixel-wise, and to construct high-quality NDVI time-series datasets for further analysis. The Savitsky-Golay filter is a moving filter that fits local polynomial regression to replace negatively biased NDVI values with an upper envelope of NDVI time series.

Phenological Metrics
In order to investigate spatial aspects of various phenological metrics over the study area, these metrics were mainly derived from rainfall data. Rainfall data were used because common phenology algorithms based on vegetation indices such as TIMESAT or singular spectrum analysis [31], which uses adjustable thresholds on the seasonal amplitude of the NDVI curve, proved unsuitable for use in this semi-desert region with highly variable rainfall. For example, the commonly used thresholds of 10% from the base of the NDVI curve to define the start of season (SOS) and 35% for the end of season (EOS) frequently returned a nine-month length of growing season (LOS) into the following year (see maps of EOS for wet and dry years (Supplementary Material S1)). This may be because, at low NDVI, or due to persistence of greenness within image pixels long after the last rains, vegetation signals are often confused with soil reflectance [32] in semi-desert regions. In fact, the maximum time between planting and harvest in Tharparkar is only 4-5 months, from June to October. Indeed Zhang et al.'s [33] phenological study of different vegetation indices compared to PhenoCam ground measurements reported significant differences among different vegetation indices, with EOS prediction showing the greatest uncertainty. Furthermore, the uncertainty was greatest for EOS in dry zone ecosystems, with up to 29 days of difference for savannas compared to seven days for forests. This high level of uncertainty was observed in spite of using daily data. Since our study uses 16-day interval MODIS composites, even greater uncertainty would be expected.
Thus, the seasonal rainfall variable for SOS was defined based on daily rainfall from the CHIRPS satellite rainfall product. We used the definition of Zhang et al. [13] in the Sahel zone and modified it according to local patterns of rainfall in Mithi and surrounding weather stations, as described below.
Rainy season onset was defined as the first occurrence of at least 20 mm of cumulative rainfall within seven consecutive days after 1 June, followed by at least 20 mm of rainfall in the next 20 days to avoid "false starts". The intensity of rainfall was defined as the total amount of rainfall during the rainy season divided by the number of rainy days with rainfall ≥1 mm. NDVI max represents the maximum NDVI values reached during a growing season, and the variables, SOS, and timing of NDVI max are expressed in Julian days. The variable NDVI max was used to represent biomass productivity, as, in this region of rain-fed agriculture, farmers cultivate a single crop per season, and the amount of offtake at harvest is likely to correspond to the maximum biomass achieved in a season. Thus, other measures such as greater or lesser integral, which also consider biomass development over the whole season, were considered less relevant for this study.

Correlation
The Pearson product moment correlation was calculated between pairs of variables (e.g., annual rainfall and NDVI max ) using the 19 grid layers (grids for each year 2000-2018) for each variable over Tharparkar. For each pixel in a dataset, a vector of values from 19-year time series is correlated between the two datasets [34]. To avoid spatial dependence, the CHIRPS dataset at 5-km resolution was firstly resampled to match the 250-m pixel size of the NDVI data. A p-value greater than 0.05 was used to mask pixels where correlation was insignificant at the 95% confidence level.

Dekadal Comparison
For 10-day (dekadal) comparisons, a good relationship was shown between satellite-based rainfall estimates and gauge rainfall, with correlation coefficient R values between 0.60 and 0.85 for individual weather stations. Correlation values for both TRMM and CHIRPS were similar, with R values of 0.75 and 0.74, respectively, while CHIRPS showed slightly higher RMSE value of 22.45 mm/dekad compared to TRMM's RMSE value of 21.93 mm/dekad. CHIRPS showed some underestimation of dekadal rainfall with a bias of 0.73 and negative mean error (ME) value of −2.08 mm/dekad, while TRMM showed overestimation, with a bias of 1.12 and positive mean error value of 0.94 mm/dekad.

Monthly Comparison
Monthly satellite rainfall estimates showed better performance than daily or dekadal timeframes, with correlation values between 0.80 and 0.97 for individual weather stations. This improvement was because aggregation of daily or dekadal data into monthly values canceled out short-term errors at daily or dekadal scales. TRMM, with a higher correlation coefficient (R = 0.90) and lower RMSE value of 28.6 mm/month, performed better overall than CHIRPS (R = 0.84 and RMSE of 38.6 mm/month). Similar to daily and dekadal scales, CHIRPS showed underestimation of monthly rainfall with a bias of 0.77 and negative mean error (ME) value of −5.05 mm/month, while TRMM showed overestimation of monthly rainfall with a bias of 1.16 and positive mean error of 3.61 mm/month. For drought monitoring studies, overestimation of satellite rainfall (ME > 0) should be avoided, whereas underestimation is undesirable for hydrological and flood forecasting studies [6]. Since the CHIRPS dataset showed a tendency to underestimate rather than overestimate rainfall, it was selected for further analysis in this study, comparing spatial distributions of phenological growth parameters with rainfall. Figure 5a shows the spatial distribution of annual rainfall over Tharparkar from CHIRPS data, with a steep decline from over 400 mm in the east to less than half of this (below 200 mm) in the west. The strong positive correlation between annual rainfall and NDVI max , with R values mainly from 0.6 to 1.0 (Figure 5b), indicates the strong dependence of agricultural productivity on rainfall, except for the northeast and along canals (masked areas in white).

Spatial Relationships between Annual Rainfall and Phenology Metrics
However, when exceptionally wet and dry years were examined separately, the spatial patterns were different (Figures 6 and 7). In dry years, there appeared to be little spatial relationship between annual rainfall (Figure 6a,d) and biomass productivity, represented by NDVI max (Figure 6e-h). Highest rainfall occurred in the southeast, but highest NDVI max values generally occurred in the west, northwest, and east. In wet years, although there still appeared to be little spatial relationship between annual rainfall (Figure 7a-d) and biomass productivity (Figure 7e-h), highest rainfall occurred in the northeast, but NDVI max was highest in the northwest, southeast, and south-central. In the northeast, with very high rainfall above 600 mm in 2006 and 2011, NDVI max values remained as low as 0.3, but reached NDVI max > 0.6 in some areas with only 400 mm of rainfall. This may suggest that farming/settlement patterns are not taking advantage of available rainfall, which is highest in the east. Figure 8 shows the time of year when NDVI max was reached, for the four dry and four wet years, classified in two-week intervals. No spatial or temporal pattern emerged, indicating that the crop harvest in any particular year, whether a dry or wet year, could be early or late, depending on the unpredictability of the seasonal rainfall pattern. Moreover, in a dry year, if there was at least some rainfall, exceeding 150 mm as in 2004 and 2014, harvest could be as late as November. and east. In wet years, although there still appeared to be little spatial relationship between annual rainfall (Figure 7a-d) and biomass productivity (Figure 7e-h), highest rainfall occurred in the northeast, but NDVImax was highest in the northwest, southeast, and south-central. In the northeast, with very high rainfall above 600 mm in 2006 and 2011, NDVImax values remained as low as 0.3, but reached NDVImax > 0.6 in some areas with only 400 mm of rainfall. This may suggest that farming/settlement patterns are not taking advantage of available rainfall, which is highest in the east.           Figure 8 shows the time of year when NDVImax was reached, for the four dry and four wet years, classified in two-week intervals. No spatial or temporal pattern emerged, indicating that the crop harvest in any particular year, whether a dry or wet year, could be early or late, depending on the unpredictability of the seasonal rainfall pattern. Moreover, in a dry year, if there was at least some rainfall, exceeding 150 mm as in 2004 and 2014, harvest could be as late as November. Figure 8. Timing of NDVI max for four dry years (a-d) and four wet years (e-h). White areas represent not enough rainfall from July to December; thus, there is no maximum NDVI value.

Correlations
The moderate, positive correlation between annual rainfall and SOS, with R = 0.47 (Table 3), suggests that an early start to the growing season (low SOS value in days of the year (DOY)) often occurred in years with low total rainfall; thus, an early SOS with early rains should not be used as indicative of a high-rainfall year. However, since the relationship between SOS and NDVI max was negative (R = −0.51), this indicates that early onset of season (low SOS value) tended to result in high crop yields. Thus, an early start to growth followed by continued even distribution of rainfall may be more important than total annual rainfall for good crop yields in a year. A fairly strong correlation (R = 0.77) between mid-to-late rainy season (August) rainfall and annual rainfall was observed, and August rainfall also showed the highest correlation (R = 0.62) with NDVI max . Unfortunately, this late rainy season observation would not be useful for enabling early prediction of drought conditions.
The metric time of NDVI max correlated negatively (R = −0.59) with NDVI max , suggesting that highest biomass was achieved when crops ripened early. This was especially true for years with adequate rainfall, as can be seen from the mapped spatial distributions. Thus, comparisons between Figures 6 and 8, and between Figures 7 and 8 also show that, spatially, high NDVI max was achieved in areas with early NDVI max , especially for wet years (compare Figure 7f,h with Figure 8f,h). This agrees with the observation (above) that the highest biomass in wet years was achieved in areas with early SOS. However, in dry years (compare Figure 6f,g with Figure 8b,c), highest biomass was apparently achieved late in the growing season, with harvesting as late as November, as in the drought years of 2004 and 2014. Table 3. Correlation coefficients (R) for mapped phenology metrics (significant pixels only (p = 0.05)). Annual rainfall is derived from amalgamation of CHIRPS monthly estimates, and start of the growing season (SOS) is derived from daily estimates.  Figure 9 shows rainfall intensity for (a-d) dry and (e-h) wet years (representing the total amount of rainfall during the rainy season divided by the number of rainy days with rainfall ≥1 mm). indicate that, in wet years, there was a tendency for more intense rainfall toward the north and east of the study area, and, in dry years, this tendency was toward the southeast. These are the areas where relatively low NDVI max occurred, corresponding to relatively high annual rainfall. As it is known that crop damage can occur from intense rainstorms preceding the harvest in Tharparkar [35], this may partly explain the apparent low productivity of these areas.
Remote Sens. 2020, 12, x FOR PEER REVIEW 12 of 16 Figure 9. Rainfall intensity in four dry years (a-d) and four wet years (e-h).

Discussion
Overall, CHIRPS showed better estimates for daily rainfall, while TRMM showed better monthly estimates, and the results were similar for dekadal estimates. The lower accuracy observed for daily estimates than for dekadal and monthly estimates is presumed to be partly due to the coarse spatial resolution of the sensors compared to the rain gauge point location. In areas of very low rainfall, spatial distribution of rain is also likely to be patchy. The longer time intervals (i.e., dekdal and monthly) would tend to average out the uncertainty of the daily data. CHIRPS tended to underestimate while TRMM tended to overestimate measures, but overestimation is considered more serious as this may overlook a drought situation. Therefore, the study used CHIRPS data for the spatial representation of rainfall metrics and their spatial correlation with biomass in the form of NDVImax and timing of NDVImax.
The very strong spatial relationship observed between CHIRPS-derived annual rainfall and NDVImax confirmed that rainfall is fundamental to farming in Tharparkar. The study observed a negative correlation between SOS and NDVImax, meaning that early start to the season tended to result in high biomass and vice versa, which would normally be expected. However, the observation of positive correlation between SOS and annual rainfall means that years with early rainfall may have low total seasonal rainfall and vice versa, which is counter-intuitive. These two points together suggest that the total rainfall in a season is not the most important factor in biomass productivity. They suggest that sustained and consistent distribution of rainfall may be more important in this region of low but highly variable rainfall. Thus, high biomass may well be achieved when SOS is early, but this is not necessarily dependent upon very high total seasonal rainfall. Rather, subsequent sustained rainfall is required to achieve good crop yields. Thus, even distribution of rainfall may be more important than total annual rainfall for good crop yields in a year. Indeed, although pearl millet can survive low annual rainfall of 300 mm, the seasonal distribution is more important. Rainfall is required before sowing, followed by steady water supply after germination, then again during flowering and fruiting 40-60 days after germination [36]. However, damage may occur with heavy rain before harvest, which may cause sprouting, breaking of the stalk, and bird feeding. Mung bean requires slightly higher rainfall, of 400 mm, with similar seasonal distribution to millet, i.e., a continued water supply after germination but more in the July-August flowering/fruiting stage [35]. Thus, the maps of rainfall intensity, combined with those of annual rainfall, enable the monitoring of evenly distributed and sustained rainfall, which appears more important than total seasonal rainfall in crop production in Tharparkar. This finding from satellite image data supports the conclusions of Memon et al. [1], who compared rain gauge data over 38 years with farm questionnaires about response to drought in Tharparkar. They reported the main cause of crop failure to be the increasingly erratic pattern of rain, especially decline in the July mid-growing season, in spite of a slight increase

Discussion
Overall, CHIRPS showed better estimates for daily rainfall, while TRMM showed better monthly estimates, and the results were similar for dekadal estimates. The lower accuracy observed for daily estimates than for dekadal and monthly estimates is presumed to be partly due to the coarse spatial resolution of the sensors compared to the rain gauge point location. In areas of very low rainfall, spatial distribution of rain is also likely to be patchy. The longer time intervals (i.e., dekdal and monthly) would tend to average out the uncertainty of the daily data. CHIRPS tended to underestimate while TRMM tended to overestimate measures, but overestimation is considered more serious as this may overlook a drought situation. Therefore, the study used CHIRPS data for the spatial representation of rainfall metrics and their spatial correlation with biomass in the form of NDVI max and timing of NDVI max .
The very strong spatial relationship observed between CHIRPS-derived annual rainfall and NDVI max confirmed that rainfall is fundamental to farming in Tharparkar. The study observed a negative correlation between SOS and NDVI max , meaning that early start to the season tended to result in high biomass and vice versa, which would normally be expected. However, the observation of positive correlation between SOS and annual rainfall means that years with early rainfall may have low total seasonal rainfall and vice versa, which is counter-intuitive. These two points together suggest that the total rainfall in a season is not the most important factor in biomass productivity. They suggest that sustained and consistent distribution of rainfall may be more important in this region of low but highly variable rainfall. Thus, high biomass may well be achieved when SOS is early, but this is not necessarily dependent upon very high total seasonal rainfall. Rather, subsequent sustained rainfall is required to achieve good crop yields. Thus, even distribution of rainfall may be more important than total annual rainfall for good crop yields in a year. Indeed, although pearl millet can survive low annual rainfall of 300 mm, the seasonal distribution is more important. Rainfall is required before sowing, followed by steady water supply after germination, then again during flowering and fruiting 40-60 days after germination [36]. However, damage may occur with heavy rain before harvest, which may cause sprouting, breaking of the stalk, and bird feeding. Mung bean requires slightly higher rainfall, of 400 mm, with similar seasonal distribution to millet, i.e., a continued water supply after germination but more in the July-August flowering/fruiting stage [35]. Thus, the maps of rainfall intensity, combined with those of annual rainfall, enable the monitoring of evenly distributed and sustained rainfall, which appears more important than total seasonal rainfall in crop production in Tharparkar. This finding from satellite image data supports the conclusions of Memon et al. [1], who compared rain gauge data over 38 years with farm questionnaires about response to drought in Tharparkar. They reported the main cause of crop failure to be the increasingly erratic pattern of rain, especially decline in the July mid-growing season, in spite of a slight increase in overall annual rainfall over the last four decades.
In wet years, such as, for example, 2006 and 2011, with above 700 mm of annual rainfall, the apparent lack of spatial correspondence between the mapped annual rainfall and NDVI max may be due to susceptibility of arid zone crops, such as millet and mung beans, to heavy and intense rain. The mapped distributions show that heaviest rainfall in these wet years occurred in areas which are usually driest, with mean annual rainfall of 250-300 mm; thus, farmers would not expect such heavy rain. This was confirmed by Figure 9e-h, showing more intense rainfall in these areas, i.e., the northeast of the study area in wet years and southeast in dry years. More intense rainfall also suggests that a considerable proportion of the annual rainfall would occur on fewer days, with uneven seasonal distribution. This would not satisfy the water requirement of millet and mung bean, which require evenly distributed water supply.
In dry years, areas showing an NDVI max response (areas with no NDVI response in dry years were not mapped) were not the areas of highest rainfall; thus, dry years also lack spatial correspondence with annual rainfall. Timing of maximum biomass in such dry years was also observed to be delayed until October/November in these areas. For example, in these areas, under 200 mm of rain fell in 2004 and 2014, and ripening/harvest extended into late October, possibly due to farmers hoping for late rain.
In terms of overall ability to predict drought based on the mapped distributions, the only moderate correlations between SOS and annual rainfall, and between SOS and NDVI max make early rain and early start to the growing season fairly poor predictors of seasonal crop productivity. Additionally, the moderately good correlations observed between mid-to-late rainy season (August) rainfall and both annual rainfall (R = 0.77) and maximum biomass (NDVI max ) (R = 0.62) would unfortunately be too late to predict drought conditions and/or low crop yields. Low August rainfall may be the first indicator of imminent low crop yields and economic hardships among farming families. However, according to Pasha et al. [2], drought in Thar, as well as starvation and death resulting from it, was reported in February of 2014, which indicates that greatest distress occurs late in the dry season following low rainfall the previous year. Indeed, Mithi station recorded only 190 mm in 2013. Therefore, knowledge of the location of failing rains in August would indeed give an early warning of imminent severe distress in the following dry season among farming families in those areas affected.

Conclusions
The study demonstrated the effective utilization of satellite-based rainfall products for spatial analysis of phenology over the Tharparkar desert region, as well as implications and policy indicators arising from the analysis. Although both CHIRPS and TRMM generally gave reliable rainfall estimates compared to ground stations, CHIRPS daily estimates were more reliable than those for TRMM and the tendency for underestimation by CHIRPS, compared to overestimation by TRMM, made CHIRPS better suited for use in a drought-prone region where overestimation may overlook serious drought periods. In this arid region where drought is common and rainfall unpredictable, determination of phenological thresholds based on vegetation indices proved unreliable. Therefore, the use of satellite rainfall products from CHIRPS to derive spatial estimates of SOS and rainfall criteria for comparison with NDVI-based biomass productivity was essential for this study.
The counter-intuitive observation of high rainfall associated with low NDVI max is unlikely to be due to soil differences, as the Thar desert is dominated by compound parabolic sand dunes, arranged in a repetitive pattern across the whole landscape. Thus, although local soil differences do exist, these are at a dune scale, rather than a regional scale. A more likely explanation based on the data presented here is damage to ripening crops by intense rainfall, as indicated by the spatial correspondence between maps depicting areas of low NDVI max , high annual rainfall, and more intense rainfall. Rainfall of high intensity in an arid region also implies uneven rainfall distribution, which would be damaging to major crops grown in Tharparkar, i.e., millet and various types of beans. A false start to the rainy season encouraging farmers to plant may see them without sufficient seed for a second planting, even if sufficient rain were to return, resulting in low yield for such areas, even given adequate rain. However, if a false start does occur in any region (rainfall over 20 mm followed by a dry spell of no rain for 20 days), this can be determined using CHIRPS dekadal data, with a time lag of three dekads, i.e., 30 days, at most. This would allow potential relief measures in areas affected, such as construction of water storage reservoirs and food distribution and storage facilities.
Overall, the study demonstrated the great difficulty in predicting seasonal rainfall due to the very high variability of rainfall in Tharparkar; thus, the occurrence of drought cannot be known even one month in advance. The ability to map the spatial distribution of rainfall and biomass using satellite products can, however, permit closer and timely analysis of the whole seasonal farming calendar. This may include locating areas of potential economic hardship, malnutrition, and hunger in the long dry season resulting from low seasonal rainfall the previous year. The data also enable recommendations for the redistribution of settlements to areas of higher and more reliable rainfall.