Methods for Characterizing Groundwater Resources with Sparse In-Situ Data

Methods for Characterizing Groundwater Resources With Sparse In-Situ Data Ren Nishimura Department of Civil and Construction Engineering, BYU Master of Science Groundwater water resources must be accurately characterized in order to be managed sustainably. Due to the cost to install monitoring wells and challenges in collecting and managing in-situ data, groundwater data is sparse in space and time especially in developing countries. In this study we analyzed long-term groundwater storage changes with limited timesseries data where each well had only one groundwater measurement in time. We developed methods to synthetically create times-series groundwater table elevation (WTE) by clustering wells with uniform grid and k-means-constrained clustering and creating pseudo wells. Pseudo wells with the WTE values from the cluster-member wells were temporally and spatially interpolated to analyze groundwater changes. We used the methods for the Beryl-Enterprise aquifer in Utah where other researchers quantified the groundwater storage depletion rate in the past, and the methods yielded a similar storage depletion rate. The method was then applied to the southern region in Niger and the result showed a ground water storage change that partially matched with the trend calculated by the GRACE data. With a limited data set that regressions or machine learning did not work, our method captured the groundwater storage trend correctly and can be used for the area where in-situ data is highly limited in time and space.


INTRODUCTION
Groundwater is a precious resource and is increasingly in demand. Groundwater is used for multiple purposes including but not limited to drinking water, agriculture, mining, or industry [1][2][3]. Approximately 30% of the earth's freshwater exists as groundwater while only 1.2% is surface water such as lakes, rivers, and streams [4]. Because of climate changes and severe droughts, water managers are shifting from surface water to groundwater [5]. The average groundwater pumping rate for the entire world has gradually increased and is estimated to be 600 -1000 km3 per year and expected to increase even more due to population growth [6].
Understanding long-term groundwater level changes is essential for sustainable groundwater management; however, groundwater is often pumped out of aquifers without consideration of sustainability. In the past half century, the groundwater storage in various places has been depleted [7]. Theis [8] explained that groundwater extraction is initially derived from removal of water in storage, but in the long-term, the source of groundwater extraction shifts from the storage to decrease in discharge and/or increase in recharge. No additional water is removed from storage when a new equilibrium is reached. When groundwater extraction exceeds natural recharge, the system cannot reach a new equilibrium resulting in perpetual groundwater mining [9]. For example, Central Valley, California experienced 80 km 3 (65,000,000 ac-ft) of groundwater storage depletion due to lack of sustainable groundwater management and over pumping since the 1960's [10]. The Beryl-Enterprise aquifer in Escalante Valley, Utah was also studied and it was estimated that the net groundwater depletion rate is approximately 65,000 acft per year, resulting in substantial water level declines [11].
Understanding groundwater is more important than ever due to increasing pressure for sustainable water resource management. In recent years, there have been many regulatory requirements initiated to protect the groundwater resources in the United States. Several states require water agencies to characterize groundwater resources and implement management plans for sustainable groundwater development and management. For example, California issued the Sustainable Groundwater Management Act (SGMA) to help protect groundwater resources over the long term [12,13]. The SGMA requires the formation of groundwater sustainability agencies (GSAs) for the high and medium priority basins and for the GSAs to develop and implement groundwater sustainability plans to avoid over usage of groundwater resources and to mitigate the over-usage [14]. In Utah, groundwater management plans were created to promote wise use of the groundwater, protect existing water rights, and address water quality issues and overpumping of groundwater.
Because of higher demand, increased pumping, and the regulatory requirements of groundwater, understanding and characterizing long-term groundwater storage changes is critical to sustainably manage aquifers, but it is technically challenging. First, groundwater is not directly visible and it has to be physically measured at discrete locations to determine water table elevations. Groundwater generally requires wells or piezometers to measure water table elevations but they need to be spatially distributed throughout aquifers. Not only is it expensive to physically sample wells on a consistent basis but it also requires manpower to log long periods of records to accurately characterize the condition of aquifers [15,16]. Because of these difficulties, groundwater level data is spatially or temporally sporadic, especially in developing countries [17].

Previous Work
Several researchers have developed techniques to address the issue of sparse data.
Researchers have used in-situ data driven models to predict groundwater levels using generic programming [18], principal component regression [19], and artificial neural networks [20].
Other researchers have incorporated remote sensing data to address the temporal data scarcity.
Evans et al. [21] used soil moisture data from the Global Land Data Assimilation System (GLDAS) developed by the National Oceanic and Atmospheric Administration along with an Extreme Learning Machine algorithm to impute missing groundwater measurement at a given well. Seo and Lee [22] used multiple satellite data together with long short term memory (LSTM) and convolutional neural network LSTM for groundwater data imputation.
Spatial data scarcity is often addressed by spatial interpolation techniques such as Kriging [23]. Gundogdu and Guney [24] analyzed multiple semivariogram models of kriging for the groundwater table estimation and comparison in Mustafakemalpasa in Turkey. Ruybal et al. suggested that three-dimensional interpolation with space and time, called spatiotemporal interpolation, yields more realistic temporal and spatial changes in water levels relative to spatial interpolation [25]. These methods are useful but require a significant amount of in-situ groundwater measurements to temporally impute time-series data or spatially interpolate the groundwater elevations at a given time. Challenges remain for developing countries where available groundwater measurements are sparser than what is necessary to perform the temporal and spatial interpolation methods described above.

Study Location and Background
Niger is located between latitude 11 to 24 degree North and longitude 0 to 16 degree East as shown in Figure 1-1. It is a landlocked country in West Africa and covers an area of 1,270,000 km2, which makes it the largest country in West Africa, and 80% of its area is part of the Sahara desert [26]. Niger has distinct wet and dry periods where the northern area has longer dry periods. The wet period starts around June to October for five months with the average precipitation of 700-1100 mm in the southern part and July to August in the northern desert area of Niger [27]. The Niger River is the only permanent river in the country whereas the other rivers are ephemeral, called "wadi" in Niger. There are also internal drainage ponds and dammed reservoirs for surface water storage [28].

Figure 1-1: The location of Niger and its adjacent countries in the African continent
Niger experienced a severe drought in the early 1980's together with desertification and needed alternative reliable water sources [29]. Aquifers in Niger were intensively studied starting around 1980 by various donors and groundwater wells were developed [30]. There are studies in the past that found increases in groundwater storage even though the area experienced severe drought. Leduc et al. [31] concluded that the increases were due to changes of land use and increase in the groundwater recharge patterns. Although groundwater has been periodically studied in Niger, a challenge still exists to sustainably manage the groundwater resources over the long term because the legal and organizational groundwater management frameworks are fragmented, inconsistent, and sometimes incomplete in the region [32,33]. As a result, there are very few historical groundwater measurements to accurately characterize the long-term changes.
The SERVIR program operated by the National Aeronautics and Space Administration (NASA) assists developing countries addressing sustainability and environmental issues incorporating Earth Observations to assess, analyze, and manage natural resources, including groundwater, to improve human lives [34]. This research was funded by NASA SERVIR and we worked with the West Africa hub that was one of the five NASA SERVIR regional hubs. The West Africa hub is partnered with an organization called AGRHYMET located in Niamey, Niger and who had provided the data for this research. The objective of our research was to assist water managers and local stakeholders to analyze and characterize local groundwater conditions.
We use two approaches to assess groundwater storage changes in Niger. The GRACE Groundwater Subsetting Tool (GGST) calculates groundwater storage anomalies for a userdefined area using data from the NASA Gravity Recovery and Climate Experiment (GRACE) missions [35]. We derive groundwater storage anomalies by subtracting terrestrial water storage anomalies simulated by GLDAS from the total water storage anomalies provided by the GRACE mission. This method does not require in-situ data and can be used to analyze the changes of groundwater storage in data-poor areas from 2002 to the present time. Barbosa et al. [36] studied the groundwater storage changes in two regions in Niger using the GGST and identified a significant increase in groundwater storage volume from 2011 to present. While this approach is relatively simple and doesn't require in-situ data, the native resolution of the GRACE grids is 3 degrees x 3 degrees so it can only be applied to large regions.
The Groundwater Data Mapping Tool (GWDM) displays historical groundwater level measurements at monitoring wells on an interactive map [37]. It can also be used to analyze groundwater storage changes through a three-step process. First, a temporal interpolation is performed to sample the water level time series for each well at selected intervals (2 years, 5 years, etc.). Prior to this step, the GWDM uses a machine learning algorithm to impute gaps in the water level times series using correlations with Earth observations. In the second step, the water levels from the wells are interpolated spatially at each time interval to create a raster of the water level that is clipped to the aquifer boundary. Finally, the volume between each subsequent pair of rasters is computed and multiplied by a storage coefficient to derive a curve of groundwater storage change vs time. Unlike the GRACE method, this approach can be applied to aquifers of any size as long as sufficient water level data is available.
To assess groundwater resources in Niger using the GWDM approach, we obtained groundwater data from the AGRHYMET Regional Centre in Niamey, Niger [38].
AGRHYMENT is a research center associated with the Permanent Interstate Committee for Drought Control in the Sahel (CILSS), and organization dedicated to achieving food security and increased agricultural production in the region. Unfortunately, while there are some wells in Niger with historical measurements, they are extremely rare. Most wells are not sampled after they are constructed. We were able to obtain data for three aquifers but only for a small set of wells and the measurements only covered a short period of time. However, we did obtain the well inventory containing 15,688 wells scattered throughout southern Niger, and the water level at the time of construction was measured for many of these wells. The GWDM method requires a time series of measurements at each well, so we can't directly apply this method, but the data set does provide information on how water levels have changed over time.

Research Objective
The objective of this study is finding a way to harness these single-value well records and calculate groundwater storage change over selected regions. were then used to analyze the groundwater storage changes in selected regions in Southern Niger.

DATA
The data used in this research was provided by AGRHYMET and included a well file that contained 15,688 rows of well information with dates ranging from 1947 to 2010, well locations in degree latitude and longitude, depth to water table in meters, and several other columns. We filtered out the wells missing critical information, such as date or location, and 4,835 wells out of the original 15,688 wells remained for further analysis.
The well records had three date-related columns; the date that the well construction started, the date the well construction finished, and the date the well was rehabilitated. For a given well, one or more of these dates could be empty. For wells with multiple dates, we had to determine which date to associate with the corresponding depth to water table value. In consultation with an AGRHYMET representative, priority was given first to the date rehabilitated, then to the date construction finished, and finally to the date construction started.
The well locations and dates are shown in   To obtain well data for the Beryl-Enterprise Aquifer, we used the USGS National Water Information System (NWIS) web service. NWIS web service provides online public access to water-resources data collected in the United States. Similar to Niger's data set, groundwater measurements were reported as depth to groundwater. We also calculated the WTE by the same method. Unlike the Niger data set, most of the wells in the Beryl-Enterprise Aquifer have historical groundwater measurements. The locations of the wells and a histogram of water level sampling frequency for the aquifer is shown Figure 2-2 To replicate the temporal data scarcity characteristic of Niger, a single groundwater measurement was randomly sampled (n=1) from each well in the Beryl-Enterprise Aquifer, and we repeated this 500 times to create 500 different random sample data sets. We used each of these data sets to calculate the annual storage depletion rate using our proposed algorithm and then compared with the known depletion rate in order to validate our methods. The date range of a random sample varies among the 500 data sets depending on which date is picked by random sampling. The spatial distribution of the wells does not change among the 500 data sets, and each well only has one groundwater measurement. One random sample for each well (n=1) is an ideal data set to replicate the temporal data scarcity in Niger. However, in order to evaluate the impact of data frequency on the spread or uncertainty of the method results, we created three additional data sets where we randomly selected more than one sample per well. We decided to use 5, 15, and 30 random sample sizes (n=5, 15,30) and created 500 additional random sample data sets per sample size to analyze the variations of resulting estimated depletion rates.

METHODS
Our proposed method for dealing with data sets with one historical WTE per well is to group wells from a common geographic location represented by a "pseudo-well" at the centroid of the location. The WTE values for each of the wells associated with the pseudo-well are from a variety of dates, so we aggregate them into a single synthetic WTE time series as shown in

K-Mean-Constrained Clustering
To overcome the issue of uneven numbers of the individual wells contributing to the pseudo-wells, we used a modified K-means clustering algorithm to set the minimum number of the individual wells to the pseudo-wells. A K-means algorithm iteratively partitions the data set into a user-defined number of distinct non-overlapping groups (K) where each data point belongs to only one group. The algorithm iteratively tries to minimize the sum of the squared distance between the data points and the cluster's centroid until the assignment of the data points to the clusters no longer changes [41]. This algorithm may still result in an uneven number of data points in each cluster and therefore does not fully resolve our issues of the grid clustering method. To ensure a more uniform distribution of wells to clusters, we used a modified version of the K-means algorithm from a Python library called k-means-constrained which is a modified version of scikit-learn's KMeans method whereby a minimum and/or maximum cluster member size can be specified [42].
While a K-means clustering takes only a user input of K, the modified K-means algorithm also takes the minimum size for each cluster and adjusts the number of clusters accordingly. As a result, each pseudo-well at least has the specified number of individual wells.
After applying the k-means-constrained method we placed pseudo wells at the centroids of the clusters and aggregated the WTE measurements from individual wells in each cluster to form WTE time series as shown in Figure 3  We tested a minimum cluster size of 13 and 25 wells for the Beryl-Enterprise aquifer and 28 and 14 clusters/pseudo-wells were created as shown in Figure 3-4. One advantage of this method is because each pseudo-well has the minimum number of individual wells, the resulting synthetic WTE time series are more likely to exhibit a similar temporal data distribution.
However, unlike the grid clustering method, the spatial size (area) of each cluster is different in this method. In regions where wells are sparser, the cluster must cover more area in order to include the minimum number of wells. On the other hand, in the area where the wells are densely located, this method creates clusters that are spatially smaller as shown in the middle and southern region of the Beryl-Enterprise aquifer in Figure 3   Prior to generating the synthetic time series for imputation, we identified outlier WTE values at each pseudo-well. Because the WTE values for each pseudo-well are derived from individual wells from various spatial locations and due to data errors, the pseudo-wells often had WTE values that seemed to be outliers. Thus, we removed WTE values at each pseudo-well that were outside of the 3rd standard deviation.
After removing outliers, we then fit a univariate spline curve to the WTE values at each pseudo-well using the Scipy univariate spline, which is a 1D smoothing spline [43]. This method fits a spline equation of degree k (k=1-5) to the provided x and y data by specifying a smoothing condition. A smoothing factor s is used to choose the number of knots, and this number is adjusted until the smoothing condition equation ( 3-1 ) is satisfactory.
We used the univariate spline because it generates a smooth curve that approximates the trend of the data points. The knot parameters were adjusted to 10,000,000 after visually inspecting the spline fit so that the fitted curves did not create an extreme overshoot and an undershoot. The fitted curve values were then resampled to at monthly intervals. Examples of the univariate spline curves are shown in Figure 3-6.
For time series that have a long gap between clusters of data points as shown in

Spatial Interpolation with Kriging
After the temporal interpolation step, we have WTE data at the pseudo-well locations scattered throughout an aquifer available for every month. We can now use these point data to spatially interpolate WTE levels throughout the entire aquifer at any or all of the time step using kriging. To perform the kriging we used the used the GSTools Python package [45] and used the  The kriging interpolation was then carried out at a selected spatial resolution throughout the aquifer. We set the spatial resolution to be 0.1 degree for both latitude and longitude for both the Utah and Niger regions. The interpolation results were placed in a time-varying NetCDF raster grid [46].

Calculation of Aquifer Storage Change
Once the NetCDF rasters of water levels were generated, the final step is to calculate groundwater storage changes vs time. This is accomplished by performing mathematical operations on series of raster data set layers of WTE at specific times produced during the spatial interpolation. The first data set, which is corresponding to the earliest time step in the series of our interest, is known as 0 and serves as the baseline from which changes in aquifer For the Niger data set that is in metric units, the mean radius of the Earth is set 6,371,000 meters [47]. The aquifer storage change is calculated in cubic kilometers (km 3 ), and also WTE and the drawdown i were calculated in meters. For the Utah data set, the changes in aquifer storage volume were calculated in acre-feet, which is the volumetric unit used in prior studies of this aquifer by the other researchers. We use a mean Radius of the Earth of 3,959 miles. In this case, the result of Equation  is in square mile-ft, which is converted to acre-ft by multiplying by 640.

Utah
The objective of the Utah data set was to use a site with a known storage depletion curve to determine if the pseudo well technique can reasonably predict storage changes in regions with sparse data. To simulate a data-poor environment, we created a sample data set by randomly selecting a single WTE value from each well. We applied the three-step method: 1. The average storage depletion rate over the designated time range was calculated and compared to the rate found by the other two studies. Table 4   were then compared to the Groundwater Storage Anomaly (GWSa) found using GRACE data via the GGST application.

Storage Change Analysis with Grid Clustering
The uniform grid clustering with 0.05-degree square grids created sparse spatial distribution of the pseudo-wells as displayed in Figure 4-2 (a). This is mainly due to the limited spatial coverage of each grid that the grid did not have more than five unique monthly time-     Although the groundwater storage trends from our methods did not match with the trend derived from GRACE for the range in time that the data sets overlapped, there are a number of possible explanations for this discrepancy. One issue could be the well-documented "leakage" effect often encountered when applying GRACE data to relatively small regions such as the one used in this study. The native GRACE grid cells are 3 x 3 degrees and the GRACE cells that overlap the region being analyzed extend over into adjacent regions, skewing the results. Our study regions may have experienced a decrease in groundwater storage as shown by our methods, but neighboring regions may have experienced gradual increase in the storage over time. Another possible explanation is that many of the wells used in this study were shallow while the GRACE results measure storage change at a deeper level. Finally, another source of uncertainty is quality of Niger's well inventory, especially in the date columns. We determined that the date when the well was rehabilitated corresponded to the groundwater elevation measurement if there are multiple date values in the well inventory. If the wells were not resampled when rehabilitated and if the water table was going up, a low WTE value from an earlier date when the well was constructed would have been artificially shifted to the later date, potentially contributing to a downward trend of the groundwater as opposed to the actual upward trend.

CONCLUSIONS
In this research, we developed a method to characterize groundwater resources with extremely sparse in-situ data, in which there was only one groundwater elevation measurement per well. The methods were shown to yield a reasonably accurate estimate but also carry a high level of uncertainty.
Various data imputation methods have developed for temporally sparse groundwater table measurements; however, most of the techniques still require a handful of in-situ measurements to apply regressions or machine learning techniques. Using GRACE data, water managers can analyze the groundwater storage changes without in-situ data; however, the spatial resolution is limited and only large aquifers can be analyzed; furthermore, the data is only available from 2002. With our method, water managers can use well data which may containing as little as one historical groundwater measurement and analyze long-term groundwater storage change in small aquifers before and after 2002. Our method can be an additional tool for water managers to characterize the groundwater storage change over time where in-situ data is highly limited as is frequently the case in developing countries.