Abstract
Detailed knowledge of the uppermost water table representing the shallow groundwater system is critical in order to address societal challenges that relate to the mitigation and adaptation to climate change and enhancing climate resilience in general. Machine learning (ML) allows for high resolution modeling of the water table depth beyond the capabilities of conventional numerical physically-based hydrological models with respect to spatial resolution and overall accuracy. For this, in-situ well and proxy observations are used as training data in combination with high resolution covariates. The objective of this study is to model the depth of the uppermost water table for a typical summer and winter condition at 10 m spatial resolution over entire Denmark (43,000 km2). CatBoost, a state of the art implementation of gradient boosting decision trees, is employed in this study to model the water table depth and the associated uncertainties. The groundwater domain has not been the most prominent field of applications of recent hydrological ML advances due to the lack of big data. This study brings forward a novel knowledge-guided ML framework to overcome this limitation by integrating simulation results from a physically-based groundwater flow model. The simulation data are utilized to (1) identify wells that represent the uppermost water table, (2) augment missing training data by accounting for simulated water level seasonality, and (3) expand the list of covariates. The curated training dataset contains around 13,000 wells, 19,000 groundwater proxy observations at lakes, streams and coastline as well as 15 covariates. Cross validation attests that the ML model generalizes well with a mean absolute error of around 115 cm considering solely well observations and a MAE of <50 cm taking also the proxy observations into consideration. Quantile regression is applied to estimate confidence intervals and the estimated uncertainty is largest for moraine clay soils that are characterized with a distinct geological heterogeneity. This study highlights a novel research avenue of knowledge-guided ML for the groundwater domain by efficiently supporting a ML model with a physically-based hydrological model to predict the depth of the water table at unprecedented spatial detail and accuracy.
Introduction
A key state variable of the hydrological cycle is the depth of the uppermost water table, i.e., shallow groundwater, with a broad range of crucial societal, and environmental implications such as securing infrastructure, food production and sustaining ecosystems (Gleeson et al., ). Following the global analysis of water table patterns by Fan et al. (), up to one-third of the land area is influenced by shallow groundwater, being either directly groundwater-fed or having the water table or capillary fringe within plant rooting depths. More concretely, the shallow groundwater system plays a key role in building mitigation and adaptation measures to address climate change, as the uppermost water table controls greenhouse gas emissions from wetlands (Tiemeyer et al., , ) and as rising groundwater exacerbates flooding and falling groundwater intensifies droughts (Taylor et al., ). Floods can either be directly induced or intensified by the shallow groundwater system which has a special relevance for urban areas (MacDonald et al., ; Bricker et al., ). In an agronomical context, the shallow groundwater system is vital to meet crop water requirements in many agricultural settings while having adverse consequences for crop yield when the water table is too close to the surface (Kahlown et al., ; Zipper et al., ). Moreover, the uppermost water table constitutes a link between subsurface and the land-surface by affecting the energy balance and near-surface climatic conditions (Larsen et al., ; Maxwell and Condon, ).
The above mentioned relevancy of the shallow groundwater system requires versatile modeling systems to meet the demands of decision makers with respect to accuracy and spatial resolution. Accuracy and spatial resolution can be considered main bottlenecks in the advancement of physically-based numerical hydrological models. Increasing the spatial resolution does not necessarily go along with an improvement of the model accuracy, as established process descriptions and parametrizations are not always scalable (Beven and Cloke, ; Clark et al., ). At the same time, increasing the number of computational units results in a computational burden limiting the possibility to conduct thorough parameter calibration and sensitivity analysis. Due to the stringent parametrization and rigid model structure, conventional hydrological models cannot fully harness the wealth of readily available environmental big data. Nevertheless, physically-based models integrate decades of hydrological knowledge and are indispensable for integrated assessments of the hydrological cycle and simulating hydrological response under non-stationarity, i.e., climate change. As outlined by Shen (), Reichstein et al. (), and Nearing et al. (), machine learning (ML) is gaining increasing attention in the hydrological science to overcome some of the previously mentioned constraints of conventional physically-based models. ML has the advantage of being data flexible with respect to optimally utilizing the wealth of environmental big data while providing accurate predictions at low computational costs. However, ML lacks process descriptions which commonly restricts trained ML models to deliver predictions within observed ranges of the training dataset. In order to reconcile advantages of both modeling perceptions, knowledge-guided ML forms a promising new research avenue (Rajaee et al., ; Kraft et al., ). Knowledge-guided ML has the aim to integrate physical consistency into ML improve model performance and robustness. There exist multiple approaches to design knowledge-guided ML models as outlined by Read et al. (), Reichstein et al. (), Konapala et al. (), and others, which are also referred to as physics- or process-guided. A clear formal definition of these modeling approaches is still lacking, but they generally aim at integrating aspects of scientific knowledge into a ML model.
This study aims at modeling the depth of the uppermost water table at 10 m spatial resolution over entire Denmark for a typical summer and winter condition by implementing a knowledge-guided ML model. The physically-based information are obtained from the Danish national water resources model, that integrates groundwater and surface water processes (Højberg et al., ; Stisen et al., ). The physically-based model is used three-fold, (1) to derive threshold depths to select wells that reflect the shallow groundwater, (2) to augment training data at wells with incomplete pairs of summer-winter observations, and (3) to inform the ML model with the typical summer and winter water table depth using simulation results at 100 m resolution.
Decision tree based ML models are popular tools for geospatial modeling of environmental variables (Hengl et al., ; Tyralis et al., ). In this context a target variable, available as point data, is used in conjunction with maps of explanatory variables to curate a training dataset. Relationships between the target variable and the explanatory variables are established via the decision trees which, once trained, can be generalized to make predictions of the target variable at all grids. This framework has been successfully applied across the geosciences to model water chemistry indicators (Tesoriero et al., ; Erickson et al., ) soil properties (Møller et al., ; Hengl et al., ), subsurface redox conditions (Close et al., ; Koch et al., ), water table depth (Bechtold et al., ; Koch et al., ), and other variables. In such modeling frameworks, uncertainty can be quantified via quantile regression (López López et al., ; Tyralis et al., ). Gradient boosting is among the state of the art techniques to build decision tree models and, for this study we have employed the CatBoost implementation of gradient boosting decision trees (Dorogush et al., ; Prokhorenkova et al., ). Numerous studies apply ML to model the temporal water table dynamics by using various ML techniques (Sun et al., ; Guzman et al., ; Wunsch et al., ). These studies highlight successful applications, but are always limited to time series modeling, at a few selected sites. Despite these efforts, the spatial dimension is often neglected in the published studies. Fienen et al. (), Bechtold et al. (), and Koch et al. () are among the few studies that model the spatial variability of water table depth using ML. Knowledge-guided ML applications in the groundwater domain are for example the physics informed neural networks model developed by Guo et al. () that allows to solve partial differential equations with less calculational time.
The three main objectives of the paper are as follows: (1) to train a ML model to predict the uppermost water table depth at 10 m spatial resolution over entire Denmark, (2) to quantify uncertainty using quantile regression, and (3) to formalize a knowledge-guided ML framework that builds upon a physically-based hydrological model.
Materials and Methods
Study Area
This study is carried out for the entire land phase of Denmark, located in Northern Europe and covering an area of around 43,000 km2 (Figure 1). Denmark is generally flat with a maximum elevation of 170 m.a.s.l. and agriculture is the main land cover with around 70%. The landscape of Denmark was formed by a sequence of Pleistocene glaciations and postglacial processes. The soils of eastern Denmark are dominated by Weichselian moraine sediments with a moderate clay content, whereas western Denmark is characterized by older moraine sediments originating from the Saalian age intertwined by sandy Weichselian outwash plains.
Figure 1
In Denmark, the groundwater system is under pressure as a consequence of climate change and abstractions, revealed by quantitative modeling assessments conducted by Henriksen et al. (
Data
In order to make seasonal estimates (typical winter and typical summer) of the shallow water table at high resolution, an initial comprehensive data processing has been conducted to curate a high-quality dataset containing water level observations at shallow wells, additional groundwater observations as well as national maps of explanatory variables.
Physically-Based Model
The national water resources model of Denmark (DK-model), which has previously been refined from the original 500 m resolution to an updated 100 m resolution, was employed in this study (Højberg et al.,
Figure 2

Simulation results for the depth of the uppermost water table from a national physically-based model (PBM) at 100 m spatial resolution: The left panel shows the simulated median summer condition (JJA). The center panel shows the simulated median winter condition (DJF). The right panel depicts the median seasonality, calculated as summer minus winter.
Groundwater Observations
The open-access Danish well database (Jupiter) contains around 100,000 wells with at least a single water table observation in the selected 30-year period between 1990 and 2019. Based on the national well dataset, we first identified the wells that represent the shallow groundwater system, i.e., characterize the uppermost water table. Given the geological complexity of Denmark, this can vary from just a few meters for locations with a thick surficial clay layer to several tens of meters for the sandy outwash plains. In order to find suitable threshold depths, that classify a well as being either shallow or not, we combined the national soil map (Figure 1) and the simulated uppermost water table (Figure 2). For this analysis, the midpoint intake depth of each well was used, which is relative to the top and the bottom of the intake. The 95th percentile of the simulated water table depth was calculated for each soil type which guided the definition of the presented threshold depths applied to the midpoint intake depths of the wells (Table 1). Wells were only selected if the intake depth was lower than the soil type dependent threshold depth. A spatially soil type distributed threshold depth, guided by a physically-based model, has the advantage to reflect the natural conditions best possibly. Alternative, a constant threshold depth of e.g., 10 m would result in the selection of wells that do not reflect the uppermost water table in clayey settings where a well may be placed in a sand unit below a 6 m surficial clay layer containing the uppermost water table.
Table 1
| Soil Type | Max depth (m.b.g.l.) |
|---|---|
| Sand | 15 |
| Moraine Sand | 20 |
| Peat | 3 |
| Chalk | 20 |
| Unclassified | 3 |
| Moraine Clay | 3 |
| Sand and Gravel | 10 |
| Marine Sand | 3 |
| Silt | 3 |
Threshold depths that define a shallow well with respect to soil types (Figure 1).
Values of max depth were derived from the 95th percentile of the simulated water table depths for each of the given soil types and were applied to the midpoint intake depths of the wells.
After applying the soil types dependent threshold depths, summer, and winter depths to the uppermost groundwater were calculated at each well. For this task, summer refers to observations from the months June, July, and August (JJA) whereas winter refers to the months December, January, and February (DJF). Inter-annual variation was not considered during this processing step and in case a well contained several summer or winter observations a median was calculated. This resulted in 13,047 well, as shown in Figure 3, of which 1,378 wells had both a summer and a winter observation, 5,651 only a winter observation and 6,018 only a summer observation. For the wells that were missing either a winter or a summer observation the simulated median seasonality (Figure 2) was used to augment the missing season, i.e., winter = summer—seasonality, and vice versa. Further, the training dataset was extended by additional observations that reflect proxy observations for surficial groundwater levels with a depth of zero. In total, 19,074 groundwater connected lakes with an aerial extent of at least 100 m2 were used for this purpose. Additionally, 1,000 points, randomly placed along each, the stream network and along the coastline, were added to better represent these under sampled settings where the water table is expected to be at the surface year-round. This resulted in a total number of 34,061 groundwater observations spread across Denmark, yielding a density of around 0.8 observations per km2. Based on the well data only, the average summer depth is 3.6 and 3.1 m for the winter season with a standard deviation of 2.8 m for both.
Figure 3

Training dataset containing wells (left panel) and additional observations (right panel). The dataset contains 13,047 shallow wells with both a winter and a summer observation. The additional observations comprise 19,074 groundwater connected lakes, 1,000 stream points and 1,000 coastal points. All additional observations were defined with zero depth of the uppermost water table.
Explanatory Variables
In Table 2, an overview of the covariates used to model the depth of the uppermost water table is presented. In total, 15 covariates were assembled as input to the ML model. This list comprises information on soil texture, geology, topography-based characteristics, water body proximity, land cover, and outputs from a hydrological simulation with the DK-model. The native spatial resolution of the covariates varied, but all covariates were resampled to 10 m to be in agreement with the defined output resolution. For the resampling we used a bilinear interpolation method for the continuous variables. The water body proximity was expressed as both the vertical and horizontal distance to the nearest water body, which contained rivers, lakes, and the coastline.
Table 2
| Variable | Abbreviation | Description | Source |
|---|---|---|---|
| Clay content 0–30 cm | ClayA | Adhikari et al., | |
| Clay content 30–60 cm | ClayB | Clay content in percentage for four soil layers at | |
| Clay content 60–100 cm | ClayC | 30 m resolution | |
| Clay content 100–200 cm | ClayD | ||
| Thickness of top clay | ClayThick | Thickness of the uppermost clay layer at 100 m resolution | DK-model |
| Landscape typology* | LType | Geomorphological classification in 13 classes as polygon shape file | Breuning-Madsen and Jensen ( |
| Land Use* | LUse | 7 land use classes at 100 m spatial resolution | Levin et al. ( |
| Degree of urbanization | Urban | Percentage of grid cell that is paved at 10 m spatial resolution | |
| Water Body* | WBody | Binary water layer containing rivers, lakes and coastline at 10 m spatial resolution | |
| Elevation model | DEM | Digital elevation model at 10 m spatial resolution | The Danish Agency for Data Supply and Efficiency (SDFE) |
| Terrain slope | Slope | Rise and fall of the terrain surface in degree at 10 m spatial resolution | |
| Horizontal distance to WBody | HDis | Horizontal distance to nearest water body at 10 m spatial resolution | |
| Vertical distance to WBody | VDis | Vertical distance to nearest water body at 10 m spatial resolution | |
| PBM—winter condition | PBMw | Depth of water table for median winter condition simulated by a PBM at 100 m spatial resolution | DK-model |
| PBM—summer condition | PBMs | Depth of water table for median summer condition simulated by a PBM at 100 m spatial resolution |
Overview of the explanatory variables used to model the uppermost water table.
Categorical variables are indicated with an asterisk.
Gradient Boosting Decision Trees
We applied a new implementation of the gradient boosting decision tree (GBDT) algorithm, i.e., Cat Boost that was first developed by Yandex engineers in 2017 (Dorogush et al.,
For the purpose of simulating the depth of the uppermost water table using CatBoost we employ two different objective functions. First, the mean absolute error (MAE) is used to train the best estimate of the winter and summer condition. One model is trained for each season. The MAE is expressed as follows:
where sim is the simulated groundwater depth and obs the observed for a total of n training data. Besides simulating the best estimate, we are also interested in quantifying the uncertainty of the model. For this, we utilized quantile regression to define objective functions targeting specific quantiles of the distribution. This yields a probabilistic model that is not trained to estimate the conditional mean, but to estimate a defined quantile q of the distribution instead:
Setting q to 0.5 yields the same result as the MAE, but setting q to other values will give asymmetric weights to the residual depending on q and the overall sign of the error. Setting q to 0.1 will estimate the 10th percentile, by associating a weight of 0.9 to over predictions and a weight of 0.1 to under predictions. Thereby the 10th percentile can be approximated, meaning that the model will be trained to over predict 90% of the times. With the uncertainty analysis we intend to estimate the 68 and 95% confidence intervals which can be achieved by training four models, each model targeting a single quantile (e.g., 0.16 and 0.84 for the 68% confidence interval). The four models have to be trained individually for both, summer, and winter condition. We used CPU on a Windows machine (4 2.2 GHz Intel Xeon processors with 56 cores in total, 256 GB RAM) for training and prediction.
Results
High Resolution Groundwater Model
CatBoost regression was used to simulate the depth of the uppermost water table for a typical winter and summer condition at 10 m resolution over entire Denmark. Hyper parameters were first manually calibrated and afterwards automatically fine-tuned using a randomized search approach. Table 3 contains the nine hyper parameters included in the randomized search, a short description, and the tested values. The randomized hyper parameter search was conducted for the summer and the winter model using a 3-fold cross validation approach using 75% of the data. The 3-fold cross validation is the default of CatBoost's randomized hyper parameter search algorithm. The MAE, based on the 2,500 hyper parameter combination derived from the randomized search algorithm, varied just 4 cm. A common hyper parameter set that was among the top 2% for both, summer, and winter model, was selected for further modeling (Table 3). We found that the performance did not deteriorate for the test data (25%).
Table 3
| Hyperparameter | Description | Tested values | Optimized value |
|---|---|---|---|
| learning_rate | Reduces the gradient step during training | 0.05, 0.075, 0.1, 0.125, 0.15 | 0.05 |
| Depth | Number of levels in the decision trees | 8, 9, 10, 11, 12, 13 | 13 |
| l2_leaf_reg | The coefficient of the L2 regularization term of the loss function | 0, 2, 4, 6, 8, 10, 12 | 4 |
| Subsample | Random selection of training data for defining splits | 0.5, 0.6, 0.7, 0.8, 0.9, 1 | 1 |
| rsm | Random selection of covariates for defining splits | 0.5, 0.6, 0.7, 0.8, 0.9, 1 | 0.6 |
| random_strength | Randomness for selecting the optimal split | 0.5, 0.75, 1, 1.25, 1.5 | 0.5 |
| min_data_in_leaf | Minimum data in each leaf | 1, 5, 9, 13, 17, 21, 25 | 25 |
| bagging_temperature | Random weights to training data | 0, 0.5, 1, 1.5 | 1.5 |
CatBoost hyper parameter used in the randomized search (2,500 combinations).
The maximum number of decision trees, i.e., iterations in the gradient boosting, was set to 1,000.
Table 4 shows results obtained from a 4-fold cross validation test for the summer and winter models. For this, four CatBoost regression models were trained using 75% of the data for training and 25% of the data was held back for validation. Overall, little variation was found across the cross validation models, which indicates an overall robustness. The MAE, averaged across the four cross validation test, is 47 cm for both, summer, and winter model. This calculation is based on the entire training dataset. Only considering well data yields an average MAE of around 1.15 m. The increase in MAE is due to the fact that the additional observations, all with a depth of zero, are well-captured by the model and generally yield lower residuals. Based on the coefficient of determination over half of the variance in the well data is accounted for by the models and over 70% of the variance of the entire training dataset. The 19,000 lakes are represented with a MAE of below 5 cm, whereas additional observations along the stream network and coastline possess a MAE of around 10 cm.
Table 4
| Performance | Summer model | Winter model | |||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|
| cv1 | cv2 | cv3 | cv4 | Mean | cv1 | cv2 | cv3 | cv4 | Mean | ||
| All data | MAE | 0.48 | 0.48 | 0.48 | 0.46 | 0.47 | 0.47 | 0.47 | 0.47 | 0.45 | 0.47 |
| RMSE | 1.24 | 1.21 | 1.23 | 1.19 | 1.22 | 1.23 | 1.21 | 1.21 | 1.18 | 1.20 | |
| R2 | 0.75 | 0.76 | 0.75 | 0.77 | 0.76 | 0.72 | 0.72 | 0.71 | 0.73 | 0.72 | |
| Only well observations | MAE | 1.17 | 1.18 | 1.19 | 1.15 | 1.17 | 1.16 | 1.17 | 1.18 | 1.13 | 1.16 |
| RMSE | 1.96 | 1.92 | 1.95 | 1.91 | 1.93 | 1.94 | 1.92 | 1.93 | 1.89 | 1.92 | |
| R2 | 0.51 | 0.53 | 0.51 | 0.54 | 0.52 | 0.51 | 0.52 | 0.51 | 0.53 | 0.51 | |
| Only lakes | MAE | 0.03 | 0.04 | 0.04 | 0.03 | 0.04 | 0.03 | 0.03 | 0.03 | 0.03 | 0.03 |
| RMSE | 0.30 | 0.27 | 0.33 | 0.27 | 0.29 | 0.26 | 0.25 | 0.28 | 0.23 | 0.25 | |
| R2 | |||||||||||
| Only river and coastline | MAE | 0.13 | 0.11 | 0.13 | 0.12 | 0.12 | 0.09 | 0.08 | 0.11 | 0.10 | 0.10 |
| RMSE | 0.31 | 0.31 | 0.35 | 0.32 | 0.32 | 0.27 | 0.22 | 0.28 | 0.27 | 0.26 | |
| R2 | |||||||||||
Results for a 4-fold cross validation (cv) applied on the summer and winter model.
Performance is quantified by means of a mean absolute error (MAE), root mean squared error (RMSE) and coefficient of determination (R2). Results are given for different subset of the data (all data, only wells, only lakes and only river and coastline).
Figure 4 depicts the simulated national maps of the depth of the uppermost water table at 10 m spatial resolution for summer and winter. The maps offer much detail revealing the interplay of geology, topography, and waterbody proximity. The overall regional patterns show resemblance with the simulations results of the physically-based model at 100 m, which were used as covariate in the 10 m ML models (Figure 2). However, disagreements between the 100 m physically-based model and the 10 m ML model are also present. The former simulates a shallower water table in the Eastern part of Demark where moraine clay is the dominant lithology. This is especially the case for the winter condition. Further, the 100 m physically-based model simulates a homogeneously deep water table in the sandy areas whereas the ML model results in more heterogeneity. The disagreements may be explained by the differences in spatial resolution, but also with the vertical discretization of the computational layers in the physically-based model. Sandy layers quickly run dry which results in a deep water table and clayey layers hold water which yields a very shallow water table. The ML based median seasonality is calculated as the difference between the summer and the winter maps. The amplitude is generally lowest close to streams and lakes and highest in the center of Denmark where topography is high and a permeable subsurface is present. Similar to the physically-based simulated amplitude, the ML derived amplitude has intermediate values around 0.5 m (yellow category) in the Western part of Denmark for the sandy outwash plains. Based on both approaches, the amplitude for the moraine clay settings is around 1 m (green category) which is mainly driven by the very shallow water table during winter.
Figure 4

Simulation results for the depth of the uppermost water table obtained from the summer and winter model at 10 m spatial resolution, left and center panel, respectively. The right panel depicts the simulated amplitude calculated as summer depth minus winter depth.
Figure 5 presents the same data as shown in Figure 4, just for a zoom section of ~15 km2. Here, the imprint of the stream network as well as the many lakes, become apparent as areas with a very shallow depth of the uppermost water table (< 0.5 m). The difference between the drier summer and the wetter winter is clearly notable and the difference between the two models is plotted as the seasonality. The largest amplitude is present in areas with high topography and the lowest amplitude is found along the river network and lakes.
Figure 5

Same data as presented in Figure 4, but zoomed to an ~15 km2 area located on the island of Fyn. The zoom location is indicated in the overview map (bottom left).
Covariate Importance
The quantified importance of each input feature for the summer and winter model is presented in Figure 6 for the entire training dataset and for a subset containing only the well data. Based on a trained model, CatBoost calculates the “prediction value change” to quantify how much on average the prediction changes if the covariate value changes. The average changes that represent the importance of a given covariate are normalized to add up to 100. The vertical distance to the closest water body (VDis) clearly stands out as the most importance covariate in both models. At locations with VDis close to zero, the water table is typically also close to the surface. However, large VDis values, which indicate small scale topographical variations often result in a deeper water level, as the shallow groundwater does not follow the topography in such settings. The most important categorical covariate is the landscape type classification (LType) which contains 13 landscape classes, such as moraine, marine plains, outwash plains, and others. There is little difference between winter and summer model, which indicates robustness between the two models. Differences between covariate importance with respect to the entire dataset and only well data conveys that the horizontal distance to the closest water body (HDis) is mostly relevant to the additional observations, namely lakes, rivers, and coast, which is expectable as these are characterized with a distance of zero. The physically-based model (PBM) stands out as second most important covariate when only considering the well-training dataset. This underlines that the physically-based model can provide the ML model with a meaningful information, despite the differences in spatial resolution, since both model the same variable.
Figure 6

Covariate importance quantified as prediction change in % calculated for summer and winter model for two subsets of training data (all data and well data only). The covariate abbreviations are explained in Table 2.
Uncertainty Analysis
Quantile regression has been applied in order to estimate the uncertainties associated with the summer and winter models. For both models, four quantile models have been trained using the same hyper parameters as applied previously. The quantile models were set with q = 0.16 and 0.84 for the 68% confidence interval and q = 0.025 and 0.975 for the 95% confidence interval. Results are shown in Figure 7 and were calculated on the basis of the same 4-fold cross validation test as presented in Table 4. Figure 7 only contains wells and data are sorted with respect to the simulated groundwater depth and increase alongside an increasing x-axis. The observations, 13,047 in total, are plotted as a density plot and an overall good agreement between model and observations can be attested to both, summer, and winter. The blue envelope plots indicate the two confidence intervals and following the quantile regression definition, 5% of the observations are expected to be outside the light blue envelope (95% confidence) and 32% are expected to be outside the dark blue envelope (68% confidence). For both, the summer, and winter model, the uncertainty increases with depth, which gets supported by the large spread of observations for larger depths. Groundwater depths above 6 m have an uncertainty of around 4 and 8 m following the 68 and 95% confidence intervals.
Figure 7

Results of the 4-fold cross validation test for summer model (top) and winter model (bottom). The simulated data are sorted, and the observations are plotted as point density. Only well observations are shown in the plot. Uncertainty bands are included for two confidence intervals.
Based on the 95% confidence interval, which reflects the spread of ~±2 standard deviations around the mean, a map of the standard deviation at 10 m spatial resolution has been calculated. Overall, the relationship between the standard deviation and the water table depth is nearly linear and the coefficient of determination, expressed as the average standard deviation over the average water table depth is 0.58. This indicates that the uncertainty is roughly half of its water table depth value. Figure 8 exemplifies the relationship between the simulated depth of water table and the associated uncertainty (standard deviation) for two soil types. In moraine clay soils, the depth of the uppermost water table is typically in the top few meters, whereas moraine sand soils are characterized with deeper water tables. The coefficient of variation for the two soil types is 0.74 for moraine clay and 0.44 for moraine sand which underlines that uncertainty is larger for moraine clay soils as opposed to moraine sand soils. This can be expected given the distinct geological heterogeneity in the moraine clay soils.
Figure 8

Density scatter plot showing the simulated depth of the groundwater for moraine clay and sand soil against the associated standard deviation (std.) for the winter model. Results are shown for two of the nine soil types (Figure 1).
Discussion
Training Dataset
The winter trainings dataset comprises all DJF water level observations and JJA represents summer conditions. The median was calculated in the case of multiple observations per well per season. This approach introduces uncertainties since it rules out inter-annual variability and further, variability within the summer and winter seasons is also ignored. This compromise was accepted in order to obtain a large trainings dataset of groundwater observations. Large-scale groundwater datasets are typically very heterogenous, which is especially the case for the temporal dimension. The spatial density of wells may be high, but the temporal resolution of water level observations poses challenges to ML applications. Heterogeneity originates from varying frequencies and periods of observations. As an example, around 66% of the 104,000 wells in Denmark have only a single observation in the period of 1990–2019. In order to capitalize on under sampled shallow wells in a big data context, this study brings forward a knowledge-guided ML framework that employs the simulated seasonality from a physically-based groundwater flow model. With this augmentation strategy, the training dataset could be expanded significantly from 1,378 shallow wells with both, summer and winter observations, to 13,047 with either gap-filled summer or winter observation. Previously, Koch et al. (
This study utilized a comprehensive set of 15 covariates to predict water table variability. This selection has been guided by a previous Danish study by Koch et al. (
Temporal Resolution
This study implements a simplified temporal dimension of the groundwater table, by simulating two seasons, namely a typical winter and summer condition. These temporal snapshots are of course a prude simplification of a variable that is known to possess a distinct temporal variability. However, previous studies that apply ML to model the spatial variability of the water table have focused on a single time step; i.e., Koch et al. (
Machine Learning Model
Overall, the accuracy of the trained ML model, which was quantified by means of a 4-fold cross validation test, was very satisfying and in the range of what is generally considered very acceptable in groundwater flow modeling (Henriksen et al.,
Results from a physically-based groundwater model were incorporated in the ML model as explanatory variable and it was shown that the importance of the 100 m groundwater model was the second highest with respect to the well observations. This is satisfying as it issues consistency between the two modeling approaches and underpins that the physically-based model can guide the ML model. It has to be noted that the wells used for training the ML model were also used to calibrate the groundwater model. However, we believe that reusing of data does not have any implications and in the end, a resemblance between the groundwater model and the ML model is also desirable. The most important explanatory variable is the vertical distance to the nearest waterbody, which indicates that the shallow groundwater system is decoupled from small scale topographical variability. The applied knowledge-guided ML framework is novel for the groundwater domain and may be applicable to other regions where results from a physically-based groundwater model are available. Adding simulation results from hydrological models into ML models has also been successfully applied for the surface water system, i.e., predicting streamflow (Konapala et al.,
We envision that the developed water table map can support the planning of climate resilient infrastructure design, with special focus on flooding or the planning of rewetting of lowlands to reduce greenhouse gas emissions from drained peatland soils. In contrast to conventional PBMs, the ML proposed herein cannot be used to run modeling scenarios (e.g., climate change or water management), but our results can be used as a national scale screening tool to identify areas where a PBM can subsequently be applied to test relevant scenarios at high spatial resolution. The ML based results reflect the current climate conditions and if the impact of climate change of the uppermost water table requires investigation, a PBM should be applied instead.
Uncertainty Analysis
The quantile regression technique was used to estimate uncertainty bounds. Uncertainty was quantified corresponding to the 68 and 95% confidence intervals of the simulated water table depth. We found that uncertainty generally increases with depth and that the average coefficient of variation is 0.58, indicating that the standard deviation is more than half of the water table depth, but depends on the geological setting with higher uncertainty for the complex moraine clay soils. One known drawback of quantile regression is that it requires individual training of each selected quantile, which can result in an invalid distribution, meaning that the estimated quantiles are not monotonically increasing (Bondell et al.,
Following the discussion provided by Vaysse and Lagacherie (
Conclusions
The study features a knowledge-guided ML model of the uppermost water level at 10 m spatial resolution at national scale for Denmark. We have applied the gradient bosting decision tree implementation of CatBoost to model a typical summer and winter condition. The associated uncertainties were estimated using quantile regression techniques. Predicting water levels at 430 million grids is unfeasible with conventional dynamic physically-based groundwater flow models, which highlights the benefits of using alternative ML modeling approaches instead, to reach unprecedented spatial detail. We draw the following main conclusions from our work:
The applied high resolution ML model could predict water table variability with high accuracy. The MAE of well observations was around 115 and 50 cm taking also the groundwater proxy observations (lakes, rivers, and coastline) into consideration.
A physically-based groundwater flow model was successfully incorporated into the model building to (1) select wells that are representative for the shallow groundwater system, (2) augment training data by accounting for simulated water level seasonality, and (3) extend the list of explanatory variables. This forms a novel application of knowledge guided ML for the shallow groundwater domain.
The water table depth was simulated for two temporal snapshots (typical summer and winter). Future ML research in the groundwater domain must focus on modeling the full spatio-temporal variability of water table depth.
Publisher's Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Statements
Data availability statement
Publicly available datasets were analyzed in this study. The water table observations are available via the Jupiter database: http://data.geus.dk/geusmap. The raw data and scripts supporting the results and conclusions of this article will be made available by the corresponding author, without undue reservation. The final summer and winter water table depth maps are visualized and freely available viahttps://hipdata.dk/.
Author contributions
JG: data curation. JG, RS, LT, SS, and HH: conceptualization, result interpretation, and quality control. JK: code development, study design, writing original draft, and visualization. All authors have read, edited, and agreed to the published version of the manuscript.
Funding
The study has been funded by the Danish Digitalization Strategy (FODS) through the HIP (Hydrological Information- and Prognosis-system) project related to FODS 6.1.
Acknowledgments
The authors wish to acknowledge contributions from GEUS colleagues, Søren Kragh, Maria Ondracek, Michael van Til, and Annesofie Jakobsen who all contributed to the development of the 100 m implementation of the DK-model. Further, the authors want to thank the Danish Agency for Data Supply and Efficiency (SDFE) for a good collaboration and support throughout the project.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
References
1
AdhikariK.KheirR. B.GreveM. B.BøcherP. K.MaloneB. P.MinasnyB.et al. (2013). High-resolution 3-D mapping of soil texture in Denmark. Soil Sci. Soc. Am. J.77, 860–876. 10.2136/sssaj2012.0275
2
BechtoldM.SchlafferS.TiemeyerB.De LannoyG. (2018). Inferring water table depth dynamics from ENVISAT-ASAR C-band backscatter over a range of peatlands from deeply-drained to natural conditions. Remote Sens.10:536. 10.3390/rs10040536
3
BechtoldM.TiemeyerB.LaggnerA.LeppeltT.FrahmE.BeltingS. (2014). Large-scale regionalization of water table depth in peatlands optimized for greenhouse gas emission upscaling. Hydrol. Earth Syst. Sci.18, 3319–3339. 10.5194/hess-18-3319-2014
4
BevenK. J.ClokeH. L. (2012). Comment on hyperresolution global land surface modeling: meeting a grand challenge for monitoring Earth's terrestrial water. Water Resour. Res.48:52. 10.1029/2011WR010982
5
BondellH. D.ReichB. J.WangH. (2010). Noncrossing quantile regression curve estimation. Biometrika97, 825–838. 10.1093/biomet/asq048
6
Breuning-MadsenH.JensenN. H. (1992). Pedological regional variations in well-drained soils, Denmark. Geogr. Tidsskr. J. Geogr.92, 61–69. 10.1080/00167223.1992.10649316
7
BrickerS. H.BanksV. J.GalikG.TapeteD.JonesR. (2017). Accounting for groundwater in future city visions. Land Use Policy69, 618–630. 10.1016/j.landusepol.2017.09.018
8
ClarkM. P.BierkensM. F. P.SamaniegoL.WoodsR. A.UijlenhoetR.BennettK. E.et al. (2017). The evolution of process-based hydrologic models: historical challenges and the collective quest for physical realism. Hydrol. Earth Syst. Sci.21, 3427–3440. 10.5194/hess-21-3427-2017
9
CloseM. E.AbrahamP.HumphriesB.LilburneL.CuthillT.WilsonS. (2016). Predicting groundwater redox status on a regional scale using linear discriminant analysis. J. Contam. Hydrol.191, 19–32. 10.1016/j.jconhyd.2016.04.006
10
DorogushA. V.ErshovV.GulinA. (2018). CatBoost: gradient boosting with categorical features support. arXiv arxiv: 1810.11363.
11
EricksonM. L.ElliottS. M.BrownC. J.StackelbergP. E.RansomK. M.ReddyJ. E. (2021). Machine learning predicted redox conditions in the glacial aquifer system, northern continental United States. Water Resour. Res.57:e2020WR028207. 10.1029/2020WR028207
12
FanJ.YueW.WuL.ZhangF.CaiH.WangX.et al. (2018). Evaluation of SVM, ELM and four tree-based ensemble models for predicting daily reference evapotranspiration using limited meteorological data in different climates of China. Agric. For. Meteorol.263, 225–241. 10.1016/j.agrformet.2018.08.019
13
FanY.LiH.Miguez-MachoG. (2013). Global patterns of groundwater table depth. Science339, 940–943. 10.1126/science.1229881
14
FienenM. N.MastersonJ. P.PlantN. G.GutierrezB. T.ThielerE. R. (2013). Bridging groundwater models and decision support with a Bayesian network. Water Resour. Res.49, 6459–6473. 10.1002/wrcr.20496
15
FriedmanJ. H. (2001). Greedy function approximation: a gradient boosting machine. Ann. Stat.29, 1189–1232. 10.1214/aos/1013203451
16
GeorganosS.GrippaT.VanhuysseS.LennertM.ShimoniM.WolffE. (2018). Very high resolution object-based land use-land cover urban classification using extreme gradient boosting. IEEE Geosci. Remote Sens. Lett.15, 607–611. 10.1109/LGRS.2018.2803259
17
GleesonT.BefusK. M.JasechkoS.LuijendijkE.CardenasM. B. (2016). The global volume and distribution of modern groundwater. Nat. Geosci.9, 161–167. 10.1038/ngeo2590
18
GuoH.ZhuangX.RabczukT. (2020). Stochastic Analysis of Heterogeneous Porous Material with Modified Neural Architecture Search (NAS) Based Physics-Informed Neural Networks Using Transfer Learning. Available online at: http://arxiv.org/abs/2010.12344 (accessed June 8, 2021).
19
GuzmanS. M.PazJ. O.TagertM. L. M. (2017). The use of NARX neural networks to forecast daily groundwater levels. Water Resour. Manag.31, 1591–1603. 10.1007/s11269-017-1598-5
20
HancockJ. T.KhoshgoftaarT. M. (2020). CatBoost for big data: an interdisciplinary review. J. Big Data. 7, 1–45. 10.1186/s40537-020-00369-8
21
HenglT.MillerM. A. E.KriŽanJ.ShepherdK. D.SilaA.KilibardaM.et al. (2021). African soil properties and nutrients mapped at 30 m spatial resolution using two-scale ensemble machine learning. Sci. Rep.11, 1–18. 10.1038/s41598-021-85639-y
22
HenglT.NussbaumM.WrightM. N.HeuvelinkG. B. M.GrälerB. (2018). Random forest as a generic framework for predictive modeling of spatial and spatio-temporal variables. PeerJ.6:e5518. 10.7717/peerj.5518
23
HenriksenH. J.KraghS. J.GodtfredsenJ.OndracekM.van ThilM. J.JakobsenA.et al. (2020). Dokumentationsrapport vedr. modelleverancer til Hydrologisk Informations- og Prognosesystem (in Danish). Copenhagen: GEUS.
24
HenriksenH. J.TroldborgL.HøjbergA. L.RefsgaardJ. C. (2008). Assessment of exploitable groundwater resources of Denmark by use of ensemble resource indicators and a numerical groundwater-surface water model. J. Hydrol.348, 224–240. 10.1016/j.jhydrol.2007.09.056
25
HenriksenH. J.TroldborgL.NyegaardP.SonnenborgT. O.RefsgaardJ. C.MadsenB. (2003). Methodology for construction, calibration and validation of a national hydrological model for Denmark. J. Hydrol.280, 52–71. 10.1016/S0022-1694(03)00186-0
26
HøjbergA. L.TroldborgL.StisenS.ChristensenB. B. S.HenriksenH. J. (2013). Stakeholder driven update and improvement of a national water resources model. Environ. Model. Softw.40, 202–213. 10.1016/j.envsoft.2012.09.010
27
HuangG.WuL.MaX.ZhangW.FanJ.YuX.et al. (2019). Evaluation of CatBoost method for prediction of reference evapotranspiration in humid regions. J. Hydrol.574, 1029–1041. 10.1016/j.jhydrol.2019.04.085
28
KahlownM. A.AshrafM.Zia-Ul-Haq (2005). Effect of shallow groundwater table on crop water requirements and crop yields. Agric. Water Manag.76, 24–35. 10.1016/j.agwat.2005.01.005
29
KarlssonI. B.SonnenborgT. O.RefsgaardJ. C.TrolleD.BørgesenC. D.OlesenJ. E.et al. (2016). Combined effects of climate models, hydrological model structures and land use scenarios on hydrological impacts of climate change. J. Hydrol.535, 301–317. 10.1016/j.jhydrol.2016.01.069
30
KidmoseJ.RefsgaardJ. C.TroldborgL.SeabyL. P.EscrivàM. M. (2013). Climate change impact on groundwater levels: Ensemble modelling of extreme values. Hydrol. Earth Syst. Sci.17, 1619–1634. 10.5194/hess-17-1619-2013
31
KochJ.BergerH.HenriksenH. J.SonnenborgT. O. (2019a). Modelling of the shallow water table at high spatial resolution using random forests. Hydrol. Earth Syst. Sci.23, 4603–4619. 10.5194/hess-23-4603-2019
32
KochJ.StisenS.RefsgaardJ. C.ErnstsenV.JakobsenP. R.HøjbergA. L. (2019b). Modeling depth of the redox interface at high resolution at national scale using random forest and residual gaussian simulation. Water Resour. Res.55, 1451–1469. 10.1029/2018WR023939
33
KonapalaG.KaoS. C.PainterS. L.LuD. (2020). Machine learning assisted hybrid models can improve streamflow simulation in diverse catchments across the conterminous US. Environ. Res. Lett.15:104022. 10.1088/1748-9326/aba927
34
KraftB.JungM.KörnerM.ReichsteinM. (2020). Hybrid modeling: Fusion of a deep approach and physics-based model for global hydrological modeling. Int. Archiv. Photogram. Rem. Sens. Spat. Inform. Sci.43, 1537–1544. 10.5194/isprs-archives-XLIII-B2-2020-1537-2020
35
LarsenM. A. D.ChristensenJ. H.DrewsM.ButtsM. B.RefsgaardJ. C. (2016). Local control on precipitation in a fully coupled climate-hydrology model. Sci. Rep.6, 1–9. 10.1038/srep22927
36
LevinG.BlemmerM. K.NielsenM. R. (2012). Basemap: Technical Documentation of a Model for Elaboration of a Land-Use and Land-Cover Map for Denmark. Aarhus: Aarhus University.
37
López LópezP.VerkadeJ. S.WeertsA. H.SolomatineD. P. (2014). Alternative configurations of quantile regression for estimating predictive uncertainty in water level forecasts for the upper Severn River: a comparison. Hydrol. Earth Syst. Sci.18, 3411–3428. 10.5194/hess-18-3411-2014
38
MacDonaldD.DixonA.NewellA.HallawaysA. (2012). Groundwater flooding within an urbanised flood plain. J. Flood Risk Manag.5, 68–80. 10.1111/j.1753-318X.2011.01127.x
39
MaxwellR. M.CondonL. E. (2016). Connections between groundwater flow and transpiration partitioning. Science353, 377–380. 10.1126/science.aaf7891
40
MøllerA. B.BeucherA.IversenB. V.GreveM. H. (2018). Predicting artificially drained areas by means of a selective model ensemble. Geoderma320, 30–42. 10.1016/j.geoderma.2018.01.018
41
MøllerA. B.IversenB. V.BeucherA.GreveM. H. (2017). Prediction of soil drainage classes in Denmark by means of decision tree classification. Geoderma352, 314–329. 10.1016/j.geoderma.2017.10.015
42
NearingG. S.KratzertF.SampsonA. K.PelissierC. S.KlotzD.FrameJ. M.et al. (2020). What role does hydrological science play in the age of machine learning?Water Resour. Res.57:e2020WR028091. 10.1029/2020WR028091
43
ProkhorenkovaL.GusevG.VorobevA.DorogushA. V.GulinA. (2018). Catboost: Unbiased boosting with categorical features. arXiv [Preprint] arXiv:1706.09516.
44
RajaeeT.EbrahimiH.NouraniV. (2019). A review of the artificial intelligence methods in groundwater level modeling. J. Hydrol.572, 336–351. 10.1016/j.jhydrol.2018.12.037
45
ReadJ. S.JiaX.WillardJ.ApplingA. P.ZwartJ. A.OliverS. K.et al. (2019). Process-guided deep learning predictions of lake water temperature. Water Resour. Res.55, 9173–9190. 10.1029/2019WR024922
46
ReichsteinM.Camps-VallsG.StevensB.JungM.DenzlerJ.CarvalhaisN.et al. (2019). Deep learning and process understanding for data-driven Earth system science. Nature566, 195–204. 10.1038/s41586-019-0912-1
47
ShenC. (2018). A transdisciplinary review of deep learning research and its relevance for water resources scientists. Water Resour. Res.54, 8558–8593. 10.1029/2018WR022643
48
StisenS.OndracekM.TroldborgL.SchneiderR. J. M.van ThilM. J. (2019). National vandressource model (in Danish). Modelopstilling Og Kalibrering Af DK-model 2019. GEUS Rapp. 2019/31, Copenhagen.
49
SunY.WendiD.KimD. E.LiongS. Y. (2016). Technical note: application of artificial neural networks in groundwater table forecasting-a case study in a Singapore swamp forest. Hydrol. Earth Syst. Sci.20, 1405–1412. 10.5194/hess-20-1405-2016
50
TaylorR. G.ScanlonB.DöllP.RodellM.Van BeekR.WadaY.et al. (2013). Ground water and climate change. Nat. Clim. Chang.3, 322–329. 10.1038/nclimate1744
51
TesorieroA. J.TerziottiS.AbramsD. B. (2015). Predicting redox conditions in groundwater at a regional scale. Environ. Sci. Technol.49, 9657–9664. 10.1021/acs.est.5b01869
52
TiemeyerB.Albiac BorrazE.AugustinJ.BechtoldM.BeetzS.BeyerC.et al. (2016). High emissions of greenhouse gases from grasslands on peat and other organic soils. Glob. Chang. Biol.22, 4134–4149. 10.1111/gcb.13303
53
TiemeyerB.FreibauerA.BorrazE. A.AugustinJ.BechtoldM.BeetzS.et al. (2020). A new methodology for organic soils in national greenhouse gas inventories: data synthesis, derivation and application. Ecol. Indic.109:105838. 10.1016/j.ecolind.2019.105838
54
TyralisH.PapacharalampousG.BurnetasA.LangousisA. (2019a). Hydrological post-processing using stacked generalization of quantile regression algorithms: large-scale application over CONUS. J. Hydrol.577:123957. 10.1016/j.jhydrol.2019.123957
55
TyralisH.PapacharalampousG.LangousisA. (2019b). A brief review of random forests for water scientists and practitioners and their recent history in water resources. Water11:910. 10.3390/w11050910
56
VaysseK.LagacherieP. (2017). Using quantile regression forest to estimate uncertainty of digital soil mapping products. Geoderma291, 55–64. 10.1016/j.geoderma.2016.12.017
57
WunschA.LieschT.BrodaS. (2021). Groundwater level forecasting with artificial neural networks: a comparison of long short-term memory (LSTM), convolutional neural networks (CNNs), and non-linear autoregressive networks with exogenous input (NARX). Hydrol. Earth Syst. Sci.25, 1671–1687. 10.5194/hess-25-1671-2021
58
ZipperS. C.SoyluM. E.BoothE. G.LoheideS. P. (2015). Untangling the effects of shallow groundwater and soil texture as drivers of subfield-scale yield variability. Water Resour. Res.51, 6338–6358. 10.1002/2015WR017522
Summary
Keywords
water table depth, machine learning, high resolution, CatBoost, quantile regression
Citation
Koch J, Gotfredsen J, Schneider R, Troldborg L, Stisen S and Henriksen HJ (2021) High Resolution Water Table Modeling of the Shallow Groundwater Using a Knowledge-Guided Gradient Boosting Decision Tree Model. Front. Water 3:701726. doi: 10.3389/frwa.2021.701726
Received
28 April 2021
Accepted
28 June 2021
Published
01 September 2021
Volume
3 - 2021
Edited by
Yoram Rubin, University of California, Berkeley, United States
Reviewed by
Jie Niu, Jinan University, China; Aldo Fiori, Roma Tre University, Italy; Jinsong Chen, Lawrence Berkeley National Laboratory, United States
Updates

Check for updates
Copyright
© 2021 Koch, Gotfredsen, Schneider, Troldborg, Stisen and Henriksen.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Julian Koch juko@geus.dk
This article was submitted to Water and Hydrocomplexity, a section of the journal Frontiers in Water
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.