Journal topic
Hydrol. Earth Syst. Sci., 22, 4213–4228, 2018
https://doi.org/10.5194/hess-22-4213-2018

Special issue: The changing water cycle of the Indo-Gangetic Plain

Hydrol. Earth Syst. Sci., 22, 4213–4228, 2018
https://doi.org/10.5194/hess-22-4213-2018

Research article 09 Aug 2018

Research article | 09 Aug 2018

Modeling the changes in water balance components of the highly irrigated western part of Bangladesh

Modeling the changes in water balance components of the highly irrigated western part of Bangladesh
A. T. M. Sakiur Rahman1, M. Shakil Ahmed2, Hasnat Mohammad Adnan3, Mohammad Kamruzzaman4, M. Abdul Khalek2, Quamrul Hasan Mazumder3, and Chowdhury Sarwar Jahan3 A. T. M. Sakiur Rahman et al.
• 1Hydrology Lab, Department of Earth and Environmental Sciences, Graduate School of Science and Technology, Kumamoto University, 2-40-1 Kurokami, Kumamoto, Japan
• 2Department of Statistics, University of Rajshahi, Rajshahi 6205, Bangladesh
• 3Department of Geology and Mining, University of Rajshahi, Rajshahi 6205, Bangladesh

Correspondence: A. T. M. Sakiur Rahman (shakigeo@gmail.com)

Abstract

The objectives of the present study were to explore the changes in the water balance components (WBCs) by co-utilizing the discrete wavelet transform (DWT) and different forms of the Mann–Kendall (MK) test and develop a wavelet denoise autoregressive integrated moving average (WD-ARIMA) model for forecasting the WBCs. The results revealed that most of the potential evapotranspiration (PET) trends (approximately 73 %) had a decreasing tendency from 1981–1982 to 2012–2013 in the western part of Bangladesh. However, most of the trends (approximately 82 %) were not statistically significant at a 5 % significance level. The actual evapotranspiration (AET), annual deficit, and annual surplus also exhibited a similar tendency. The rainfall and temperature exhibited increasing trends. However, the WBCs exhibited an inverse trend, which suggested that the PET changes associated with temperature changes could not explain the change in the WBCs. Moreover, the 8-year (D3) and 16-year (D4) periodic components were generally responsible for the trends found in the original WBC data for the study area. The actual data was affected by noise, which resulted in the ARIMA model exhibiting an unsatisfactory performance. Therefore, wavelet denoising of the WBC time series was conducted to improve the performance of the ARIMA model. The quality of the denoising time series data was ensured using relevant statistical analysis. The performance of the WD-ARIMA model was assessed using the Nash–Sutcliffe efficiency (NSE) coefficient and coefficient of determination (R2). The WD-ARIMA model exhibited very good performance, which clearly demonstrated the advantages of denoising the time series data for forecasting the WBCs. The validation results of the model revealed that the forecasted values were very close to actual values, with an acceptable mean percentage error. The residuals also followed a normal distribution. The performance and validation results indicated that models can be used for the short-term forecasting of WBCs. Further studies on different combinations of wavelet analysis are required to develop a superior model for the hydrological forecasting in the context of climate change. The findings of this study can be used to improve water resource management in the highly irrigated western part of Bangladesh.

1 Introduction

The water balance model is considerably important for water resource management, irrigation scheduling, and crop pattern designing (Kang et al., 2003; Valipour, 2012). The model can also be used for the reconstruction of catchment hydrology, climate change impact assessment, and streamflow forecasting (Alley, 1985; Arnall, 1992; Xu and Halldin, 1996; Molden and Sakthivadivel, 1999; Boughton, 2004; Anderson et al., 2006; Healy et al., 2007; Moriarty et al., 2007; Karimi et al., 2013). Therefore, accurately forecasting the water balance components (WBCs) and detecting the changes in them is important for achieving sustainable water resource management. However, hydrometeorological time series are contaminated by noises from hydrophysical processes. This affects the accuracy of the analysis, simulation, and forecasting (Sang et al., 2013; Wang et al., 2014). Hence, denoising the time series is essential for improving the accuracy of the obtained results. In this study, the wavelet denoising technique was coupled with the ARIMA (autoregressive integrated moving average) model for forecasting the WBCs after detecting the changes in them by using different forms of the Mann–Kendall (MK) test. Moreover, the time period responsible for the trends in the WBC time series was identified using discrete wavelet transform (DWT) time series data.

2 Study area, data, and methods

2.1 Study area

The climate of Bangladesh is humid, warm, and tropical. The western part of Bangladesh covers approximately 41 % or 60 165 km2 of the country. The geographic coordinates of the study area extend between a latitude of 2136–2638 N and longitude of 8819–9101 E. The annual rainfall and average temperature in the study area vary from 1492 to 2766 mm, with an average of 1925 mm, and 24.18 to 26.17 C, with an average of 25.44 C, respectively (Kamruzzaman et al., 2016a). Bangladesh is the fourth-largest producer of rice in the world (Scott and Sharma, 2009), and the livelihood of a majority of the people (approximately 75 %) (Shahid and Behrawan, 2008; Kamruzzaman et al., 2016b) depends on agricultural practices. The crop calendar of Bangladesh is related to the climatic seasons. Rice is grown during three seasons (Aus, Aman, and Boro) in Bangladesh. Almost 73.94 % of the cultivable area in the country is used to cultivate Boro rice (Banglapedia, 2003). The Aus and Aman rice varieties are mainly rain-fed crops. However, Boro rice is almost completely groundwater-fed (Ravenscroft et al., 2005) and requires approximately 1 m of water per square meter in Bangladesh (Harvey et al., 2006; Michael and Voss, 2009).

2.2 Data

The national climate database of Bangladesh prepared by the Bangladesh Agricultural Research Council (BARC) was used for this study. The database is available for research and can be obtained from the BARC website (http://climate.barcapps.gov.bd/, last access: 27 July 2018). The database has been prepared from the data recorded by the Bangladesh Meteorological Division and contains long-term monthly climate data, such as rainfall, minimum, maximum, and average temperatures, humidity, sunshine hours, wind speed, and cloud cover. The locations of the meteorological stations in the study area are displayed in Fig. 1. The data are rearranged according to the hydrological year for the period from 1981–1982 to 2012–2013. The hydrological year in Bangladesh begins in April and ends in March.

Figure 1Study area in the western part of Bangladesh with locations of meteorological stations.

2.3 Methods

In this study, the WBCs were calculated and their trends were identified using the MK or Modified MK (MMK) test for evaluating the long-term water balance of the highly irrigated western part of Bangladesh. The DWT data of the WBC time series were analyzed for identifying the time period responsible for the trend in the data. The WBCs were forecasted using the ARIMA model, whose performance was statistically evaluated. If the performance of the model was unsatisfactory for forecasting the WBCs, denoising of the original time series was conducted using DWT techniques to improve the performance of the model. The descriptions of the methods are presented in the following sections.

2.3.1 Calculation of the potential evapotranspiration and water balance components

The potential evapotranspiration (PET) is a key parameter to estimate the WBCs. In this study, the potential evapotranspiration was calculated using the Penman–Monteith equation (Allen et al., 1998). The soil–water balance concept proposed by Thornthwaite and Mather (1955) is one of the most widely used methods for estimating the WBCs. This method is suitable for assessing the effectiveness of agricultural water resource management practices and regional water balance studies because it allows the actual evapotranspiration (AET), water deficit, and water surplus to be estimated (Chapman and Brown, 1966; Bakundukize et al., 2011; Karim et al., 2012; Viaroli et al., 2017). The actual evapotranspiration (AET) is the amount of water removed from the surface due to evaporation and transpiration. The amount by which the PET exceeds the AET is termed as the deficit. The surplus is the excess rainfall received after the soil has reached its water-holding capacity (de Jong and Bootsma, 1997). Calculating the field capacity of the soil is essential for estimating the WBCs. The field capacity of the soil in the study area was calculated using the soil texture map of Bangladesh prepared by the Soil Resource Development Institute, Bangladesh (SRDI, 1998), where the description of the soils was presented by Huq and Shoaib (2013). The values suggested by Thornthwaite and Mather (1957) for the water-holding capacity of the soil and rooting depth of the plants were used for estimating the WBCs in this study. The first step of the calculation involves subtracting 5 % rainfall from the monthly rainfall data because this amount of water is lost due to direct runoff (Wolock and McCabe, 1999; Karim et al., 2012; Kanoua and Merkel, 2015). The remaining rainfall amount is included in the calculation. The WBCs, such as the AET, surplus, and deficit, were estimated using the formulas presented in Table 1. The details of the WBC calculations are available in the Supplement.

Table 1Calculations of water balance components (Thornthwaite and Mather, 1957).

P is the rainfall (mm), R0 is the direct runoff (mm), PET is the potential evapotranspiration (mm), AET is the actual evapotranspiration (mm), and ΔSB is the changes in soil moisture storage (mm).

2.3.2 Trend test

In this study, the trends in the WBCs were detected using the nonparametric MK test (Mann, 1945; Kendall, 1975) because it exhibits a better performance than the parametric test (Nalley et al., 2012) for identifying trends in hydrological variables, such as rainfall (Shahid, 2010), temperature (Kamruzzaman et al., 2016a), PET (Kumar et al., 2016), soil moisture (Tabari and Talaee, 2013), runoff (Pathak et al., 2016), groundwater level (Rahman et al., 2016), and water quality (Lutz et al., 2016). The MK test cannot be used to accurately calculate the test statistic (Z) if there exists a significant serial correlation at lag 1 in the time series data (Yue et al., 2002) because the variance is underestimated (Hamed and Rao, 1998). The autocorrelation at lag 1 was checked before analyzing the time series data. If there existed a significant lag-1 autocorrelation at the 5 % level, the MMK test (Hamed and Rao, 1998) was applied instead of the MK test. The estimated Z statistic from the MK or MMK test was evaluated for the direction of the trend (a positive Z statistic indicated an increasing trend and vice versa). Moreover, the Z statistic indicated the level of significance of the obtained trend. If the calculated Z statistic is equal to or higher than the tabulated value of the Z statistic (+1.96), it indicates a significant positive trend at the 95 % confidence level. If the calculated Z statistic is equal to or less than 1.96, it indicates a significant decreasing trend. Moreover, the sequential values of the u(t) statistic derived from the sequential MK (SMK) test (Sneyers, 1990) are used for detecting the change point. The u(t) statistic is similar to the Z statistic (Partal and Küçük, 2006). The magnitude of the change was calculated using Sen's slope estimator (Sen, 1968). Numerous studies have already been conducted (notably Nalley et al., 2012) using the methods described in this section. Further details regarding these methods can be obtained from Mann (1945), Sen (1968), Kendall (1971), Hamed and Rao (1998), Sneyers (1990), and Yue et al. (2002).

2.3.3 Wavelet transform (WT) and periodicity

Wavelet analysis has been used in different parts of the world to identify the periodicity in hydroclimatic time series data (Smith et al., 1998; Azad et al., 2015; Nalley et al., 2012; Araghi et al., 2015; Pathak et al., 2016). WT, a multiresolution analytical approach, can be applied to analyze time series data because it offers flexible window functions that can be changed over time (Nievergelt, 2001; Percival and Walden, 2000). WT can be applied to detect the periodicity in hydroclimatic time series data (Smith et al., 1998; Pišoft et al., 2004; Sang, 2012; Torrence and Compo, 1998; Araghi et al., 2015; Pathak et al., 2016) and exhibits better a performance than traditional approaches (Sang, 2013). There exist two main types of WT, namely continuous WT (CWT) and DWT. Applying the CWT is complex because it produces numerous coefficients (Torrence and Compo, 1998; Araghi et al., 2015), whereas DWT is simple and useful for hydroclimatic analysis (Partal and Küçük, 2006; Nalley et al., 2012). The wavelet coefficients of the DWT with a dyadic format can be calculated as follows (Mallat, 1989):

$\begin{array}{}\text{(1)}& {\mathit{\psi }}_{m,n}\left(\frac{t-\mathit{\tau }}{s}\right)={s}_{o}^{-m/\mathrm{2}}\mathit{\psi }\left(\frac{t-n{\mathit{\tau }}_{o}{s}_{o}^{m}}{{s}_{o}^{m}}\right),\end{array}$

where ψ is the mother wavelet, m is the wavelet dilation, and n is the wavelet translation. The specified fixed dilation step (so) is larger than 1, and τo is the location parameter. For practical application, the values of so and τo are considered as 2 and 1, respectively (Partal and Küçük, 2006; Pathak, 2016). After substituting these values in Eq. (1), the DWT for a time series xi becomes the following:

$\begin{array}{}\text{(2)}& {W}_{m,n}={\mathrm{2}}^{-m/\mathrm{2}}\sum _{i=\mathrm{0}}^{N-\mathrm{1}}{x}_{i}\mathit{\psi }\left({\mathrm{2}}^{-m}\phantom{\rule{0.125em}{0ex}}i-n\right),\end{array}$

where W indicates the wavelet coefficient at a scale s=2m and location τ=2mn

In the DWT, details (D) and approximations (A) of the time series can emerge from the original time series after passing through low-pass and high-pass filters, respectively. When approximations are the high-scale and low-frequency components, details are the low-scale and high-frequency components. Successive iterations are performed to decompose the time series into its several low-resolution components (Mallat, 1989; Misiti et al., 1997). In this study, four levels (D1–D4) of decomposition were performed following the dyadic scales. The decompositions are referred to as D1, D2, D3, and D4, which correspond to a 2-, 4-, 8-, and 16-year periodicity, respectively. The Daubechies wavelet was used because of its superior performance in hydrometeorological studies (Nalley et al., 2012, 2013; Ramana et al., 2013; Araghi et al., 2015). To confirm the periodicity present in the time series, the correlation coefficient (Co) between u(t) of the original data, u(t) of the decomposition (D) time series data, and different models (D1 + A … D4 + D3 + A) of the time series data were calculated and the obtained results were compared (Partal and Küçük, 2006; Partal, 2009).

2.3.4 ARIMA models

ARIMA models (Box and Jenkins, 1976) are used in hydrological science to identify the complex patterns in data and project future scenarios (Adamowski and Chan, 2011; Valipour et al., 2013; Nury et al., 2017; Khalek and Ali, 2016). ARIMA models include (1) an autoregressive process (AR) represented by order p, (2) nonseasonal differences for nonstationary data termed as order d, and (3) a moving average (MA) process represented by order q. An ARIMA model of order $\left(p,d,q\right)$ can be written as follows:

$\begin{array}{}\text{(3)}& {\mathrm{\varnothing }}_{p}\left(L\right){\left(\mathrm{1}-L\right)}^{d}{Y}_{t}={\mathit{\theta }}_{\mathrm{0}}+{\mathit{\theta }}_{q}\left(L\right){U}_{t},\end{array}$

where θ0 is the intercept with a mean of 0, Ut is the white process with constant variance, p(L) represents the AR term $\left(\mathrm{1}-{\mathrm{\varnothing }}_{\mathrm{1}}L-\phantom{\rule{0.25em}{0ex}}\mathrm{\dots }\phantom{\rule{0.25em}{0ex}}-{\mathrm{\varnothing }}_{p}{L}^{p}\right)$, and θq(L) represents the MA term $\left(\mathrm{1}-{\mathit{\theta }}_{\mathrm{1}}L-\phantom{\rule{0.25em}{0ex}}\mathrm{\dots }\phantom{\rule{0.25em}{0ex}}-{\mathit{\theta }}_{p}{L}^{p}\right)$.

2.3.5 Wavelet denoising

Wavelet denoising based on the thresholds introduced by Donoho et al. (1995) has been applied to hydrometeorological analysis (Wang et al., 2005, 2014; Chou, 2011). In this study, the following three analysis steps were performed for denoising the time series data.

1. Decomposing the time series data x(t) into M resolution levels for obtaining the detail coefficients (Wj,k) and approximation coefficients using the DWT.

2. The detail coefficients obtained from the DWT (1 to M levels) were treated using threshold (Tj) selection. A soft or hard threshold can be used to deal with detail coefficients and obtain the decomposed coefficient. In this study, a soft threshold was selected because it performed better than a hard threshold (Wang et al., 2014; Chou, 2011).

3. Detail coefficients from levels 1 to M and approximate coefficients at level M were reconstructed to obtain denoising time series data.

Selecting the threshold value is essential for denoising the data. In this study, the universal threshold (UT) method (Donoho and Johnstone, 1994) was used for estimating the threshold value because it exhibited satisfactory performance in analyzing hydrometeorological data (Wang et al., 2005; Chou, 2011).

2.3.6 Assessment of model performance

There exist several indicators to assess the performance of the models. The Nash–Sutcliffe efficiency (NSE) (Nash and Sutcliffe, 1970) coefficient, a normalized goodness-of-fit statistic, is the most powerful and popular method for measuring the performance of hydrological models (McCuen et al., 2006; Moussa, 2010; Ritter and Muñoz-Carpena, 2013). The NSE coefficient was used in this study to evaluate and compare the ARIMA and WD-ARIMA models. The NSE is calculated as follows (Nash and Sutcliffe, 1970):

$\begin{array}{}\text{(4)}& \mathrm{NSE}=\mathrm{1}-\frac{{\sum }_{i=\mathrm{1}}^{N}\left({O}_{i}-{P}_{i}{\right)}^{\mathrm{2}}}{{\sum }_{i=\mathrm{1}}^{N}\left({O}_{i}-\stackrel{\mathrm{‾}}{O}{\right)}^{\mathrm{2}}}=\mathrm{1}-{\left(\frac{\mathrm{RMSE}}{\mathrm{SD}}\right)}^{\mathrm{2}},\end{array}$

where N, Oi, Pi, $\stackrel{\mathrm{‾}}{O}$, and SD are the sample size, number of observations, model estimates, mean, and standard deviation of the observed values, respectively. The performance of a model can be evaluated according to its NSE value as very good (NSE≥0.9), good (NSE = 0.8–0.9), acceptable (NSE≥0.65), and unsatisfactory (NSE < 0.65) (Ritter and Muñoz-Carpena, 2013). ERMS is the root-mean-square error and can be calculated as follows:

$\begin{array}{}\text{(5)}& {E}_{\mathrm{RMS}}=\sqrt{\frac{{\sum }_{i=\mathrm{1}}^{N}\left({O}_{i}-{P}_{i}{\right)}^{\mathrm{2}}}{N}}.\end{array}$

The coefficient of determination (R2) is another goodness-of-fit test to measure the performance of models. The perfect fit of the model draws a line between the actual values and fitted values, where R2 is 1. If yi is the observation data, ${\stackrel{\mathrm{^}}{y}}_{i}$ represents the model-forecasted values of yi and N is the number of data points used. R2 is given as follows (Sreekanth et al., 2009):

$\begin{array}{}\text{(6)}& {R}^{\mathrm{2}}=\mathrm{1}-\frac{{\sum }_{i=\mathrm{1}}^{N}{\left({y}_{i}-{\stackrel{\mathrm{^}}{y}}_{i}\right)}^{\mathrm{2}}}{{\sum }_{i=\mathrm{1}}^{N}{\left({y}_{i}\right)}^{\mathrm{2}}-\frac{\left({\sum }_{i=\mathrm{1}}^{N}{y}_{i}{\right)}^{\mathrm{2}}}{N}}.\end{array}$

Moreover, the mean percentage error (EMP) and mean error (EM) were also calculated to evaluate the validation of the model for forecasting. EMP indicates the percentage of bias (large or small) between the forecasted and actual data (Khalek and Ali, 2016). EMP and EM can be calculated as follows:

$\begin{array}{}\text{(7)}& & {E}_{\mathrm{MP}}=\left(\frac{\mathrm{1}}{n}\phantom{\rule{0.125em}{0ex}}\sum _{t=\mathrm{1}}^{n}\frac{{Y}_{t}\left(\mathrm{actual}\right)-{Y}_{t}\left(\mathrm{forecasted}\right)}{{Y}_{t}\left(\mathrm{actual}\right)}\right)×\mathrm{100}\phantom{\rule{0.125em}{0ex}}\mathit{%},\text{(8)}& & {E}_{\mathrm{M}}=\frac{\mathrm{1}}{n}\sum _{t=\mathrm{1}}^{n}\left[{Y}_{t}\left(\mathrm{actual}\right)-{Y}_{t}\left(\mathrm{forecasted}\right){\right]}^{\mathrm{2}}.\end{array}$
3 Results

3.1 Exploratory statistics of the water balance components

The mean annual PET in the study area between 1981–1982 and 2012–2013 varied from 1228 to 1460 mm (Fig. 2a), with an average of 1338 mm. High PET values were observed in the central part of the area, where the annual rainfall was low, but the temperature was high (Kamruzzaman et al., 2016a). The standard deviations of the PET varied from 205 (Jessore station) to 41 mm (Bhola station). The AET (Fig. 2b) (average = 925 mm) was almost 31 % less than the PET because during the dry months (December–May), the soil moisture condition reached a critical stage. The annual surplus of water varied from 515 to 1277 mm (Fig. 2c), with an average of 838 mm. According to Wolock and McCabe (1999), 50 % of the surplus water can be considered as runoff for the major parts of the world. A high amount of surplus water was found in the northern part of the study area and along the coastal area. The annual deficit of water, which mainly occurred during the dry season (December–May), varied from 329 to 556 mm, with an average of 416 mm (Fig. 2d). The highest annual deficit of water was observed in Rajshahi, which is located in the central–western part of the study area, where the depth of groundwater below the surface increases rapidly (Shamsudduha et al., 2009; Rahman et al., 2016).

Figure 2Distribution of mean annual (a) PET, (b) AET, (c) surplus, and (d) deficit of water in the study area during the hydrologic year 1981–1982 to 2012–2013.

3.2 Trend and periodicity of the water balance components

3.2.1 Potential evapotranspiration

The MK or MMK test based on lag-1 autocorrelation was applied to detect the trend in the PET. Table 2 represents the Z statistic of the MK or MMK test for the original PET time series data and the Z statistic of the decomposition (D1–D4), approximation (A), and model (D1 + A … D3 + D4 + A) time series. The estimated Z statistic of the original data ranged from 2.07 (Satkhira station) to 2.37 (Bhola station). The Satkhira and Bhola stations exhibited significant PET trends. The plots of the sequential u(t) statistic obtained from the SMK test for these two stations are displayed in Fig. 3, where the dashed lines correspond to a 5 % significance level (±1.96). The decreasing PET trend for the Satkhira station began in 1985–1886, and a significant decreasing trend occurred in 1993–1994. The trend reversed after 2007–2008. However, the significant increasing PET trend of the Bhola station began very recently (2010–2011) after some fluctuation.

Table 2Z statistic of MK or MMK of original time series, approximation, and different models PET of DWT (the dominant components are shown in bold and the asterisks denote significance at a 5 % level).

MSE, total mean square error; Co, correlation between original data and DWT models.

Figure 3Sequential values of the u(t) statistics of (a) Satkhira station and (b) Bhola station.

Most of the trends (73 %) observed in the PET time series data of the study area were negative and statistically insignificant at the 95 % confidence level or 5 % significance level. Moreover, the Z statistic of the approximation (A) time series obtained using the DWT indicated decreasing PET trends for all the stations. The calculated Z statistic of the approximation (A) time series was approximately 1.8 after rounding the figures for all the stations. The approximation time series data of all the stations exhibited a similar pattern (Fig. S1 of the Supplement) over time. The magnitude of PET changes ranged from 10.89 mm yr−1 for the Satkhira station to 1.67 mm yr−1 for the Bhola station (Fig. 4a). The MK or MMK test was also applied to the decomposition time series and model time series generated from the combination of the approximation and decomposition time series data. Table 2 represents the results for four stations arranged in alphabetical order, and the complete results can be found in Table S1 of the Supplement. To determine the dominant periodicity affecting the PET trends, a two-step analysis was performed. First, the value closest to the Z statistic of the original time series data was obtained from the Z-statistic values of different model and decomposition time series data. Second, the correlation coefficients (Co) of pairs of data (such as the Co between the u(t) statistics obtained from the SMK test for the original and decomposition time series data) were estimated, and the highest Co was determined from the estimated Co values for different pairs (Table 2). The Z statistic of the D4 time series data for the Barisal station was 0.76, which was the closest to the Z statistic (0.72) of the original time series data (Table 2). Moreover, the Z statistic of the model (D3 + D4 + A) time series data was 0.56, which is the second-nearest value to the Z statistic of the original time series and has the highest correlation coefficient (Co = 0.85). The D4 (16-year) component was the dominant periodic component in the trend of the original data. However, D3 also affected the trend of the data. The Z-statistic value (2.47) of the original time series for the Bhola station was the closest to that (2.36) of the model (D2 + D4 + A) time series data. However, the Z-statistic values of the D2, D4, D2 + A, and D4 + A time series were 0.61, 1.2, 0.48, and 0.9, respectively. These values were not close to the Z statistic of the original time series data. Hence, in this case, the Z statistic was unable to determine which periodic component (D2D4) was the basic periodic component for the significant trend in the original data. To determine the dominant periodic component, the values of Co were analyzed. The correlation coefficient (Co) between the u(t) statistic of the SMK test for the original and D4 time series data was higher than the correlation coefficient between the u(t) statistic of the SMK test for the original and D2 time series data (Table 2). Moreover, the values of the Z statistic for time series with the D4 components, such as the D4 and D4 + A model time series, were higher than those for time series with the D2 component (D2 and D2 + A) (Table 2). Therefore, D4 was the main periodic component responsible for the PET trend of the Bhola station. However, the Z-statistic values of D4 and D4 + A were not close to the Z statistic of the original data (Table 2). Moreover, there existed a statistically significant positive trend in the original PET data of the Bhola station, whereas the trends of the D4 and D4 + A model time series data were not statistically significant. When the D2 time series was added to the D4 + A model time series data, the Z statistic of the resultant (D2 + D4 + A) model time series data was very close to that of the original time series data. The trend of the D2 + D4 + A model time series was statistically significant, similar to the trend in the original time series data (Table 2). Hence, D2 affected the trend of the original time series data. Station-wise analysis indicated that almost half of the stations exhibited harmony between the Z-statistic values of the D3 + D4 + A model and original time series data. Individual analysis of the D3 and D4 time series data indicated that a higher relationship existed between the D4 and original time series data. Three stations (Dinajpur, Ishurdi, and Jessore) exhibited similar Z-statistic values for the original and D1 + D4 + A model time series data, with higher Co values of the u(t) statistic for the SMK test on the D4 time series data than that for the SMK test on the original data (except for the Ishurdi station). Moreover, two stations (Bhola and Satkhira) exhibited significant trends in the original data. The closest Z statistic was found between the original and D2 + D4 + A time series data for both of the stations. D4 (16-year periodicity) was the dominant periodic component according to the Co values for both these stations. Therefore, 16-year periodicity was the main periodic component responsible for the trends in the PET data over the study area. Moreover, D3 (8-year) periodicity also had an effect on the trends for some stations (Tables 2 and S1). D4 (16-year) periodicity dominates the annual rainfall trend for the Marmara region in Turkey (Partal and Küçük, 2006). Araghi et al. (2015) determined that 8–16-year (D3 to D4) periodicity is responsible for the trends in the annual temperature in Iran.

Table 3Comparison of performance of ARIMA model and WD-ARIMA model.

3.2.2 Actual evapotranspiration

All the stations except the Bogra station exhibited decreasing trends in the AET. The calculated Z statistic ranged from 2.90 for the Bogra station to 0.31 for the Ishurdi station. Similar to the PET trends, the AET trends were also insignificant at a 5 % significance level. However, the Ishurdi station exhibited a significant (at a 5 % significance level) decreasing trend. The magnitudes of the trends of the original AET data varied from 5 mm yr−1 for the Faridpur station to 0.75 mm yr−1 for the Bogra station. The distribution of the trend magnitude is displayed in Fig. 4b. The periodicity in the AET was marginally different from that in the PET (Table S2). For almost half of the stations (five), D2 (4-year) was the main periodic component. D4 (16-year) also affected the trend because the Z statistic of the D2 + D4 + A model time series was the nearest to that of the original series for the Khulna and Ishurdi stations. Moreover, D4 (16-year) was the main periodic component for the Rangpur and Rajshahi stations. D1 (2-year) was the dominant periodic component for the Barisal, Bhola, and Bogra stations. The AET value depends on climatic factors, such as the PET, rainfall, and soil moisture conditions. The variations in the periodicities of the AET and PET were mainly related to the soil moisture conditions of the area.

Figure 4Distribution of rate of changes of WBCs during the period of 1981–1982 to 2012–2013.

3.2.3 Surplus

Almost 82 % of the stations exhibited insignificant decreasing trends for the annual surplus of water. The magnitude of the trends of the original annual surplus data ranged from 11.63 to 6.71 mm yr−1 (Fig. 4c). The periodicity characteristics of the PET and surplus were similar (Table S3). D4 (16-year) was the main periodic component present in seven stations. In most cases, D2 was also present (D2 + D4 + A), except in Rajshahi. D3 (8-year) was mainly responsible for the surplus trend of three stations. Surplus mainly occurred during the rainy season (June–October) in the study area, when the soil pores were almost completely filled with water and the AET was equal to the PET. Surplus mainly depends on rainfall and hence provides insight regarding the periodicity in rainfall.

3.2.4 Deficit

Approximately 73 % of the stations exhibited increasing trends for the annual deficit of water. The increasing trends were significant for two stations at the 95 % confidence level (Table S4). However, the Satkhira station exhibited a significant decreasing trend (Z $=-$2.08) in the annual deficit of water. The magnitude of the trends of the original annual deficit data ranged from 8.1 to 7.7 mm yr−1 (Fig. 4b). Periodicity analysis revealed that D4 was mainly responsible for the trends in the annual deficit of water. The Z statistic of the (D2 + D4 + A) model time series data was close to the Z statistic of the original time series data (Table S4). D3 (8-year periodicity) was also responsible for the trends in the data of the two stations.

3.3 Model selection and forecasting ability

The ARIMA model was selected for forecasting the WBC time series. A four-step analysis was performed during time series modeling. (1) First, the stationarity of the data was checked using the Augmented Dickey–Fuller (ADF) test. (2) Then, the autocorrelation function (ACF) was used for selecting the order of the MA process (Figs. S2–S5). (3) The partial autocorrelation function (PACF) was then used for selecting the order of the AR process (Figs. S2–S5). (4) Finally, the appropriate model was selected based on several trials and model selection criteria, such as Akaike information criterion (AIC) and Bayesian information criterion (BIC). In addition to the manual model selection based on the ACF, PACF, AIC, and BIC, the auto ARIMA function of the “forecast” package (Hyndman et al., 2017) of R (R 3.4.0 language developed by R Core Team, 2016) was used during the trails for model selection to obtain information regarding the nature of the data for modeling. The model with the lowest AIC and BIC values and highest R2 value was selected. The QQ plot was prepared to examine the normality of the residuals. The performance of the ARIMA model (parameters are given in Table S5) was evaluated using the NSE coefficient and R2 values (Table 3). The estimated NSE coefficient of the ARIMA model for the PET time series varied from 0.6 for the Bhola station to 0.81 for the Jessore station (Table 3). The ARIMA model exhibited an unsatisfactory performance for almost all the stations. The average NSE coefficient of the 11 stations was 0.38, and the R2 values ranged from 0.1 to 0.81, with an average of 0.38. Moreover, the NSE coefficient of the Bhola station indicated that the ARIMA model was unsuitable for forecasting the PET. The ARIMA model was also applied to the AET, surplus, and deficit time series data. There existed no significant spikes in the ACF and PACF of the AET (Fig. S3). Moreover, the results obtained from the auto ARIMA functions exhibited similar results. Therefore, the ARIMA model was unsatisfactory for forecasting the variability in the AET. For WBCs such as surplus and deficit, the performance of the ARIMA model was similar to that of the AET, except for a few cases. Because hydrometeorological data are affected by noises from different hydrophysical processes (Wang et al., 2014), the results obtained using the ARIMA models were unsatisfactory. To improve model performance, noise must be removed from the data. In this study, DWT denoising was applied to the WBC data and the quality of the denoising time series data was examined before further processing. When selecting a method for denoising the time series using WT, the mean of the original and denoising time series data should be close and the standard deviation of the denoising time series should be less than that of the original time series (Wang et al., 2014). Figure 5a displays the means of the actual and wavelet denoising PET time series. No visible difference was observed between the mean of the original and DWT wavelet denoising time series data. Moreover, the standard deviation of the PET for the wavelet denoising time series was lower than that for the original time series (Fig. 5b). The AET, surplus, and deficit time series also exhibited similar results. Furthermore, the lag-1 autocorrelation of the wavelet denoising time series data must be higher than that of the original time series (Wang et al., 2014). Under this condition, the absolute lag-1 value of autocorrelation for the wavelet denoising time series was higher than that for the original series (Figs. S2b, S3b, S4b, and S5b). The performance of the WD-ARIMA model is represented in Table 3. After denoising the data, the performance of the ARIMA model was satisfactory for all the WBC time series data (Table 3). The average NSE coefficient of the WD-ARIMA model for the PET time series of the 11 stations located in the western part of Bangladesh was 0.76, with an average R2 value of 0.67. The R2 and NSE coefficient values indicated that the performance of the WD-ARIMA model was better than that of the classical ARIMA model for the modeling of PET (Table 3). Moreover, the average NSE value of the WD-ARIMA model for the AET time series of the 11 stations was 0.92, which indicated that the performance of the model was very good. The average R2 value was 0.89, which indicated that the model could explain almost 89 % of the variance in the data (Table 3). The WD-ARIMA model also exhibited a very good forecasting performance for the annual surplus and deficit (Table 3). The average NSE coefficient of the WD-ARIMA model for the annual surplus of the 11 stations was approximately 0.92, and the average R2 value was 0.9. The WD-ARIMA model exhibited a good performance in forecasting the annual deficit (average NSE = 0.88). The performance of the WD-ARIMA model was good or very good for forecasting the AET, annual surplus, and annual deficit. However, the performance was acceptable for forecasting the PET. This deviation may have arisen because the variability of the PET was higher than that of the other WBCs or the deviation may be related to the variability of climatic variables.

Figure 5Comparison between actual and wavelet denoise PET time series (a) mean and (b) standard deviation.

The WD-ARIMA models were validated to explore their forecasting ability. The mean percentage error (EMP) of the forecasted values for the 4-year period from 2008–2009 to 2012–2013 was calculated to determine the percentage bias of the forecasted data (Table 4). The average EMP of the WD-ARIMA model for the PET values of the 11 stations was 0.6 (ranging from 0.75 to 3.34), which indicated that the forecasted values were marginally lower than the actual values. The typical plots of the actual time series data versus the fitted model data, normal QQ plots of the residuals of the models, and actual and observed values for the WBCs (plots for all the stations are displayed in Figs. S6–S9) are illustrated in Fig. 6. The plot of the actual values versus the forecasted values (Fig. 6) indicates that the actual and forecasted values were very close for the hydrologic years 2009–2010 and 2010–2011. The normal QQ plots revealed that the residuals of the models were near normal. However, the differences in the values increased after these two hydrologic years for all the WBCs (Figs. S6–S9). The EMP values of WD-ARIMA models for the AET ranged from 0.7 to 0.2, with an average of 0.09, which indicated that the forecasted AET values were marginally lower than the actual AET values. The EMP values for the annual surplus (average $=-$0.75) and annual deficit (average $=-$0.12) were similar to that for the AET and PET. The average EMP values for all the WBCs were negative, which indicated that the forecasted values for the WBCs were marginally lower than the actual values for most of the stations.

Table 4Accuracy of WD-ARIMA models of WBCs for validation of the model's predictive ability for the period of 2009–2010 to 2012–2013.

Figure 6Plot of best WD-ARIMA model first panel represents actual versus fitted values for the period of 1981–1982 to 2012–2013, the second panel is normal QQ plot of residuals of the model, and the third panel shows actual, fitted, and forecasted values for 2009–2010 to 2012–2013. (a) PET of Rangpur station located in the north, (b) AET of Ishurdi station located in the central part, (c) deficit of Rajshahi station located in NW Bangladesh, and (d) surplus of Bhola station located in the south of the study area.

3.4 Discussion

This study indicated that a decreasing PET trend dominated the study area. However, positive trends in the rainfall and temperature dominated the western part of Bangladesh (Shahid and Khairulmaini, 2009; Kamruzzaman et al., 2016a). Moreover, a recent study found a negative trend in the evapotranspiration for four stations located in northwest Bangladesh (Acharjee et al., 2017). Although the annual rainfall and temperature of the Satkhira station exhibited positive trends (Kamruzzaman et al., 2016a), its PET exhibited a significant decreasing trend. Increasing temperature and decreasing PET trends were observed in the Yunnan Province of South China (Fan and Thomas, 2012). McVicar et al. (2012) also found decreasing PET trends in different parts of the world. Therefore, although the temperature is the primary factor driving changes in the PET (IPCC, 2007), temperature-based models cannot suitably explain the causes of PET changes. To obtain a detailed insight regarding the mechanisms underlying the PET changes, a detailed analysis must be conducted of all climatic variables, such as rainfall, temperature, sunshine hours, wind speed, and humidity, and climate-controlling phenomena, such as El Niño–Southern Oscillation.

The WD-ARIMA model was used in this study for forecasting the WBCs. The performance of the model indicated the benefit of denoising hydrological time series data, such as the PET, AET, surplus, and deficit. However, the NSE coefficient indicated that the performance of the model was acceptable for PET forecasting (NSE  0.65). The deviation between the forecasted values and actual values increased with increasing time steps. Therefore, the WD-ARIMA model was unsuitable for long-term forecasting. The WD-ARIMA model was developed by coupling the discrete wavelet denoising time series data and ARIMA model. The soft threshold method was selected for denoising the time series data, and the UT method was used for determining the threshold value. However, there exist other approaches, such as SURE (Stein, 1981) and MINMAX (Donoho and Johnstone, 1998), for determining the threshold value. Moreover, Wang et al. (2014) developed a hybrid method called the adaptive wavelet denoising approach using sample entropy (AWDA-SE) for denoising hydrometeorological time series data, such as rainfall and streamflow data. The study (Wang et al., 2014) indicated that the performance of the developed denoising method was better than that of conventional methods for denoising rainfall and streamflow data. The aforementioned approaches may be used to increase the performance of the ARIMA model for forecasting hydrological variables, such as the PET. Moreover, there exist several mother wavelet families, such as Daubechies, Haar, Coiflets, Morlet, and Mexican hat (Sang, 2013). In this study, only Daubechies 6 from the Daubechies wavelet family was applied as the mother wavelet for the DWT. The WD-ARIMA model exhibited very good performance for forecasting the AET, surplus, and deficit, whereas the classical ARIMA model exhibited poor performance or was unable to forecast the WBCs. Moreover, studies (Chou, 2011; Kisi, 2008; Partal, 2009; Santos and da Silva, 2014; Rahman and Hasan, 2014; Nury et al., 2016; Adamowski and Chan, 2011; Khalek and Ali, 2016) have indicated that the performance of wavelet-aided models is better than that of the classical ARIMA and ANN models for forecasting nonstationary hydrometeorological variables. Because traditional methods such as Wiener filtering, Kalman filtering, and Fourier transform are unsuitable for nonstationary hydrological time series data (Adamowski and Chan, 2011; Sang, 2013), wavelet denoising can be used to improve the performance of the classical ARIMA model for forecasting hydrological variables.

4 Summary and conclusions

In this study, the changes in the WBCs were explored using various forms of the wavelet-aided MK test. Moreover, a wavelet-aided ARIMA model was used for forecasting the WBCs. The results obtained from trend analysis indicated that decreasing trends were dominant in all the WBCs in the western part of Bangladesh during the period from 1982–1983 to 2012–2013. However, most of the trends were insignificant at the 95 % confidence level. One significant positive and one significant negative PET trend was found for the Satkhira and Bhola stations, respectively. Different combinations of the D and A (i.e., D + A and D + A + A) components of the DWT were analyzed using the Co value of the u(t) statistic from the SMK test, which provides detailed information regarding the dominant periodicity and time period affecting the trend of the original data (see the Trend and periodicity section or the example of the Bhola station). The findings of this study revealed that to obtain details regarding the time period responsible for the trends in the data, different combinations of components (D + A and D + A + A) must be analyzed rather than only the details (D) or approximation (A) components of the WT data. Moreover, this study indicated that the changes in temperature and rainfall were not only associated with the changes in the PET. To determine the attributes of PET changes, a detailed analysis must be conducted of all the relevant climatic variables. In the western part of Bangladesh, the D3 (8-year) and D4 (16-year) components had a dominant effect on the trends in the original WBC time series data. D2 (4-year) periodicity was also present in some cases, especially for the AET. Because surplus occurs during the monsoon season and most of the rainfall occurs during this season, the rainfall pattern may have a similar periodicity (D3 to D4).

Modeling of the study revealed that the WBC time series data was affected by noises from different hydrophysical interactions. As a result, the classic ARIMA model exhibited unsatisfactory performance in most of the cases (e.g., PET) or was unable to model the variability and changes in the AET, surplus, and deficit. This study indicated that the ARIMA model can be used to model the time series data of WBCs after denoising the data using DWT with a UT. The quality of the wavelet denoising time series data was evaluated, and satisfactory results were obtained for WBC data denoising. The performance of the fitted WD-ARIMA model was evaluated using the NSE and R2 values. The average NSE and R2 values of the 11 stations located in the western part of Bangladesh were 0.76 and 0.67, respectively, for the PET; 0.92 and 0.89, respectively, for the AET; 0.92 and 0.9, respectively, for the annual surplus; and 0.88 each for the annual deficit. The validation of the WD-ARIMA model for the period of 2009–2010 to 2012–2013 provided an acceptable EMP value. Thus, the WD-ARIMA model had an acceptable to very good performance for the short-term forecasting of WBCs. However, the gap between the actual and forecasted data increased with increasing time. The obtained results encourage further studies to determine a realistic model for real-world application under changing climate. The results of this study can be incorporated into water resource management plans for the highly irrigated western part of Bangladesh, where the groundwater resource is at a critical stage. Further studies regarding the denoising of hydrological time series data using different mother wavelets, such as Haar and Coiflet, and the determination of thresholds by using the MINMAX, SURE, or entropy-based adaptive denoising approaches would enable the development of superior models for forecasting hydroclimatic time series in the context of climate change and be beneficial for sustainably managing water resources.

Data availability
Data availability.

The national meteorological database of Bangladesh prepared by the Bangladesh Agricultural Research Council (BARC) was used to accomplish this study. Data are available for research and can be obtained from the BARC website (http://climate.barcapps.gov.bd/).

Supplement
Supplement.

Author contributions
Author contributions.

ATMSR designed and wrote the manuscript with input from all co-authors. MSA, MK, MAK, and ATMSR prepared the R code and ATMSR, MK, MAK, and MSA performed the statistical analysis. HMA and ATMSR performed the water balance analysis. QHM and CSJ supervised the whole work.

Competing interests
Competing interests.

The authors declare that they have no conflict of interest.

Special issue statement
Special issue statement.

Acknowledgements
Acknowledgements.

We thank the two anonymous reviewers for their constructive comments that greatly improved the manuscript. We would like to thank editor Ana Mijic of the special issue for her comments and support for publication in Hydrology and Earth System Sciences.

Edited by: Ana Mijic
Reviewed by: two anonymous referees

References

Acharjee, T. K., Halsema, G., Ludwig, F., and Hellegers, P.: Declining trends of water requirements of dry season Boro rice in the north-west Bangladesh, Agr. Water Manage., 180, 148–159, 2017.

Adamowski, J.: Development of a short-term river flood forecasting method based on wavelet analysis, Warsaw Polish Academy of Sciences, Monograph, 172, 2007.

Adamowski, J. and Chan, H. F.: A wavelet neural network conjunction model for groundwater level forecasting, J. Hydrol., 407, 28–40, 2011.

Ahasan, M. N., Chowdhary, M. A. M., and Quadir, D. A.: Variability and trends of summer monsoon rainfall over Bangladesh, J. Hydrometeorol., 7, 1–17, 2010.

Allen, R. G., Pereira, L. S., Raes, D., and Smith, M.: Crop evapotranspiration: guidelines for computing crop water requirements, FAO Irrigation and Drainage Paper, No. 56, Rome, Italy, 328 pp., 1998.

Alley, W. M.: Water balance models in one-month-ahead streamflow forecasting, Water Resour. Res., 21, 597–606, 1985.

Anderson, R., Hansen, J., Kukuk, K., and Powell, B. Development of watershed-based water balance tool for water supply alternative evaluations, Proceedings of the Water Environment Federation, (WEF`06), Water Environment Federation, 2817–2830, 2006.

Araghi, A., Baygi, M. M., Adamowski, J., Malard, J., Nalley, D., and Hasheminia, S. M. Using wavelet transforms to estimate surface temperature trends and dominant periodicities in Iran based on gridded reanalysis data, Atmos. Res., 155, 52–72, https://doi.org/10.1016/j.atmosres.2014.11.016, 2015.

Arnall, N. W.: Factors controlling the effects of climate change on river flow regimes in a humid temperate environment, Journal Hydrogeology, 132, 321–342, 1992.

Azad, S., Debnath, S., and Rajeevan, M. Analysing predictability in Indian monsoon rainfall: a data analytic approach, Environ. Process., 2, 717–727, 2015.

Bakundukize, C., Camp, M. V., and Walraevens, K.: Estimation of groundwater recharge in Bugesera region (Burundi) using soil moisture budget approach, Geol. Belg., 14, 85–102, 2011.

Boughton, W.: Catchment water balance modelling in Australia 1960–2004, Agr. Water Manage., 71, 91–116, 2004.

Box, G. E. P. and Jenkins, G. M.: Time Series Analysis: Forecasting and Control (Revised edition), San Francisco, Holden Day, 1976.

Chapman, L. C. and Brown, D. M.: The climates of Canada for agriculture. Canada Land Inventory Report No. 3, Environment Canada, Lands Directorate, 24, 1966.

Chou, C.: A threshold based wavelet denoising method for hydrological data modelling, Water Resour. Manag., 25, 1809–1830, https://doi.org/10.1007/s11269-011-9776-3, 2011.

de Jong, R. and Bootsma, A.: Estimates of water deficits and surpluses during the growing season in Ontario using the SWATRE model, Can. J. Soil Sci., 77, 285–294, 1997.

Donoho, D. L.: De-noising by soft-thresholding, IEEE T. Inform. Theory, 41, 613–627, 1995.

Donoho, D. L. and Johnstone I. M.: Ideal Denoising in an Orthonormal Basis Chosen from a Library of Bases, CR Acad. Sci. I-Math., 319, 1317–1322, 1994.

Donoho, D. L. and Johnstone, I. M.: Minimax estimation via wavelet shrinkage, Ann. Stat., 26, 879–921, 1998.

Fan, Z. and Thomas, A.: Spatiotemporal variability of reference evapotranspiration and its contributing climatic factors in Yunnan Province, SW China, 1961–2004, Climatic Change, 116, 309–325, https://doi.org/10.1007/s10584-012-0479-4, 2012.

Fulton, J. W., Risser, D. W., Regan, R. S., Walker, J. F., Hunt, R. J., Niswonger, R. G., Hoffman, S. A., and Markstrom, S. L.: Water-budgets and recharge-area simulations for the Spring Creek and Nittany Creek Basins and parts of the Spruce Creek Basin, Centre and Huntingdon Counties, Pennsylvania, Water Years 2000–06: U.S. Geological Scientific Investigations Report 2015–5073, 86 pp., https://doi.org/10.3133/sir20155073, 2015.

Hamed, K. H. and Rao, A. R.: A modified Mann–Kendall trend test for autocorrelated data, J. Hydrol., 204, 182–196, 1998.

Harvey, C. F., Ashfaque, K. N., Yu, W., Badruzzaman, A. B. M., Ali, M. A., Oates, P. M., Michael, H. A., Neumann, R. B., Beckie, R., Islam, S., and Ahmed, M. F.: Groundwater dynamics and arsenic contamination in Bangladesh, Chem. Geol., 228, 112–136, 2006.

Hasan, M. A., Islam, A. K. M. S., and Bokhtiar, S. M.: Changes of reference evapotranspiration ETo in recent decades over Bangladesh, 2nd International Conference on Advances in Civil Engineering, 26–28 December 2014 CUET, Chittagong, Bangladesh, 2014.

Healy, R. W., Winter, T. C., LaBaugh, J. W., and Franke, O. L.: Water Budgets: Foundations for Effective Water Resources and Environmental Management, United States Geological Survey, Reston, Virginia, 2007.

Huq, S. M. I. and Shoaib, J. U.: The Soils of Bangladesh, 1–172, https://doi.org/10.1007/978-94-007-1128-0, 2013.

Hyndman, R., Mitchell O'Hara, Wild, M., Bergmeir, C., Razbash, S., and Wang, E. R.: Language Forecast Package Development team, 131 pp., 2017.

IPCC (Inter-governmental Panel on Climate Change): Technical summary of climate change 2007: the physical science basis, Contribution of working group I to the fourth assessment report of the intergovernmental panel on climate change, edited by: Solomon, S., Qin, D., Manning, M., Marquis, M., Averyt, K., Tignor, M. M. B., Miller Jr., H. L., and Chen, Z., Cambridge University Press, Cambridge, 2007.

Kamruzzaman, M., Rahman, A. T. M. S., Kabir, M. E., Jahan, C. S., Mazumder, Q. H., and Rahman, M. S.: Spatio-temporal Analysis of Climatic Variables in the Western Part of Bangladesh, Environ. Dev. Sustain., 18, 89–108, https://doi.org/10.1007/s10668-016-9872-x, 2016a.

Kamruzzaman, M., Kabir, M. E., Rahman, A. T. M. S., Mazumder, Q. H., Rahman, M. S., and Jahan, C. S.: Modeling of Agricultural Drought Risk Pattern using Markov Chain and GIS in the Western Part of Bangladesh, Environ. Dev. Sustain., 18, 569–588, doi 10.1007/s10668-016-9898-0, 2016b.

Kang, S., Gu, B., Du, T., and Zhang, J.: Crop coefficient and ratio of transpiration to evapotranspiration of winter wheat and maize in a semi-humid region, Agr. Water Manage., 59, 239–254, 2003.

Kanoua, W. and Merkel, B. J.: Groundwater recharge in Titas Upazila in Bangladesh, Arab. J. Geosci., 8, 1361–1371, https://doi.org/10.1007/s12517-014-1305-2, 2015.

Karim, M. R., Ishikawa, M., and Ikeda, M.: Modeling of seasonal water balance for crop production in Bangladesh with implications for future projection, Ital. J. Agron., 7, e21, https://doi.org/10.4081/ija.2012.e21 2012.

Karimi, P., Bastiaanssen, W. G. M., and Molden, D.: Water Accounting Plus (WA+) – a water accounting procedure for complex river basins based on satellite measurements, Hydrol. Earth Syst. Sci., 17, 2459–2472, https://doi.org/10.5194/hess-17-2459-2013, 2013.

Kendall, M. G.: Rank Correlation Methods, Griffin, London, 1975.

Khalek, M. A. and Ali, M. A.: Comparative Study of Wavelet-SARIMA and Wavelet-NNAR Models for Groundwater Level in Rajshahi District, IOSR Journal of Environmental Science, Toxicology and Food Technology (IOSR-JESTFT), 10, 1–15, 2016.

Kisi, O.: Stream flow forecasting using neuro-wavelet technique, Hydrol. Process., 22, 4142–4152, 2008.

Kumar, M., Denis, D. M., and Suryavanshi, S.: Long-term climatic trend analysis of Giridih district, Jharkhand (India) using statistical approach, Modeling Earth Systems and Environment, 2, 116, https://doi.org/10.1007/s40808-016-0162-2, 2016.

Leta, O. T., El-Kadi, A. I., Dulai, H., and Ghazal, K. A.: Assessment of climate change impacts on water balance components of Heeia watershed in Hawaii, J. Hydrol.-Regional Studies 8, 182–197. https://doi.org/10.1016/j.ejrh.2016.09.006, 2016.

Lutz, S. R., Mallucci, S., Diamantini, E., Majone, B., Bellin, A., and Merz, R.: Hydroclimatic and water quality trends across three Mediterranean river basins, Sci. Total Environ., 571, 1392–1406, 2016.

Mallat, S. G.: A theory for multi-resolution signal decomposition: the wavelet representation, IEEE T. Pattern Anal., 11, 674–693, 1989.

Mann, H. B.: Nonparametric tests against trend, Econometrica, 13, 245–259, 1945.

McCuen, R. H., Knight, Z., and Cutter, A. G.: Evaluation of the Nash–Sutcliffe Efficiency Index, J. Hydrol. Eng., 11, 597–602, 2006.

McVicar, T. R., Roderick, M. L., Donohue, R. J., Li, L. T., Van Niel, T. G., Thomas, A., Grieser, J., Jhajharia, D., Himri, Y., Mahowald, N. M., Mescherskaya, A. V., Kruger, A. C., Rehman, S., and Dinpashoh, Y.: Global review and synthesis of trends in observed terrestrial near-surface wind speeds: Implications for evaporation, J. Hydrol., 416–417, 182–205, https://doi.org/10.1016/j.jhydrol.2011.10.024, 2012.

Michael, H. A. and Voss, C. I.: Controls on groundwater flow in the Bengal Basin of India and Bangladesh: regional modeling analysis, Hydrogeol. J., 17, 1561–577, 2009.

Misiti, M., Misiti, Y., Oppenheim, G., and Poggi, J.: Wavelet Toolbox User's Guide, The Math Works, Inc., 1997.

Molden, D. and Sakthivadivel, R.: Water accounting to assess use and productivity of water, Int. J. Water Resour. D., 15, 55–71, https://doi.org/10.1080/07900629948934, 1999.

Moriarty, P., Batchelor, C., Abd-Alhadi, F., Laban, P., and Fahmy, H.: The Empowers Approach to Water Governance Guidelines, Methods and Tools. Jordan: Inter-Islamic Network on Water Resources Development and Management (INWRDAM), 2007.

Moussa, R.: When monstrosity can be beautiful while normality can be ugly: assessing the performance of event–based flood models, J. Hydrolog. Sci., 55, 1074–1084, 2010.

Nalley, D., Adamowski, J., and Khalil, B.: Using discrete wavelet transforms to analyze trends in streamflow and precipitation in Quebec and Ontario (1954–2008), J. Hydrol., 475, 204–228, 2012.

Nalley, D., Adamowski, J., Khalil, B., and Ozga-Zielinski, B.: Trend detection in surface air temperature in Ontario and Quebec, Canada during 1967–2006 using the discrete wavelet transform, Atmos. Res., 132–133, 375–398, 2013.

Nash, J. E. and Sutcliffe, J. V.: River flow forecasting through conceptual models, part I: a discussion of principles, J. Hydrol., 10, 282–290. https://doi.org/10.1016/0022-1694(70)90255-6, 1970.

Nasher, N. M. R. and Uddin, M. N.: Maximum and Minimum Temperature Trends Variation over Northern and Southern Part of Bangladesh, Journal of Environmental Science and Natural Resources, 6, 83–88, 2013.

Nievergelt, Y.: Wavelets Made Easy, Birkhäuser, 297 pp., 2001.

Nury, A. H., Hasan, K., and Alam, J. B.: Comparative study of wavelet-ARIMA and wavelet-ANN models for temperature time series data in northeastern Bangladesh, Journal of King Saud University – Science, 29, 47–61, 2017.

Nury, A. H., Hasan, K., Erfan, K. M., and Dey, D. C.: Analysis of Spatially and Temporally Varying Precipitation in Bangladesh, Asian Journal of Water, Environment and Pollution, 13, 15–27, https://doi.org/10.3233/AJW-160023, 2016.

Partal, T.: Modelling evapotranspiration using discrete wavelet transform and neural networks, Hydrol. Process., 23, 3545–3555, 2009.

Partal, T.: Wavelet transform-based analysis of periodicities and trends of Sakarya basin (Turkey) streamflow data, River Res. Appl., 26, 695–711, 2010.

Partal, T. and Küçük, M.: Long-term trend analysis using discrete wavelet components of annual precipitations measurements in Marmara region (Turkey), Phys. Chem. Earth, 31, 1189–1200, 2006.

Pathak, P., Kalra, A., and Ahmed, S.: Wavelet-Aided Analysis to Estimate Seasonal Variability and Dominant Periodicities in Temperature, Precipitation, and Streamflow in the Midwestern United States, Water Resour. Manag., 30, 4649–4665, https://doi.org/10.1007/s11269-016-1445-0, 2016.

Percival, D. B. and Walden, A. T.: Wavelet Methods for Time Series Analysis, 1106 Cambridge University Press, New York, 594 pp., 2000.

Pišoft, P., Kalvová, J., and Brázdil, R.: Cycles and trends in the Czech temperature series using wavelet transforms, Int. J. Climatol., 24, 1661–1670, 2004.

Rahman, A. T. M. S.: Sustainable Groundwater Management in the Context of Climate Change in Drought Prone Barind Area, NW Bangladesh, Unpublished Master of Philosophy Thesis, Institute of Environmental Science, University of Rajshahi, Bangladesh, 2016.

Rahman, A. T. M. S., Jahan, C. S., Mazumder, Q. H., Kamruzzaman, M., and Hossain, A. Evaluation of spatio-temporal dynamics of water table in NW Bangladesh: An integrated approach of GIS and Statistics, Sustainable Water Resource Management, 2, 297–312, https://doi.org/10.1007/s40899-016-0057-4, 2016.

Rahman M. A., Yunsheng, L., and Sultana, N.: Analysis and prediction of rainfall trends over Bangladesh using Mann–Kendall, Spearman's rho tests and ARIMA model, Meteorol. Atmos. Phys., 129, 409–424, https://doi.org/10.1007/s00703-016-0479-4, 2016.

Rahman, M. J. and Hasan, M. A. M.: Performance of Wavelet Transform on Models in Forecasting Climatic Variables, in: Computational Intelligence Techniques in Earth and Environmental Sciences, edited by: Islam, T., Srivastava, P., Gupta, M., Zhu, X., and Mukherjee, S., Springer, Dordrecht, 2014.

Rahman, M. R. and Lateh, H.: Spatio-temporal analysis of warming in Bangladesh using recent observed temperature data and GIS, Clim. Dynam., 46, 2943–2960, 2016.

Ramana, R. V., Krishna, B., Kumar, S. R., and Pandey, N. G.: Monthly rainfall prediction using wavelet neural network analysis, Water Resour. Manag., 27, 3697–3711, 2013.

Ravenscroft, P., Burgess, W. G., Ahmed, K. M., Burren, M., and Perrin, J.: Arsenic in groundwater of the Bengal Basin, Bangladesh: Distribution, field relations, and hydrogeological setting, Hydrogeol. J., 13, 727–751, 2005.

R Core Team: A language and environment for statistical computing, R Foundation for Statistical Computing, Vienna, Austria, 2016.

Ritter, A. and Muñoz-Carpena, R.: Performance evaluation of hydrological models: Statistical significance for reducing subjectivity in goodness-of-fit assessments, J. Hydrol., 480, 33–45, https://doi.org/10.1016/j.jhydrol.2012.12.004, 2013.

Sang, Y.-F.: A practical guide to discrete wavelet decomposition of hydrologic time series, Water Resour. Manag., 26, 3345–3365, 2012.

Sang, Y. F.: A review on the applications of wavelet transform in hydrology time series analysis, Atmos. Res., 122, 8–15, https://doi.org/10.1016/j.atmosres.2012.11.003, 2013.

Sang, Y.-F., Wang, Z., and Liu, C.: Discrete wavelet-based trend identification in hydrologic time series, Hydrol. Process., 27, 2021–2031, 2013.

Santos, C. A. G and da Silva, G. B. L.: Daily streamflow forecasting using a wavelet transform and artificial neural network hybrid models, Hydrolog. Sci. J., 59, 312–324, 2014.

Scott, C. A. and Sharma, B.: Energy supply and the expansion of groundwater irrigation in the Indus-Ganges Basin, International Journal of River Basin Management, 7, 1–6, 2009.

Sen, P. K.: Estimates of the regression coefficient based on Kendall's tau, J. Am. Stat. Assoc., 63, 1379–1389, 1968.

Shahid, S. and Behrawan, H.: Drought risk assessment in the western part of Bangladesh, Nat. Hazards, 46, 391–413, 2008.

Shahid, S. and Khairulmaini, O. S.: Spatial and temporal variability of rainfall in Bangladesh, Asia-Pac. J. Atmos. Sci., 45, 375–389, 2009.

Shahid, S.: Recent trends in the climate of Bangladesh, Clim. Res., 42, 185–193, 2010.

Shamsudduha, M., Chandler, R. E., Taylor, R. G., and Ahmed, K. M.: Recent trends in groundwater levels in a highly seasonal hydrological system: the Ganges-Brahmaputra-Meghna Delta, Hydrol. Earth Syst. Sci., 13, 2373–2385, https://doi.org/10.5194/hess-13-2373-2009, 2009.

Smith, L. C., Turcotte, D. L., and Isacks, B. L.: Stream flow characterization and feature detection using a discrete wavelet transform, Hydrol. Process., 12, 233–249, 1998.

Sneyers, R.: On the Statistical Analysis of Series of Observations, Secretariat of the World Meteorological Organization, 192 pp., 1990.

SRDI (Soil Resources Development Institute): Soil map of Bangladesh, Soil Resources Development Institute, 1998.

Sreekanth, P., Geethanjali, D. N., Sreedevi, P. D., Ahmed, S., Kumar, N. R., and Jayanthi, P. D. K.: Forecasting groundwater level using artificial neural networks, Curr. Sci., 96, 933–939, 2009.

Stein, C. M.: Estimation of the Mean of a Multivariate Normal-Distribution, Ann. Stat., 9, 1317–1322, 1981.

Syed, M. A. and Al Amin, M.: Geospatial Modeling for Investigating Spatial Pattern and Change Trend of Temperature and Rainfall, Climate, 4, 21, https://doi.org/10.3390/cli4020021, 2016.

Tabari, H. and Talaee, P. H.: Moisture index for Iran: spatial and temporal analyses, Global Planet. Change, 100, 11–19, 2013.

Thornthwaite, C. W. and Mather, J. R.: The Water Balance, Publications in Climatology, VIII(1), 1–104, Drexel Institute of Climatology, Centerton, New Jersey, 1955.

Thornthwaite, C. W. and Mather, J. R.: Instructions and tables for computing potential evapotranspiration and the water balance, Publications in Climatology, Laboratory of Climatology, Drexel Institute of Technology, Centerton, New Jersey, USA, 10, 183–311, 1957.

Tiwari, M. K. and Chatterjee, C.: Development of an accurate and reliable hourly flood forecasting model using wavelet–bootstrap–ANN (WBANN) hybrid approach, J. Hydrol., 1, 458–470, 2010.

Torrence, C. and Compo, G. P.: A practical guide to wavelet analysis, B. Am. Meteorol. Soc., 79, 61–78, 1998.

Valipour, M.: Ability of Box-Jenkins Models to Estimate of Reference Potential Evapotranspiration (A Case Study: Mehrabad Synoptic Station, Tehran, Iran), IOSR Journal of Agriculture and Veterinary Science, 1, 1–11, 2012.

Valipour, M., Banihabib, M. B., and Behbahani, M. R.: Comparison of the ARMA, ARIMA, and the autoregressive artificial neural network models in forecasting the monthly inflow of Dez dam reservoir, J. Hydrol., 476, 433–441, 2013.

Viaroli, S., Mastrorillo, L., Lotti, F., Paolucci, V., and Mazza, R.: The groundwater budget: a tool for preliminary estimation of the hydraulic connection between neighboring aquifers, J. Hydrol., 556, 72–86, doi.org/10.1016/j.jhydrol.2017.10.066, 2017.

Wang, W., Ding, J., and Li, Y.: Hydrologic Wavelet Analysis, Chemical Industry Press, Beijing, China, 2005 (in Chinese).

Wang, D., Singh, V. P., Shang, X., Ding, H., Wu, J., Wang, L., Zou, X. Chen, Y., Chen, X., Wang, S., and Wang, Z.: Sample entropy based adaptive wavelet de-noising approach for meteorological and hydrologic time series, J. Geophys. Res.-Atmos., 119, 8726–8740, https://doi.org/10.1002/2014JD021869, 2014.

Wolock, D. M. and McCabe, G. J.: Effects of potential climatic change on annual runoff in the conterminous United States, J. Am. Water Resour. As., 35, 1341–1350, 1999.

Xu, C.-Y. and Halldin, S.: The effect of climate change in river flow and snow cover in the NOPEX area simulated by a simple water balance model, Proc. of Nordic Hydrological Conference, Alkureyri, Iceland, 1, 436–445, 1996.

Yoon, H., Jun, S. C., Hyun, Y., Bae, G., and Lee, K. K.: A comparative study of artificial neural networks and support vector machines for predicting groundwater levels in a coastal aquifer, J. Hydrol., 396, 128–138, 2011.

Young, P. C.: Nonstationary time series analysis and forecasting, Progress in Environmental Science, 1, 3–48, 1999.

Yue, S., Pilon, P., Phinney, B., and Cavadias, G.: The influence of autocorrelation on the ability to detect trend in hydrological series, Hydrol. Process., 16, 1807–1829, 2002.

Special issue