Minimizing errors in the prediction of water levels using kriging technique in residuals of the groundwater model (MODFLOW)

Groundwater monitoring and water level predictions have been a challenging issue due to the complexity of groundwater movement. Simplified numerical simulation models have been used to represent the groundwater system; these models however only provide the conservative approximation of the system and may not always capture the local variations. Several other efforts such as coupling groundwater models with hydrological models and using geostatistical methods are being practiced to accurately predict the groundwater levels. In this study, we present a novel application of geostatistical tool on residuals of groundwater model. Kriging method was applied on the residuals of the numerical model (MODFLOW) generated by the TWDB (Texas Water Development Board) for the Edwards-Trinity (Plateau) aquifer. The study was done for the years 1995 through 2000 where 90% of the observation data was used for model simulation followed by cross-validation with the remaining 10% of the observations. The kriging method reduced the average absolute error of approximately 31 m (for MODFLOW simulation) to less than 5 m. Also, the residuals’ average standard error reduced from 9.7 to 4.7. This implies that the mean value of residuals over entire period can be a good estimation for each year separately. The use of kriging technique thus can provide with improved monitoring of groundwater levels resulting in more accurate potentiometric surface maps.


Introduction
The exponential growth of population, rapid socio-economic development, increasing food demand, and changing climatic factors have led to decline in both the quality and quantity of freshwater resources. The decreasing available resources have posed serious challenges in the agricultural sector with limited water available for irrigation. As an alternative resource, the reuse of wastewater in irrigation has been increasingly recognized as an essential, and economical strategy (Adhikari & Fedler, 2020a, 2020bMehan, 2019;). However, only a small fraction of wastewater with less than 6% in the US (FAO, 2008) and less than 3% globally (FAO, 2016) is reclaimed; and the irrigation largely relies on groundwater sources. Groundwater is one of the primary sustainable water resources, especially during the high demand seasons, due to its lower susceptibility towards sudden changes. About 70% of groundwater withdrawal worldwide is used for agriculture while irrigating nearly 38% of irrigated lands (NGWA, 2020). Likewise, nearly 50% of irrigated lands in the United States are based on groundwater sources (NGWA, 2020).
Over the past centuries, extreme drought events have significantly affected both surface and groundwater resources (Fendekova et al., 2011;Lange et al., 2017). While low surface water level might be an immediate indicator of drought, changes to groundwater levels indicate long-term water scarcity. Further, it is straightforward to monitor and assess surface water changes while measuring variations in groundwater resources is very challenging and time-consuming. Also, surface water resources can be replenished by adequate precipitation, groundwater sources might be nonreplenishable once exploited or may take several years. While these conditions make groundwater resources extremely valuable, groundwater monitoring and measurements have been inconsistent geographically and over time (Ziolkowska & Reyes, 2017).
The over-exploitation from aquifers to address the irrigation needs along with frequent drought events have caused a severe drop in the groundwater table in many areas (Ahmadi & Sedghamiz, 2007). Serious concerns have already been raised on the knowledge and understanding of actual levels of groundwater resources in developing possible mitigation/adaptation measures to the current and anticipated droughts in the future. Groundwater monitoring plays a significant role in water resources management and is fundamental to understanding system dynamics, trends in storage, and the long-term sustainability of an aquifer (Reghunath et al., 2005;Taylor & Alley, 2002). Additionally, knowledge of groundwater levels provides crucial information for the sustainable use of groundwater resources and efficient and effective conjunctive management of surface and groundwater resources (Oikonomou et al., 2018).
Simplified numerical simulation models (Chakraborty et al., 2020;Hutchison et al., 2011;Xue et al., 2018) have been widely used in the last decades to represent the groundwater system; these models provide the conservative approximation of the system, however, may not always capture the local variations. These modeling techniques use mathematical equations to describe the physical model and thus the model accuracy depends on how precisely the conceptual models describe the real-world system. These models largely rely on the available water-level data ). However, groundwater-level data are often irregularly sampled, leading to temporal gaps in the record, and are not adequately distributed spatially across an aquifer (Marchant & Bloomfield, 2018;Oikonomou et al., 2018;Varouchakis & Hristopulos, 2013). The spatial sparseness of data presents challenges when spatially interpolating potentiometric surfaces and creating groundwater maps (Marchant & Bloomfield, 2018;Ruybal et al., 2019;Ziolkowska & Reyes, 2017).
Several research efforts, including coupling of groundwater models with hydrological models and use of geostatistical tools, have been done to improve our understanding of groundwater complexity. Such efforts have also been used in addressing the issues dealing with utilizing and developing water resources in global and local scales (Condon et al., 2015;Huo et al., 2016;Li et al., 2016;Wang et al., 2016). Kim et al. (2008) linked the quasidistributed watershed model (SWAT), with the fully distributed groundwater model (MODFLOW). The integrated model was applied to the Musimcheon Basin in Korea where the Spatio-temporal distribution of the groundwater level was simulated while improving the recharge and evapotranspiration calculated by SWAT. Rodriguez and Cello et al. (2008) used the coupled HEC-RAS-MODFLOW model to adjust the water surface profiles along open drain channels by upgrading the Drain module available within the MODFLOW model. Baily et al. (2017) presented SWAT-MODFLOW using a graphical user interface.
Additionally, GIS-based models are also being applied to address the recharge volume of aquifers from the hydrogeology perspective or groundwater potential zones (Gachari et al., 2011;Minor et al., 2007;Singh et al., 2013). The most vivid method applied in the spatially auto-correlated variable like groundwater level is the kriging method namely, Ordinary Kriging technique. Aboufirassi et al. (1983) employed Universal kriging to estimate the water table for Souss aquifer in central Morocco. Pucci and Murashige (1987) used kriging for optimizing data collection and utility in a regional groundwater investigation in central New Jersey and confirmed that kriging is a useful tool especially in areas lacking enough data for developing a water table management network. Hoeksema et al. (1989) applied the cokriging method to estimate the groundwater level at unknown points. A similar study was done by other researchers where the kriging method was used to estimate the water levels at wells (Desbarats et al., 2002;Olea & Davis, 1999).  used the kriging method on 31 wells in evaluating and optimizing the groundwater level observation networks. Ahmadi and Sedghamiz (2007) evaluated the spatial and temporal variations of groundwater level of 39 observation wells using kriging. In their later article (Ahmadi & Sedghamiz, 2008), the kriging and co-kriging methods were applied for groundwater depth mapping in southern Iran.  used Artificial Neural Networks (ANNs) to estimate the temporal prediction of the water level and applied the kriging method to spatial parts. Ruybal et al. (2019) used spatiotemporal kriging in evaluating groundwater levels in the Arapahoe aquifer.
Conforming the former researchers, this study combines the geo-spatial kriging tool with a groundwater model. In this study, the Kriging method was applied on the residuals of the numerical model (MODFLOW) generated by the TWDB (Texas Water Development Board) for the Edwards-Trinity (Plateau) aquifer to improve the estimation of the water table spatially. The study was done for the years 1995 through 2000 where 90% of the observation data was used for model simulation followed by cross-validation with the remaining 10% of the observation data. To the authors' knowledge, no prior efforts have been done using the technique adopted in this paper. Additionally, there are limited studies that have resulted in such a significant improvement in groundwater level predictions that we were able to achieve through this research.

Study area
The Edwards-Trinity (Plateau) Aquifer, as shown in Figure 1, expands over west-central Texas between 97° and 105° west longitudes and between 29° and 33° north latitudes. The topology of the aquifer is known as a plateau, lightly leveling from about 610 m (2000 ft) above sea level in the southeast to about 915 m (3000 ft) in the northwest. The precipitation in the region ranges between 86 cm (34 in) in the east to 30 cm (12 in) in the west. The maximum average annual temperature for the study area ranges between 23 °C (73° F) in the Trans-Pecos uplands to 26 ° C (79° F) in southern Val Verde County. The study area had numerous drought events during the past hundred years (Anaya & Jones, 2009). Yet, the drought events are expected to have minimal impacts for the study period (1995 through 2000) as the last drought was during the 1950s.

Data
Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 16 September 2021 doi:10.20944/preprints202109.0293.v1 Data on groundwater levels, climate data, boundaries shapefiles, data on hydraulic properties of the aquifer, elevation data, and MODFLOW simulated heads were required. All these data were obtained from the TWDB website (TWDB, 2020). The water levels were obtained from the Groundwater Database (GWDB) consisting of shapefiles with geospatial information on water level and quality for management, monitoring, and characterization of the water in the Edwards-Trinity aquifers. The precipitation data was gathered using raster data for the years 1995 through 2001. The boundaries and hydraulic properties were gathered from separate shapefiles consisting of aquifer boundaries, model boundaries, Texas county boundaries, hydraulic conductivity for each one sq. km. The elevation data was obtained from the raster files consisting of DEM (Digital Elevation Model), top and bottom elevation of the Edwards and Trinity aquifers. Lastly, MODFLOW generated heads were obtained as publicly accessible binary files produced from the MODFLOW model developed by TWDB. Figure 2 shows the location of observation wells for all the years (1995)(1996)(1997)(1998)(1999)(2000).  Figure 3 represents a schematic diagram of the research methodology adopted. As shown in the figure, the general scheme of this study consists of three major components -1) Data Preparation, 2) Calibration, and 3) Cross-Validation and includes a series of steps used interactively as listed below:
Mapping MODFLOW simulated groundwater heads (model imported from TWDB) into their corresponding coordinates and overlap with observation data to find the MODFLOW estimated values in observation point.
2. Subtracting the observed groundwater head with MODFLOW simulated head and consider as the model residuals.

3.
Dividing the residuals into two separate dataset, 90 percent of data for fitting kriging methods (calibrating residuals) and 10 percent for validating part (validating residuals), in a random selection.
4. Pre-evaluating the calibrating residuals and fit kriging method to generate the estimated residual map for study domain. 5.
Comparing the validating residuals with estimated one to evaluate the accuracy of kriging method.

Data preparation
To minimize the distortion of the area and appropriate evaluation of the data, the NAD 1983 Texas Centric Mapping System Albers was used. First, the data were categorized into two groups -Observed data and MODFLOW generated output. The observed data consists of the mean values for the water tables during the winter season (winter months were selected as the water tables are expected to be more stable during winter and thus minimizes the anthropogenic effects due to uncertainties in demand seasons). The observed data were then provided with coordinates and transformed into a shapefile. For the MODFLOW generated data, the model grid shapefile produced by the MODFLOW program was employed in the study area and saved based on the needed attributes for future analysis. The binary files consisting of groundwater head data were converted to a text file using Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 16 September 2021 doi:10.20944/preprints202109.0293.v1 python programming. The extracted head data were matched with the corresponding observed data by locating the observation data in the grid shapefile. Finally, the shapefile of observation data was overlaid on the MODFLOW grid shapefile to determine the value of the MODFLOW generated head for each observation point. The MODFLOW generated heads were compared with the observed data for the corresponding locations thus obtaining the residuals.

Kriging method
To find the spatial distribution of the model residuals, the interpolation method, kriging, was applied. Kriging is a regionalization technique that incorporates the autocorrelation between known data values in its estimation of values at unmeasured sites. The kriging technique accepts irregularly spaced data and uses the measured values for closely surrounding points and is reproducible. Also, kriging is a robust interpolator at measured data points and calculates an error of estimate (Dunlap & Spinazola, 1984). Equation 1 shows the stochastic function that is considered in the simple kriging method: where Z(s) is the estimating function -which is simulated groundwater head elevation in our study, µ is a constant value representing the mean value of the groundwater head and ɛ(s) is a random error function regarding the model error from the observed records.
Estimating the value of Z(s) relies on two primary assumptions of considering that variables are randomly distributed, and they preserve second-order stationary condition respect to the location (Wackernagel, 2013), which means: (2) where h is a vector that connect point s to point s+h. Equation (2) implies that the expected value [ ( )] is constant in all domains as represented in equation (1) by the constant µ. However, in reality this value fluctuates from place to place due to inherent trends and variability in data. To deal with this problem, it is usually assumed that the estimating variable (groundwater head) comprise two components ( ) = ( ) + ( ) (4) which m(s) presents the dynamic part and ( ) is the statistical component of the estimating variable which includes the specially correlated random variable. The deterministic part of the variable Z consists of auxiliary variables known as drifts. Since m(s) is treated as a deterministic part, it can be determined in separate process and summed up with residuals random function component e(s). With this perspective, to come up with the nonstationary problem, the auxiliary variables are usually generalized by using the linear regression followed by a simple kriging in statistical part of variables on the regression residuals (Wackernagel, 2013). This method is called universal kriging (UK). The other way to solve the problem is that deterministic part estimated by external drifts generated by auxiliary models such as numerical models considering that there is a linear correlation between external drifts and auxiliary variables (Chiles & Delfiner, 2009). This type of models is called kriging with external drift (KED) models. KED is usually used not only to better understand the behaviour of the predicting variables but also to minimize the sources of errors resulted by uncertain parameters in deterministic models (Bourennane et al., 2000;Bourennane & King, 2003;Hudson & Wackernagel, 1994). The main idea behind this approach is that each two parts of estimating variable (deterministic and stochastic parts) can be stimulated separately, assuming there is no cross correlation between variables to skip some computational steps.
In this study, the deterministic part of the head variable is estimated by existing models. We then applied the ordinary kriging (OK) to the residuals of the numerical model for statistical part assuming the mean value of error is not known and needs to be obtained over optimization process by minimizing the variance of errors. However, we add intermediate checking and correcting stage before fitting OK on the residuals to remove the trend and anisotropy that may exist due to uncertainties in simulation process of modelling.

Data investigation
In standard statistical problems, the first step before starting to model the phenomena is to examine the data. In this study, a commonly used statistical tool, histogram, was used to see if the MODFLOW residuals are normally distributed. As shown in Figure 4, the residuals in general, follow the normality pattern. However, some level of negative skewness and distortion of normality was observed in the data which could be attributed to the complexity of nature.  Next, the residuals were checked for spatial correlation as the data tends to share some information with their neighbors. A variogram method was used to estimate how strongly data are related to each other. The variogram in space is usually implemented to check two major assumptions in application of kriging method -stationary in space (variance is independent of location) and isotropy (variance is independent of direction). Thus, a semivariogram cloud or plot was used to provide better visual understanding of the data distribution and to detect any possible trends or geometric anisotropic behaviour. Figure 4 shows the semi-variogram clouds of water level residuals calculated from subtraction between observation points and the MODFLOW outputs for the years 1995 to 2000. The Difference Squared (γ) in Figure 5 represent the MODFLOW residuals dissimilarity which is defined as: where s is sample location and Z is the residual value. The reddish circles in the figure indicate that there is a strong gradient observed at short distances, as the value of the semi-variogram for each pair in these areas is significantly high. This significant change can be an indication of non-stationarity in higher ranks caused by neighboring drainage areas or rivers. However, to get a proper conclusion over the sources that resulted in non-stationary, more investigation needs to be considered which is beyond the scope of this study. The semi-variogram was further used for interpolation where semivariance was used to represent the expected value of the residuals' dissimilarity. A theoretical model was fit into the sample data. Among others, the spherical model produced the best representation of the sharp increase and bending over of semi-variance for most years except year 2000. As observed from the sample variogram in Figure 6, the Nugget variance of 50 and range of 50000 was defined for fitted theoretical semi-variogram.

Model simulation
Ordinary Kriging method was applied on the MODFLOW residuals. Not all years have observation data available for the same locations thus Kriging was applied for each year separately with randomly selected observation points. To cope with the problem of existing trend, which is raised from uncertainties of simulated model and reflected to the residuals, the first-order trend was assumed in the data and the exponential Kernel function was applied. The Spherical model was chosen to fit on a semivariogram, and the maximum number of neighbors affecting the predicted data was limited to ten points. The kriging when applied as a single model resulted in higher mean root square error (not presented in this paper), thus the area was divided into four quadrants with 45 offsets to minimize the distortion due to anisotropy as observed in Figure 5. Figure 7 shows the prediction for Ordinary Kriging applied to Edwards -Trinity aquifer between the years 1995 to 2000. High residuals (as shown by dark red and blue color) were observed for the areas outside the boundary while the residual values are minimal inside the study area. Also, similar observations were observed in deviations as shown in Figure 8. The high residuals or larger deviation (as shown by purple color) outside the boundary is attributed to the missing observation data. Also, some areas inside the study area showed relatively higher residuals with larger deviations which could also be due to the missing observation data.

Model Validation
Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 16 September 2021 doi:10.20944/preprints202109.0293.v1 As mentioned before, the entire dataset was divided into two sets -a simulation dataset (90% data) and a validation dataset (remaining 10% data). First, kriging was applied on 90% of the data where the validation was performed automatically using the ArcGIS Geostatistical Analyst toolbar (ArcGIS Desktop, 2021). The Geostatistical Analyst toolbar provides the measured and predicted values and the standard error for each point. During the process, the software keeps any individual points separate from other data (referred to as dataset) for parameter estimation. The parameters are then estimated using every such datasets and the process continues until optimal parameters that best-fit the entire datasets are obtained. Table 1 presents the average water head for observation points along with predicted water heads using MODFLOW and Kriging, while Table 2 provides a detailed summary of kriging application on the MODFLOW residuals. As observed in Table 1, the residuals obtained from the MODFLOW simulation are much higher as compared to the residuals obtained after the application of kriging. The residuals significantly reduced from an average value of approximately 37 m to less than 1 m, with a standard error of approximately 0.5 m.
As a next step, a manual cross validation was adopted where the calibrated model was applied to the remaining 10% of the data. Tables 3 presents the average water heads for observation points along with predicted water heads using MODFLOW and Kriging. Similar to observations in Table  1, higher residuals were obtained for MODFLOW simulation for all the years. However, after application of kriging on MODFLOW residuals, the average absolute error of approximately 31m (from MODFLOW simulation) reduced to less than 5 m. Similar reduction was observed in the residuals' average standard error where the values reduced from 9.7 to 4.7. The observed reduction in standard error thus indicates that average value of residuals over entire period can be a good estimation for each year separately. Table 4 provides a detailed summary of the kriging application on the residuals.  For all the years, the water head for observation points obtained after kriging is much closer to the observed values. However, for the year 2000, three observation points (11, 12, and 15) showed higher residuals after the application of kriging as shown in Table 5. This could possibly be due to the higher prediction accuracy of MODFLOW simulation in those points. The table with water heads for all the years (1995 through 2000) has been provided as a supplementary file in the Appendix 1 (Table A1).  The kriging method was applied to improve spatial confidence in groundwater-level predictions at unsampled locations. Kriging was applied on the MODFLOW residuals for the groundwater levels in the Edwards-Trinity aquifer in Texas. The study was done for the years 1995 through 2000 where 90% of the observation data was used for model simulation followed by validation with the remaining 10% of the observations. The Kriging method significantly improved water level predictions. The average absolute error for of approximately 31m (from MODFLOW simulation) reduced to less than 5 m after the application of kriging on MODFLOW residuals. Also, the average residuals' standard error decreased from 9.7 to 4.7, which indicates that average value of residuals over entire period can be a good estimation for each year separately. With improved water level predictions, geostatistical tools such as kriging can be used to produce more accurate potentiometric surface maps. Such improved results and accurate monitoring of groundwater resources will lead to the sustainable use of groundwater resources while also aiding in efficient and effective conjunctive management of surface and groundwater resources.