Mapping of the Water Table Levels of Unconfined Aquifers Using Two Interpolation Methods ()
Received 13 June 2016; accepted 31 July 2016; published 3 August 2016
1. Introduction
Boreholes in situ investigations performed with Standard Penetration Test (SPT) represent the most common type of subsurface investigation in the Brazilian geotechnics. Information obtained from this subsurface investigation is widely used in foundation engineering (Nspt values, soil texture, water table depth, compactness or consistency of the soils). The information can also subsidize the planning, design and execution of others types of civil works (e.g., exacavations, slope stability and earth works) and provide parameters to urban and environmental management.
The SPT test is used worldwide as an important geotechnical investigation and the first attempt at normalization was made in 1958 by American Society for Testing and Materials (ASTM). The North American standards are used frequently in South America. Specific standards established by ABNT [1] are adopted in Brazil.
The use of deterministic and geostatistical interpolation methods to predict the spatial distribution of water table levels and standard penetration tests values (Nspt) has been studied in engineering geology and geotechnics research in the past twenty years, some using Geographic Information System (GIS)-based investigation. The advance of computational tools is providing the intensified use the geostatistical methods to predict the spatial distribution of geological and geotechnical variables. But for some of these variables, deterministic methods can produce results as or more reliable than geostatistical, with the advantage that they demand less specific theoretical knowledge [2] - [8] .
This work was part of a PhD research project aimed at the development of a Decision Support System using a geological-geotechnical database to subsidize foundation engineering activities and other activities related to urban and environmental management [9] . This paper focuses on the prediction of the water table levels (unconfined aquifers) using deterministic and geostatistics interpolation methods. The study area, totaling 59.8 km2, comprises 40 districts and important geological formations of João Pessoa City, Paraiba State, Brazil (Figure 1).
The geological-geotechnical database contains the following main information layers: the Digital Elevation Model (DEM) obtained from topographic maps at 1:10,000 scale (contours of 5 meters); a map of engineering geological units adjusted at 1:10,000 scale, but mainly based on geological mapping at 1:50,000 scale; and the data from soundings boring logs with Standard Penetration Test (SPT) and deep wells.
Geostatistical interpolation methods consider the spatial characteristics of autocorrelation regionalized variables, whereas the deterministic methods consider only the amounts involved in the interpolation. Both methods rely on the similarity of nearby sample points to create the surface. Deterministic techniques use mathematical functions for interpolation. Geostatistics relies on both statistical and mathematical methods, which can be used to create surfaces and to assess the uncertainty of the predictions.
The main methods of deterministic interpolation available in GIS software include triangulation with linear interpolation, Inverse Distance Weighted (IDW), spline and spline modified. This study used a spline modified method named Topo to Raster command in the raster interpolation tools of the ArcMap module of the Arcgis Desktop software.
This routine uses an interpolation technique specifically designed to create a surface that more closely represents a natural drainage surface and better preserves both ridgelines and stream networks from the input contour data. The method is essentially a discretized thin plate spline technique for which the roughness penalty has been modified to allow the fitted Digital Elevation Model (DEM) to follow abrupt changes in terrain, such as streams, ridges and cliffs [10] . ESRI [11] says that the Topo to Raster routine is optimized to have the computational efficiency of local interpolation methods, such as IDW interpolation, without losing the surface continuity of global interpolation methods, such as spline (a deterministic method) and kriging (a geostatistical method).
Geostatistical methods have the capability of producing a prediction surface, but they can also provide some measure of the certainty or accuracy of the predictions. Ordinary kriging is a geostatistical method that is common in GIS analysis tools. Kriging is similar to IDW (deterministic method) in that it weights the surrounding measured values to derive a prediction for each location. However, the weights are based not only on the distance between the measured points and the prediction location but also on the overall spatial arrangement among the measured points.
In the Kriging method, the empirical semivariogram is a means to explore the assumption that things that are close to one another are more alike than those farther away (quantified as spatial autocorrelation). The following models of empirical semivariogram are available: circular, spherical, tetraspherical, pentaspherical, exponential, Gaussian, rational quadratic, hole effect, K-Bessel, J-Bessel and stable. The selected model influences the prediction of the unknown values, particularly when the shape of the curve near the origin differs significantly.
2. Method and Materials
The ArcMap module of the Arcgis Desktop software (versions 9.1 to 10.2) was used to combine information from surface investigations, subsurface investigations (percussion drilling with SPT essays) and interpolation of the water table level. The ArcCatalog module was used to manage the database files. The ArcScene module was used for tridimensional visualizations of some thematic maps. The main steps, acitivities and results of the study are presented in Figure 2.
The first step included data collection and preparation, from the digital topographic base to the information from the boring logs in the city of João Pessoa (capital of Paraiba State, Brazil, Figure 1). Information from the boring logs was obtained by selection between 2000 geotechnical profiles provided by a local foundation enterprise that works within the study area.
The topographic information was supplied by several institutes and boards from the city council, such as contour lines of 5 meters at 1:10,000 scale. Elevation points and all drainage nets and streams were obtained and digitalized after a georeferencing process with geographical coordinates. A map of the engineering geological units was adjusted to this digital topographic base at 1:10,000 scale, based on geological mapping at 1:50.000 scale available to all study area and field research.
The DEM modeling process was then undertaken to provide a continuous topographic surface of the study area to adjust the start elevations of all boreholes considered in the analysis. The Topo to Raster ArcMap routine was used. Acell size of 4 meters was defined as the interpolation grid, which is 20% lower the maximum permissible cartographic error at 1:10,000 scale (0.5 mm or 5 m). The filter and the fill tools were applied to obtain
![]()
Figure 2. Main steps, activities and results of the study.
the final hydrologically correct interpolated raster surface.
The final DEM was qualitatively and quantitatively evaluated to ensure that the data and parameters supplied to the program resulted in a realistic representation of the topographic surface. The qualitative evaluation created contours from the new surface with the Contour tool and compared them to the input contour data. The quantitative evaluation required to withholding a percentage of the input data from the interpolation process. After generating the surface, the height of these known points was subtracted from the generated surface to examine how closely the new surface represented the true surface. These differences were used to calculate the root mean squared error (RMSE) using the following equation [10] :
(1)
where n is the total number of sample points, and zi and z1,i are, respectively, the measured and estimated values at the ith location.
The second stage consisted of field activities and desk work to select the boreholes, to check their locations and to insert the geotechnical information in the database. The structure of the database of the boreholes is a vector file (shapefile) of point feature type with two tables of attributes attached.
The main attribute table contains fields that are unique to each borehole and do not vary with depth: identifier borehole; location coordinates; starting elevation of the hole (extracted from the final DEM); date of survey; depth of water level and total borehole depth. The secondary attribute table contains the identifier borehole field and the fields that vary with depth: Nspt; soil type; soil texture; soil color and compactness or consistency index. The borehole identifier field allows the joining of the two tables needed to conduct the interpolations of the different variables for deep levels.
The third stage comprised interpolation of the water table data using deterministic (Topo to Raster) and geostatistical (Ordinary Kriging) methods. The analysis started in the northeast portion of João Pessoa City, where more boreholes data were available. This area comprises alluvial deposits (Quaternary) and recent marine sediments with fluvial contribution (Quaternary) located along the Atlantic coast and the mudstones and sandstones of the Barreiras Group (Tertiary).
For interpolation of the water table levels, the unconfined aquifer systems associated with engineering geological units were considered. The compiled geodata were separated according to the period of the year and the associated rain season. In the study area, the rainy season is from March to August and the dry season is from September to February.
The algorithm used for the topo to raster method is essentially a discretized thin plate spline technique [12] for which the roughness penalty has been modified to allow the fitted interpolated surface to follow abrupt changes. The method uses the following general formula:
(2)
where: j = 1, 2, ・・・, N; N is the number of points; λj are coefficients found by the solution of a system of linear equations; rj is the distance from the point (x, y) to the jth point; T(x, y) and R(r) are defined differently, depending on the selected option.
The modeling with the kriging method, in addition to the interpolation of the water table, also produced a standard deviations map to provide an evaluation of the accuracy of the results. According Oliver [13] , the general formula for the kriging interpolator is formed as a weighted sum of data:
(3)
where Z(Si) = the measured value at the ith location; λi = an unknown weight for the measured value at the ith location; S0 = the prediction location and N = the number of measured values.
After interpolation processes, a validation step was conducted out by using data of water table levels that were not previously used for the estimates. The estimated values were compared to calculate the error between the two situations.
Maps of the water table depths were also produced to further analyze the quality of the water table surfaces interpolated by both methods. These maps were obtained by subtracting the Digital Elevation Model (DEM) of the water table surfaces obtained by both methods, on a cell-by-cell basis math operation.
3. Results and Discussion
3.1. Database, Geological and Geotechnical Aspects
The DEM is an important layer in any geotechnical database because all subsequent spatial analyses are dependent on the accuracy of it. In the present research, it was fundamental to provide the starting elevations of all boreholes considered in the interpolation analysis. The DEM also provided other auxiliary maps and layers (Slope map, Hypsometric map, Aspect map, Hill shade map) that were used in the study. Both validations (qualitative and quantitative) indicated good accuracy of the DEM. The RMSE obtained in the quantitative validation was 0.67 meters, or approximately 13% of the maximum permissible cartographic error to the 1:10,000 scale. This RMSE can be considered adequate to scale and objectives of the study.
The marine sediments with fluvial contribution located along the Atlantic coast (Quaternary) and the Barreiras Group (Tertiary) are the main engineering geological units in the study area considering their extensions and demand for new civil works. The first comprises a narrow coastal strip, with higher expression in the extreme northeast of the study area (Figure 3) and consists of well-sorted sands, fine to medium grained and remains of marine animals (shells). In the estuaries of the major drainages, these sediments become finer material (silt and clay) and organic matter.
The Barreiras Group consists of mudstones, sandstones and silty clay and sandy-conglomeratic beds, without the presence of fossils. The sediments are often poorly selected, with predominance of sand and clay, with red colors and variegated whitish horizons.
These units also define two different terrain compartments. The marine sediments occur at elevations from zero to 10 meters, and the Barreiras Group occurs at elevations from 10 to 50 meters. The transition between these two compartments occurs abruptly, forming a steep cliff along the coast. The recent alluvial deposits (Quaternary) occur in the valleys of the main drainages and rivers along different elevations (Table 1 and Figure 3).
3.2. Interpolation of the Water Table
For the interpolations of the water table levels (unconfined aquifers), it was considered the alluvial deposits, the recent marine sediments and the Barreiras Group engineering geological units comprised an area of 59.8 km2. Of the 311 boreholes that reached the water table, 170 were made in the dry season and 141 during the rainy period of the year. The interpolations used the boreholes from the dry season because they outnumbered the ones carried out at the rainy season (Figure 3).
In addition to the sounding boring logs, additional information regarding water table level was provide from 16 groundwater wells (data from the dry season), which that made it possible to obtain more data from drilled wells located in terrains where the depth of the water table was higher than 20 meters, the maximum average depth of the percussion boreholes.
Additionally, 127 elevation points over natural drainages were used to model the water table surface. This procedure was used to obtain a water table surface adjusted to the natural drainage system present in the
![]()
Figure 3. Engineering geological units, elevations, SPT boreholes and weels.
![]()
Table 1. Percentage areas of the terrain elevations by engineering geological units.
interpolation area. The average density of all data considered in the interpolations of water table for both methods was 5.2 points/km2 (170 boreholes + 16 wells + 127 points over the drainages/59.8 km2).
The site disposition demonstrates that most of sounding data belong to the marine sediments engineering geological unit, with the terrain elevation varying from 0 to 5 meters high. At this unit, water table levels were measured from −1 (zero elevation equal to seal level) to 5 meters high. As expected, the water table levels obtained in the dry seasons throughout the years were lower than those measured during rainy periods.
3.2.1. Topo to Raster Method
The topo to raster method considers several parameters to perform the interpolation. The input feature containing the water table levels to be interpolated was a point feature class (170 boreholes, 16 wells and 127 points over drainages). The others input features were a line feature class of stream locations (natural drainages and rivers) and a feature class containing a polygon that represents the outer boundary of the interpolated surface. A cell size of 4 meters was used for the interpolation of the water table, the same used for the digital elevation model. The others parameters were defined considering the characteristics of the interpolated data and the recommendations of the topo to raster command.
As illustrated in Figure 3, the water table surface (dry season) interpolated using the topo to raster method was reclassified into six intervals of elevations (meters): −9.4 (minimum elevation obtained) to 0, 0.1 to 2, 2.1 to 5; 5.1 to 10, 10.1 to 20 and 20.1 to 43.7 (maximum elevation obtained). Table 2 presents the percentage areas of the water table levels by engineering geological units, detailing the negative values.
The interpolated water table surface presents a coherent spatial arrangement when considering the general topography of the area studied. The topo to raster model also shows adequate behavior at valley floors.
The lowest water table level obtained was −9.4 meters. This anomalous value occurred in a very restricted area of the unit of the Barreiras Group and is associated with interpolation deficiency (less than 0.1% of the total area, Table 2). The water table levels obtained tend to be −1.9 to zero meters near the coastal contour line (15% of the marine sediments total area), gradually increasing towards the continent, and reaching values 0 to 5 meters over all marine sediments (85% of the total area, Table 2 and Figure 4).
There is a sharp increase in these values in the contact between the marine sediments and the mudstones/ sandstones of the Barreiras Group, which is characterized topographically by a steep cliff along the coast (See Figure 2 and Figure 3). In the Barreiras Group, the estimated water table predominantly has values greater than 5 meters (85.5% of its total area, see Table 2). The largest sector of this unit, with water levels between 2 to 5 meters, is in the south portion of the interpolated area, near the valleys of the Jaguaribe and Timbo Rivers. In stretches of higher ground of the Barreiras Group, the interpolated water table reached elevations between 10 to almost 44 meters (almost 64% of its total area)). The alluvial deposits presented water levels between 0 and 20 meters, but 51.5% of the total area had water levels between 2 and 5 meters. Only 0.1% of this unit had negative water levels (see Table 2 and Figure 3).
3.2.2. Kriging Method
Initially, variograms and covariance functions were prepared to estimate the statistical dependence values that depend on the model of autocorrelation (Angle direction = 0˚; Angle tolerance = 90˚). Estimates by ordinary kriging from the obtained variogram proved satisfactory in the analysis by cross-validation, especially in the stretch where there is a greater number of samples. The best fit was obtained for the spherical model with lag (distance) of 300 meters and maximum distance of 2016 m.
The water table level was interpolated using the ordinary kriging method, the variogram model defined above, and the same cell size and data used in the topo to raster interpolation. Unlike the topo to raster method, where it is possible to set a specific area to perform the interpolation, in the kriging method, this limit depends on the distribution of the sampling points and the definitions used in the interpolation model (lag and maximum distance). To facilitate comparison of the results, a cut of this interpolated area was made, considering the same limit adopted in the top to raster method.
Figure 5 presents a map with the water table surface interpolated using the kriging method. To facilitate comparison between the two methods, this water table surface was sliced in the same intervals used in the topo table to raster interpolation but considering the extremes values of water levels obtained in the kriging interpolation (−9.3 and 43.5 meters).
![]()
Figure 4. Water table levels interpolated using the deterministic method topo to raster (dry season).
The ordinary kriging produced a water table surface similar to the topo to raster method (see Figure 3 and Figure 4). The extreme values obtained with the kriging method were slightly lower than those obtained by the topo to raster method, but this method showed a slight general tendency to be less affected by local values in relation to the kriging method.
The lowest water table level obtained was −9.3 meters using the kriging method. This anomalous value occurred in a restricted area of the unit of the Barreiras Group (0.2% of its total area, Table 3). The water table levels obtained were −1.9 to zero meters near the coastal contour line (14.6% of the marine sediments total area), gradually increasing towards the continent, and reaching values approximately 0 to 5 meters over the entire marine sediments unit (85.3% of the total area, Table 3 and Figure 5).
The main difference between the two interpolation methods occurred in the south sector of the study area, on both sides of the Timbó Valley River, where the kriging method resulted in negative to zero values of water
![]()
Table 2. Percentage areas of water table levels by engineering geological units (topo to raster method).
![]()
Table 3. Percentage areas of the water table level by engineering geological units (kriging method).
level in the Barreiras Group (2.2% of the total area ranging from −1.9 to zero meters, Figure 4 and Figure 5 and Table 3).
A standard deviations map was also prepared to provide some measure of the accuracy of the interpolated water table using the kriging method. This map was obtained from the application of the square root to the predicted variance values generated by the ordinary kriging interpolation.
The values of standard deviations varied from 0.1 to 5.8 meters with an average value of 2.5 meters. These values were sliced into five intervals: 0.1 to 1, 1.1 to 2, 2.1 to 3, 3.1 to 4 and 4.1 to 5.8 meters (Figure 6). Table 4 presents the percentage area of these five classes of standard deviations in the three engineering geological units present in the interpolation area.
The marine sediments showed the best level reliability for the water level interpolated with almost 61% of its total area with standard deviation ranging from 1.1 to 2 meters. The Barreiras Group presented the worst reliability for the interpolated water level, with almost 46% of its total area having standard deviations ranging from 3.1 to 4 meters and 11% of its total area having standard deviations ranging from 4.1 to 5.8 meters. All engineering geological units presented low percentages of total areas with standard deviations less than 1 meter (Table 4). The results of the standard deviations associated with the interpolated water table surface using the kriging method can be explained by the distribution of SPT boreholes in these units. The marine sediments have a much higher number of SPT boreholes than other engineering geological units.
3.3. Validation of the Interpolation Results
Table 5 presents a comparison between the measured and estimated values for the water table level obtained in the boreholes that were not included in the modeling process to validate results. The two methods presented similar individual errors. Considering these individual errors, the RMSE computed with equation 1 was 1.8 meters for both methods.
![]()
Figure 5. Water table levels interpolated using the probabilistic method ordinary kriging (dry season).
![]()
Table 4. Percentage areas of the standard deviations by engineering geological units (kriging method).
![]()
Table 5. Comparison between the measured water table levels and the estimated values using both methods.
*See Figure 6.
![]()
Figure 6. Standard deviations of the water table levels interpolated with the kriging method.
The individual errors and the RMSE obtained with the validation presented values compatible with the standard deviations obtained from the kriging interpolation method. The main limitations of this validation are the reduced number of SPT boreholes that could be used (only five) and the fact that all the boreholes are located in the marine sediments unit. The positive aspect to note is that the holes are well-distributed throughout the greater extent of this unit (Figure 6).
The minimum individual error (0.2 m) occurred in SPT borehole number 84, located in the extreme south region of the marine sediments unit for the topo to raster method. For the kriging method, the minimum error was found in SPT borehole number 83, located almost 1.7 km to the northwest of the previous borehole. The maxi-
mum individual errors for both methods were found in SPT borehole number 17, located in the northeast region of the marine sediments unit (Table 5 and Figure 6).
Maps of the water table depth were also prepared to further analyze the quality of the water table surfaces interpolated by both methods. Negative depths indicate where the interpolated water table surface failed, or where the water table level obtained with the interpolation was above the terrain surface. This situation should occur only where there is some type water body (e.g., river, lake, reservoir or ocean).
Figure 6 presents the map of the water table depth obtained using the topo to raster method, which was reclassified into six ranges (meters): −17.1 to 0, 0.1 to 2, 2.1 to 5, 5.1 to 10, 10.1 to 20 and 20.1 to 58.7. Figure 8 presents the map of the water table depth obtained using the kriging method. To facilitate comparison between the two methods, the water table depth was sliced in the same intervals used in the topo to raster interpolation, using the minimum water table depth obtained in the kriging interpolation (−16.8 meters). Table 6 presents the percentage areas of the water table depth classes by engineering geological units, considering both methods and detailing the negative depths (−17.1/16.9 to −10, −9.9 to −5, −4.9 to −2 and −1.9 to 0).
The results were similar for the two interpolation methods. The negative water table depths occurred mainly in the alluvial deposits, comprising 43.2% and 45.2% of the total area for the topo to raster and kriging methods, respectively. Of these totals, more than 25% ranging from −1.9 to 0 meters, for both methods. The Barreiras Group had 4.0% and 4.8% and the marine sediments had 1.0% and 0.8% of their total areas with negative depths, considering the topo to raster and kriging methods, respectively. Additionally, 97% of the total area of the marine sediments unit had water table depths between 0.1 to 5 meters for both interpolations methods (Figure 7 and Figure 8 and Table 6).
The negative water table depths or the failures of the interpolation methods are associated mainly with alluvial deposits units because sudden changes in the land conformation occur in the drainage valleys, making adjustments of interpolated phreatic surfaces difficult. These results can be enhanced by considering more sampling points on the drainages and rivers.
![]()
Table 6. Percentage areas of the water table depth classes by engineering geological units (MDE―water table surfaces obtained with both methods).
*1: Alluvial Deposits; 2: Marine Sediments; 3: Barreiras Group.
![]()
Figure 7. Water table depth (MDE minus the water table surface interpolated using the topo to raster method).
4. Conclusions
Phreatic surface interpolation (unconfined aquifers) was guided by comparison of two methods: deterministic topo to raster and geostatistical ordinary kriging. Both methods provide satisfactory and similar results.
The topo to raster method (deterministic) had a slight general tendency to be less affected by local values in relation to the kriging method. The first method also had the advantage of integrating the drainage flow system during the modeling process, which is a relevant aspect when addressing the water table levels of unconfined aquifers. The kriging method had the advantage of providing some measure of the accuracy of the interpolated water table levels.
The validation processes indicated that reliable water tables surfaces were obtained by both methods (RMSE = 1.8 meters). The negative water table depths or the failures of the interpolation methods are associated mainly
![]()
Figure 8. Water table depth (MDE minus the water table surface interpolated using the topo to raster method).
with alluvial deposits units. Interpolation tests were made using only the borehole and well data elevation. The use of 127 elevation points over natural drainages improved the results for both interpolation methods. These results can be enhanced by considering more sampling points on the drainages and rivers. These points should be carefully selected and with field control to ensure good correlation with the elevation of the water table.
The deterministic methods can provide as or more reliable results than the geostatistical demanding less theoretical knowledge. The geostatistical methods have the main advantage of producing a prediction surface and also provide some measure of the certainty or accuracy of the predictions, which is fundamental for probabilistic approaches in geotechnical engineering.
The maps of water table depth obtained by the difference between DEM and the interpolated water table maps provided, in map form, the performance of each interpolation method compared to each other, with the advantage that they also supplied preliminary hydrogeological models of the study area for engineering geotechnical works.
The study provided the possibility of making reliable previsions about water table levels and depths (unconfined aquifers) for 40 districts of João Pessoa City, which can be useful in different applications related to foundation works and other urban and environmental management activities. GIS proved to be an important tool during realization of all research steps, including database generation and structuring to spatial treatment of stored information.
Acknowledgements
The authors gratefully acknowledge the financial support of the Brazilian research funding agencies FAPESP and CAPES.