Assessing and forecasting of groundwater level fluctuation in Joypurhat district, northwest Bangladesh, using wavelet analysis and ARIMA modeling

Groundwater resource plays a crucial role for agricultural crop production and socio-economic development in some parts of the world including Bangladesh. Joypurhat district, the northwest part of Bangladesh, a crop production hub, is entirely dependent on groundwater irrigation. A precise assessment and prediction of groundwater level (GWL) can assist long-term groundwater resources (GWR) management, especially in drought-prone agricultural regions. Therefore, this study was carried out to identify trends and magnitude of GWL fluctuation (1980–2019) using the modified Mann–Kendall test, Pettitt’s test, and Sen slope estimators in the drought-prone Joypurhat district, northwest Bangladesh. Time-series data analysis was performed to forecast GWL from 2020 to 2050 using the Auto-Regressive Integrated Moving Average (ARIMA) model. The findings of the MMK test revealed a significant declining trend of GWL, and the trend turning points were identified in the years 1991, 1993, 1997, and 2004, respectively. Results also indicate that the declining rate of GWL varied from 0.104 to 0.159 m/year and the average rate of GWL declination was 0.136 m/year during 1980–2019. The outcomes of wavelet spectrum analysis depicted two significant periods of the declining trend in Khetlal and Akkelpur upazilas. The results obtained from the optimal identified model ARIMA (2, 1, 0) indicate that GWL will decline at a depth of 13.76 m in 2050, and the average declination rate of GWL will be 0.143 m/year in the study area. The predicted results showed a similar declining tendency of GWL from 2020 to 2050, suggesting a disquieting condition, particularly for Khetlal upazila. This research would provide a practical approach for GWL assessment and prediction that could help decision-makers implement long-term GWR management in the study area.


Introduction
Groundwater is one of the most valuable natural resources in the world. It plays a pivotal role in the economic, social, and environmental sustainability context. The potential use of groundwater, particularly in agricultural-based developing countries like Bangladesh, is of paramount concern to the policy-makers on groundwater resources (GWR) management perspective (Morris et al. 2003;Mackay et al. 2015;Salam et al. 2020a;Islam et al. 2021a). In recent times, climate variability (e.g., rainfall infiltration rate, surface runoff, evaporation, and increase in temperature), rapid population increase, and over-exploitation have had a detrimental effect on this valuable GWR in many drought-prone regions worldwide (Shahid and Hazarika 2010a, b;Kalhor and Emaminejad 2019;Yadav et al. 2020;Ajibade et al. 2021). This is primarily true for the northwest region of Bangladesh, where high population density and increasing agricultural crop production raise water demand (Shahid 2011a, b;Husna et al. 2016;Islam et al. 2022). As a result, GW levels (GWL) are declining, and there is rising water scarcity in the droughtprone Barind tract region. However, groundwater scarcity is another challenge around the globe. The investigation for sustainable GWR management includes GWL assessment and prediction. In addition to management concerns, the long-term prediction will be helpful to close the missing dataset too. So, accurate GWL assessment and prediction are important for effective and sustainable GWR management in water-stressed areas, especially in the drought-prone northwest of Bangladesh.
The use of groundwater in the irrigation sector for crop production is increasing significantly in Bangladesh. In such a case, the northwest region of Bangladesh is facing a critical condition of groundwater due to the impact of overexploitation. Abstraction of groundwater from deep tube wells has resulted in substantial increases in agricultural production and significant falls in aquifers' groundwater heads (Patle et al. 2015;Rushton et al. 2020;Mallick et al. 2021). More than 75% of irrigation water is managed by groundwater withdrawal in the northwest part of Bangladesh (Shahid and Hazarika 2010a, b). During the dry season, most surface water resources like rivers, canals, beels (static water bodies), ponds, and other reservoirs get dried out because of limited rainfall and other adverse climatic factors that compel people to use groundwater in the northwest part of Bangladesh. When groundwater withdrawal exceeds the natural recharge to groundwater storage for a long time in a considerable area, the declination of GWL occurs. The falling trends of groundwater levels indicate that groundwater exploitation exceeds recharge, indicating unsustainable withdrawal (Zannat et al. 2019). The GWL decline rate is higher in the northern region than in the southern part (Sumiya and Khatun 2016). The groundwater table has reduced noticeably in the northwest part of Bangladesh due to the overexploitation of groundwater (Dey et al. 2017).
Joypurhat district is the country's agricultural hub. In this area, farmers use a lot of irrigation equipment to withdraw groundwater for irrigation purposes to boost cropping intensity (Islam et al. 2018). Rice and potatoes are the main crops in this area, which are both highly water-consuming crops, and their cultivation relies on groundwater irrigation. The groundwater table in most of the Kalai, Khetlal, Panchbibi, and Akkelpur upazilas goes below the suction limit, causing hand tube-wells and shallow tube-wells to be partially or fully operable (Hossain et al. 2015). Intensive agricultural, rapid industrial development, and domestic demand have put more pressure on groundwater use and accelerated the rate of water table declination in this area. These conditions demand proper monitoring and investigations into GWL depletion. The trend analysis will help in understanding the variation, changing direction, and trend of GWL and GWL projection for future management.
Few researchers in northern Bangladesh have studied GWL depletion at the natural and regional levels Jahan et al. 2010;Rahman et al. 2016;Dey et al. 2017;Salam et al. 2020a;Elbeltagi et al. 2022;Pham et al. 2022;Di Nunno et al. 2022). Many research scholars have investigated the reasons and magnitudes of GWL depletion in NW Bangladesh. For example, Zafor et al. (2017) concentrated on the falling trend direction in GWL, which implies an unbalanced groundwater setting in northern Bangladesh. Zahid (2015) reported that the tendency of GWL fluctuation is elevated in Bangladesh's northwest Barind area. Rahman et al. (2016) stated that there was a rising tendency of GWL in the Barind tract region. Besides, Dey et al. (2017) reported that the GWL has decreased noticeably in northwest Bangladesh due to over-withdrawn groundwater, increasing Boro rice farming, decreasing river water levels, lessening in wetland area, and scarcity of water reservoirs. Although there has been a rapid rise in the irrigated cropping area during the last two decades, poorer performance and lack of groundwater efficacy and surface water utilization are still hindrances to substantial GWR management. Recently, Pham et al. (2022) used machine learning algorithms to predict groundwater levels in a drought-prone region. Di  tested the predictive capabilities of non-linear autoregressive with exogenous inputs (NARX) and extreme learning machine (ELM) for GWL forecasting in northern Bangladesh. Elbeltagi et al. (2022) applied the hybrid ensemble modeling framework to estimate GWL in the northern region of Bangladesh. However, its long-term viability is jeopardized due to the low quality and quantity of groundwater in northwest Bangladesh (Shahid 2011a, b). Yet, irregular trend pattern and changing characteristics of GWL are less understood in the vital drought-prone region. Therefore, it is essential to have a better insight into the unusual behavior, changing patterns in GWL, and future situations of groundwater fluctuation in northwest Bangladesh.
When data is limited and producing an accurate forecast is more essential than understanding underlying processes, statistical models are a viable option (Takafuji et al. 2019). Bidwell (2005), for example, created an autoregressive moving average (ARMA) model to forecast GWL in New Zealand. Using an ARMA model, Lee et al. (2009) forecasted the GWL in Changwon, Korea. Salam et al. (2020a) used an ARIMA model to generate monthly GWL at a local scale in Bangladesh's northwestern area. However, the ARIMA model typically has difficulties with nonstationary data and cannot handle nonstationary data without preprocessing the input data (Zhang et al. 2021;Islam et al. 2021b). Nonstationary data handling techniques are not as sophisticated as those for stationary data, and further study is required to figure out ways that can handle nonstationary data more successfully. Wavelet analysis may be used in this context. Previous research (Jerin et al., 2021;Ghose et al., 2021) shows that wavelet analysis works best when studying nonstationary time series.
No such comprehensive study has been performed before in the drought-prone region, part of the Barind Tract, Joypurhat district, northwest Bangladesh. This type of study is urgently required due to severe water stress. Based on earlier literature, the two goals have been devised to overcome the knowledge gaps. First, the modified Mann-Kendall test, Sen slope estimator, Pettitt's test, and Wavelet power spectrum analyses have been applied to assess the spatial and temporal trend patterns of GWL variation in Joypurhat district, northwest Bangladesh. On the other hand, statistical analyses require a more in-depth understanding of future GWL characteristics. Due to the oversimplification of complicated hydrological mechanisms and the scarcity of groundwater datasets, achieving precise forecasting through statistical methods, such as the Auto-Regressive Integrated Moving Average (ARIMA), is a promising task Akhter et al. 2019;Salam et al. 2020a). Second, this study aims to forecast the GWL from 2020 to 2050 in Joypurhat district using ARIMA modeling.
However, based on the earlier discussion, the following may be the most significant novel aspect of the current research: Firstly, the finding contributes to knowledge by developing and executing GWL forecasting methods in a droughtprone region of northwest Bangladesh with high groundwater usage.
Secondly, improved knowledge of groundwater forecasting in Bangladesh's drought-prone northwest. This endeavor would provide a strong basis for earth scientists, government officials, and stakeholders to improve water resource management and agricultural sustainability.
Finally, an ARIMA model for GWL forecasting was presented, which combines multi-statistical approaches to evaluate GWL. These approaches, to the best of the author's knowledge, have never been applied in prior research on GWL evaluation and forecasting.

Study area description
The study area has been selected for the study area located in the Barind tract in the northwest region of Bangladesh. It is situated between 24°46'N to 25°22'N latitude and 88°56' E to 89°14'E longitude, as shown in Fig. 1. Joypurhat district is bounded on the north by Gaibandha district, Dinajpur district, and the India border; on the south by Bogra district and Naogaon district; on the east by Bogra district and Gaibandha district; and on the west by Naogaon district and the India border. The total area of the district is 965.88 km 2 (Islam et al. 2018). It consists of five upazilas, namely, Joypurhat Sadar, Akkelpur, Khetlal, Kalai, and Panchbibi. The economy of Joypurhat district is entirely dependent on agriculture and is famous for its granaries. There are 1978 nos. deep tube-wells (DTW) and 8050 nos. shallow tubewells (STW) used for irrigation purposes in agricultural production (BADC 2019). The study region shows a rich agricultural hub of about 80,270 ha of net cultivatable site, which depends on groundwater irrigation. Rice and potatoes are the main cash crops in this area. Both crops are highly water-consuming in an irrigated field. Water scarcity is a common problem during the dry winter period in the study area.

Geology and hydrological conditions
The Joypurhat district is the portion of the Tista floodplain that belonged to the shelf region (Bogra slope), encompassing the Himalayan Foredeep region of the Bengal Basin (Reimann and Hiller 1993). The study area lies in the subhumid area bounded by geological formations deposited from the Dupi Tila Sandstone Formation. The sediment is deposited as an older alluvial fan where different rivers flow through the study areas, like Chiri River and Haraboti River, which cross via Panchbibi upazila. The small Jamuna River crosses via Joypurhat Sadar and Panchbibi upazilas. The Tulshiganga River passes via Joypurhat Sadar, Khetlal, and Akkelpur upazilas. The Sree River crosses via Joypurhat Sadar upazila.
The topography of the study area varies from 24 m in Panchbibi to 14.4 m in Akkelpur. The study area is relatively flat, sloping from north to south. The topsoil of the study area is mainly clay and the aquifer properties like transmissivity vary from 1240 to 1700 m2/day (IWM 2009). The principal aquifers of the study area consist of fine sands, medium sands, and coarse sands with mild stone particles, as shown in Fig. 2(a-e). The study area is covered in a thick clay sheet with Plio-Pleistocene and Holocene ages . The clay layer serves as an aquitard, increasing surface runoff and restricting groundwater recharge. The soil texture in Joypurhat district is mainly silty loam. There are 2 types of aquifers found in the study area, which are as follows: first, a semi-confined shallow aquifer located immediately below the top soil and ranging from 5 to 25 m, second, a lower shallow aquifer extending from 15 to 60 m and composed of coarse-medium grained sand, silt, and fine material (Islam et al. 2018). Considering the lithological characteristics of five selected wells in the study area, the average thickness of the aquifers in the Joypurhat district varies from about 40 to 51 m. The screenable thickness of the study area ranges from 29.39 to 35.06 m, with a bore depth of up to 65 m from the ground surface. Joypurhat is a district with a sub-tropical climate. The average annual temperature is 24.4 °C, and the average annual rainfall is 1673 mm. More than 60% of the rain occurs in July-August, with an average of 364 mm (Das and Islam 2021). The pre-monsoon period occurs from March to May, and the monsoon period occurs from June to October. The warmest month of the year is May, with a mean temperature of 28.9 °C (Das 2021). The coldest month of the year is January, when the average temperature is 18 °C and the average annual evaporation of the study area is around 1060 mm (Salam et al. 2020b).

Data source and quality control check
The monthly GWL data (1980-2019) of the various Joypurhat district stations were obtained from the Bangladesh Water Development Board (BWDB), Dhaka. Because there are no BWDB stations in Akkelpur and Kalai upazilas for measuring GWL data , the nearest stations in the study area, Parbotipur and Nischinta, are used for Akkelpur and Kalai upazilas, respectively. The study area is based on monthly rainfall and temperature data from the Bangladesh Meteorological Department (BMD). In all cases, systematic quality control techniques are applied to evaluate data quality. In the Buishand Range test, Normality Test, and Homogeneity Test, all of the data used in the study area is significant at 99% confidence level (p-value < 0.001). These tests also revealed that the data remains of good quality.

Modified Mann-Kendall test
The modified Mann-Kendall test is the modified version of the non-parametric Mann-Kendall test, which is a robust tool in the presence of autocorrelation. This test may be used to identify trend detection in hydrometeorological data series when datasets are non-random and impacted by autocorrelation (Mann 1945;Kendall 1975;). Hamed and Rao (1998) first developed a variance correction method to address serial correlation in trend assessment. Praveen et al. (2020) have also suggested the same approach in this respect. Datasets are primarily detrended, and the adequate sample size is computed using the ranks of significant serial correlation analysis, which is then used to correct the test statistic's inflated (or deflated) variance. The Mann-Kendall test statistics S is given as: where, S is Mann-Kendall statistic and sgn is the signum function. Each of the data point xi is taken as a reference point which is compared with the rest of the data points xj as follows: The value S shows the direction of trend. A positive value of S indicates an upward trend and negative value indicates downward trend (Islam et al., 2021a;Jerin et al. 2021). The test statistics S is expressed as: where, t i is the number of ties up to sample i. The modified variance V * (S) is given by the following equation proposed by Yue and Wang (2004): where n n * is considered as correction factor and it is calculated by the following equation derived by Bayley and Hammersley (1946): In time series analysis, the presence of a statistically significant trend is calculated using Zc value.
In this method, a positive value indicates an increasing trend and a negative value indicates decreasing trend. The trend is considered as statistically significant in this study when p ≤ 0.05 at 5% significance level.

Sen's slope estimator
Sen's slope estimator test is a non-parametric method (Sen, 1968) which is used to estimate the magnitude of trends (change per unit time) in time series data. The magnitude of the trend is expressed as Q calculated by Eq. (7): where, X j and X k are denoted by data values at times j and k (j > k), respectively. A positive value of Q indicates an increasing trend and negative value indicates decreasing trend in time series data.

Pettitt's test
The Pettitt's test is a common non-parametric technique which was used to detect an abrupt change point in hydrological time series datasets (Pettitt 1979). The non-parametric statistics is defined as follows: K T is the location point where the abrupt change point of time series data occurred. The probability of significance level of K T is approximately lied on p ≤ 0.05 followed by Eq. (10).

Wavelet transform analysis
Wavelet analysis is a common tool to analyze local spectral and temporal variations of power within a time series for a particular scale and location. It contains non-stationary power at various frequencies (Jerin et al. 2021;Rahman et al. 2021). Depending on time and non-dimensional frequency, the Morlet wavelet consists of a plane wave modulated by Gaussian following Eq. (11).
where Ψ 0 ( ) is denoted by wavelet value at non-dimentional time and 0 is the non-dimensional frequency, here taken to be (Eq. 6) in this study in order to satisfy the admissibility condition (Farge 1992). According to the function of basic wavelet, the scaled wavelets are ascertained as follows: where S is dilation parameter used to change the scale and n is the translation parameter used to slide in time. The wavelet transform W n (S) is expressed as the convolution of the wavelet function which is defined by the following Eq. (13): where Ψ * indicates the complex conjugate and n is the localized time index. The wavelet-based approach reduced the majority of the "noisy" data and rendered variations more reliable.

ARIMA modeling
The Auto-Regressive Integrated Moving Average (ARIMA) model is the popular model for modeling time series data and predicting future data in a series (Box and Jenkins, 1976;Rahman et al., 2017;Akhter et al., 2019;Islam et al., 2021a). In this study, a non-seasonal ARIMA model is subjected to understanding past data for forecasting. It consists of three statistical parts, such as autoregressive (AR) terms, integrated (I) terms, and moving average (MA) terms. These terms are classified as ARIMA (p, d, q) in which p is an autoregressive part that forecasts future values based on past values, d is an integrated part used for removing non-stationarity in forecasting, and q is the moving average part used for maintaining lagged forecast errors in the prediction system. In this study, the time series of GWL data (denoted by Y) was used to assign appropriate models and forecasts. In terms of Y, the general forecasting equation could be expressed as follows with the convention introduced by Box and Jenkins (1976).
where Ŷ t is the forecasted value at time t, Y t-1 is the previous forecasted value, is constant term, φ is autoregressive parameter, and is the moving average parameter. The Box and Jenkins method of ARIMA modeling follows four steps to choose the best model and forecast. These steps are model identification, parameter estimation, diagnostic checking, and forecasting. In the first step, the time series must be standardized and normalized by proper order differencing. The most popular method for identifying the appropriate demands of an ARIMA model is the adjustment behavior of the autocorrelation function (ACF) and partial autocorrelation function (PACF). In the second step, model parameters are estimated by the maximum likelihood method. Model adequacy is performed with diagnostic checking if the errors of model assumptions are satisfied. When the model seems inadequate, a new tentative model is assigned. The steps described above are repeated until the model assumption associated with the error is satisfied, and a satisfactory model is finalized. Finally, the finally selected model is used for ARIMA fitting and forecasting.
In the present study, SPSS software was used for ARIMA modeling. GWL time series data were partitioned in a 30:10 ratio. The first part of time series data from 1980-2009 was used for model identification, and (14) the second part of the time series from 2010-2019 was used for model validation and model performance and predictability testing. First, the Bayesian Information Criteria (BIC) is first considered for the efficient evaluation of model performance. Lastly, the most appropriate forecasting model was finally chosen based on model fit statistics values of R2, RMSE, MAPE, MaxAPE, MaxAE, and Ljung-Box statistics. The main justification for employing the ARIMA model is that it is the most common and broadly used stochastic model, and it has been effective in groundwater resource studies (Gibrilla et al. 2017;Zhang et al. 2021) and other cited works Akhter et al. 2019;Islam et al., 2021a). The ARIMA model has received a lot of interest because of its ease of learning and use in characterizing many kinds of time-series datasets (Lee et al. 2009).

Variation in climatic parameters
Rainfall, temperature, and evaporation distributions play an essential role in GWL fluctuating in a specific area. Figure 3(a-c) shows the monthly maximum and minimum temperature and rainfall distributions of the northwest region, Joypurhat district. The hottest months are identified as April to May, and the coldest month is January. The highest rainfall occurred in July, and the lowest rainfall occurred in January. Heavy rainfall is observed at temperatures of 24-35 °C in July. It is also noticeable that GWL started depletion in the hottest month due to high temperatures and evaporation differences. When temperature and evaporation decreased at a favorable condition for increasing rainfall, groundwater level began to rise toward the ground surface once more.

Apparent variation in groundwater level
Consecutive records of GWL are essential for assessing groundwater fluctuation and its trend, as well as groundwater potential. The observed GWL data generally represents the apparent variation of the groundwater table in a study area. Figure 4(a-e) is generated to explain the noticeable variation of GWL. It is observed that the groundwater table is highly depleted in April due to groundwater pumping for irrigation, high temperature, evapotranspiration, and other related climatic factors. But GWL is raised towards the ground surface in September for high monsoon rainfall, which enhances the recharging of groundwater storage. The annual groundwater depth variation limit is set at a maximum of 7.07 m in Khetlal and a minimum of 4.6 m in Panchbibi upazila, according to Fig. 4(a-e). It was also observed that the highest GWL declination from the ground surface is 12.82 m in Khetlal and the lowest GWL declination is 2.35 m in Akkelpur upazila. It can be mentioned that the amount of groundwater pumping, temperature, evapotranspiration, and rainfall distribution might affect the apparent variation of GWL. During the dry season, the adverse climate and the rate of groundwater pumping for irrigation reach an extreme level, and rainfall is limited; thus, GWL declination is increased.
On the other hand, favorable climatic conditions, limited groundwater irrigation, and heavy rainfall during the wet season create an auspicious situation during the wet season, so GWL declination started decreasing, and GWL rose towards the ground surface successively. Therefore, GWL is highly depleted in April and increases towards the ground surface in September and October. The highest decreasing rate of GWL was found in Khetlal upazila. The rate of decline was the slowest in Panchbibi upazila, which is in the Joypurhat district in the northwest of Bangladesh.

Trend of groundwater level fluctuation
The modified Mann-Kendal test is carried out on time-series data (1980-2019) on a monthly scale. The modified Mann-Kendal test is carried out on time-series data (1980-2019) on the monthly GWL in selected stations of different upazilas of Joypurhat district. To find out the trend in time series GWL, MMK "Zc" statistics were determined and are shown  Table 2. In all cases, "Zc" values are at a 5% level of significance. So, it can be said that the null hypothesis was rejected, and there was found a significant trend in time series data of groundwater levels. The positive sign "Zc" statistics suggested that GWL was declining for the ground surface.
The "Zc" statistics values obtained from the MMK test varied from 7.37 to 17.78 in different stations, respectively. It is observed that the lowest declining trend (Zc = 17.78) of GWL level is found at Panchbibi upazila. The comparatively highest declining trend (Zc = 7.37) is found at Khetlal upazila of the study area, Fig. 5 a GWL declining trend analysis based on MMK test; b GWL trend slope analysis; and c monthly GWL trend analysis  which is shown in Fig. 5a. After determining the trend, Sen's slope estimator was applied to calculate the slope of trend/rate of groundwater level declination (Q) in m/decade (Table 1) Figure 5b shows that the high rate of GWL level declination (Sen's slope, Q) is found at Ketal and Kalai upazilas; a medium changing rate at Joypurhat Sadar and Akkelpur; and a low changing rate at Panchbibi upazila, which is illustrated in Fig. 5c. After determining trend status and trend slope, the change point of GWL was identified at the times of Dec 1991, Dec 1993, Jan 1997, Feb 1997, and Dec 2004 for Akkelpur, Joypurhat Sadar, Khetlal, Kalai, and Panchbibi upazila, respectively, following the result of Pettitt's test, which is shown in Fig. 6. The first GWL level declination trend was started at Akkelpur in 1991 and then at Joypurhat Sadar in 1993, Khetlal and Kalai in 1997, and Panchbibi in 2004, respectively. It can be said that the groundwater pumping in Khetlal and Kalai upazilas during the dry season is greater than in the other upazilas. Two upazilas are more vulnerable to groundwater resources than others. Figure 7(a-f) depicts the Morlet wavelet power spectrum of the GWL level. The black contour line identifies the significance level, and the black line marks the cone of influence. An analysis of the wavelet power spectrum shows that two significant periods of increasing trend were detected in Khetlal upazila in a scale range of 4.0-6.0 m from 1980 to 1992 and a scale range of 11.4-16 m from 1992 to 2019, and Akkelpur upazila in a scale range of 6.9-11.2 m from 1980 to 1993 and a scale range of 11.4-16 m from 1994 to 2019. In Joypurhat Sadar, Panchbibi, and Kalai upazilas, there is only one significant period of increased GWL level. However, one significant period of the GWL level increasing trend started from 1991 to 1999 in Joypurhat district.

ARIMA model selection and validation
Plots of the autocorrelation function (ACF) and partial autocorrelation (PACF) are crucial for ARIMA model selection. Based on plots of ACF and PACF, the orders p and q were decided to develop the tentative ARIMA models for the 1980-2019 time series of groundwater levels. For each station, ten tentative ARIMA models were selected with different orders of p, d, and q with logical and reasonable ranges. The tested models were as follows: (0, 0, 0), (0, 1, 1), (1, 1, 0), (1, 1, 1), (1, 1, 2), (2, 1, 1), (2, 1, 2), (0, 1, 2), (2, 1,0). To analyze the ACF and PACF residual plots within the confidence limit, some tentative suitable models for each station were selected primarily based on the lowest normalized Bayesian Information Criteria (BIC  Schwarz et al., 1978). It was also considered to test if the ACF residuals at different lags existed at different lags from zero, where not being different from zero was anticipated. In ARIMA Model selection, the higher the R squared value, the better, and the lower the rest of the values. R squared values typically range from 0 to 1. In most cases, the standard value of R squared is 0.80 or above, showing high-level correlation. A root mean squared value is accepted up to a maximum of 0.80. MAPE value is excellent when it is less than 10%, but 10-50% is acceptable. There is no correct value for MaxAPE, MAE, MaxAE, and Ljung-Box Staistics. But the lower these values, the better the model. These parameters are relative to the data set used. On considering the higher R squared value and the least value of RMSE, MAPE, MaxAPE, MAE, MaxAE, and Ljung-Box Statistics in selecting the best ARIMA model for each station, ARIMA (1, 1, 0), ARIMA (2, 1, 0), ARIMA (1, 1, 0), ARIMA (0, 1, 1), and ARIMA (0, 1, 0) were finally chosen respectively for Joypurhat Sadar, Panchbibi, Akkelpur, Kalai, and Khetlal shown in Table 2. The R squared, RMSE, and MAPE values were all within the allowable range in all cases of the finally chosen model, and other values were lowest in comparison to the nearest model. Similarly, the best model for the study area was identified as the ARIMA model (2, 1, 0) by evaluation of ACF and PACF plots and model goodness of fit statistics. This   Fig. 8(a-c). A comparison of observed GWL during 2010-2019 was accomplished for further validation of the selected ARIMA model for each station and Joypurhat district. The comparison plot of observed GWL vs. forecasted GWL during 2010-2019 of the Joypurhat district average is in good consistency with the obtained R squared value of 0.9672 for the Joypurhat district shown in Fig. 9. So, the selected ARIMA model for the study area was considered for forecasting GWL.

Forecasting of groundwater levels
The selected best models for each station and Joypurhat district were applied to forecast GWL as shown in Fig. 10(a-f). Beginning in 2020, GWL is expected to fall linearly in Joypurhat Sadar, Panchbibi, Akkelpur, Kalai, and Khetlal upazilas. These forecasting results reveal a sharp trend of increasing from the ground surface in Khetlal, Kalai, and Akkelpur upazilas. Joypurhat Sadar and Panchbibi upazila have a moderate trend of GWL going down from the ground surface.
Overall, a moderate increasing trend of GWL levels from the ground surface was found in the Joypurhat district. As shown in Fig. 10(a-c) In Joypurhat district, the observed average rate of declination in GWL level was found to be 0.136 m/year from 1980 to 2019 as per Sen's estimator. The forecasted rate of declination during 2020-2050 was determined at 0.143 m/year, which shows a moderate increasing trend. The climate change effect, groundwater irrigation, and rainfall distribution were evaluated to examine the predicted results. There was a slight adverse effect on rainfall distribution and climate. However, the pressing question for increased groundwater irrigation in the study area's agricultural development may be the major factor influencing groundwater level declination trend. The main crops in this area are rice and potatoes, which can only be grown with irrigation from groundwater.

Discussion
This study explores the GWL fluctuation in a drought-prone Barint tract, namely, the Joypurhat district of northwest Bangladesh. GWL in the Joypurhat district of northwest Bangladesh is declining gradually. It is observed that the declining trend of GWL is more significant in Khetlal and Kalai upazilas than in the other upazilas. The main reason is that the flow of the Atria River occurs in the upper part of the Joypurhat district. The base flow happens in the SW direction because of the average slope, triggering the minimum aquifer recharge from river water in this region. Besides, the soil texture in this region is dominated by silty clay soil, which has relatively lower hydraulic conductivity characteristics (Zannat et al. 2019). Consequently, groundwater recharge from rainfall is very low.
The rates of GWL change ranged from 1.04 to 1.59 m/ decade during 1980-2019 in Joypurhat district. Salam et al. (2020a) found that the rates of GWL change varied from 0.1 to 1.3 m/decade during 1981-2017 in the northwest part of Bangladesh. The outcome is in line with Zahid (2015) and Rahman et al. (2016). The highest change point was detected at Khetlal upazila in 1997. After the 2000s, the GWL declined rapidly due to climatic effects (Islam et al. 2021b). This may have happened due to huge Boro rice and potato cultivation dependent on groundwater irrigation. Wavelet analysis exhibits two significant periods of GWL in Khetlal and Akkelpur, whereas Pahnchbibi shows no considerable period of GWL. The Khetlal upazila reveals a drier condition than the wet condition of groundwater. The results of the ARIMA model reveal a rising trend of declining GWL in Khetlal, and the trend of declining GWL is lower in Panchbibi than in Kalai. The rate of GWL decline from 2020 to 2050 will be 0.143 m/decade compared to the present period (0.136 m/decade). This finding coincided with the outcome of Dey et al. (2017), where they observed significant variations in GWL across the northwest region of Bangladesh. Salam et al. (2020a) also found that the GWL goes down before the monsoon and goes up after the monsoon in the northwestern part of Bangladesh.
According to Adham et al. (2010), approximately 8.6% of annual rainfall (1670 mm) infiltrates into the aquifer and eventually contributes to groundwater recharge. Although there has been plenty of surface water available for irrigation of agricultural crop varieties, the surface water supply has become limited during the dry season due to the changing climate and erratic rainfall patterns Pham et al. 2022). In recent decades, annual rainfall in the Joypurhat district has never exceeded 1700 mm, which is 52% less than the national average of 2550 mm. The rate of GWL in this drought-prone area has dropped faster than in other districts during the previous decade, as supported by the findings of Dey et al. (2017), Salam et al. (2020a), andDi Nunno et al. (2022). Furthermore, elevation is a factor that influences GWL, since the Joypurhat district is about 10 to 20 m above mean sea level. Hence, the possibility for groundwater recharge from rainfall and surface water sources is lower here, and the wetland area is less, resulting in a considerable declination of GWL in the study area.
Wavelet analysis is very important for assessing any region's long-term wet and dry phases. Wavelet analysis reveals two significant periods of groundwater level in Khetlal and Akkelpur upazilas, but only one significant period in the other three upazilas in Joypurhat district. Khetlal exhibits a dry state as opposed to the wet condition of groundwater. The findings are comparable to those of Salam et al. (2020a), who reported a considerable period of GWL in the Rangpur area of northern Bangladesh.
The ARIMA model, as reported by other researchers (e.g., Lee et al., 2009;Gibrilla et al., 2017;Takafuji et al. 2019;Salam et al. 2020a;Zhang et al. 2021;Islam et al., 2021a), is usually well adapted to predict GWLs for diverse geological settings globally. Bidwell (2005) and Lee et al. (2009) found the applicability of ARIMA for forecasting GWL in permeable aquifers, which is consistent with our findings. However, additional data like pumping rates from nearby wells or water levels of connected observed wells may considerably improve forecast accuracy. The main advantage of the ARIMA model is that it may utilize widely accessible information and is therefore very simple to implement on uninterruptible wells in various hydrogeological scenarios. They have a high degree of precision and accuracy. This model may be automatically modeled using AIC; also, its findings match the data's trend and seasonality patterns with acceptable precision and accuracy. In comparison to the findings of other studies that used machine learning models other than ARIMA (e.g., Elbeltagi et al. 2022;Pham et al. 2022;Di Nunno et al. 2022), the findings of this study show that the ARIMA model is remarkably well apt to do the GWL prediction for five observation wells. Furthermore, due to its simple input variable, the ARIMA is easily applied. When compared to the standard ANN technique, this model is especially effective at predicting a time-series dataset with a monthly scale (Takafuji et al. 2019).
In recent times, the drought-prone humid region of northwest Bangladesh has experienced severe scarcity of GWL due to the over-extraction of GW. Over-exploitation of aquifers to offset the balance between demand and supply and the variations in surface water flow rate is considered the major cause of GWL declination. The long-term GWL fluctuation is a function of many-fold parameters, including annual rainfall, annual groundwater withdrawal, and local subsurface geological conditions. The increase or decrease in precipitation can significantly affect GWL as it is the primary source of groundwater recharge (Park et al. 2011). Consistent with a widespread and prevailing decreasing trend in GWL over the study area, the decreasing trend of precipitation cannot be ignored as a responsible factor. Local geology essentially controls the timing and pathways of groundwater recharge to shallow aquifers in Bangladesh (WARPO 2000). The declining shallow groundwater storage in northwestern Bangladesh relates to the intensity of abstraction of GW and areas of the high thickness of surface clay where rainfall-fed recharge rates are low due to the low hydraulic conductivity of this subsurface geology (Shamsudduha et al. 2011). Hence, earlier evidence indicates that rigorous abstraction of groundwater for irrigation, which is one of the impacts of climate change and the local geology of the study region, can be attributed to rapidly declining GWL in the drought-prone urban area.
The decreasing rate of GWTs may be due to declining recharge resulting from low rainfall, an expanding plow pan triggered by enhancing conservation agriculture practices, and a thick silty clay surface in the Khetlal upazila. In this upazila, the rainfall has decreased by 12-80 mm/year, resulting in a decreased natural recharge over the past 10 years . The abuse of existing GWL, primarily being extracted for paddy farming, rice mill operations, and other industrial purposes, is an alarming issue that affects the Khetlal upazila, where groundwater provides 80% of the water needed for rice irrigation. Here, random use of water in industry and agriculture is causing GWL declination. Inefficient use of water for domestic and irrigation purposes leads to over-extraction of groundwater and also contributes, to some extent, to declining GWTs. On average, GWL varies between 1 and 6 m in the northern region, consistent with the Master Plan Organization (MPO 1987) claim, which states that groundwater is mainly available within 5 m of ground within alluvial aquifers. Evidence suggests that aquifer replenishment slowly ends during the dry season due to low rainfall in the short monsoon season and less available soil moisture Zinat et al. 2020). This is primarily because of more groundwater mining than aquifer recharge (Salam et al. 2020b). The over-exploitation of groundwater for irrigation without any rise in rainfall triggered a decline in GWLs to some degree. They cannot be replenished naturally, which has caused a groundwater deficit in this drought-prone region of Bangladesh (Hasanuzzaman et al., 2017).
Our study has some implications for sustainable groundwater development. Utilizing underutilized rivers, measuring aquifers' safe yields, recharging the aquifers by natural and artificial ways, harvesting rainwater, and using wastewater are the possible choices for this alternative Dey et al. 2017;. Suppose rainfall collected in sustainable ways during the monsoon season is used for domestic purposes after the low-cost filtering process. In that case, this will substantially lessen the water shortage issue in the northwest region in the long run. To aid this, rainwater harvesting facilities should be properly planned and designed based on our findings , with a substantial performance in the Khetlal and Kalai upazilas that need collaborative efforts among governments, non-government organizations, and hydro-geologists. Also, practices such as mulching can increase the waterholding capacity of field soils and decrease evaporation. Organic mulching can also aid groundwater recharge by increasing soil organic matter, especially in drought-prone upazilas such as Khetlal and Kalai. Rapid urbanization, climate change, and poor water management policies threaten groundwater resources in the study area, perplexing their sustainable management. So, strategies like pumping water back into aquifers, collecting rainwater, using technologies to save water, and managing water resources as a whole must be used.
Groundwater study must be coupled with other geoscience and hydro-engineering research. Multi-disciplinary tools should be managed in case of GWL depletion in a sustainable manner. Groundwater should be used as a supplement to surface water. Deep insight and protection of groundwater for future communities, particularly overexploitation, should be monitored from time to time and a strengthened legal approach should be taken (Bhattacharjee et al. 2019). Future groundwater issues cannot be overlooked if GWR management is consistent (Dey et al., 2017). Thus, sustainable GWR management depends on precisely evaluating groundwater resources' current and future trend patterns (Richter et al., 2006). This study concentrates on assessing the present and forthcoming trends of GWL in the Joypurhat district of northwest Bangladesh, which covers how decision-makers take the right initiative for combined sustainable GWR management. The increased irrigation coverage coupled with increased groundwater extraction from the aquifers over the recent past has been attributed to a decrease in GWLs. Hence, the northwest region, especially the Barind area, is most concerned about dropping GWL, which leads to a lack of access to water for drinking and irrigation in several regions. GWL declination can have long-term consequences for drinking water and, in particular, human health. For GWR management, it is important to keep an eye on the changes in GWL using a simple but useful tool. Although the findings are reasonable when using the GWL dataset for future predictions, however, because of the associated operational costs of groundwater monitoring, long time series for GWL forecasts are now very scarce in various places throughout the globe.
Some of the vital indicators of GWL fluctuations, including bore-well depths, closeness of well to sink or source, and hydrogeological features, e.g., aquifer porosity and permeability, are ignored in this study due to unavailability of data. Further investigation is essential by incorporating these critical hydrogeological parameters into data-driven tools to enhance short-and mid-term prediction accuracy and confirm that these hybrid methods agree well with the theoretical and numerical-based approaches. Such an understanding could give further numerical insights and would probably lead to a higher acceptance of soft computing methods in the sustainable management of GWR. Further study should consider all the wells for predicting regional and national levels of GWL. Future research should concentrate on finding a way to combine expert knowledge from different groundwater hydrology fields. In addition, with novel machine learning models that would enable us to effectively work and achieve satisfactory outcomes. This study was only confined to the GWL of one district in northwest Bangladesh. It has been suggested that more research be done on other districts so that more accurate forecasts can be made for other drought-prone districts in the northwest of Bangladesh.

Conclusion
In recent years, groundwater declination has become a major global and regional concern. This study represents how alarming the GWL depletion is in the study area. This research proves the magnitudes of GWL fluctuation scenarios obtained from the modified Mann-Kendall test, Sen's slope estimator, Pettitt's test, and wavelet analysis. The ARIMA modeling outcomes are realistic and consistent with the present condition. The study shows a decreasing trend of GWL, and the declination rate of GWL in Joypurhat district is 0.136 m/year from 1980 to 2019. Out of five upazilas, the groundwater declination rates of Khetlal and Kalai are higher than the district average, which is liable to excessive groundwater abstraction in high waterconsuming crop production like rice and potato. The predicted results revealed that the GWL will decrease at a rate of 0.143 m/year between 2020 and 2050, which is very concerning for the future sustainability of groundwater resources. Proper crop water management should be emphasized to avoid such adverse impacts in this study area. Though rainwater harvesting is recommended as a long-term remedy to cope with drought effects, simultaneous monitoring can be a short-term remedy. Future research should concentrate on GWL fluctuation as it is tele-connected with large-scale atmospheric oscillation in the study region. Our study provides sustainable groundwater potential under current and future climate change scenarios in the northwest region of Bangladesh. The findings of this study can help policy-makers make sound decisions about the future installation of irrigation equipment and the sustainable use of groundwater resources in this area and other parts of the world.