Spatial-Temporal Patterns of Agricultural Drought in Upper Progo Watershed Based on Remote Sensing and Land Physical Characteristics

Agricultural drought is alarmed by meteorological drought characterized by lower year-to-year rainfall. Under long period and continuous water deficits, plants may demonstrate stress symptoms and wilt or die. Furthermore, agricultural drought leads to crop failures and threaten the food security of an area. Progo Hulu sub-watershed is a major agricultural area in Temanggung Regency. Spatial-temporal pattern-based information about agricultural drought can be a basis for decision making in drought mitigation. This study aims to analyze spatial and temporal distribution patterns of drought, analyze the physical characteristics of land and their influence on drought pattern, and establish a prediction model of drought distribution patterns based on four physical characteristics of the land. Landsat 8 imagery is used to determine the spatial and temporal patterns of agricultural drought in Upper Progo watershed using an improved Temperature vegetation Dryness Index (iTVDI). Slope, land use, landform, and soil texture are the physical characteristics of land as the variables to determine the most influential factor of drought pattern. They are analyzed using multiple regression analysis techniques. Pixel samples are obtained through purposive sampling method based on land units. The results reveal that the spatial-temporal distribution of agricultural drought occurs rapidly on the slopes and foothills of Sumbing and Sindoro. These areas have the highest average value of the iTVDI index. Agricultural drought extends gradually in line with the number of days without rainfall. Landform is a physical characteristic that most influences the distribution of agricultural drought. The established model by utilizing four variables of physical characteristics generates an average value which almost similar to the iTVDI value produced by remote sensing data. The model can be useful to estimate drought distribution based on the number of days without rainfall. Keywords— agricultural drought; spatial-temporal analysis; spatial modeling.


I. INTRODUCTION
Drought occurs in almost all climatic zones regardless the precipitation level of a region that is associated with a decrease in the precipitation level over some time, such as in a season or a year [1]. Nevertheless, the characteristics of drought may vary significantly from region to region [2]. Drought is a natural phenomenon that occurs when rainfall decreases significantly below the normal level leading to serious hydrological imbalances and lower land production system (UN Convention to Combat Drought and Desertification, 1994) [1]. Drought is a type of complex hazard caused by disruption of water balance and impacts on agriculture, ecology, and socio-economic [3]. Such a hazard has typical characteristics as it progresses gradually and undergoes propagation. The level of drought preparedness and mitigation relies on information on the beginnings and propagation of drought in space and time [4]. Repeated occurrence of drought will become more frequent and severer due to an increase in water demand in line with population growth and limitations and uncertainties of water supply due to climate change [5].
Drought propagation varies depending on the physical conditions of the site. These conditions can be determined through a geomorphological survey. The drought susceptibility survey is a concept departed from the consideration that geomorphological and geological-related factors, pedology and vegetation characteristics have a significant effect on the large-scale drought distribution of the climatological factor itself. Potential drought-prone areas can be estimated and an action plan can be formulated and to be implemented when the hazard takes place [6]. Hao et al. in 2015 [7] also suggested that drought is complicated because it is linked to the variation of landform, climate, and vegetation types.
Agricultural drought is a type of drought hazard that involves the meteorological, hydrological, and socioeconomic field at once [1]. It is associated with a reduction of the moisture content in the soil; hence it fails to fulfill water needs for plants in a certain period. It occurs after the symptoms of meteorological drought [8]. During an agricultural drought, the moisture content of the soil will continually decrease steadily due to drainage, plant needs, and evapotranspiration to the point of the plant roots are unable to absorb moisture, and eventually, they will wilt [9]. Areas affected by agricultural drought expand gradually as a symptom of soil moisture stress on plants and propagate slowly [10]. Improper management towards agricultural drought will lead to hydrological and socioeconomic drought. Moreover, agricultural drought menaces the food security of a region [11].
The agricultural sector is the foremost sector suffered by the onset of drought due to crop dependence on water resources and soil moisture during different stages in plant growth [12]. Drought affects vegetation physiologically and morphologically, but there is a response pause in vegetation since drought is the accumulation of gradual processes. The vegetation growth situation during observation is a combined effect of the previous conditions [13]. The low precipitation value is exacerbated by the high evapotranspiration resulting in low soil moisture [1].
Agricultural drought can be identified using remote sensing technology. The methods and techniques in remote sensing are widely utilized in many aspects of drought, such as assessment invulnerability, damage, and warning systems [14]. Remote sensing has been used rapidly and extensively over the past three decades because it provides a less expensive and faster technique to calculate agricultural statistics [21]. Accuracy in temporal monitoring, spatial distribution, and severity of drought provide an effective means of reducing farmer losses [15]. Remote sensing methods provide high spatial and temporal resolutions that are capable of monitoring water stress with relatively low cost and time [16]. Indicators used to detect the symptoms of drought are withered foliage or rise in ground surface temperatures [12]. This indicator is estimated through an index, one of which is a remote sensing image-based index. A Normalized Difference Vegetation Index (NDVI) is the most commonly used index, but it is frequently linked with the Green index than the Humidity index. A Land Surface Temperature (LST) is a more sensitive index than water stress. Information about the distribution of LST can be useful, particularly those of the surface temperature and climate that are important in evapotranspiration, climate change, hydrological cycle, vegetation monitoring, urban climate (urban heat island/UHI) and environmental studies [22]. Remote sensing has played a significant role in monitoring land surface temperature, especially where instrumental data is not available [23]. In general, evapotranspiration can extensively control the surface temperature through the balance of energy on the surface. The lower evapotranspiration causes a higher energy to increase surface temperature [17]. The combination of NDVI and LST can provide more accurate information on soil moisture information [18].
Furthermore, an improved Temperature Vegetation Dryness Index (iTVDI) is an index that utilizes a combination of NDVI and LST and air temperature data to detect an agricultural drought [16]. Rahimzadeh-Bajgiran et al. in 2012 [16] researched the use of iTVDI indicators to identify agricultural drought in Iran using MODIS imagery. It demonstrates that iTVDI is strongly affected by actual rain during the dry season and can be used to calculate the moisture content. It also indicates that iTVDI has been successfully used for monitoring water deficit in vegetation or soil. The iTVDI index has moderate accuracy with a correlation of 0.8 to the soil moisture, hence this index is appropriate to detect agricultural drought. Geospatial techniques have a major role to make and update relevant databases at relatively low cost. They are useful for analyzing and modeling the spatial distribution and dynamics of biophysical objects [24] and also for the spatialtemporal change detection analysis [25].
The Progo upstream watershed is a critical agricultural area in Temanggung Regency. The 2015 dry season brought severe drought on many villages in this regency since rain did not occur for a long period. Consequently, the Local Government of Temanggung declared a drought emergency based on the Decree of Temanggung Regent No. 337/2015 on a 90-day drought emergency dated from August 1st, 2015 [19].
The basic determination of drought emergency is due to many villages experiencing drought. As a result, people endured clean water shortage and unproductive agricultural land [20]. This study aims to build a prediction of agricultural drought propagation using remote sensing data and land physical condition in the Upper Progo watershed.

II. MATERIAL AND METHOD
Several measures were conducted for achieving the research objectives. First, the analysis of the duration of the rainfall deficit is used as a basis for starting and ending the analysis through the remote sensing imagery. The duration of the rainfall deficit was analyzed by using the water balance chart between the monthly rainfall and the evapotranspiration rate calculated by using the Thornwaite formula. Drought occurs when the evapotranspiration rate is higher than the amount of rainfall and ends if precipitation exceeds the evapotranspiration rate. In this case, precipitation was measured based on monthly rainfall in 2015.
Second, the iTVDI map was done based on the time series data of NDVI, LST, and air temperature by drought duration. The air temperature was determined by utilizing SRTM image using the lapse rate method with Mock formula, in which a 100m increase in elevation typically produces a 0.6oC drop in air temperature. Furthermore, air temperature becomes the surface temperature calibration; hence iTVDI has the equation as follows: (1) Ts is the land surface temperature (LST) observed at certain pixel, while Ta is the observed air temperature calibrated with the DEM data. (Ts -Ta)obs is the land surface temperature (LST) minus the air temperature during the observation that reflects the evapotranspiration rate. (Ts -Ta)max is the dry edge, while (Ts -Ta)min is the wet edge. These edges are obtained from the scatter plot between Ts -Ta and NDVI. This scatter plot is presented on Fig. 1.
Third, the analysis of the physical condition of the land was done using the physical parameters of the land to drought. The parameters estimated to have substantial contribution toward the soil moisture deficit have a higher score than those with a less significant contribution. Residual hills D2 2 A landform map was made using contour data from Landsat 8 Satellite Imagery and geological map. The landform classification is presented in Table 1. Slope map was also created using SRTM image with 30 m spatial resolution. The slope class is presented in Table 2. Information on the spatial distribution of soil texture was obtained from the findings reported by Euis Etty Alhakim (2013). The soil texture in the study site was divided into three classes as presented in Table 3. Subsequently, the land use map was made using high-resolution imagery from Google Earth, Landsat satellite imagery, and RBI maps. Land use classification according to RBI map is presented in Table 4.   The determination of physical condition that has the most influential impact on soil moisture was analyzed using multiple regression analysis techniques. The Standardized Coefficient reflects the influence of each independent variable. Its effect on the soil moisture scored each physical condition and then converted into raster data. Purposive sampling technique was done to obtain a physical land condition in the land unit map. Subsequently, information about the score of the four land physical conditions and the iTVDI value of each recording date was determined. The samples were analyzed by using multiple regression analysis on SPSS 17 with iTVDI value as the dependent variable and the physical conditions of land as independent variables.
Finally, the prediction of soil moisture propagation was created. It is a method to predict the condition of soil moisture on then-day after the occurrence of rainfall deficit. Prediction is done by utilizing the regression analysis statistical result on the physical conditions of the land. Regression analysis is done on every recording date of Dry edge Wet edge

III. RESULT AND DISCUSSION
The Progo upstream watershed is located in Temanggung Regency, Central Java. It lies between 109o59'45" E -110o19'49" E and 7o11'42" S -7o24'3" S. It has an area of approximately 56,842 ha (568.42 km2) or 65.29% of the total area of Temanggung Regency. The map of the study site is presented in Figure 2. A.

B. Land Physical Condition of the UpperProgo Watershed
A slope map of the Upper Progo watershed is presented in Figure 3. The west and east part of the study site are classified as steep slope class. The west part of the stratovolcanic slope is steeper towards the peak. The east part is a high hill area with a moderately steep slope. A relatively flat area is located in the middle between the hills and the stratovolcanic formation.
The volcanic landform is located in the west part of the study site. It is located precisely on the west of Progo Main River. Volcanic denudational is an inactive volcano that experiences denudation but the volcanic body is still visible through the appearance of a rounded contour line. The shape of the diluted hill terrain is characterized by a dendritic river flow pattern which indicates that the constituent rocks of this landform are homogeneous. The remaining hillside is in the form of small hills located between the plains of the foot of the mountain with annual plant vegetation cover. The map of the distribution of existing landforms in the Progo upstream watershed is presented in Figure 4. The composition of land use in the Progo upstream watershed consists of irrigated paddy field, rainfed paddy field, dryland agriculture, and garden. The total area of agricultural land in the Progo upstream watershed is approximately 82% of the total area of the study site. The large percentage of agricultural land in the Progo upstream watershed provides information that most of the local people are farmers. The distribution of land use in the Progo upstream watershed is presented in Figure 5.
Characteristics of soil texture in the Progo upstream watershed were derived from data reported by Euis Etty Alhakim (2013). The soil texture is divided into three types, namely sandy loam, sandy clay, and clay. Spatial distribution of soil texture in the Progo upstream watershed is demonstrated in Figure 6.
The beginning of meteorological droughts in the Progo upstream watershed is divided into two parts. There were five weather stations detected the onset of meteorological droughts in May, while three weather stations detected the meteorological droughts in June. The beginning and end of the meteorological drought were used as the basis for initiating Landsat 8 Satellite Imagery to identify soil moisture deficits as an indicator of agricultural drought in the Progo upstream watershed. Landsat 8 has a temporal resolution of 16 days thus the agricultural drought analysis was conducted every 16 days during May-November 2015.

C. The Green Index in the Upper Progo watershed
The green index shows the condition of chlorophyll in the foliage that represents the condition of the plant. This study used the NDVI method to determine the green index in the study site. The average NDVI value increased in four images of the beginning of the recording date. The increase in NDVI values is possibly caused by the adequacy of water supply for plant growth. It is a result of the small precipitation in June. Rainfall reaches 0 mm for all weather stations from July to October. It caused a decline in the NDVI value from the images of the recording date of August 1st to October 20th. On those dates, there is no average increase in the NDVI values that indicates the disruption in plant growth.
The decline in NDVI values indicates a reducing amount of chlorophyll that can be attributed to the availability of water quantities for plant survival. The non-productive withered plants will hinder the agricultural activities in the study site. The declining rate in the NDVI value on each recording date indicates the wilting rate of plants during the groundwater deficit.

D. The Land Surface Temperature (LST) in the Upper Progo Watershed
The land surface temperature is inversely proportional to the soil moisture of the land surface. Low land surface temperature occurs because the solar radiation is absorbed by the moisture content in the soil (soil moisture). The soil moisture-temperature will rise and turn into water vapor, and evaporation occurs. As a result, the land surface temperature will be relatively stable. Precipitation that does not occur within a certain period will lead to a water supply deficit. A decrease in soil moisture will lead to the reduction of radiation absorption by water and evaporation occurs, yet solar radiation will raise the land surface temperature. The high land surface temperature in the satellite imagery will interpret the amount of soil moisture in the study site; hence the physical land conditions can be identified in the dry condition.

E. Air Temperature in the Progo upstream Watershed
Daily temperature in the recording date was obtained from the Temanggung PHP Climatology Laboratory Station at an elevation of 660 m from sea level. It was then generalized for the entire study site using the Mock formula. The results of calculation show that the air temperature will decrease in line with the increase in elevation.

F. Difference between LST and Air Temperature
The combination of LST and air temperature is expected to improve the sharpness in the iTVDI calculation. The use of air temperature data can also represent the condition of the study site that has hilly and mountainous topography. The combination of LST and air temperature used different techniques. Since the LST value is usually higher than the air temperature value, the LST value is reduced by the air temperature value. The LST-temperature value represents the evapotranspiration rate in the study site.

G. Spatial-Temporal Soil Moisture Content
The moisture content in the study site was identified using the improved Temperature Vegetation Index (iTVDI) method. This method combines two indicators to detect the soil moisture in an area. The green index and land surface temperature corrected by air temperature data became two factors for identifying the soil moisture. The NDVI and LST -air temperature data were then analyzed using scatter plot diagrams to determine the relationship of both at each recording date. The result of the scatter plot between NDVI and LST difference -air temperature is presented in Figure 7. The NDVI value on the diagram is represented by the X-axis while the LST-air temperature difference is represented by the Y-axis.
The results of the scatter plot at each recording date have a triangular or trapezoidal shape and on the top or bottom edge of the trapezoid can be drawn lines called wet boundary and dry boundary. The wet boundary and dry boundary are used to determine the equation to be included in the general iTVDI equation. The shape of the scatter plot will affect the gradient of wet boundary and dry boundary. The dry boundary is the upper edge of the trapezoid, and the wet boundary is the bottom edge of the trapezoid. Figure 7 shows that the blue color indicates a denser point. On the diagram of May 29, 2015, to August 17, 2015, the blue color is identified at a fairly high NDVI value of about 0.8 and then shifts left as the time passes. It suggests that the green index decreases due to the absence of precipitation. The shape of the triangle in the scatter plot also increasingly evident along with the increase in duration. The set of points on the left of the diagram (low NDVI values) is widened, and the set of points on the right of the diagram (high NDVI values) is narrowed. The decrease in NDVI values is inversely proportional to the rise in land surface temperatures, which triggers the sharper shape of the triangle. The points of the diagram ultimately end at a certain value and form a pattern; hence a line can be drawn at the stop point, and an equation can be formulated.
The result of the calculation using the edge equation was applied to each pixel. Hence the value of each pixel would be different. The average iTVDI value of each image was calculated using raster statistics. The water availability is estimated to decline in the absence of precipitation for the next several months. It is indicated by the increase in the value of iTVDI from the recording date on August 1, 2015, to October 20, 2015. The average value of iTVDI increased significantly on August 1 because the east part of the study site was covered by cloud. The west part of the study site was not covered; hence it had a relatively high iTVDI value and high average value as well. The iTVDI value slightly decreased on the recording date of September 18, 2015. The highest average iTVDI value occurred on October 20, 2015. On that date, the evapotranspiration process was relatively low thus the soil moisture is also less. It was due to the rainfall duration was short, causing water availability was increasingly limited.
The green color shows a low iTVDI value, implying the evapotranspiration rate is quite high and the moisture content is also high. The yellow color indicates the moderate iTVDI values, while orange color indicates a high iTVDI value indicating low moisture content. The orange color on each recording date dominates the area around Mount Sindoro and Mount Sumbing. It indicates the very low evapotranspiration rate in those areas and also the low moisture content. Low moisture content is possibly due to intensive agricultural activity on the steep slopes.
Soil moisture propagation in the study site can be seen by observing the pixel color at each recording date in Figure 8. In general, the change in pixel color into yellowish ones starting from the image date of August 1, 2015. Previously, the green color dominates the images. The yellow color continues to spread from August 1, 2015, to October 20, 2015. It originally dominates the west part of the study site and begins spreading to the east part since August 17, 2015, and dominates the east part on October 20, 2015.

H. The Effect of Land Physical Conditions on Soil Moisture
Content Physical conditions of land as the variables in this study were slope, landform, land use, and soil texture. The four variables of the physical condition of the land were then scored according to their effects on the soil moisture condition. Subsequently, the variables were analyzed using regression analysis to find out which physical condition had the most influence on soil moisture condition. The coefficient value is directly proportional to the effect of the physical condition of the land. The higher the coefficient value in a variable indicates the greater the influence of independent variables on the dependent variable.
Constant and variable values always change according to the duration of the rainfall deficit in 2015. The tendency of constants and variables is presented in Figure 9, 10, 11, 12, and 13. The coefficient of landform typically has a higher value than other variables. The coefficient of landform has a value > 0.1 except at the data on September 18, 2015, and October 20, 2015. The slope is the second variable that has a significant role in the soil moisture conditions. Soil texture is the third variable that affects the iTVDI value, while land use has an insignificant role than other variables.

I. Soil Moisture Propagation Model based on Precipitation Deficit Duration
The soil moisture propagation model is used to predict the extent of moisture content deficit-affected area as observed from the evapotranspiration rate in the study site. The model was constructed from data of iTVDI and physical condition variables to determine soil moisture conditions. It also has the advantage, which is the capacity to predict the cloudcovered areas of the image. The iTVDI pixel value is recreated by exploiting the constant value and the coefficient of each variable for each recording date of the results of multiple regression calculation.
The blue color indicates a wet soil moisture condition (low iTVDI) and light brown indicates dry soil moisture (high iTVDI). Distribution of soil moisture classes fluctuated from the first to the 145th day of rainfall deficit. The volcanic areas of Sindoro and Sumbing have almost consistent dry soil moisture conditions, and several times in very dry and normal class. In general, changes in soil moisture conditions during the 145 days of rainfall deficit fluctuated. The first to the 49th day show a decline average iTVDI value. It can also be seen in the spatial pattern that shows the blue color is increasingly widespread to reach the highest area on the 49th day. Furthermore, in the 49th to 145th, the iTVDI values fluctuated but tended to rise, and the highest average value occurred on the 145th day of rainfall deficit. The spatial pattern also fluctuated but in general, the study site has a drier class shift as observed from the color change in Fig. 14. The agricultural drought has been identified through the spatial-temporal iTVDI index, which indicates the rapid agricultural drought occurs mainly on the slopes and at the foot of Mount Sindoro and Mount Sumbing. The soil moisture deficit extends gradually in line with the longer the days of rainfall deficit. High iTVDI values clustered in the west part of the study site while the low and moderate iTVDI values spread in the middle and east part of the study site.
Landform is the variable of physical land condition with the highest effect on the soil moisture deficit propagation. Furthermore, the second variable is the slope of an area.
The modeling of soil moisture deficit using land physical condition variables generates average value that is almost the same with the iTVDI value produced from pixel image calculation. Predictive modeling has the advantage of eliminating cloud elements in the calculation.