Pages 45 - 55
Keywords Running average cross-correlation, Lag-normalized correlation, VI anomalies, Drought-prone areas, Hazard mapping
Anjillyn Mae C. Perez * and Ariel C. Blanco
The authors are with the Department of Geodetic Engineering, University of the Philippines, Diliman, Quezon City, Philippines 1101
* Corresponding author: firstname.lastname@example.org
The Philippines is experiencing recurring drought events accentuated by the increasing incidence of the El Niño phenomenon, particularly due to its position in the equatorial region. To address the negative effects of unpredicted drought events, especially in light of perceptible climate change, the present research developed a methodology for drought assessment in the Cagayan River Basin (CRB) by combining remotely-sensed data including 16-day MODIS Vegetation Indices (VIs) (250 m resolution), monthly Terrestrial Water Storage Changes (TWSC) in 1° grids from NASA’s Gravity Recovery and Climate Experiment (GRACE), and daily in-situ measurements of rainfall and water levels. Time series data of these parameters from 2002 to 2011, a period when major and minor drought events occurred in the country, were statistically analyzed. To smooth out short-term and random variations in the individual time series, a 3-period running average was first applied to each time series. The resulting time series datasets were then subjected to cross-correlation analysis to examine whether the VI anomalies and TWSC values were statistically related to rainfall and water level. The cross-correlation results showed that VI anomalies and TWSC values exhibited strong correlation with in-situ rainfall and water level measurements, having coefficients greater than 0.90 and low time lags ranging from 1 to 3 months during drought events. From this observation, lag-normalized correlation analysis between TWSC and VI anomalies was developed to characterize the onset of drought events and identify drought-prone areas. The majority of the areas identified as susceptible to drought are located in the provinces of Cagayan and Isabela, with areas of 2,596 ha or 68.21 % of the total land area of Cagayan, and 3,751 ha or 56.89 % of the total land area of Isabela. The findings in this research can serve as a basis for proper water resource allocation, especially in the identified drought-prone provinces in the country.
Drought is defined by the United Nations Convention to Combat Desertification (UNCCD) as “a naturally occurring phenomenon that exists when precipitation has been significantly below normal recorded levels, causing serious hydrological imbalances that adversely affect land resource production systems” . It is a normal and recurring climatic event that is considered as “a creeping phenomenon” since its effects often take weeks or months to manifest. In addition, the beginning and ending of drought events are difficult to determine  as well as the quantification of thier severity . According to the Intergovernmental Panel on Climate Change (IPCC) , drought events of different spatial and temporal scales are expected to increase due to the variability in precipitation patterns.
In the Philippines, several El Niño Southern Oscillation (ENSO) episodes in the central and eastern equatorial Pacific have caused minor to major droughts in vulnerable provinces, impacting mostly agriculture, livestock, fisheries, and forestry . The increasing frequency of drought can result in loss of land productivity, crop failure, loss of biodiversity and fragile ecosystems, and even loss of human lives and properties . Drought can affect not only the natural environment, but also socio-economic aspects, such as health, power generation, and food security. However, decision-makers, stakeholders and local people in drought-prone areas are not sufficiently equipped with efficient drought monitoring systems and preparedness plans to minimize the adverse effects of such phenomena.
Drought events in the Philippines are monitored by PAGASA Climatology and Agrometeorology Branch (CAB) based on downscaled predictions from global climate centers in the USA, Australia, Great Britain, and Japan . CAB also uses different drought indices such as Moisture Availability Index (MAI), Yield Moisture Index (YMI), Rainfall Extreme Index (REI), Generalized Monsoon Index (GMI), and Percent of Normal , which are mainly derived from precipitation data. A recent research combined remotely-sensed Normalized Difference Vegetation Index (NDVI) and Land Surface Temperature (LST) from MODIS to derive the Standardized Vegetation-Temperature Ratio (SVTR) to forecast and monitor agricultural drought in the Philippines .
Generally, drought monitoring schemes commonly rely heavily on precipitation data and do not incorporate deep soil moisture and groundwater storage conditions . Due to intensive logistical requirements needed to continuously monitor actual soil moisture and groundwater levels through field measurements, the need to incorporate remotely-sensed satellite images to supplement ground-based weather data is needed. This is necessary to detect the onset and identify the location and severity of drought events, for a spatially extensive, close to real-time, and comprehensive assessment of current drought conditions , . With the advancement of Remote Sensing (RS) systems, drought metrics with improved spatial detail, without the need for extensive meteorological precipitation networks, have become one of the most promising tools for large-scale drought monitoring .
The relationship between RS-derived Terrestrial Water Storage (TWS) anomalies and Vegetation Indices (VIs) has been explored to test their potential use in drought monitoring in Western and Central Europe . It was observed that the correlation between ground water storage anomalies and VIs becomes stronger when warm seasons start. However, considering the spatial resolution of available remotely-sensed TWS and VI data (e.g. 250-m pixel size and 1° grid size), their use for drought studies are mostly confined to large river basins and national-scale analysis. Thus, there is still a need to assess the accuracy and applicability of these datasets for local-scale drought monitoring in smaller river basins.
The present research aims to assess the applicability of combining remotely-sensed data and in-situ measurements of rainfall and water levels for the identification of drought-prone areas in a local river basin in the Philippines. Specifically, it intends to develop a method of drought assessment by utilizing freely-available remotely sensed data.
2. Material and Methods
The methodology used in this study is presented in Figure 1. The general procedures include the selection of the study area, data acquisition, pre-processing, statistical analyses, identification of drought-prone areas, characterization of drought onset and progression, and validation of the results.
2.1 Study area
The criteria in selecting the study area for this research include: (1) an area where periodic drought occurs with varying levels of severity; (2) an area of considerably large size considering the resolution of the datasets that will be used; and (3) an area with existing monitoring stations that have historical rainfall and water level data. The Cagayan River Basin (CRB) was selected as the study area for this research since it is the largest river basin in the country, it is frequently affected by drought events, and has the highest number of rainfall and water level stations.
The CRB is located in the north-eastern part of Luzon, Philippines, within latitudes 15° 52’ N to 18° 23’ N and longitudes 120° 51’ E to 122° 19’ E. It is the largest river basin in the country with 33 sub-basins and approximately 505 km of streams, covering approximately 27,280 km2 of drainage area , . The extent of the CRB, as shown in Figure 2, encompasses parts of nine (9) provinces. It is bounded in the north by the coastline of the Babuyan Channel where inland waters are drained, in the south by the Caraballo-Maparang Mountains, in the east by the Sierra Madre Mountains, and in the west by the Cordillera Mountains . It is relatively flat in the central plains, and hilly to mountainous in the western, southern, and eastern sides . The majority of land cover types in the CRB are croplands, range lands, and forests. It is considered vulnerable to recurring drought since it is composed of dry sub-humid agricultural areas with seasonal aridity .
2.2 Data Acquisition and Pre-processing
2.2.1 Watershed, sub-basins and HRU delineation
The watershed delineation and its sub-basins were extracted using the 90-meter Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM) of the study area using the ArcSWAT extension tool in ArcGIS®. Hydrologic response units (HRUs) were generated by combining the soil, slope, and land cover maps of the CRB. The slope map of the CRB was derived from SRTM DEM, while the most recent soil type map published in 2003 was obtained from the Bureau of Agricultural Research Spatial Analysis and Information Systems Laboratory for Research and Development (BARSAIL for R&D) of the Department of Agriculture (DA). The general land cover map of the CRB generated by J. Principe  was produced by performing Maximum Likelihood Classification (MLC) to mosaicked Landsat 7 ETM+ images acquired in February and March, 2009. The HRUs were considered as the minimum mapping unit, defining the smallest level of analysis performed in this study. Since this research deals with the use of Vegetation Indices (VIs), 453 agricultural HRUs, 444 range land HRUs, and 354 forest HRUs, totaling 1,251 HRUs were retained for analysis.
2.2.2 MODIS Vegetation Indices (VIs) and VI anomalies
Vegetation Indices (VIs) are used in many applications, including drought monitoring and assessment. Specifically, MODIS Level 3 VI products containing Normalized Difference Vegetation Index (NDVI) and Enhanced Vegetation Index (EVI) images can be used to assess the spatial and temporal dynamics of vegetation , which may be associated with the occurrence of drought events. These indices exhibit strong relationships with rainfall data, as presented in various studies , . Table 1 shows the descriptive parameters of the MODIS VI datasets used.
The NDVI is computed as the ratio of the difference and sum of near infrared (N) and red (R) reflectance values  for each pixel, i:
EVI incorporates the blue reflectance band to reduce atmospheric effects and correct soil background :
where C1 and C2 are the coefficients of the aerosol resistance term, which uses the blue band (B) to correct aerosol influence on the red band, L is the canopy background brightness correction factor, and G is the gain factor. For MODIS, the values used are L=1, C1=6, C2=7.5, and G=2.5 .
In order to better relate the VIs with drought occurrences, VI anomaly data were calculated and used as effective indicators in assessing drought conditions and mapping of drought-risk areas -. The NDVI and EVI anomalies were calculated in this study using the equation below:
pixel on the jth day of the year, respectively; and VImean,i,j is the average of the VI values for the ith pixel on the jth day of the year, applied for the whole time series . Using the online USGS Global Visualization Viewer (GloVis) tool, a total of 262 MOD13Q1 16-day VI images from February 2000 to June 2011 were downloaded. The NDVI and EVI layers were extracted and transformed from sinusoidal projection to WGS-84 geographic coordinate system using the MODIS Conversion Toolkit of ENVI®. The layers of each VI type were then stacked and the Savitzky-Golay filter was applied to improve the quality of the images , . The means of the NDVI, NDVI anomaly (NDVIA), EVI and EVI anomaly (EVIA) were calculated for each HRU and the values were tabulated in spreadsheets to produce the time series datasets. Homogeneity tests were performed in order to ensure that no major shift in the values with respect to time occurred for the whole time series and only individual trends were present in the data .
2.2.3 GRACE Terrestrial Water Storage Changes (TWSC)
Terrestrial Water Storage (TWS) refers to the total amount of water stored in the surface and sub-surface  of land regions. It is a vertically integrated measure of groundwater, soil moisture, snow and ice, permafrost, surface water, and wet biomass or vegetation water , , over a given period of time. The most dominant factors in TWS variability are seasonal groundwater changes (ΔGW) and soil moisture changes (ΔSM) . Thus, the simplified equation for Terrestrial Water Storage Changes (TWSC) is:
The GRACE twin satellites, a joint mission of NASA and the German Aerospace Centre (DLR) launched in March 2002, aims to measure small monthly variations in the Earth’s gravity field . The processing of raw spherical harmonic coefficients into Level 3 monthly equivalent water thickness enables reliable detection of global mass distribution and spatio‐temporal variations in TWS . Various studies verified the correspondence of GRACE TWSC with real groundwater data, as well as with model-simulated hydrologic parameters in large basins , -. For the present study, monthly GRACE TWSC data processed by Sean Swenson with the support of NASA’s MEaSUREs Program , , were downloaded from the Jet Propulsion Laboratory’s TELLUS website. The descriptive parameters of the data are shown in Table 2. The TWSC values were expressed as the difference between the masses for each month and the average from January 2003 to December 2007 based on the documentations of the processing conducted by Swenson , . The values for the two (2) tiles covering CRB were scaled and the monthly TWSC time series data was generated (Figure 3). Homogeneity tests were also performed on the TWSC time series.
2.2.4 Ground-based rainfall and water level data
In-situ rainfall (RF) and water level (WL) measurements were gathered from the Climate Data Section of the Philippine Atmospheric, Geophysical and Astronomical Services Administration (PAGASA-CDS) and the Bureau of Research Standards of the Department of Public Works and Highways (DPWH-BRS), respectively. Lists of drought events and dry spells from 1950 to 2011 were also obtained from PAGASA-CDS. Daily rainfall data in millimeters were available for three (3) weather stations from 2000 to 2009, while daily water level data in meters from three (3) gauging stations were provided from 2001 to 2006. The data were averaged to obtain 16-day and monthly RF and WL values to form to the temporal resolution of the MODIS VI and GRACE TWSC images used for the study.
2.3 Statistical Analyses
2.3.1 Running averages cross-correlation analysis
To smooth out short-term and random variations in the individual time series, 3-period running averages were first applied to each time series. This process involves the averaging of the values per 3 successive months to highlight the temporal trends in the data. The resulting time series datasets were then subjected to cross-correlation analysis to examine whether the VI anomalies and TWSC values are statistically related to rainfall and water level, and can therefore be used in assessing drought events. The cross-correlation between two (2) time series having the same time step and duration is the product-moment correlation, which is a function of time lag :
where cuu(0) and cyy(0) are the sample variances of the two (2) time series, ut and yt, and cuy(k) is the cross-covariance function. The cross-correlation value is simply the correlation coefficient with the second parameter lagging in time. The values range from -1 to 1, where -1 indicates a strong inverse relationship, 1 means a strong direct relationship, and 0 signifies that the parameters are not related.
All the running-averaged time series pairs were subjected to cross-correlation analysis implemented for 1-year, 2-year, and whole time ranges to check for short-term and long-term correlation trends. The correlation coefficients and corresponding time lags of the first positive peak in the values were analyzed to check for indications of drought events.
2.3.2 Lag-normalized correlation analysis between TWSC and VI anomalies
To be able to integrate VI anomalies and TWSC per HRU into a single index that may be used in drought analysis, lag-normalized correlation was developed, which can be computed by dividing the correlation coefficients (ruy) of VI anomalies and TWSC (i.e. time series pair, k) by the time lag (laguy) plus 1. The lag should be plus 1 to avoid indeterminate solutions, as shown in Equation 6:
When the correlation coefficients are high and the time lags are low, the values that will be computed for the lag-normalized correlation will also be high, which indicates higher probability of an incoming or current drought event. Otherwise, the lag-normalized correlations will be moderate to low corresponding to lesser chance of drought occurrence for that particular period.
2.4 Identification of Drought-prone Areas
The calculated 1-year and 2-year TWSC-NDVIA and TWSC-EVIA lag-normalized correlation values for the HRUs were expressed in deciles to consistently scale the data in all maps based on the frequency of occurrence of each lag-normalized correlation sum. The HRUs within the 8th, 9th, and 10th decile were considered as drought-prone. In total, four (4) sets of drought-prone areas were identified by this technique based on different VI anomaly types (i.e. NDVI, EVI) and time series intervals (i.e. 1-year, 2-year). The corresponding drought-prone barangays were then mapped based on the identified drought-prone HRUs.
2.5 Characterization of Drought Onset and Progression
Drought onset and progression may be characterized by performing lag-normalized correlation analysis between TWSC and VI anomalies using running 12-month time series. However, it is necessary to assess the applicability of this method in characterizing the onset and progression of drought events to identify whether it has a potential value for early warning applications.
The cross-correlation analysis between the 3-month running using a sliding time range of 12 months instead of the fixed, non-overlapping 1-year and 2-year ranges. The 12-month sliding time series covering years 2002 to 2011 employed 1-month time steps to assess whether the method would be effective in understanding the behavior of the lag-normalized correlation values before, during, and after the actual drought periods in the CRB. Using two (2) test HRUs (i.e. HRU 8-1 and HRU 28-1) which are both drought-sensitive, proximate to rainfall and water level stations, and within each of the two (2) TWSC grids, the correlation coefficients and time lags of the first positive peaks in the cross-correlation results were identified and lag-normalized correlations were computed.
Aside from the correspondence between observed and reported drought occurrences, supplementary datasets on the quarterly volume of crop production available from 2000 to 2012 were used for further validation of the results. The crop production data in the identified drought-prone provinces were evaluated to assess any matches to the trends in the lag-normalized correlations calculated previously.
3. Results and Discussion
3.1 Statistical Analyses
3.1.1 Running averages cross-correlation analysis
The results of the running averages cross-correlation analysis between RF-TWSC, WL-TWSC, RF-VI anomalies, WL-VI anomalies, and TWSC-VI anomalies are presented and discussed below.
Rainfall, water level, and TWSC
From the results of the analysis, it was evident that RF and TWSC exhibit high positive correlation for all weather stations since excessive or deficient RF can lead to a corresponding increase or decrease in TWSC. For the 1-year time series, the year 2002 obtained the highest correlation coefficient of 0.932 with a 1-month time lag, as seen in Table 3. The 2002-2003 time series has the highest correlation of 0.826 with a 1-month time lag amongst all 2-year intervals. From historical records, severe droughts from May 2002 to March 2003 and June 2009 to June 2010, moderate events from October 2006 to July 2007 , , and localized dry spells in 2005  occurred in the Philippines, all of which specifically affected the provinces covering the CRB. The findings were all prominent for the recorded drought periods since, as rainfall deficiency persists, the decline in groundwater storage follows, therefore increasing the correlation coefficients and lowering the time lags.
Similarly, WL and TWSC gave very strong correlations for all stations, mostly for years 2002 and 2003, with the highest value of 0.983 for year 2002. The time lags were also better than RF-TWSC, mostly having only a 0 to 1 month response delay. Years with high coefficients and low time lags corresponded to drought periods since the relationship between groundwater and water level becomes stronger as extended periods of insufficient rainfall occur.
Rainfall, water level, and VI anomalies
Since the vegetation types and conditions vary across the HRUs in the CRB, the result of the cross-correlations between rainfall/water level and VIs gave inconsistent trends. This can be attributed to the insensitivity of VIs in detecting water shortage when vegetation is senescent or when the vegetation coverage is low, as well as other limitations related to the seasonality of vegetation . The periods of low rainfall and water levels positively correspond to the gradual decrease in vegetation health; and while this extends for a longer period, the correlation becomes stronger. Based on , lag times between rainfall and VI anomalies range from 0 to 2 months since rainfall deficiency can only affect plant growth after 1 to 12 weeks. This trend was consistent for both NDVIA and EVIA, having 0 to 2 time lags and high correlation coefficients as can be seen in Figure 4, highlighting their potential in the identification of drought events.
Similarly for the WL and VIs, the peak of correlation coefficients with 0- to 1-month time lags were also detected for the years 2002, 2003, and 2006, with 2002 having the highest coefficient of 0.802. 2009 was also identified as a period with high correlation (0.783) but is not consistent for all stations, unlike the other years, which may be due to the strong typhoon ‘Pepeng’ that occurred in October 2009  that brought excessive rains in the study area and disrupted the patterns.
TWSC and VI anomalies
It was distinctly observed that when the first positive peak in the cross-correlation between TWSC and VI anomalies consistently has a lower time lag, there is a recorded drought occurrence. Time lags become shorter when extended periods of dryness occur because the gradual decrease in the TWSC correlates better with the slow decline of vegetation health expressed as VI anomalies. This trend was observed between TWSC and VI anomalies for the years 2002, 2003, 2006, 2007, and 2010, which are all periods of drought. Figure 5 shows the 1-year TWSC-EVIA correlation and time lag maps in the CRB.
From the maps, it can be seen that higher correlation coefficients above the critical value 0.576 (i.e. 95 % confidence interval from the Pearson’s correlation coefficient table) and lower time lags (i.e. < 3 months) were associated with the droughts that occurred especially in 2002-2003 and 2010. High coefficients with low time lags primarily occurred in the central part of the basin during severe droughts, where the terrain is flat and the land cover types are mostly agricultural and range lands. During localized and moderate droughts, some areas covered by existing irrigation systems obtained low correlations.
It was evident that both the coefficients and time lag values should be considered in monitoring drought conditions since there are periods and locations where the coefficients alone cannot detect droughts. Given that TWS is highly associated with ground water storage levels and the amount of rainfall, higher correlations between TWSC and VI anomalies can be an indicator of an impending hydrological drought. It was also observed that even though both NDVIA and EVIA successfully identified the drought periods, TWSC-EVIA exhibited better correspondence with known drought events, temporally and spatially.
3.1.2 Lag-normalized correlation analysis between TWSC and VI anomalies
In order to integrate the correlations and lag times of TWSC and VI anomalies obtained from Section 3.1.1, the lag-normalized correlations of TWSC and VI anomalies were assessed for each HRU. Since the trend shows that high correlation coefficients with low time lags between TWSC and VI anomalies correspond to drought periods, their ratio can serve as an effective indicator of drought events and a basis for mapping drought-prone areas in the CRB. The results significantly improved for TWSC-NDVIA, although TWSC-EVIA still gave better results. For the 1-year and 2-year time series of TWSC-EVIA, the drought-affected areas became more pronounced not only in severe occurrences, but even for years with reported localized and moderate droughts. The progression of drought events can also be seen from the 2-year time series trends based on the increase of drought-associated areas shown in Figure 6, since most drought events encompass two successive years.
The drought that peaked in 2002 was proven to be the most extensive and severe occurrence in the CRB based on coverage of the affected areas, followed by the event in 2010. It was also apparent that the peak of droughts in 2006-2007 occurred in 2007. In fact, a second, separate drought event happened from June to July 2007, because of an abnormality in the transition between El Niño and La Niña . Similarly for 2010, the 2010-2011 lag-normalized correlations showed larger drought-affected areas compared to 2009-2010, due to the occurrence of typhoons during the last quarter of 2009 which most likely affected the trend in the correlations.
After analyzing the lag-normalized correlation trends for all HRUs separated per land cover type, per slope, and per soil type, it was found that range lands, located in areas with 0-17 % slope composed of loam, silty loam, or clay loam soils are the most sensitive to drought. This means that these types of HRUs respond to changes in the amount of rainfall with the lowest time lags and with higher correlation coefficients. The high correlation and low time lags between TWSC and VI anomalies can be attributed to the land cover’s dependence on natural sources of water , , warmness and dryness of the air and land , and the water-holding capacity of the soil .
3.2 Identification of Drought-prone areas
Drought-prone areas were identified using the sum of 1-year and 2-year lag-normalized correlation values between TWSC and VI anomalies. The areas identified using TWSC-NDVIA were spatially contiguous and mostly in the central portion of the river basin; while those identified using EVIA consist of additional small and distributed HRUs mostly on the western side of the CRB, covering range and agricultural lands with slopes of 0-17 %. It was observed that EVI is more responsive to drought than NDVI for range lands due to its sensitivity to soil background.
Between NDVI and EVI anomalies, the latter generally gave more consistent and realistic results based on correct association with different land cover types, especially for the 2-year time series. A higher percentage of drought areas detected using EVI anomalies matched the previously identified drought-sensitive areas. The summary of the identified drought-prone barangays in the CRB using 2-year lag-normalized correlations of TWSC and EVIA is presented in Table 4 and the corresponding map is shown in Figure 7. Most of the barangays identified as drought-prone are located in the provinces of Cagayan and Isabela. These provinces were the most affected by drought events based on current historical records , .
3.3 Characterization of Drought Onset and Progression
The statistical analyses performed identified historical drought events in the CRB from 2000 to 2011, both spatially and temporally. The cross-correlation analysis of paired parameters exhibited strong relationships between each other, which intensified during extended dry periods. Lag-normalized correlations of TWSC-VI anomalies depicted trends in the values associated with spatial characteristics and distribution of the affected areas. The test HRUs used, HRU 8-1 and HRU 28-1, were analyzed to confirm the potential of the method for drought onset and progression monitoring.
After analyzing the plots of TWSC-NDVIA (Figure 8) and TWSC-EVIA correlation coefficients, time lags, and running 12-month lag-normalized correlations for the test HRUs, it was observed that the correlations and time lags of severe droughts that peaked in the 3rd quarter of 2002 and first quarter of 2010 obtained high correlations and low time lags consistently for both VI types. Although the correlations were also high during the moderate drought in 2006 to 2007, the time lags observed were also high because of the shorter temporal extent of the drought event that occurred. All other time intervals showed variable peaks and troughs in the correlations and time lags, which can mainly be attributed to the occurrence of excessive rains in short intervals due to typhoons.
The drought event that peaked in January to February 2007 started to see a gradual increase in the lag-normalized correlations of TWSC-NDVIA 1 month after the identified beginning of the drought period (i.e. 2 months before the peak) for both test HRUs. For TWSC-EVIA, test HRU 8-1 in Cagayan obtained a 2-month delay (i.e. 1 month before the peak), while no delay was observed for the test HRU 28-1 in Isabela (i.e. 3 months before the peak). Similarly, the drought event that peaked in January to March 2010 exhibited a gradually increasing trend for both test HRUs, but with several intermediate fluctuations, 2 months after the specified start of the drought period for the TWSC-NDVIA (i.e. 5 months before the peak), and 4 months for the TWSC-EVIA (i.e. 3 months before the peak).
Even though there were time delays on the detection of the drought’s onset, the potential use of the running lag-normalized correlation method for early warning is still highly relevant. This is because all observations showed that drought can be detected by looking into the trend of the lag-normalized correlations months before the peak period of the phenomenon; specifically about 1 to 5 months before the peak of the drought event, depending on its severity and temporal extent. As also observed from the plots, the progressing drought is indicated by the continuous, gradual increase in the trend or approximately constant level of the running lag-normalized correlation values, disregarding the minimal troughs and peaks in between. When the values start to decrease, the amount of rainfall can be said to be going back to normal, thus indicating the end of the drought phenomenon.
The trend lines can also be used to detect the potential occurrence of droughts earlier; nevertheless, the lag-normalized correlation plots must be filtered and processed first, especially EVIA, to ensure that the changes in the values can only be associated to drought episodes. However, the patterns present in the lag-normalized correlations of TWSC-VIA can already be used to characterize the onset and progression of drought events in the CRB.
The quarterly volume of production of palay (i.e. rice) and corn in drought-prone provinces were assessed to show whether the trends correspond with the actual drought periods and other detected fluctuations in the results of running lag-normalized correlations between TWSC and VI anomalies. The plot of the crop production for the province of Isabela is shown in Figure 9.
In most instances, relatively lower than normal volumes of palay and corn production for each province coincided with the periods of moderate to severe drought. There were declines in the volume of production of both palay and corn, not only during the actual and identified drought events in the CRB (i.e. 2002 to 2003, 2006 to 2007, 2009 to 2010), but also for the other intervals associated with local droughts (i.e. 2005) and typhoons/floods, as suggested by the high lag-normalized correlations. In addition, it was observed that in some cases, the decline in crop production was reflected only in the next quarter because of the differences in the type and growth status of the crops during the event.
The drought in 2009 to 2010 led to the greatest decrease in crop production in all provinces, followed by the one that occurred in 2002 to 2003; while the other drops in values were relatively minimal. It was also confirmed that Isabela and Cagayan are the most drought-susceptible provinces, since aside from having the largest drought-prone areas, the relative decrease in the volume of crop production during the drought events were also the highest. With this additional validation, it can be concluded that running lag-normalized correlation analysis between GRACE-derived TWSC and MODIS VI anomalies can be used as an effective indicator to characterize water-related events such as droughts, which in turn affect crop production.
This research aimed to assess the applicability of integrating remotely sensed data and in-situ measurements of rainfall and water levels for the identification of drought-prone areas in a local river basin in the Philippines. Specifically, it intends to develop a method of drought assessment by utilizing freely-available remotely sensed data.
The findings showed that combining GRACE-derived TWSC and MODIS VIs using the methodology presented in this research can provide a means for the characterization and assessment of drought occurrences in the Cagayan River Basin. The method developed using lag-normalized correlation analysis between TWSC and VI anomalies can be used in monitoring and characterizing the onset and progression of drought events, as well as in the identification of drought-prone areas. The reliability of the method to assess drought events was verified through the strong correlation of TWSC and VI anomalies with in-situ rainfall and water level measurements. There were observed response delays between variations in TWSC and VIs from precipitation deficiencies, which were found to range from 0 to 3 months during drought periods. The parts of the CRB located in flat areas (0-17 % slope), having natural vegetation cover (e.g. range lands), and with loam, silty loam, or clay loam soil types are the most drought-sensitive areas and can potentially be used in monitoring drought. In addition, it was identified that most of the drought-prone barangays derived from the lag-normalized correlations between TWSC-EVIA are located in Cagayan and Isabela, with areas of 2,596 ha or 68.21 % of the total land area of Cagayan, and 3,751 ha or 56.89 % of the total land area of Isabela. These findings are important for various stakeholders in regard to proper resource allocation, planning, and decision-making related to drought.
The combined use of GRACE TWSC values and MODIS VI images can be considered as a promising tool in assessing the temporal and spatial extent of droughts, especially when incorporated with other ground-measured hydro-meteorological parameters such as rainfall and water level. Since monitoring schemes usually depend on precipitation-derived indices alone, the use of RS-derived indicators demonstrated in this study is considered a practical means of drought assessment and monitoring, considering that the datasets used are free, relatively up-to-date, and have large spatial coverage.
This study was carried out with the financial support of the Department of Science and Technology (DOST), through the HRDP Scholarship Grant under the Philippine Council for Industry, Energy and Emerging Technology Research and Development (PCIEERD). The authors would also like to acknowledge the work of Sean Swenson, through the NASA MEaSUREs Program for the processed GRACE land data, as well as the MODIS Land group through the LPDAAC at the EROS Data Center for the MODIS data products used in this research.
 DA, DENR, DOST and DAR. (2012, Jan. 06). “The Updated Philippine National Action Plan to Combat Desertification, Land Degradation, and Drought FY 2010-2020,” 2010. [Online]. Available: http://www.bswm.da.gov.ph/ladaphilippines/pdf/The%20Updated%20Philippine%20NAP%20to%20Combat%20Desertification,%20Land%20Degradation%20and%20Drought%20(DLDD).pdf
 D. Wilhite, “Drought Monitoring as a Component of Drought Preparedness Planning,” in Coping with Drought Risk in Agriculture and Water Supply Systems, Dordrecht, Springer Netherlands, vol. 26, A. Iglesias, A. Cancelliere, D. A. Wilhite, L. Garrote and F. Cubillo, Eds., 2009, pp. 3-19.
 V. Alexandrov. (2010, May. 06). “Monitoring Soil Drought (A Review),” October 2006. [Online]. Available: http://www.meteo.bg/
 IPCC, “Climate Change 2007: Synthesis Report,” Cambridge University Press, New York, 2007.
 R.N. Concepcion. (2010, May. 05). “Summary Report: National Country Report on the UNCCCD Implementation (Philippines),” 2000. [Online]. Available: http://ow.ly/F9n130fl5Vv
 CAB-PAGASA,. (2012, Feb. 28). “Coping with Climate,” 2012. [Online]. Available: http://www1.pagasa.dost.gov.ph/
 R. G. De Guzman. (2010, May. 06). “Status of Drought Monitoring and Preparedness in the Philippines,” 2009. [Online]. Available: http://www.wamis.org/agm/meetings/etdret09/ETDRET_Guzman.pdf
 G. J. Perez, M. Macapagal, R. Olivaresa, E. M. Macapagal and J. C. Comiso, “Forecasting and Monitoring Agricultural Drought in the Philippines,” The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, vol. XLI-B8, pp. 1263-1269, 2016. Doi: https://doi.org/10.5194/isprsarchives-xli-b8-1263-2016
 R. Houborg, M. Rodell, B. Li, R. Reichle, R. Heim and J. Lawrimore. (2010, May. 06). “Toward Integrating Enhanced GRACE Terrestrial Water Storage Data into the U.S. and North American Drought Monitors,” 2010. [Online]. Available: http://ams.confex.com/ams/pdfpapers/161080.pdf
 G. Berhan, S. Hill, T. Tadesse and S. Atnafu, “Using Satellite Images for Drought Monitoring: A Knowledge Discovery Approach,” Journal of Strategic Innovation and Sustainability, vol. 7, no. 1, pp. 135-153, 2011.
 G. Hu, W. Ying and C. Weijia, “Drought Monitoring Based on Remotely Sensed Data is the Key Growing Period of Winter Wheat: A Case Study In Hebei Province, China,” The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, vol. XXXVII, pp. 403-408, 2008.
 A. Ghulam, Z. Li, Q. Qin and Q. Tong, “Designing of the Perpendicular Drought Index,” Environ. Geol., vol. 52, pp. 1045–1052, 2007. Doi: https://doi.org/10.1007/s00254-006-0544-2
 B. Li, M. Rodell, B. Zaitchik, R. Reichle, R. Koster and T. van Damm, “Assimilation of GRACE Terrestrial Water Storage into a Land Surface Model: Evaluation and Potential Value for Drought Monitoring in Western and Central Europe,” Journal of Hydrology, pp. 1-15, 2012. Doi: https://doi.org/10.1016/j.jhydrol.2012.04.035
 RBCO (2011, Oct. 15). “Plans and Programs of River Basin Control Office Relative to Water Resources Management and River Basin Management,” October 2010. [Online]. Available: http://www.wepa-db.net/pdf/0710philippines/8_RBCO.pdf
 NWRB, “Principal River Basins of the Philippines,” National Water Resources Board, Manila, 1976.
 ABS-CBN News (2011, Oct. 21). “The Cagayan River Basin,” October 2009. [Online]. Available: http://news.abs-cbn.com/research/10/23/09/cagayan-river-basin
 BRS-DPWH, “Water Resources: Region No. 2 Cagayan Valley,” DPWH, Quezon City, 2002.
 J. A. Principe, “Exploring Climate Change Effects on Watershed Sediment Yield and Land Cover-based Mitigation Measures Using SWAT Model, RS and GIS: Case of Cagayan River Basin,” M.S. thesis, University of the Philippines - Diliman, Quezon City, 2012.
 F. Ndayisaba, H. Guo, A. Bao, H. Guo, F. Karamage and A. Kayiranga, “Understanding the Spatial Temporal Vegetation Dynamics in Rwanda,” Remote Sensing, vol. 8, no. 2, 2016. Doi: https://doi.org/10.3390/rs8020129
 P. Chopra, “Drought Risk Assessment Using Remote Sensing and GIS: A Case Study of Gujarat,” M.S. thesis, International Institute for Geo-information Science and Earth Observation, Enschede, The Netherlands, 2006.
 D. Fitch, D. Stow, A. Hope and S. Rey, “MODIS Vegetation Metrics as Indicators of Hydrological Response in Watersheds of California Mediterranean-type Climate Zones,” Remote Sensing of Environment, vol. 114, p. 2513–2523, 2010. Doi: https://doi.org/10.1016/j.rse.2010.05.026
 A. Huete, C. Justice and W. Van Leeuwen. (2010, May. 05). “MODIS Vegetation Index (MOD 13) Algorithm Theoretical Basis Document Version 3,” 1999. [Online]. Available: https://modis.gsfc.nasa.gov/data/atbd/atbd_mod13.pdf
 P. Thenkabail, M. Gamage and V. Smakhtin, “The Use of Remote Sensing Data for Drought Assessment and Monitoring in Southwest Asia,” Research Report 85, International Water Management Institute, Colombo, Sri Lanka, 2004.
 N. Hammouri and A. El-Naqa, “Drought Assessment Using GIS and Remote Sensing in Amman-Zarqa Basin,” Jordan Journal of Civil Engineering, vol. 1, no. 2, pp. 142-152, 2007.
 H. Murad and A. Saiful Islam, “Drought Assessment using Remote Sensing and GIS in North-west Region of Bangladesh,” in 3rd International Conference on Water & Flood Management, Dhaka, 2011.
 J. Bojanowski, W. Kowalik and Z. Bochenek, “Noise Reduction of NDVI Time-Series: A Robust Method Based on Savitzky-Golay Filter,” Polish Association for Spatial Information: Annals of Geomatics 2009, vol. VII, no. 2(32), pp. 14-23, 2009.
 Y. Huang, D. Jiang, D. Zhuang, H. Ren and Z. Yao, “Filling Gaps in Vegetation Index Measurements for Crop Growth Monitoring,” African Journal of Agricultural Research, vol. 6, no. 12, pp. 2920-2930, 2011.
 X. Addinsoft. (2013, Jan. 16). “Homogeneity tests for time series,” 2013. [Online]. Available: http://www.xlstat.com/en/products-solutions/feature/homogeneity-tests-for-time-series.html
 K. Haile, Estimation of Terrestrial Water Storage in the Upper Reach of Yellow River, “ M.S. thesis, University of Twente, 2011.
 G. Strassberg, B. Scanlon and M. Rodell, “Comparison of Seasonal Terrestrial Water Storage Variations from GRACE with Groundwater-level Measurements from the High Plains Aquifer (USA),” Geophysical Research Letters, vol. 34, no. 14, 2007. Doi: https://doi.org/10.1029/2007gl030139
 M. Rodell. (2010, May. 05). “Remote Sensing of Terrestrial Water Storage and Application to Drought Monitoring,” [Online]. Available: https://www.drought.gov/drought/sites/drought.gov.drought/files/media/resources/workshops/20080206_NIDIS_Knowledge_Remote_Sensing_CO/matt_rodell.pdf
 R. Schmidt, P. Schwintzer, F. Flechtner, C. Reigber, A. Guntner, P. Doll, G. Ramillien, A. Cazenave, S. Petrovic, H. Jochmann and J. Wunsch, “GRACE Observations of Changes in Continental Water Storage,” Global and Planetary Change, vol. 50, pp. 112-126, 2006. Doi: https://doi.org/10.1016/j.gloplacha.2004.11.018
 M. Rodell, J. Chen, H. Kato, J. Famiglietti, J. Nigro and C. Wilson, “Estimating Groundwater Storage Changes in the Mississippi River Basin (USA) using GRACE,” Hydrogeology Journal, vol. 15, no. 1, pp. 159-166, 2007. Doi: https://doi.org/10.1007/s10040-006-0103-7
 S. Swenson, P. Yeh, J. Wahr and J. Famiglietti, “A Comparison of Terrestrial Water Storage Variations from GRACE with In Situ Measurements from Illinois,” Geophysical Research Letters, vol. 33, pp. 1-5, 2006. Doi: https://doi.org/10.1029/2006gl026962
 B. Zaitchik, M. Rodell and R. Reichle, “Assimilation of GRACE Terrestrial Water Storage Data into a Land Surface Model: Results for the Mississippi River Basin,” Journal of Hydrometeorology, vol. 9, no. 3, p. 535–548, 2008. Doi: https://doi.org/10.1175/2007jhm951.1
 M. Leblanc, P. Tregoning, G. Ramillien, S. Tweed and A. Fakes, “Basin-scale, Integrated Observations of the Early 21st Century Multiyear Drought in Southeast Australia,” Water Resources Research, vol. 45, pp. 1-10, 2009. Doi: https://doi.org/10.1029/2008wr007333
 J. Chen, C. Wilson, B. Tapley, L. Longuevergne, Z. Yang and B. Scanlon, “Recent La Plata Basin Drought Conditions Observed by Satellite Gravimetry,” Journal of Geophysical Research, vol. 115, pp. 1-12, 2009. Doi: https://doi.org/10.1029/2010jd014689
 F. Landerer and S. Swenson, “Accuracy of Scaled GRACE Terrestrial Water Storage Estimates,” Water Resources Research, vol. 48, pp. 1-11, 2012. Doi: https://doi.org/10.1029/2011wr011453
 C. Chatfield, “The Analysis of Time Series, An Introduction” (6th Ed.), New York: Chapman & Hall/CRC, 2004.
 G. P. Yumul Jr., C. B. Dimalanta, N. T. Servando and F. D. Hilario, “The 2009-2010 El Niño Southern Oscillation in the Context of Climate Uncertainty: The Philippine Setting,” Philippine Journal of Science, vol. 139, no. 1, pp. 119-126, 2010.
 C. M. Reyes, S. N. Domingo, C. D. Mina and K. G. Gonzales, “Climate Variability, SCF, and Corn Farming in Isabela, Philippines: a Farm and Household Level Analysis,” Philippine Institute for Development Studies, Makati, 2009.
 R. Lasco, F. Pulhin, P. Jaranilla-Sanchez, K. Garcia and R. Gerpacio, “Mainstreaming Climate Change in the Philippines. Working Paper Nr 62,” World Agroforestry Centre, Los Banos, 2008.
 T. Udelhoven, M. Stellmes, G. del Barrio and J. Hill, “Assessment of Rainfall and NDVI Anomalies in Spain (1989–1999) Using Distributed Lag Models,” International Journal of Remote Sensing, vol. 30, no. 8, p. 1961–1976, 2009. Doi: https://doi.org/10.1080/01431160802546829
 C. Bhuiyan, “Desert Vegetation during Droughts: Response and Sensitivity,” Beijing, 2008.
 C. R. Townsend, M. Begon and J. L. Harper, “Essentials of Ecology”, New York: John Wiley and Sons, 2009.
 J. Rahn, “Making the Weather Work for You: A Practical Guide for Gardener and Farmer”, Charlotte, Vermont: Garden Way Publishing, 1979.