Analysis of groundwater table variability and trend using ordinary kriging: the case study of Sylhet, Bangladesh

Many developing countries experience widespread groundwater declination. Sustainable management actions include generation of an accurate groundwater distribution based on an extensive groundwater monitoring network which is often cost prohibiting in the context of a developing country such as Bangladesh. Further, such knowledge is lacking for the Sylhet region where groundwater was documented to be under tremendous pressure. Specifically, the gap in the current literature exists regarding groundwater trends and its areal extent for this region. This paper bridges the gap in research by focusing on trends and spatial and temporal variation of groundwater level changes for this area. This study addresses this problem by creating groundwater level predictions at the ungauged areas using geostatistical methods applied to a detailed set of data. In this study, the spatial variability of annual-average depth to the water table at 46 observation wells in the Sylhet division in Bangladesh is analyzed for 2000, 2005, 2010, and 2015. The geostatistical analysis applies the ordinary kriging method with cross-validation to create the water table maps for the study area. The results indicate a substantial increase in groundwater depths during the studied period from 2000 to 2015 in some locations in the study area. Importantly, this work identifies the vulnerable zones in the area due to the groundwater lowering trend. The study adds to the groundwater management research in developing countries and focuses on the spatial and temporal groundwater variation. The findings from the modeling exercise contribute to identification of the vulnerable areas and therefore help policymakers in making informed decisions to manage groundwater resources in this sensitive region sustainably.


Introduction
The groundwater system, which is very dynamic by nature, constitutes about 98% percent of the Earth's freshwater (Zaporozec 2004). The increasing demand due to population growth and rapid industrialization has rendered it vulnerable which especially concerning in developing countries where resources for sustainable groundwater management are scarce. Over-abstraction, poor irrigation management, and reduced groundwater recharge are the primary causes of groundwater lowering (Bellingeri et al. 2017). The decline of groundwater also has significant environmental impacts that include groundwater quality degradation (Foster et al. 2018). Consequently, the combined effect of groundwater depletion along with water quality degradation results in the uncertainty of future water availability (Tabari et al. 2012). Other environmental effects include the increased risk of land subsidence (Cui et al. 2018), reduction in the amount of water in streams and lakes (Perkins et al. 2017), and increased production costs due to greater energy needs for pumping (Turner et al. 2019).
Groundwater declination is a common and widespread problem in developing countries; Bangladesh is no exception. Irrigation and household water supply in the country primarily depend on groundwater. It is declining at a rate of 0.1 to 0.5 cubic-meter per year (m 3 /yr) in the country (Dey et al. 2017). Consequently, the country's groundwater resource needs proper management actions for sustainability.

Page 3 of 12 120
For effective groundwater management, the knowledge of the spatial and temporal variation of groundwater conditions over the concerned area is essential (Varouchakis et al. 2012). However, generating an accurate groundwater distribution over an area needs an extensive groundwater monitoring network. In reality, the monitoring networks are not always uniformly distributed due to inaccessibility, high installation and maintenance costs. Therefore, groundwater level predictions at the ungauged areas are warranted while making informed management decisions.
In hydrogeology, the environmental quality assessment analysis is commonly conducted to evaluate the total pollution level in a water body including Nemerow pollution index and risk assessment index (Mishra et al. 2016), and a scoring methodology is used to assess the performance of freshwater conservation plan (Singh et al. 2021). Groundwater research commonly relies on geostatistics as it produces the accurate estimates (Ma et al. 1999). For example, geostatistical interpolation methods are extensively used for generating data for the ungauged areas (Theodossiou and Latinopoulos 2006), and to derive the groundwater trends over the long term (Reghunath et al. 2005). Similarly, Kumar (2007) used geostatistics to interpolate the groundwater levels by applying universal kriging for optimal contour mapping of groundwater levels. Besides, geostatistical techniques are used in hydrogeochemistry to identify the contaminant sources and control or mitigate them (Elumalai et al. 2017). It is also applied in agricultural management to assess the spatiotemporal variability of groundwater quality parameters for irrigation (Yazdanpanah 2016).
Apart from groundwater quality evaluation, geostatistics is used for assessing aquifer vulnerability (Machiwal et al. 2018). Machiwal et al. (2012) used kriging to analyze the groundwater level data from 50 sites in Western India to find the spatial autocorrelation and variances in groundwater levels. The study found the kriging advantageous to spot critical locations where groundwater resource management strategies are required. Xiao et al. (2016) used seven discrete interpolation methods to interpolate the groundwater level in Beijing, China, and found that the simple kriging provides the best fit for the model. Furthermore, Nikroo et al. (2010) applied different kriging methods to identify the best geostatistical interpolation technique and predicted groundwater depth and elevation in Southwest Iran.
Understanding the groundwater variability and trend in the Sylhet division of Bangladesh is the primary motivation of this work. Although the area is vulnerable to groundwater lowering, understanding its trend and areal extent has mostly been ignored in the previous works (Zafor, 2017). To fill the gap, we explore the spatial and temporal groundwater variation using the collected data from 46 observation wells in the Sylhet division for 2000, 2005, 2010, and 2015. The specific aim of the study is creation of groundwater level predictions at the ungauged areas using geostatistical methods and identification of vulnerable areas in the Sylhet region of Bangladesh. The results of this study will help authorities to make informed decisions for efficient management of the groundwater resource in the study area.

Study area and data collection
The study area is the Sylhet division (24°30′N, 91°40′E), which is located at the northeast region of Bangladesh ( Fig. 1) and also addresses as the tea capital of Bangladesh (Nury et al. 2017). It is also known as the tea capital of Bangladesh (Nury et al. 2017). It borders the Chittagong and Dhaka divisions to the southwest and west, respectively, as well as Meghalaya, Assam, and Tripura states of India to the north, east, and south, respectively. The study area is vulnerable to groundwater lowering as it is the primary source for domestic and irrigation water supplies (Zafor et al. 2017). This study investigates the depth to the water table throughout the study area.
The study area, the Sylhet division, is subdivided into four major districts, namely Habiganj, Moulvibazar, Sunamganj, and Sylhet. These four districts are subdivided into 36 subdistricts called Upazila (Nury et al. 2017;Zafor et al. 2017). Its climate is humid subtropical, with a predominantly hot and humid summer and a relatively cold winter. The average annual highest and lowest temperature recorded as 23 °C (Aug-Oct) and 7 °C (Jan), respectively (Zafor et al. 2017).
The aquifer of the study area is semi-confined to confined in nature. It is made of weathered alluvial sands of the Dupi Tila formation (Zafor et al. 2017). Its thickness varies from 20 to 98 m in some areas (Ahmed et al. 2019). Groundwater is the primary drinking water source in the area because the study area does not have a central water distribution network. Most of the houses have deep tube wells for groundwater extraction. For this reason, a continuous monitoring of groundwater depth and quality is necessary.
The Bangladesh Water Development Board (BWDB) (www. bwdb. gov. bd) is one of the main responsible authorities to monitor the groundwater throughout the year. The organization maintains a distributed network of monitoring wells across the country. The BWDB staffs take water table reading weekly from the monitoring wells and send the data

Methodology
Kriging is the best linear optimum unbiased interpolation method of estimating unknown values of spatial and temporal variables with a minimum mean interpolation error (Chung et al. 2019). Generally, several types of kriging methods are available: ordinary, simple, universal, Poisson probability, and more (Böhner & Bechtel 2017). The selection of the method depends on the characteristics of the available data. For instance, universal kriging is appropriate for nonstationary data, and cokriging suits better for a group of correlated data (Chung et al. 2019). The ordinary kriging (OK) is a standardized type of kriging and known as the best linear unbiased estimator (BLUE) (Cressie 2015). This method assumes that the data sets are stationary and exploits semivariogram, which approximates the values without bias and least variance, to determine the best fit for the spatial relationship model (Cressie 2015). Considering the precedence of ordinary kriging and the characteristics of the available data, this study adopted the ordinary kriging method.
After the acquisition, the water table data (i.e., depth to the water table measured in meters) were checked for normality by observing histograms and normal Quantile-Quantile (QQ) plots. To fit the spatial map and to get a better prediction result in the map, the data should exhibit a normal distribution either in a histogram or a normal QQ plot. If the data show a skewed pattern revealed in the histogram, then the data need to be transformed.
The histogram of the mean annual groundwater depth data of 2000 ( Figure S-2) is skewed right, meaning that a small number of samples with high groundwater depths lie on the right tail. Although 2000 data appear to be positively skewed, the skewness is not severe. This is also supported by the similarity of mean and median values, and most of the points falling on the straight line seen in the Normal QQ plot ( Figure S-3). Therefore, no transformation is needed for the 2000 data.
However, the skewed pattern is observed for the water  table data Table 1 shows the required transformations for the mean annual water table data of 2000, 2005, 2010, and 2015. Following the data transformation, global trends of groundwater (i.e., a rising or deepening trend of water table along with the spatial directions) were evaluated. The semivariogram (h) is the primary tool in geostatistics that expresses the spatial dependence between the neighboring observations (Ahmadi and Sedghamiz 2008). It is defined as one-half of the variance of the difference between the attribute values at all points separated by h as given in Eq. 1. where z(x i ) and z(x i + h) are the magnitudes of the variable (e.g., depth to the water table measured from the ground surface in meter as shown in Fig. 2) at the point x i and a point of distance h from the point x i . N(h) is the total number of feature pairs (i.e., pair of wells) separated by the distance of h.
For examining and quantifying spatial autocorrelation, empirical semivariogram models were developed for 2000, 2005, 2010, and 2015 using Eq. 2. Each model was fitted to the points which were used for its development. The semivariogram modeling is identical to fitting a least-squares line in the regression analysis. (1) where value at location i and j is a depth to the water table in meters.
Following the model fitting, the differences between estimated and observed values were summarized using the cross-validation statistics described by Ahmed and De Marsily (1987). The crossvalidation statistics include the estimates of root-mean-square error (RMSE) and average standard error (SE) (Chung et al. 2019;Taye et al. 2018). A smaller RMSE and SE indicate a closer observed and predicted groundwater depths by the empirical semivariogram models. Besides, the value of RMSE close to zero suggests an unbiased prediction. Finally, the residual for the prediction map were checked for normality using a QQ plot.

Result and discussion
The trend analysis identifies the presence of trends, if any, existing in the data and the order of polynomial that fits the data. Figure Figure S-4 shows that the best fit lines using data points projected on the perpendicular planes (North-South or East-West) are curvilinear. Therefore, a second-order polynomial (quadratic) can capture the groundwater trend for 2000, 2005, 2010, and 2015. These trends need to be removed beforehand to capture the short-range local variation. After modeling the residuals, the global trends are added to the analysis to obtain a better prediction surface. Before semivariogram modeling, the spatial autocorrelation and a directional influence among the collected water table data were checked. Four pairs of locations are selected in the semivariogram within the close range of 0 to 450 m distance indicating different adjacent points (Fig. 2). Besides, the lines linking the locations fall in the area within different colors. Consequently, the spatial autocorrelation exists in the data for 2000. Similarly, results show that spatial autocorrelation exists for 2005, 2010, and 2015. The semivariogram modeling was implemented to determine the best fit for the spatial relationship models. Some values were grouped according to the separation distance (bin) to obtain the semivariogram curve. Figure 3 illustrates the directional influence on water table data for 2000. It shows pairs of the semivariogram value within an approximately same distance. Since most of the lines are scattered, there is no directional influence exists in the 2000 data. Similarly, no directional influence is detected for 2005, 2010, and 2015 data. Finally, the values of the parameters, including nugget, range, sill, and shape, were obtained. To illustrate, Fig. 4 shows the obtained parameters for 2000.
Cross-validation checks for the predictive models were performed to identify the model's accuracy before generating the prediction maps. As stated earlier, both the RMSE and average SE should be as small as possible. The results of the checks for the water table data are shown in Table 2. Besides, Fig. 5 shows the QQ plot for the residuals of the prediction model for 2000, 2005, 2010, and 2015. The results suggest that the residuals are normally distributed, with most of the points falling on the straight line of the QQ plot.
The study area is divided into six zones to analyze the groundwater variability-east, west, south, central, Depth to water table (m) Year North-East northeast, and northwest. Figure 6 shows the groundwater trend corresponding to each zone. Specifically, the temporal plots show the mean annual-average depths to the water tables in the six zones (Fig. 6). In addition, Table 3 lists all the Upazilas inside the study area with the corresponding zone and provides the mean annual depth to water table along with the temporal percent change for each Upazila. The major findings of the trend analysis are as follows:  Groundwater showed a declining trend in some of the past study areas (Zafor et al. 2017). The present study has also found a water table declining trend in the northeast zone, which is in consistent with Zafor et al. (2017). In contrast, for the other five zones, results show an overall rising trend of groundwater depth in recent years (2010)(2011)(2012)(2013)(2014)(2015). In 2015, Upazilas like Kanaighat, Sylhet Sadar, Derai, Sunamgonj, Hobigonj were flooded above their respective danger level from mid-June to mid-September, reported by the Bangladesh water development board (BWDB, 2015). Also, 32% percent of the total area of Bangladesh affected by flood in the year 2015, whereas total percentage was 18% in the year of 2010 (BWDB, 2015). This might be the reason of rising groundwater trend in the year of 2015 compared to 2010. However, local declination is also observed in some Upazilas like Madhabpur in west, Rajnagar in east, Bahubal, Chunarughat, and Sreemangal in south zones as shown in Fig. 7 and Table 3. The change in rainfall and local irrigation practices could be the primary reason behind the groundwater table fluctuation which is beyond the scope of this current study and requires further investigation.

Conclusion
The study compared the groundwater spatial distribution map of 2000,2005,2010, and 2015 to assess the groundwater trend in the study area. It was found that the water table rose to a shallower depth in 62.5 percent of the 32 Upazilas studied under the study during 2000-2015. In contrast, water table deepened in 37.5 percent Upazilas. Upazilas like Jaintapur, Gowainghat, Maulvibazar Sadar, Kanaighat, and Bahubal are the susceptible zones due to groundwater lowering. Indiscriminate groundwater withdrawals and reduced recharge due to industrialization and urbanization could be the probable causes of groundwater lowering in the study area. A relatively small sample size was the primary limitation of the study. Additional piezometric wells with higher resolution data for a long duration might generate a better prediction map and capture the trend with seasonal variation. Therefore, the future works should focus on the application of various geostatistical techniques and inclusion of seasonality into the data analysis.
Funding This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.
Data availability Data can be provided upon request.

Conflict of interest
There is no potential conflict of interest reported by the authors.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.