Rainfall Trend and Variability Over Opak River Basin, Yogyakarta, Indonesia

Rainfall intensity seems to be increasing nowadays due to climate change as presented in many studies of both global and regional scale. Consequently, cities worldwide are now more vulnerable to flooding. In Indonesia, increasing frequency of floods was reported for the past decades by The National Agency for Disaster Countermeasure (BNPB). To understand the rainfall changes, longterm trend evaluation over a specific area is then crucial due to the large variability of spatial and temporal rainfall distribution. This study investigates the homogeneity and trend of rainfall data from 20 stations over the Opak River basin, Yogyakarta, Indonesia. A long-term ground observation rainfall data whose period varies from 1979 to 2019 were analyzed. Non-parametric Mann – Kendall test was applied to assess the trend, while the magnitude was calculated using the Sen’s slope estimator. An increasing annual maximum of daily rainfall intensity was observed at four stations on a 0.95 confidence level based on the Mann – Kendall test, while the Sen’s slope estimator shows a positive trend at almost all stations. The trend of heavy rainfall frequency was also found to be significantly increased, with only one station showed a decreasing trend. Furthermore, this paper also described the spatial rainfall variability. Positive trend was mostly found during the rainy season, while the negative trend occurred during the dry season. This could pose a challenge for water resource management engineering and design, such as water supply systems or reservoir management. Understanding this phenomenon will benefit hydrologists in preparing future water resource engineering and management.


INTRODUCTION
Attention to climate change impact on our environment has been increasing recently, with rainfall trend changes as one of the most attractive issues for hydrologic studies. (IPCC 2014) had observed a clear increasing precipitation over the mid -latitude, and slight changes over the other regions. It also projected a more intense and frequent extreme precipitation over the mid -latitude land masses and wet tropical regions. Moreover, both the spatial and temporal rainfall variability is also amplified: higher intensity and more frequent rainfall have occurred during the wet season and longer drought during the dry season (Held & Soden 2006;Zhang & Fueglistaler 2019).
An extensive study using a global rainfall data series had suggested an increase in global annual maximum daily precipitation (Westra, Alexander & Zwiers 2013). Positive trends were also indicated through regional studies, such as in northern Australia (Roderick, Wasko & Sharma 2019), India (Goswami et al. 2006), and many other regions. (As-syakur et al. 2013) observed an increasing rainfall over Indonesia, especially in Kalimantan, Java, Sumatra, and Papua, from 1998 to 2010, based on the TRMM multi-satellite dataset. Although, they also reported a decreasing trend over the other islands. This conclusion was supported by (Avia 2019), which described a significant increase in rainfall intensity over Java, from 1981 to 2016, using the CMAP rainfall dataset.
On the contrary, decreasing trends in rainfall were also marked. A recent study using a regional climate model of western Maritime Continents had confidently suggested a decreasing trend of rainfall in many regions, although it also identified a slight increase of rainfall around Java and Sumatra during the rainy season (Kang, Im & Eltahir 2019). Another regional study had also shown a similar pattern, a decreasing extreme rainfall over most of South East Asia and South Pacific countries based on the data period from 1961 to 1998 (Manton et al. 2001). The negative trend of annual precipitation over Indonesia was also reported through a Mann -Kendall test on 63 rainfall datasets whose period varies from 1950 to 1997 (Aldrian 2007). Opposite ideas of rainfall trends could be understood since rainfall distribution has high variability, both temporary and spatially. Therefore, it is essential to add more evidence by investigating local rainfall variability.
A regional meteorological event, such as tropical cyclones, could also induce extreme precipitation over the affected area, as pointed out globally (Lau, Zhou & Wu 2008;Khouakhi, Villarini & Vecchi 2017) and locally (Knight & Davis 2009;Bagtasa 2017). Such an event also occurred over Java, Indonesia, when tropical cyclone Cempaka passed the southern coast of Java in late November 2017, as explained in (Samrin et al. 2019). Similar events had occurred in Bima City, West Nusa Tenggara, Indonesia, as described in (Sok 2019). The intensity of the tropical cyclones might increase in the future due to global warming (Knutson et al. 2010). Furthermore, tropical cyclone usually occurs during the rainy season in Indonesia.
Consequently, cities worldwide are now more vulnerable to urban flooding (Hassan, Nile & Al-Masody 2017;Park, Oh & Won 2020). In Indonesia, the increasing number frequency of floods was recorded for the past decades, and Java Island had the highest number of floods in 2019, as reported by The National Agency for Disaster Countermeasure (BNPB 2020). For example, an urban flood had occurred on March 11, 2020, over some regions in Yogyakarta, Indonesia, due to extreme precipitation. The heavy rainfall, whose intensity is comparable to the average monthly rainfall, lasted around two hours as recorded by an automatic rainfall recorder ( Figure 1). The flood had affected residents in the area, and they were able to capture the event (Figure 2). Local news reported inundations at some other locations around Yogyakarta as well. The flood has devastated many public facilities around the areas, although the damages could be anticipated quickly.  Despite the increasing evidence of rainfall changes around the globe, the trend over the tropical regions seems to be unclear due to the insufficient sample (Manton et al. 2001). In Indonesia, rainfall data was also scarce, even though the region would likely be affected by global warming, as it is a maritime country and located in the tropical region. Nevertheless, with the increasing available rainfall data and an improved hydrological data management system, it is now possible to investigate a long-term precipitation trend.
This work intends to clarify the rainfall variability in terms of intensity and frequency over the Opak River basin, Indonesia, using ground observation data whose period varies from 1979 to 2019. Several studies of rainfall trends on a larger scale have been carried out. However, it is still necessary to investigate and renew such research at a local scale due to a high rainfall variability both spatially and temporally. We believe that this study will be beneficial for future water resource management and design.

MATERIAL AND METHOD
This work focuses on identifying rainfall variability at the rainfall stations over the Opak River basin. The Opak River has a catchment area of about 2500 km 2 covering the southern part of Merapi Mountain to Opak River mouth at the southern coast of Java ( Figure 3). In general, the climate in this region is characterized by two major seasons. First, the rainy season, which commonly lasts from October to April, with a typical rainfall intensity around 10 -20 mm/day or 40 -60 mm/day during its peak. Second, the dry season which lasts between May to September with mostly zero rainfall. The area covers 20 rainfall stations, managed by the River Basin Management Organization of Serayu -Opak (BBWSO). The available data period varies from 1979 to 2019 for each station, as displayed in Figure 4. Most of the rainfall data were available from the early 2000s, with several stations started from 1985. However, in 2010 many records were missing. One might associate this with the eruption of Mount Merapi, which occurred in October 2010, since the missing data were located near the mountain (Station 1 to Station 9 in Figure 3). Nevertheless, there is no clear evidence of this causation.
Prior to the trend analysis, a homogeneity test was carried out to each station based on Buishand's test, the rescaled adjusted partial sums (RAPS) (Buishand 1982). Then, nonparametric tests against the trend of several rainfall variables were carried out using the Mann -Kendall test (Mann 1945;Kendall 1948).

Homogeneity test
Homogeneity testing for a time series of rainfall datasets is crucial to detect inconsistency data due to non-climatic reasons, such as instrumentation error, human error, and abrupt change in the surrounding areas of the instrument (Patakamuri, Muthiah & Sridhar 2020). This study applied Buishand's Range test (Buishand 1982) to examine the consistency of rainfall data at each station in the study area. This method tests the cumulative deviation of each data point from their mean. It was selected to treat each dataset from each station separately. Thus, the inequality of data length ( Figure 4) for each station will not have any significant impact.
Conceptually, homogeneity of a long-term data series is reached when each value ( ) in a dataset does not deviate from the mean ( ̅ ). In this case, the cumulative deviation will be zero ( 0 * = 0) or fluctuate around zero. The data deviation could be calculated by using Equation (1). * = ∑ ( − ̅ ) =1 (1) with k = 1, 2, 3, … N. Then, the inhomogeneity was examined using a 'rescaled adjusted partial sums', * * , which is formulated as in Equation (2) where is the standard deviation.
The sensitivity and significance of the deviation were identified by and values as in Equation (3) and (4).
where smaller values of and shows a better consistency. To determine whether a data series is consistent, the /√ and /√ values should be compared to the critical values (Buishand 1982). Lower values represent a better homogeneity, while values exceeding the critical value indicate inhomogeneity.
We present the homogeneity test results in Figure  3 by classifying the stations into five categories: a) Very useful = both and values pass the critical test on 95% confidence, b) Useful = one of the values falls between the critical test on 95% and 99% confidence, c) Doubtful = both values fall between 95% and 99% critical value, d) Suspect = one of the values still less than the critical value on 99% confidence, e) Ignore = both values are larger than the critical values on 99% confidence, or the data length is less than three years.
A similar categorization is seen in (Patakamuri et al. 2020) and (Wijngaard, Klein Tank & Können 2003) where they classified the homogeneity test results into three categories based on absolute and relative homogeneity tests.

Nonparametric Mann -Kendall test
Trends of time series data were identified by conducting the nonparametric Mann-Kendall test (Mann 1945;Kendall 1948). The Mann -Kendall (MK) test checks the null hypothesis ( 0 ) of no trend and alternate hypothesis ( ) of increasing or decreasing monotonic trend. The Mann-Kendal test statistic ( ) is calculated using Equation (5).
where is the number of data points, and are the data values in time series and ( > ), respectively, while and denote the data time.
The statistic of is normally distributed with the mean and the variance as in Equations (6) and (7).
where is the number of data points, is the number of tied groups, and is the number of the th tied group. The -test is then computed by following Equation (8).
The detected trend ( ) is accepted if the computed value of | | is greater than the /2 obtained from the standard normal cumulative distribution tables. An increasing trend will be indicated by the positive sign of and vice versa. For example, at 0.95 confidence level, the value is 1.645. Therefore, a data series with a calculated greater than 1.645 will indicate an increasing trend with a 0.95 confidence level.

Sen's slope estimator
Sen's slope estimator was used to predict the magnitude of the trend (Sen 1968). The slope, , is then calculated to each data series of = 1 , 2 , … , , using Equation (9).
where is the number of data points, and are the data values in time series and ( > ), respectively. The Sen's slope estimator, , is then calculated by taking the median of .

RESULT AND DISCUSSION
Results of both homogeneity and trend tests are summarized in Table 1. The homogeneity test was applied to all 20 stations in the Opak River basin. The "data length" column in the table shows the total number of years of available data with the detailed period presented in Figure 4. The data length at each station varies from 10 to 40 years, giving enough confidence to investigate the rainfall trend. The average annual rainfall at each station ranges from 1420 mm to 2700 mm with the highest record at around 4000 mm as shown in Figure 5.
Generally, most datasets show homogeneity, except for Gemawang (S07), Terong (S13), and Sanden (S15), which could not pass the consistency test as described in Table 1 and Figure 3. Both station Terong (S13) and Sanden (S15) do not have enough data for the test. The test showed inhomogeneity in the rainfall dataset of station Gemawang (S07) with a value of /√ is 1.407 and /√ of 1.589, which exceed the critical values with a 99% confidence level. This inconsistency was due to a relatively large data variation compare to the total data length. In 2009, the total annual rainfall was about 894.8 mm, which is relatively low compare to the average value (1956.53 mm). Excluding the records before 2009, the rainfall dataset in Gemawang station is consistent with the /√ and /√ values are 0.593 and 0.888 respectively. However, a nine-year dataset was considered insufficient for a climate study. Thus, we excluded Gemawang station (S07) in the trend analysis.
The rainfall trend was first identified through a preliminary graphical inspection and a simple linear trend test. Sen's slope estimator was also calculated to estimate the trend magnitude. The nonparametric Mann -Kendall trend test was then applied to 17 stations that have previously passed the homogeneity test. The test was employed for six data series at each station: total annual rainfall, annual maxima of daily rainfall, monthly rainfall, and three values of heavy rainfall ( > 50 mm, > 100 mm, and > 150 mm). Result of the trend tests indicates a tendency of increasing rainfall, especially for heavy rainfall occurrences. The detail of the analysis will be discussed in the following sections.

Trends of rainfall intensity
Rainfall trends both in the area-averaged and at each station were first analyzed. A decreasing trend is shown in the area-averaged value of annual rainfall as illustrated in Figure 6. Both the linear regression and Sen's slope estimator show a significant decrease. However, the Mann -Kendall test suggests a no trend as the negative | | value is less than the statistic (Table 1). In contrast, a slightly increasing trend is confirmed in the area-averaged value of the annual maximum of daily rainfall (Figure 7). A similar pattern is shown at each station, where the annual maximum of daily rainfall suggests a more pronounced increase while the annual rainfall shows a slight change in both positive and negative directions. This might indicate uneven time distribution of rainfalls, a higher intensity of rainfall might occur in a short duration, which could cause a more severe flood than a lower intensity with a longer duration.  Following the linear trend and Sen's slope estimator, the Mann -Kendall test was also performed on the 17 stations that have consistent rainfall datasets. The test was done to the following data series: a) Annual rainfall, b) Annual maximum of daily rainfall, c) Monthly rainfall, and d) Heavy rainfall frequency.
The test results are displayed in Figure 8 to Figure  10 consecutively. In Figure 8, we spatially represent rainfall changes at each station over the study area for both the annual precipitation ( Figure 8a) and its daily maximum (Figure 8b). We then present the seasonal change through the trend analysis of the monthly rainfall (Figure 9), which is followed by the increasing number of heavy rainfalls over the study area ( Figure 10).
A monotonic increasing trend of the annual rainfall, with 95% confidence level, was identified only at a few stations and no significant decreasing trend was reported. An increasing annual rainfall trend was observed at the Wanagama station at 95% confidence level. In addition to that, at the Karangploso and Bedugan stations, increasing trends are observed at 92.5% confidence level. Furthermore, the calculated values are dominated by positive values, while negative values resulted only at five stations as presented in Table 1. Therefore, it might be safe to assume an increasing rainfall trend in the annual rainfall at most of the stations over Opak River basin.
A similar result was obtained from the trend test of the annual maximum of daily rainfall. Most of the calculated values are positive with four stations showing an increasing trend at 95% confidence level: Kemput, Pundong, Kedungkeris, and Wanagama. Areas of the increasing and possibly decreasing maximum rainfall can be seen from Figure 8b. Increasing trends appear around Wanagama station for both the total annual rainfall and annual maxima. The negative trends were indicated at the Beji/Ngawen and Gedangan stations, although these stations are in a similar environment with the Wanagama station. This result might show the high spatial variability of the rainfall trend.
Analysis of the monthly rainfall trend at each station was also done to see the seasonal rainfall change. The bar chart in Figure 9 shows the number of stations having a positive and negative trend based on Mann -Kendall test. More stations seem to have an increasing trend during the rainy season and vice versa. A significant increasing trend was detected around march and April, while the significant decreasing trend was found around August. These results are likely aligned with (Zhang & Fueglistaler 2019) which suggested a more unevenness of rainfall pattern as well as (Held & Soden 2006) which showed evidence of the "wet gets wetter and dry gets drier" phenomenon.

Trend analysis of extreme rainfall
We further analyzed the trend by applying the Mann -Kendall test to the frequency of heavy and extreme rainfall. In this work, we used the rainfall classification from the Agency for Meteorology Climatology and Geophysics (BMKG) as follows: a) Light rain : 5 -20 mm/day, b) Moderate rain : 20 -50 mm/day, c) Heavy rain : 50 -100 mm/day, d) Very heavy rain : 100 -150 mm/day, e) Extreme rain : >150 mm/day.
We applied the trend analysis to the frequency of heavy rain ( > 50 mm/day), very heavy rain ( > 100 mm/day), and extreme rain ( > 150 mm/day). The analysis result is displayed in Figure 10. Horizontal lines in the figure refer to the statistic value at 95% confidence level. The values outside the horizontal lines indicate a monotonic trend in both increasing and decreasing directions. The graph suggests a dominant increasing trend of heavy rainfall at most of the stations. This result supports the assumption of more frequent occurrences of heavy rainfall in the future.

CONCLUSIONS
We have collected daily rainfall data from 20 rainfall stations over the Opak River basin, Yogyakarta, Indonesia, for a period of 12 to 38 years, ranging from 1979 to 2019. We aggregated the data into several data series: annual, annual maximum daily, and monthly rainfall. The frequency of heavy rainfall was also extracted and analyzed. All data series were then subjected to the trend tests based on the linear regression analysis, Sen's slope estimator, and Mann -Kendall tests. Generally, the trend tests seem to suggest increasing trends at most of the stations over the study area, aligned with the previous work (As-syakur et al. 2013). However, a minor decreasing trend is also observed. Both the linear regression and Sen's slope estimator agreed that rainfall has been increasing at most of the stations. The intensifying rainfall at several stations appears to be rather significant as suggested by the Mann -Kendall test on the annual maximum of daily rainfall. Moreover, a positive trend of the heavy rainfall frequency was observed at ten stations in a broad area of the Opak River basin. On the other hand, this study contrasts with (Aldrian 2007), which suggested a decreasing precipitation trend over Indonesia.
Their study used 63 rainfall data samples over Indonesia, with one rainfall station in Yogyakarta from 1951 to 1997. Their result suggested that the rainfall during that period was decreasing. This difference is understandable due to a different data period used in the analysis. Furthermore, it might be evidence of the high temporal variability of rainfall.
Therefore, it might be safe to expect higher and more frequent rainfall in the future, especially during the rainy season. Therefore, increasing rainfall design could also be anticipated in the water management engineering, such as drainage or flood protection system. On the contrary, less rainfall might be expected during the dry season as detected by the trend test on the monthly rainfall data series. This could pose a challenge for water resource management, such as water supply systems or reservoir management. This issue is significant not only for the present study area but also for a broad tropical region as the same trends might be shared under the similar climate system. Further studies will be needed to clarify the river discharge trends as it was also affected by land-use changes under the rapid development of the catchment area.

DISCLAIMER
The authors declare no conflict of interest.

ACKNOWLEDGMENTS
The authors would like to thank the River Basin Management Organization of Serayu-Opak (BBWSO) for providing the data necessary for this research. We also appreciate the Hydraulic Laboratory of Civil Engineering and Environmental Department of Universitas Gadjah Mada for maintaining the automatic rainfall recorder and providing the real time rainfall data.