2018 WR 023939
2018 WR 023939
2018 WR 023939
10.1029/2018WR023939
Resolution at National Scale Using Random
Special Section: Forest and Residual Gaussian Simulation
Big Data & Machine Learning
in Water Sciences: Recent Julian Koch1 , Simon Stisen1 , Jens C. Refsgaard1 , Vibeke Ernstsen2,
Progress and Their Use in
Advancing Science Peter R. Jakobsen2, and Anker L. Højberg1
1
Department of Hydrology, Geological Survey of Denmark and Greenland, Copenhagen, Denmark, 2Department of
Key Points: Geochemistry, Geological Survey of Denmark and Greenland, Copenhagen, Denmark
• A random forest model was trained
to predict subsurface redox
conditions at 100‐m resolution Abstract The management of water resources needs robust methods to efficiently reduce nitrate loads.
across Denmark
• Residuals were simulated using
Knowledge on where natural denitrification takes place in the subsurface is thereby essential. Nitrate is
geostatistics to estimate uncertainty naturally reduced in anoxic environments and high‐resolution information of the redox interface, that is, the
and to generate multiple plausible depth of the uppermost reduced zone is crucial to understand the variability of the denitrification
realizations
• A framework to quickly integrate
potential. In this study we explore the opportunity to use random forest (RF) regression to model redox
newly acquired field data into local depth across Denmark at 100‐m resolution based on ~13,000 boreholes as training data. We highlight the
redox models was outlined and importance of expert knowledge to guide the RF model in areas where our conceptual understanding is not
exemplified for a field site
represented correctly in the training data set by addition of artificial observations. We apply random
forest regression kriging in which sequential Gaussian simulation models the RF residuals. The RF model
reaches a R2 score of 0.48 for an independent validation test. Including sequential Gaussian simulation
Correspondence to: honors observations through local conditioning, and the spread of 800 realizations can be utilized to map
J. Koch,
[email protected] uncertainty. Emphasis is put on adequate handling of nonstationarities in variance and spatial correlation of
the RF residuals. The RF residuals show no spatial correlation for large parts of the modeling domain,
and a local variance scaling method is applied to account for the nonstationary variance. Moreover, we
Citation:
Koch, J., Stisen, S., Refsgaard, J. C., present and exemplify a framework where newly acquired field data can easily be integrated into random
Ernstsen, V., Jakobsen, P. R., & forest regression kriging to quickly update local models.
Højberg, A. L. (2019). Modeling depth
of the redox interface at high resolution Plain Language Summary Nutrients, such as nitrogen in form of nitrate are essential for plant
at national scale using random forest growth. Despite their benefits, nutrient surplus can cause adverse health and ecological effects. In fact,
and residual Gaussian simulation.
Water Resources Research, 55, nitrate leaching from agricultural fields is one of the major water resources management challenges in
1451–1469. https://doi.org/10.1029/ today's agricultural landscapes. During the transport through the subsurface, nitrate can potentially be
2018WR023939 degraded in anoxic layers below the redox interface. This is referred to as denitrification and its potential
varies spatially depending on the local landscape. The location of the redox interface is thus essential
Received 21 AUG 2018
Accepted 23 JAN 2019 knowledge toward a spatially differentiated water resources management. We applied a machine learning
Accepted article online 29 JAN 2019 method to model the redox depth at 100m spatial resolution across Denmark, based on redox data obtained
Published online 20 FEB 2019 from ~13,000 boreholes and environmental variables that explain redox variability. We highlight the
importance of expert knowledge in guiding the machine learning model by adding artificial observations in
underrepresented landscapes. As an add‐on, a geostatistical model is used to honor local data at grid scale
and quantify the spatial uncertainty elsewhere. Lastly, we outline and exemplify a framework that allows an
easy integration of newly acquired field data for updating local models of the depth to the redox interface.
1. Introduction
The transport of excess nitrate from agricultural fields to surface waterbodies, groundwater, and marine
waters constitutes a significant problem for environmental and human health across the globe (Fenton
et al., 2009; Hinkle & Tesoriero, 2014). The widespread nutrient pollution has been accelerated by the
twentieth‐century agricultural intensification and threatens secure supplies of drinking water and mainte-
nance of ecosystem quality (Erisman et al., 2013).
Natural removal of nitrate occurs during the transport from source areas by different natural biogeochemical
processes. These processes include denitrification, where nitrate is reduced to N2 gas through a microbially
©2019. American Geophysical Union. facilitated process, sedimentation, or the transformation of nitrate into other N‐compounds (Højberg et al.,
All Rights Reserved. 2017; Jessen et al., 2017; Thayalakumaran et al., 2008). Following the Danish Nature Agency (2016), nitrate
loads in Denmark need to be further reduced by up to 70% in some catchments in order to meet require-
ments of the EU Water Framework Directive. Previously, researchers in Denmark have investigated the
environmental and economic gains of spatially targeted measures that are, opposed to uniform regulations,
spatially differentiated with respect to the natural denitrification potential (Hashemi et al., 2018; Jacobsen &
Hansen, 2016). This is in line with efforts from other countries such as Germany (Kuhr et al., 2013) or the
United Kingdom (Bateman et al., 2013). Hansen, Gunderman, et al. (2014) detailed that in order to support
spatially differentiated measures, integrated hydrological modeling can be used as a tool to gain insights into
the spatial variability of the denitrification potential at catchment scale. In line with McMahon et al. (2008)
and Hansen et al. (2009), this potential can be estimated based on combined knowledge on subsurface flow-
paths and the depth to the redox interface, that is, the location of the uppermost reduced (anaerobic) zone,
below which denitrification takes place. Recent advances on how to model redox conditions in the subsur-
face were highlighted by Close et al. (2016), Tesoriero et al. (2015), and Hansen, Christensen, et al. (2014). In
the former, Close et al. (2016) simulated categorical redox conditions, that is, oxidized, mixed, and reduced,
for two regions in New Zealand using a multivariate linear model. They used groundwater chemistry data to
assign redox conditions following the approach of McMahon and Chapelle (2008). Tesoriero et al. (2015)
mapped the probability of encountering oxic groundwater at various depths for a regional scale watershed
in the United States using logistic regression. In the latter, Hansen, Christensen, et al. (2014) studied the dis-
tinct spatial heterogeneity of the redox depth and brought forward the concept of a physically based model
that inferred the redox variability from the spatial pattern of groundwater recharge. In their work, local data
on the redox depth were derived from boreholes where a change in sediment color indicates the transition
from oxidized (brown) to reduced (gray) conditions (Ernstsen et al., 1998; Hendry et al., 1984). Based on
the description of sediment color in boreholes, supplementary geological data, and expert knowledge, a
national map of the redox depth across Denmark has been published by Ernstsen et al. (2010). The map
was developed at 1‐km spatial resolution and identified large‐scale variations but was originally not designed
for mapping local‐scale variations. Furthermore, quantifying the uncertainty in the manually interpreted
map is not possible. To support a successful implementation of spatially differentiated measures, data are
needed at a finer spatial scale, and uncertainty related to the depth of the redox interface with its implica-
tions on nitrate transport and denitrification must be quantified.
This study highlights the applicability of machine learning, namely, the random forest (RF) algorithm
(Breiman, 2001), to model depth to the redox interface at 100‐m spatial resolution. Capturing the distinct
heterogeneity of the redox depth requires a flexible modeling technique that is facilitated by data‐driven
models, such as RF. This marks the great potential of RF to model redox conditions at high spatial resolution.
The underlying concept of RF is built on an ensemble of decision trees that can be used for both, regression
and classification tasks. A related study carried out by Tesoriero et al. (2017) used RF as a classification tool
to model the probability of groundwater having a nitrate concentration above a certain threshold through-
out a large watershed. Further, nitrate concentrations were modeled by multiple studies applying RF or clo-
sely related machine learning techniques (Nolan et al., 2015; Ransom et al., 2017; Rodriguez‐Galiano et al.,
2014). Within a broader hydrological context, the applicability of RF and similar machine learning
approaches have already been proven versatile by several studies. For instance, Bachmair and Weiler
(2012) used RF to investigate drivers of shallow subsurface flow at a hillslope and could identify a complex
interaction of soil properties and topography that explained the observed variability. Baudron et al. (2013)
studied groundwater samples by using RF to classify them based on their geochemical characteristics with
the aim of tracing them back to their origin in a complex multilayer aquifer. Soil drainage conditions at high
spatial resolution were classified by Møller et al. (2017) with the help of decision trees. Very recent advances
in machine learning in the hydrological sciences promote deep learning, a machine learning method based
on artificial neural networks, as a versatile scientific tool (Shen, 2018). While the adaption of deep learning
in hydrology is slow, it already advanced related fields, such as soil sciences (Behrens et al., 2018). The
increase in machine learning applications in hydrology and related fields is facilitated by readily available
programming packages and a continual growth in availability, detail, and wealth of environmental data.
Hybrid models, where machine learning is linked to geostatistical techniques, are popular tools for the spa-
tial modeling of environmental variables (Kanevski et al., 2004; Li et al., 2011; 2017), but its applicability in
water resources research has not yet been tested. We utilized the random forest regression kriging (RFRK)
approach, as outlined by Hengl et al. (2015), where an RF model constitutes the trend and ordinary
kriging is utilized to predict the RF residuals. The clear benefit of geostatistical techniques such as kriging is
that uncertainties are estimated at grid scale (Hengl et al., 2004). Geostatistical methods are well established
in the field of water resources research to interpolate sparse data sets and to assess uncertainty of subsurface,
surface, and atmospheric variables (Abdu et al., 2008; Bárdossy & Li, 2008; Koch et al., 2014). As an alterna-
tive, machine learning methods provide tools to estimate uncertainty as well, but their application is so far
limited to digital soil mapping (Szatmári & Pásztor, 2018; Vaysse & Lagacherie, 2017). Their breakthrough in
water resources research awaits thorough testing and, in the meantime, hybrid models combining machine
learning and geostatistics provided powerful and well‐established tools to address uncertainties of machine
learning models.
The three main objectives of this paper are as follows: (1) to apply RF to model the depth to the uppermost
redox interface across Denmark at high spatial resolution (100 m), (2) to extend the machine learning model
with geostatistics (RFRK) to assess uncertainties and to obtain conditioned realizations of the redox depth,
and (3) to present a framework that facilitates easy incorporation of new field data into the modeling process.
2. Methods
2.1. Study Area
Denmark is located in northern Europe (54.5–57.8°N and 8.1–15.2°E) with an extent of approximately
43,000 km2. Topography varies from flat to gently undulating with a maximum elevation of 170 m. The
sequence of Pleistocene glaciations and postglacial processes have shaped the Danish landscape of today.
The geology of the eastern part of Denmark is dominated by Weichselian clayey till, whereas the western
part of the Danish peninsular (Jutland) is characterized by sediments of Saalian age (“hill islands”) inter-
twined by Weichselian sandy outwash plains.
2.2. Data
This study aims at modeling the depth of the redox interface, that is, the location of the upper reduced zone
in the soil column, at 100‐m resolution in Denmark by using a machine learning algorithm. For the purpose
of training the algorithm, we made use of an extensive data set containing ~13,000 boreholes with informa-
tion on the redox depth. These data were extracted from the national borehole database, Jupiter (Hansen &
Pjetursson, 2011). The transition from oxic to anoxic conditions, that is, the location of the redox interface,
can be determined by the color change in the sediment from brown (oxidized) to gray (reduced). The color
change has successfully served as a proxy for identifying the location of the redox depth in previous investi-
gations (Ernstsen et al., 1998; Hendry et al., 1984). Figure 1 depicts a map of Denmark indicating the location
of the boreholes and their derived redox depths. An initial visual assessment of the redox depth indicates
some apparent spatial patterns that are linked to the georegions and landscape types. Hansen,
Christensen, et al. (2014) brought forward the hypothesis that the present location of the redox interface
in Denmark is governed by postglacial processes. These promote a downward movement of the interface
since the beginning of the Holocene (11.700 BP) when reduced compounds were present throughout the
entire soil column. The cumulative influx of oxygen and nitrate into the subsurface controls the degree of
depletion of reduced compounds and hence the resulting downward mitigation of the redox interface. The
rate at which the redox interface mitigates has been accelerated in the past 60–70 years where nitrate, acting
as an additional oxidant next to oxygen, has been added to the system (Højberg et al., 2017). The variability in
mitigation rates can be explained by the amount of reduced compounds in the sediment (organic carbon and
pyrite), concentration of oxygen and nitrate in the infiltrating water, and the overall vertical flux. The Danish
moraine landscape, which covers most of the EastDK georegion, is characterized by a high field capacity,
which only allows air diffusion in the root zone, and the downward transport of oxygen through advection
is slow (Refsgaard et al., 1991). This typical till setting results in a shallow redox depth with small variability
across the moraine landscape. An exception is the Himmerland georegion, where most of the deep redox
depths are located. An underlying fractured kart functions as a very effective natural drain in the
Himmerland region, promoting a large influx of oxygen.
From initial analyses, it was found that some key landscapes, with distinct hydrogeochemical properties,
were not represented adequately in the data set, due to limited boreholes in these areas. This is the case
for lowland areas such as marshes, reclaimed areas, and organic soils (Figure 1). These landscapes are char-
acterized by a high content of organic or other reducing compounds, such as iron, which facilitate bacterial
Figure 1. (left) Map of Denmark showing the redox data (~13k) used for model training. The additional data are artifi-
cially added in undersampled regions (marshes, reclaimed areas, and lowlands with organic soils) where expert knowl-
edge supports a very shallow redox depth of 0.1 m. (right) Map of Denmark showing eight delineated georegions,
landscape types, and organic lowlands (mineral lowlands not shown).
activity, a groundwater table close to the surface, and low recharge rates. Based on expert knowledge, it is
very characteristic for these landscape settings that reduced conditions are evidently very close to the
terrain (Hoffmann et al., 2012). It was argued by Blicher‐Mathiesen (2012) that peat generally facilitates a
high denitrification potential due to its high organic content. This was further supported by Petersen et al.
(2018), who conducted intensive field work in a lowland with organic soils in eastern Jutland focusing on
the fate of nitrate. Their results indicated a very shallow redox depth of not more than 10 cm below
terrain for peat soils. Therefore, we decided to add artificial observations to the training data set with a
shallow redox depth of 0.1 m. The artificial observations were placed randomly in the abovementioned
landscapes with the same sampling density as in the borehole data set (0.3/km2), which resulted in 600
additional entries in the training data set.
The environmental covariates used to model the redox depth are stated in Table 1. The covariates comprise
soil texture, topography‐related variables, lithological information, outputs from a hydrological simulation
with the Danish national water resources model (DK‐model: Højberg et al., 2013), the geographic location,
a lowland classification into organic and mineral soils, the dominating landscape types (Figure 1), and the
general geological setting. The native spatial resolution of the covariates varies, but all covariates have been
resampled to 100 m to be in agreement with the desired output resolution. Following Hengl et al. (2018), the
addition of coordinates as explanatory variables is the simplest form to add spatial correlation to an RF
model. Moreover, general spatial trends of the redox conditions can be captured. The final data set used
for training the machine learning model was built by extracting values of all covariates at the locations of
the boreholes.
In order to continuously update the model of the depth to the redox interface, we present a concept (
section 2.4) allowing a straightforward incorporation of new data into the modeling process. In order to
proof this concept, we investigated recently acquired point data using a newly developed sampling probe
that measures the vertical profile of the redox potential. These data were obtained from a field campaign
in the eastern part of the Danish peninsular (Jutland) during autumn 2017. The redox profiles at 30 locations
across three field sites were used to derive the location of the redox depth.
Table 1
Overview of the Covariates Used to Model Redox Depth Using Random Forest
Mean/min/max
Name Description Original resolution (number of classes) Source
Topo Elevation above sea level (m) 1.6 m 11.9/−7.4/167.5 Danish Agency for Data Supply
Flowacc Number of upslope cells 243.5/0.0/485.0 and Efficiency
Slope Surface slope gradient (deg) 0.4/0.0/24.4
Stream Dis Horizontal distance to the 226.3/0.0/3196.3
nearest stream (m)
Quart Thick Vertical thickness of Quaternary N/A (25‐m contour 89.0/0.0.342.3 Binzer & Stockmarr (1994)
deposits map)
Depth Clay Depth to clay occurrence 500 m 13.5/0.0/355.3 DK‐model: Højberg et al. (2013)
min GWT Minimum simulated groundwater 5.5/0.0/106.3
table (meters below ground)
mean GWT Mean simulated groundwater 2.1/0.0/73.1
table (meters below ground)
Recharge Mean simulated recharge 0.8/0.0/6.2
(mm/day)
Georegion Geological regions of Denmark Scale 1:100000 8 classes Adhikari et al. (2014)
prediction can be derived by averaging over all trees where a certain data
event was retained as oob. Typically, the oob error converges as the number
of trees increases. Furthermore, it is possible to assess the importance of
each covariate, which can be important to better comprehend the trained
RF model and to evaluate its physical meaningfulness. This can be achieved
by permuting each covariate at a time, while holding the remaining covari-
ates unchanged, and recording how the oob error increases. Covariates that
show a high increase, when being permuted, are naturally considered as
important to model the target variable. A covariate that is strongly corre-
lated to another may show an alleged low importance, when being per-
muted, because its information is still contained in the correlated
covariate. Therefore, we applied two permutation strategies: first, permut-
ing each covariate at a time and, second, permuting several covariates,
which are physically related, at the same time as groups. Determining cov-
ariate importance is common practice in most RF applications in
geoscience (Cutler et al., 2007; Ließ et al., 2012). We applied the Scikit‐learn
package in Python (Pedregosa et al., 2012) for the RF modeling.
2.4. Random Forests Regression Kriging
Machine learning can easily be integrated in the classic regression kriging
(RK) approach (Hengl et al., 2007; Odeh et al., 1995). RK is based on two
components: (1) a linear model as trend and (2) a geostatiscal model of the
residuals of the linear model. In order to simulate complex nonlinear
environmental systems, one can simply replace the linear trend model
with a more advanced machine learning technique. The hybrid model
applied in this study integrated RF into the standard RK approach, which
is referred to as RFRK. RFRK has previously been successfully applied by
several studies modeling environmental variables (Ahmed et al., 2017;
Guo et al., 2015; Hengl et al., 2015; Viscarra Rossel et al., 2014), but to
our knowledge, it has not yet been tested on water resources‐related vari-
Figure 2. Workflow diagram for the random forest regression kriging ables such as the redox depth.
framework applied in this study. The top part (light blue background) covers
the standard random forest modeling. The bottom part (light green back-
In the applied RFRK model the geostatistical model was based on a vario-
ground) shows the residual simulation. oob = out of bag; sGs = sequential gram, which described the spatial continuity of the oob residuals. RK typi-
Gaussian simulation. cally applies ordinary kriging to model the best estimate of the residuals.
Alternatively, a conditional sequential Gaussian simulation (sGs) can pro-
duce equally probable stochastic realizations of the residuals, each honoring the defined variogram and con-
ditioned to the local residuals (Leuangthong & Deutsch, 2004). The variability of such a set of realizations
can be used to express the uncertainty at grid scale. Also, the individual realizations can be used to propagate
the uncertainty through a hydrological model to assess the predictive uncertainty of relevant hydrological
variables. This makes the method attractive to support decision making (Refsgaard et al., 2014). RK extended
by RF can be expressed by
PRFRK ðs0 Þ ¼ T RF ðs0 Þ þ besGs ðs0 Þ; (1)
where TRF is the RF output at location s0. êsGs is the simulated residual obtained from a conditional sGs. The
sum of trend (TRF) and simulated residual (êsGs) is the final RFRK prediction (PRFRK). The RFRK workflow
is depicted in Figure 2.
Following the original work by Matheron (1963), the calculation of the omnidirectional empirical semivar-
iance γ of the oob residuals (e) at lag h was calculated as
1 nðhÞ
γ ðhÞ ¼ ∑ ½eðsi Þ−eðsi þ hÞ2 ; (2)
2nðhÞ i¼1
where n(h) is the total number of data pairs at a given lag h, e(si) represents the RF residual at location si, and
e(si + h) is the value separated by lag h from si. Based on γ, a variogram model was fitted to describe the spa-
tial autocorrelation structure of the oob residuals (Deutsch & Journel, 1998). A variogram model is defined
by a type, range, sill, and nugget. The type sets the general shape of the variogram. Typically, spherical or
exponential models are found suitable. The range defines the distance at which spatial correlation vanishes.
The sill represents the variance that is attained beyond the range. Lastly, the nugget attributes the variance at
distance equal to 0, and it can be understood as measurement uncertainty. As described by Pebesma and
Wesseling (1998), a conditional sGs uses a random path through all simulation grids, and at each grid, a local
conditional distribution, based on observations and previously estimated data, is calculated. For this
purpose, the predefined variogram model is utilized and the simulated value is then drawn randomly from
the obtained distribution and added to the conditioning data set. Once all simulation grids are filled, a single
realization is built, and further realizations can be estimated by changing the random seed that affects the
simulation path. The R package Gstat (Pebesma, 2004) was used for variogram modeling and sGs.
2.5. Local Variance Scaling
An RF model will typically only describe a fraction of the spatial variability that is contained in the training
data set. The resulting spatial variability of the RF residuals is caused by missing or poorly resolved covari-
ates, sub‐grid variability, and measurement uncertainty. As outlined by Kanevski et al. (2004), the variogram
analysis of the RF residuals can result in two conclusions: (1) the residuals are spatially correlated or (2) no
correlation is evident. In the first case, the available covariates used to train the RF model cannot describe
the full spatial variability of the target variable and one or more crucial covariates are most likely missing.
However, the second case underlines that RF modeled all spatial variability, as represented in the training
data set best possibly. The resulting residuals are characterized as white noise and are mere effect of
small‐scale heterogeneity and measurement uncertainty. Therefore, additional covariates will likely not
result in a better performance. This is based on the assumption that a nonrandom environmental covariate
cannot explain residuals that are white noise.
Previous research has rarely touched upon this distinction when applying RFRK. In most studies, a vario-
gram was fitted automatically to γ without further commenting (Guo et al., 2015; Hengl et al., 2015). In this
study, we highlighted the importance to investigate the spatial correlation of the residuals in more detail. We
suggested identifying areas where spatial correlation is evident and thus the standard RFRK can be applied.
For other areas where this is not the case, we proposed a local variance scaling approach that simulates the
residuals as noise while allowing for local conditioning and considering nonstationarities. In other words,
the scaling adjusted the sill of the defined variogram model locally to meet nonstationary conditions.
This was achieved by setting the range of the variogram model to approximately the size of the grid to be
simulated. Such a steep variogram model simulates grids without observations as white noise following
the sill of the variogram model, while observations are honored at coinciding grids. Following this approach,
the variance is constant over the entire domain despite the fact that the residual variance may actually vary.
This has direct implications on the uncertainty estimations and results in crude overestimation and under-
estimation. This limitation could be overcome by scaling the simulated variance to a local variance. Here, we
suggested a two‐stage transformation of the estimated RF residuals for each simulation grid. First, the z‐
score was computed for the sGs‐based set of simulated residuals:
where z(si) is the z‐score of ê(si), the predicted residual at location si using sGs. The z‐score was calculated by
subtracting the mean of the realizations at location si (μsGs(si)) and dividing by the standard deviation of the
realizations (σsGs(si)). This transformation resulted in a set of realizations with a mean of 0 and a standard
deviation of 1 for each grid. The backtransformation from z(si) was then undertaken by considering the local
standard deviation (σloc) instead of σsGs while maintaining the original mean (μsGs(si)):
σ sGs ðsi Þ
bzðsi Þ ¼ zðsi Þ*σ loc ðsi Þ* þ μsGs ðsi Þ; (4)
maxðσ sGs Þ
where σloc is the local standard deviation at location si. σloc does not consider conditioning to observations,
and therefore, it is multiplied by the standard deviation of the realizations (σsGs) normalized by the maxi-
mum standard deviation, which corresponds to the sill defined in the variogram model. This multiplier
can be interpreted as borehole proximity and varies between 0 and 1. Values close to 0 indicate vicinity to
3. Results
3.1. RF Model
For the prediction of the redox depth, an RF model was trained using the
18 selected covariates and log‐transformed redox data at ~13,000 bore-
Figure 3. Results obtained from training a random forest model to simulate
holes specifying the location of the redox interface in meters below ter-
the redox depth. The out‐of‐bag (oob) prediction represents the internal
cross‐validation of the random forest algorithm. The data labeled as “all” are rain. After log transformation, the redox depth data exhibit roughly a
the final prediction averaging over all (training and oob) decision trees. The normal distribution, and initial tests indicated that transforming the data
2
R scores (coefficient of determination) quantify the performance of the has a positive effect on the RF accuracy.
random forest model. The data are log10 transformed.
The RF model was parameterized in the following way: The number of
decision trees was set to 1,000, bootstrapping with replacement was cho-
sen to sample the training data set, 33% of the covariates were considered
when splitting the data, the trees were fully expanded and thus not pruned, the mean square error was used
as a criterion to find the best split, and regression was selected as method.
Figure 3 shows the RF estimations of the redox depth based on (1) averaging only over trees where data were
held back as an oob sample and (2) averaging over all trees. The latter is just included as reference and
should not be consulted when quantifying the accuracy of the RF model. On the contrary, the oob prediction
can be regarded as an independent validation test, and the resulting coefficient of determination (R2) of 0.48
indicates that the RF model is able to capture half of the overall variability of the redox depth.
Figure 4. Variable importance quantified for the covariates used in the random forest training by implementing the con-
2
cept of permutation accuracy. The decrease in out‐of‐bag (oob) performance (R ) was used as a measure of variable
importance. Permutation was applied not only to single covariates (red) but also to groups of covariates (blue).
Figure 5. The computed semivariances for the random forest (RF) residuals (circles). The variogram analysis was split
into two parts: the hill island domain in Western Denmark and the rest of the country.
Figure 4 shows the derived covariate importance based on the concept of permutation accuracy and stresses
that the categorical maps have the highest importance in the trained RF model as the performance drops by
approximately 60% when the entire group is permuted. The lowland covariate indicates where the redox
depth is very close to the terrain, and therefore, it is the most crucial individual covariate with a decrease
of around 40%. The second most important covariate group is the one related to the digital elevation model
(~30%) followed by the national hydrological model (~20%). The least important covariates are soil texture
related as well as the Quaternary thickness and depth of clay, all describing the hydrogeological conditions.
3.2. Residual Model
The geostatistical part of RFRK simulates the RF residuals using sGs, which requires the definition of a var-
iogram model. Semivariances for the log‐transformed RF residuals were computed omnidirectionally with a
lag distance of 200 m. Initially, the semivariances were calculated for the entire domain as well as for indi-
vidual landscape types and georegions (Figure 1) to investigate the existence of any given nonstationarities
in the RF residuals. We found that the key distinction was between the hill island landscape in the western
part of the Jutland and the rest of the country (Figure 5). For the hill island landscape, a clear correlation
length was detected, which indicates a certain degree of spatial continuity in the RF residuals. However,
for the latter, no spatial correlation was discovered, indicating that the RF residuals can be understood as
a random variable. As shown by Figure 5, spherical variogram models with a range of 200 m were fitted
manually to the calculated semivariances for both domains. The short range allows sGs to condition indivi-
dual grids to local observations. The hill island variogram is nested with an additional spherical model with a
range of 2,500 m. Both variograms are characterized by a nugget of 0.02 m2, which refers to approximately
5% of the original RF prediction when being backtransformed. In this way, we quantitatively considered the
uncertainty of the redox depth related to uncertainties in the borehole description.
The results from the conditional sGs of the RF residuals are depicted in Figure 6. The mean and the variance
are based on 800 realizations. The effect of having two independent simulation domains, namely, the hill
island and the remaining part of the country, is clearly visible. In the hill island domain, the mean of the rea-
lizations shows smooth patterns typical for variogram modeling. However, for the rest of the country, the
mean is dominantly 0 for cells without borehole observations where the individual sGs realizations cancel
each other out. Nevertheless, observations are conditioned at grid scale given the variogram range of
200 m. The variance in the hill island domain is highest at grids unaffected by any observations; thus, the
distance to the closest observation is beyond the variogram range. The initial variance outside the hill island
reflected the sill that was defined by the variogram model in Figure 5. However, the applied scaling approach
(equation (4)) adjusted the variance to be in agreement with the local conditions, that is, the variance of
Figure 6. (left) The modeled redox depth using random forest (RF). In the random forest regression kriging framework, the RF result was used as trend, and the RF
residuals were simulated using sequential Gaussian simulation (sGs). (middle) The average of 800 realizations simulating the RF residuals using sGs. The hill
islands (delineated by a black line) and the rest of Denmark were simulated with specific variograms, as indicated by Figure 5. (right) The variance of 800 reali-
zations simulating the RF residuals using sGs. The variance in cells that are not in the hill islands domain was rescaled to agree with the local variance as presented
in equation (4).
residuals within a 10‐km radius. It becomes clear that the variance varies drastically across the country and
that this variation has to be accounted for in the final set of redox depth realizations. At grids that coincide
with observations, the variance is low, which is a result from conditioning the residual simulation to
boreholes. The residuals outside the hill island landscape are simulated as white noise at grids not
coinciding with a borehole, whereas grids that contain a borehole are conditioned to the RF residual. For
the former, this results in a mean of 0 when averaging the entire set of residual realizations. Figure 6 also
contains the RF estimate of the redox depth that can be used in combination with the mean and variance
of the realizations in the RFRK technique to provide a best estimate and an uncertainty assessment.
3.3. Random Forest Regression Kriging
The RFRK results are visualized in Figure 7 for a focus area in the western part of Jutland. The RFRK esti-
mate is the sum of the trend model (RF) and the mean of the simulated realizations of the RF residuals using
sGs (Figure 6). Through RFRK, the R2 increases from 0.48 (RF oob) to 0.83. This cannot be understood as an
increase in overall performance, because the residual simulation is conditioned to all available boreholes
and did not undergo an independent cross‐validation test. The variance of the realizations can be used to
define confidence intervals to express uncertainty. For some grids, the spread between the 95% confidence
intervals can be wide, for instance, from 1 to 30 m, while for others, it is rather narrow, for instance, from
0 to 1 m. The uncertainty is generally larger in the hill island domain, where the RF residuals underlie a clear
spatial correlation with a high variance. Moreover, the individual realizations of sGs can be used to visualize
the stochastic nature of the applied RFRK framework. The realizations of the redox depth are conditioned to
local borehole data but manifest clear structural differences in their patterns while maintaining an
overall resemblance.
3.4. Adding New Data
One advantage of RFRK is that new data can easily be integrated into the modeling process. We propose a
framework where new data, acquired from field campaigns, can be utilized to build local models of the redox
depth. We suggest that the RF model does not require new training after each field campaign; instead, the
residual model can be rerun to condition sGs to the new data and to update the local residual variance.
We have investigated the number of realizations needed to reach a stable mean in order to make this frame-
work as efficient as possible. Figure 8 visualizes the convergence of the cumulative normalized mean at
2
Figure 7. The results of the random forest regression kriging (RFRK) framework are displayed for a ~2,500‐km area in
Western Jutland. The top row shows the mean RFRK in the middle, and the 95% confidence intervals are given to
its left (−2 * std) and to its right (+2 * std). The bottom row exemplifies three of the 800 simulated realizations. The hill
islands domain is delineated with a black line.
50,000 randomly selected grids across the modeling domain. It is apparent that the mean stabilizes after
approximately 500 realizations, where 98% of the grids are within ±5% relative to the mean obtained with
800 realizations. Therefore, we recommend to run 500 realizations when updating a local model of the
redox depth.
Figure 8. Analyzing convergence of the stochastic realizations of the random forest residuals generated by sequential
Gaussian simulation. The normalized cumulative mean is plotted for 50,000 randomly selected grids across Denmark.
The cumulative mean was calculated by dividing the cumulative sum by the number of realization. The values were then
normalized by the redox depth obtained using the entire population of 800 realizations. The dashed lines in red
indicate a deviation of ±5% from the modeled redox depth based on the entire population and can be understood as an
acceptable error, which is defined subjectively. The dotted lines in black represent the 1st and 99th percentiles from the
50,000 grids at each cumulative step and help visualize the spread.
Figure 9. Example of the effect of incorporating new data into the random forest regression kriging (RFRK) framework
2
for the Norsminde catchment (~100 km ) in eastern Jutland. The left panel shows the final RFRK estimate based on
the average of 500 residual simulations. The blue dots represent redox depth data from the existing training data set. The
black dots indicate new data, obtained from recent fieldwork, that were added to the residual error simulation. The middle
panel shows the relative change in the final RFRK estimate considering the mean of 500 realizations. The right panel
shows the corresponding relative change of RFRK uncertainty. The latter was computed as the range of 95% confidence
intervals as illustrated in Figure 7.
In autumn 2017, a field campaign was carried out in the Norsminde catchment located in the eastern part of
Jutland. A newly developed sampling probe, measuring vertical profiles of the redox potential, was tested at
30 locations across three sites. The interpreted redox depth was added to the residual simulation as condi-
tioning data, and 500 stochastic realizations were generated (Figure 9). In the eastern part of Denmark,
the residuals are simulated by a variogram with a range of 200 m. Therefore, the updated redox depth is only
affected at grids coinciding with the new data, where the RFRK model was changed by up to 25% by adding
new conditioning data to the residual model. Adding new data also updates the local variance of the resi-
duals and thereby the uncertainty. Hence, in contrast to the mean of the realizations, the uncertainty is
affected over a larger area (the 10 km over which the local variance is calculated). The uncertainty is repre-
sented as the range of the 95% confidence intervals as shown in Figure 7. At grids coinciding with the new
observations, the uncertainty reduction was most pronounced with up to 50% and varies between 5% and
15% for the rest of the area, which corresponds to approximately 2 m in most places.
The time it takes to retrain and rerun the RF model is around 10 min. However, we recommend not to repeat
this step of the workflow when building a local model, because the RF model was hardly affected by addition
of 30 new observations to a training data set containing over 13,000 data entries. Our tests indicated that the
R2 of the oob prediction changed in the third digit after including the 30 new observations. We propose a
framework where the RF is first updated after an extensive amount of new data has been collected.
Further investigations are required to identify how much new data are needed to have a significant effect
on the RF model.
4. Discussion
4.1. The Redox Depth
The color change in the sediment, as reported by borehole descriptions, is used as a proxy to derive the
uppermost transition from oxic to anoxic conditions in this study. The process of borehole interpretation
underlies intrinsic uncertainties that are difficult to quantify. We chose to account for this by defining a nug-
get in the variogram model used for simulating the RF residuals.
We applied RFRK to identify the uppermost redox interface, where conditions change from oxic to anoxic
conditions. Knowledge of this transition is crucial in simulating the fate of nitrate, which requires combining
the redox interface with groundwater flow modeling of transport pathways. In complex aquifer systems, oxi-
dized groundwater may be found below the redox interface determined from sediment profiles or ground-
water chemistry. This may be caused by horizontally connected systems of sand aquifers or incomplete
reduction due to a combination of limited thickness of the reducing layer and the existence of fractures,
allowing preferential flow paths (Robertson et al., 1996). This has been expressed by several studies, where
redox zonation maps were derived based on groundwater chemistry for complex aquifer systems (Close
et al., 2016; Hsu et al., 2010; Lee et al., 2008). However, the importance of such conditions, which are com-
monly related to the hydrogeochemistry of the system and the scale of application, must be considered and
reflected in the use of the redox map, which is outside the scope of the present study. Nevertheless, we would
expect a good overall agreement between the estimate redox depth and occurrence of oxic, suboxic, and
anoxic groundwater.
conditioned to all available data. Based on the spatial structure of the RF residuals, it is not possible to expect
an improvement of the RF model at unsampled locations, because cells without coinciding observations are
predominantly simulated as white noise. We investigated this with a simple cross‐validation test, where 10%
of the training data were retained for building the RF model. This test revealed that the residual kriging,
based on the remaining 90% of data, showed no measurable improvement for the retained 10%.
Nevertheless, the residual simulation is valuable, as it ensures that local redox data are honored best possibly
and, further, it allows to assess uncertainty by means of the set of realizations. Even though there is little
knowledge about the residual mean of the realizations at unsampled locations, there exists reliable knowl-
edge on the variance.
Figure 10. Comparison of two approaches of mapping redox depth at national scale in Denmark: using (left) expert
knowledge and (right) machine learning, that is, random forest regression kriging. The manual interpretation of the
redox depth was conducted at 1‐km spatial resolution and has been published by the Geological Survey of Denmark and
Greenland (Ernstsen et al., 2010). The machine learning model provides results with a resolution of 100 m. Two focus
areas were selected to better exemplify the differences between the two approaches. The rectangular extend zooms into the
Himmerland area in Northern Jutland, and the quadratic extend shows the Norsminde catchment located in Western
Jutland.
Although the final RFRK prediction may only be affected at grids that coincide with the location of the new
observations, the impact on the variance of the realizations is apparent for a larger area. This is controlled by
the 10‐km radius we chose to compute the local variance of the RF residuals used to scale the variance of the
realizations. Ideally, the radius should be further reduced, but this was limited by data availability.
Nevertheless, considering the nonstationarity of the RF residual variance across the domain (Figure 6),
we found the trade‐off quite reasonable and inevitable.
& von Platen, 2014), and in the second stage, 1,000 new boreholes were added. These additional boreholes
acted as an independent validation test, and 42% of them were placed in the correct category in the
originally map.
The natural drawback of a data‐driven model is that it does not contain any a priori conceptual understand-
ing of the governing physical processes. Only processes that are conveyed by the training data set to the
machine learning model can be represented accordingly. Certain processes that are not represented can
therefore not be reflected in the output. This underlines the importance of expert knowledge in the process
of building a data‐driven model. In our study, expert knowledge became a key asset by guiding the RF model.
This was achieved by selecting appropriate covariates as well as by adding 600 additional synthetic observa-
tions in underrepresented landscape settings to the training data set. Further dissimilarities between the 1‐
km and 100‐m maps can result from the fact that some areas may still not be represented adequately in the
training data set. The map developed by the RFRK method may thus be further improved by adding artificial
data, if other landscape features with distinct redox conditions are undersampled in the present data set. This
study has demonstrated that a geologist's expert knowledge should ideally be consulted when building a
machine learning model. This interplay will play a vital role in advancing machine learning applications
in the field of water resources.
5. Conclusions
The overall aim of this study was to model the depth to the uppermost reduced layer where denitrification
takes place, at 100‐m resolution across Denmark. This was achieved by combining the machine learning
technique RF with geostatistics. The applied RFRK framework uses RF as trend and a sGs of the RF residuals
to condition to local observations and to quantify uncertainties. We draw the following main conclusion
from our work:
1. We highlight the vital role of expert knowledge in the process of building data‐driven models. In our
study, expert knowledge guided the RF model in underrepresented landscapes by addition of artificial
observations.
2. RF was able to capture the redox depth variability in the training data set with an R2 score of 0.48. Subgrid
variability and measurement uncertainty are likely the main factors causing the remaining mismatch
between model and observations. Among the continuous covariates, the ones related to topography
and groundwater table were ranked as the most important variables.
3. The complex hill island landscape in western Denmark manifests a clear spatial correlation of the RF
residuals. In the remaining part of the country, residuals are not spatially correlated, and therefore,
two individual variograms were defined to model the residuals adequately.
4. RFRK using conditional sGs to estimate the RF residuals honors local observations, and the spread of 800
realizations is used to map uncertainties. Uncertainties vary across the domain and strongly depend on
local conditions that was captured by the proposed local variance scaling approach.
5. In order to reduce uncertainties, we designed a framework that allows for easy integration of newly
acquired redox data for the purpose of building local redox models. In an example case, where 30 new
Acknowledgments observations were integrated into the RFRK model, uncertainties were reduced by up to 50% at grids con-
The work has been carried out with taining the additional data and by up to 15% in its vicinity.
financial support granted by the
TReNDS research project funded by the
Innovation Fund Denmark (grant no.
4106‐00027B). Results will be made References
available online at https://eng.geus.dk/ Abdu, H., Robinson, D. A., Seyfried, M., & Jones, S. B. (2008). Geophysical imaging of watershed subsurface patterns and prediction of soil
products‐services‐facilities/data‐and‐ texture and water holding capacity. Water Resources Research, 44, W00D18. https://doi.org/10.1029/2008WR007043
maps/maps‐of‐denmark/. The borehole Adhikari, K., Hartemink, A. E., Minasny, B., Bou Kheir, R., Greve, M. B., & Greve, M. H. (2014). Digital mapping of soil organic carbon
data are accessible via the Jupiter data- contents and stocks in Denmark. PLoS ONE, 9(8), e105519. https://doi.org/10.1371/journal.pone.0105519
base at https://eng.geus.dk/products‐ Adhikari, K., Kheir, R. B., Greve, M. B., Bøcher, P. K., Malone, B. P., Minasny, B., et al. (2013). High‐resolution 3‐D mapping of soil texture
services‐facilities/data‐and‐maps/ in Denmark. Soil Science Society of America Journal, 77(3). https://doi.org/10.2136/sssaj2012.0275
national‐well‐database‐jupiter/. The Agency, D. N. (2016). Vandområdeplan 2015–2021 for Vandområdedistrikt Jylland og Fyn (Water management plan 2015–2021 for River
authors will gladly share covariate data Basin district Jutland and Funen).
as well as intermediate results or scripts Ahmed, Z. U., Woodbury, P. B., Sanderman, J., Hawke, B., Jauss, V., Solomon, D., & Lehmann, J. (2017). Assessing soil carbon vulnerability
upon request. We gratefully acknowl- in the Western USA by geospatial modeling of pyrogenic and particulate carbon stocks. Journal of Geophysical Research: Biogeosciences,
edge Lars Troldborg and Anders Bjørn 122, 354–369. https://doi.org/10.1002/2016JG003488
Møller who provided us with covariate Bachmair, S., & Weiler, M. (2012). Hillslope characteristics as controls of subsurface flow variability. Hydrology and Earth System Sciences,
data. 16(10), 3699–3715. https://doi.org/10.5194/hess‐16‐3699‐2012
Bárdossy, A., & Li, J. (2008). Geostatistical interpolation using copulas. Water Resources Research, 44, W07412. https://doi.org/10.1029/
2007WR006115
Bateman, I. J., Harwood, A. R., Mace, G. M., Watson, R. T., Abson, D. J., Andrews, B., et al. (2013). Bringing ecosystem services
into economic decision‐making: Land use in the United Kingdom. Science, 341(6141), 45–50. https://doi.org/10.1126/
science.1234379
Baudron, P., Alonso‐Sarría, F., García‐Aróstegui, J. L., Cánovas‐García, F., Martínez‐Vicente, D., & Moreno‐Brotóns, J. (2013). Identifying
the origin of groundwater samples in a multi‐layer aquifer system with random forest classification. Journal of Hydrology, 499, 303–315.
https://doi.org/10.1016/j.jhydrol.2013.07.009
Behrens, T., Schmidt, K., MacMillan, R. A., & Viscarra Rossel, R. A. (2018). Multi‐scale digital soil mapping with deep learning. Scientific
Reports, 8(1), 15244. https://doi.org/10.1038/s41598‐018‐33516‐6
Binzer, K., & Stockmarr, J. (1994). soil, Kortserie‐ DGU, Danmarks Geol. Undersøgelse.
Blicher‐Mathiesen, G. (2012). Notat om status for N‐udledning fra lavbundsarealer.
Breiman, L. (2001). Random forests. Machine Learning, 45(1), 5–32. https://doi.org/10.1023/A:1010933404324
Close, M. E., Abraham, P., Humphries, B., Lilburne, L., Cuthill, T., & Wilson, S. (2016). Predicting groundwater redox status on a
regional scale using linear discriminant analysis. Journal of Contaminant Hydrology, 191, 19–32. https://doi.org/10.1016/j.
jconhyd.2016.04.006
Cutler, D. R., Edwards, T. C., Beard, K. H., Cutler, A., Hess, K. T., Gibson, J., & Lawler, J. J. (2007). Random forests for classification in
ecology. Ecology, 88(11), 2783–2792. https://doi.org/10.1890/07‐0539.1
Deutsch, C. V., & Journel, A. G. (1998). GSLIB: Geostatistical software library and user's guide. [online] Available from: http://orton.catie.
ac.cr/cgi‐bin/wxis.exe/?IsisScript=orton.xis&method=post&formato=2&cantidad=1&expresion=mfn=072844
Eidem, J. M., Simpkins, W. W., & Burkart, M. R. (1999). Geology, groundwater flow, and water quality in the Walnut Creek Watershed.
Journal of Environmental Quality, 28(1), 60–69. https://doi.org/10.2134/jeq1999.00472425002800010006x
Erisman, J. W., Galloway, J. N., Seitzinger, S., Bleeker, A., Dise, N. B., Petrescu, A. M. R., et al. (2013). Consequences of human modification
of the global nitrogen cycle. Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences, 368(1621), 20130116.
https://doi.org/10.1098/rstb.2013.0116
Ernstsen, V., Binnerup, S. J., & Sørensen, J. (1998). Reduction of nitrate in clayey subsoils controlled by geochemical and microbial barriers.
Geomicrobiology Journal, 15(3), 195–207. https://doi.org/10.1080/01490459809378076
Ernstsen, V., Højberg, A. L., Jakobsen, P. R., von Platen, F., Tougaard, L., Hansen, J. R., et al. (2006). Beregning af nitratreduktionsfaktorer
for zonen mellem rodzonen og frem til vandløbet. Data og metode for 1. generationskortet, GEUS, Geological Survey of Denmark and
Greenland.
Ernstsen, V., Jakobsen, P. R., & von Platen, F. (2010). Et første landsdækkende redoxkort. Vand Jord, 4, 159–160.
Ernstsen, V., von Platen, F., & for Danmark og Grønland, D. N. G. U. (2014). Opdatering af det nationale redoxkort fra 2006: til brug for den
nationale kvælstofmodel 2015. GEUS, Geological Survey of Denmark and Greenland. Retrieved from [online] Available from: https://
books.google.dk/books?id=D‐B8AQAACAAJ
Fenton, O., Richards, K. G., Kirwan, L., Khalil, M. I., & Healy, M. G. (2009). Factors affecting nitrate distribution in shallow groundwater
under a beef farm in south eastern Ireland. Journal of Environmental Management, 90(10), 3135–3146. https://doi.org/10.1016/j.
jenvman.2009.05.024
Francke, T., López‐Tarazón, J. A., & Schröder, B. (2008). Estimation of suspended sediment concentration and yield using linear models,
random forests and quantile regression forests. Hydrological Processes, 22(25), 4892–4904. https://doi.org/10.1002/hyp.7110
Guo, P. T., Li, M. F., Luo, W., Tang, Q. F., Liu, Z. W., & Lin, Z. M. (2015). Digital mapping of soil organic matter for rubber plantation at
regional scale: An application of random forest plus residuals kriging approach. Geoderma, 237‐238, 49–59. https://doi.org/10.1016/j.
geoderma.2014.08.009
Hansen, A. L., Christensen, B. S. B., Ernstsen, V., He, X., & Refsgaard, J. C. (2014). A concept for estimating depth of the redox interface for
catchment‐scale nitrate modelling in a till area in Denmark. Hydrogeology Journal, 22(7), 1639–1655. https://doi.org/10.1007/s10040‐
014‐1152‐y
Hansen, A. L., Gunderman, D., He, X., & Refsgaard, J. C. (2014). Uncertainty assessment of spatially distributed nitrate reduction potential
in groundwater using multiple geological realizations. Journal of Hydrology, 519(PA), 225–237. https://doi.org/10.1016/j.
jhydrol.2014.07.013
Hansen, J. R., Ernstsen, V., Refsgaard, J. C., & Hansen, S. (2008). Hydrogeology Journal, 16(7), 1251–1266. https://doi.org/10.1007/s10040‐
008‐0330‐1
Hansen, J. R., Refsgaard, J. C., Ernstsen, V., Hansen, S., Styczen, M., & Poulsen, R. N. (2009). An integrated and physically based nitrogen
cycle catchment model. Hydrology Research, 40(4), 347. https://doi.org/10.2166/nh.2009.035
Hansen, M., & Pjetursson, B. Free, online Danish shallow geological data, Geol. Surv. Denmark Greenl. Bull. 2011.
Hashemi, F., Olesen, J. E., Hansen, A. L., Børgesen, C. D., & Dalgaard, T. (2018). Spatially differentiated strategies for reducing nitrate loads
from agriculture in two Danish catchments. Journal of Environmental Management, 208, 77–91. https://doi.org/10.1016/j.
jenvman.2017.12.001
Hendry, M. J., McCready, R. G. L., & Gould, W. D. (1984). Distribution, source and evolution of nitrate in a glacial till of southern Alberta,
Canada. Journal of Hydrology, 70(1‐4), 177–198. https://doi.org/10.1016/0022‐1694(84)90121‐5
Hengl, T., Heuvelink, G. B. M., Kempen, B., Leenaars, J. G. B., Walsh, M. G., Shepherd, K. D., et al. (2015). Mapping soil properties of Africa
at 250 m resolution: Random forests significantly improve current predictions. PLoS ONE, 10(6). https://doi.org/10.1371/journal.
pone.0125814
Hengl, T., Heuvelink, G. B. M., & Rossiter, D. G. (2007). About regression‐kriging: From equations to case studies. Computational
Geosciences, 33(10), 1301–1315. https://doi.org/10.1016/j.cageo.2007.05.001
Hengl, T., Heuvelink, G. B. M., & Stein, A. (2004). A generic framework for spatial prediction of soil variables based on regression‐kriging.
Geoderma, 120(1–2), 75–93. https://doi.org/10.1016/j.geoderma.2003.08.018
Hengl, T., Nussbaum, M., Wright, M. N., Heuvelink, G. B. M., & Gräler, B. (2018). Random forest as a generic framework for predictive
modeling of spatial and spatio‐temporal variables. PeerJ, 6, e5518. https://doi.org/10.7717/peerj.5518
Hinkle, S. R., & Tesoriero, A. J. (2014). Nitrogen speciation and trends, and prediction of denitrification extent, in shallow US groundwater.
Journal of Hydrology, 509, 343–353. https://doi.org/10.1016/j.jhydrol.2013.11.048
Hoffmann, C. C., Heiberg, L., Audet, J., Schønfeldt, B., Fuglsang, A., Kronvang, B., et al. (2012). Low phosphorus release but high nitrogen
removal in two restored riparian wetlands inundated with agricultural drainage water. Ecological Engineering, 46, 75–87. https://doi.org/
10.1016/j.ecoleng.2012.04.039
Højberg, A. L., Hansen, A. L., Wachniew, P., Żurek, A. J., Virtanen, S., Arustiene, J., et al. (2017). Review and assessment of nitrate
reduction in groundwater in the Baltic Sea Basin. Journal of Hydrology: Regional Studies, 12, 50–68. https://doi.org/10.1016/j.
ejrh.2017.04.001
Højberg, A. L., Troldborg, L., Stisen, S., Christensen, B. B. S., & Henriksen, H. J. (2013). Stakeholder driven update and improvement of a
national water resources model. Environmental Modelling and Software, 40, 202–213. https://doi.org/10.1016/j.envsoft.2012.09.010
Høyer, A. S., Jørgensen, F., Piotrowski, J. A., & Jakobsen, P. R. (2013). Deeply rooted glaciotectonism in western Denmark: Geological
composition, structural characteristics and the origin of Varde hill‐island. Journal of Quaternary Science, 28(7), 683–696. https://doi.org/
10.1002/jqs.2667
Hsu, C. H., Han, S. T., Kao, Y. H., & Liu, C. W. (2010). Redox characteristics and zonation of arsenic‐affected multi‐layers aquifers in the
Choushui River alluvial fan, Taiwan. Journal of Hydrology, 391(3‐4), 351–366. https://doi.org/10.1016/j.jhydrol.2010.07.037
Jacobsen, B. H., & Hansen, A. L. (2016). Economic gains from targeted measures related to non‐point pollution in agriculture based on
detailed nitrate reduction maps. Science of the Total Environment, 556, 264–275. https://doi.org/10.1016/j.scitotenv.2016.01.103
Jessen, S., Postma, D., Thorling, L., Müller, S., Leskelä, J., & Engesgaard, P. (2017). Decadal variations in groundwater quality: A legacy
from nitrate leaching and denitrification by pyrite in a sandy aquifer. Water Resources Research, 53, 184–198. https://doi.org/10.1002/
2016WR018995
Jørgensen, P. R., Urup, J., Helstrup, T., Jensen, M. B., Eiland, F., & Vinther, F. P. (2004). Transport and reduction of nitrate in
clayey till underneath forest and arable land. Journal of Contaminant Hydrology, 73(1–4), 207–226. https://doi.org/10.1016/j.
jconhyd.2004.01.005
Kanevski, M., Parkin, R., Pozdnukhov, A., Timonin, V., Maignan, M., Demyanov, V., & Canu, S. (2004). Environmental data mining and
modeling based on machine learning algorithms and geostatistics. Environmental Modelling and Software, 19(9), 845–855. https://doi.
org/10.1016/j.envsoft.2003.03.004
Keller, C. K., & Van der Kamp, G. (1988). Hydrogeology of two Saskatchewan tills, II. Occurrence of sulfate and implications for soil
salinity. Journal of Hydrology, 101(1‐4), 123–144. https://doi.org/10.1016/0022‐1694(88)90031‐5
Koch, J., He, X., Jensen, K. H., & Refsgaard, J. C. (2014). Challenges in conditioning a stochastic geological model of a heterogeneous
glacial aquifer to a comprehensive soft data set. Hydrology and Earth System Sciences, 18(8), 2907–2923. https://doi.org/10.5194/hess‐
18‐2907‐2014
Kuhr, P., Haider, J., Kreins, P., Kunkel, R., Tetzlaff, B., Vereecken, H., & Wendland, F. (2013). Model based assessment of nitrate pollution
of water resources on a federal state level for the dimensioning of agro‐environmental reduction strategies. Water Resources
Management, 27(3), 885–909. https://doi.org/10.1007/s11269‐012‐0221‐z
Lee, J. J., Jang, C. S., Wang, S. W., Liang, C. P., & Liu, C. W. (2008). Delineation of spatial redox zones using discriminant analysis and
geochemical modelling in arsenic‐affected alluvial aquifers. Hydrological Processes, 22(16), 3029–3041. https://doi.org/10.1002/hyp.6884
Leuangthong, O., & Deutsch, C. V. (2004). Transformation of residuals to avoid artifacts in geostatistical modelling with a trend.
Mathematical Geology, 36(3), 287–305. https://doi.org/10.1023/B:MATG.0000028438.48852.b0
Li, J., Alvarez, B., Siwabessy, J., Tran, M., Huang, Z., Przeslawski, R., Radke, L., et al. (2017). Application of random forest, generalised
linear model and their hybrid methods with geostatistical techniques to count data: Predicting sponge species richness. Environmental
Modelling and Software, 97, 112–129. https://doi.org/10.1016/j.envsoft.2017.07.016
Li, J., Heap, A. D., Potter, A., & Daniell, J. J. (2011). Application of machine learning methods to spatial interpolation of environmental
variables. Environmental Modelling and Software, 26(12), 1647–1659. https://doi.org/10.1016/j.envsoft.2011.07.004
Ließ, M., Glaser, B. and Huwe, B: Uncertainty in the spatial prediction of soil texture. Comparison of regression tree and Random Forest
models, Geoderma, doi:https://doi.org/10.1016/j.geoderma.2011.10.010, 2012, 170, 70, 79.
Matheron, G. (1963). Principles of geostatistics. Economic Geology, 58(8), 1246–1266. https://doi.org/10.2113/gsecongeo.58.8.1246
McMahon, P. B., Böhlke, J. K., Kauffman, L. J., Kipp, K. L., Landon, M. K., Crandall, C. A., et al. (2008). Source and transport controls on
the movement of nitrate to public supply wells in selected principal aquifers of the United States. Water Resources Research, 44, W04401.
https://doi.org/10.1029/2007WR006252
McMahon, P. B., & Chapelle, F. H. (2008). Redox processes and water quality of selected principal aquifer systems. Ground Water, 46(2),
259–271. https://doi.org/10.1111/j.1745‐6584.2007.00385.x
Meinshausen, N. (2006). Quantile regression forests. Journal of Machine Learning Research, 67(3), 701–710. https://doi.org/10.1111/j.1541‐
0420.2010.01521.x
Møller, A. B., Iversen, B. V., Beucher, A., & Greve, M. H. (2017). Prediction of soil drainage classes in Denmark by means of decision tree
classification. Geoderma. https://doi.org/10.1016/j.geoderma.2017.10.015
Nolan, B. T., Fienen, M. N., & Lorenz, D. L. (2015). A statistical learning framework for groundwater nitrate models of the Central Valley,
California, USA. Journal of Hydrology, 531, 902–911. https://doi.org/10.1016/j.jhydrol.2015.10.025
Odeh, I. O. A., McBratney, A. B., & Chittleborough, D. J. (1995). Further results on prediction of soil properties from terrain attributes:
Heterotopic cokriging and regression‐kriging. Geoderma, 67(3–4), 215–226. https://doi.org/10.1016/0016‐7061(95)00007‐B
Pebesma, E. J. (2004). Multivariable geostatistics in S: The gstat package. Computational Geosciences, 30(7), 683–691. https://doi.org/
10.1016/j.cageo.2004.03.012
Pebesma, E. J., & Wesseling, C. G. (1998). Gstat: A program for geostatistical modelling, prediction and simulation. Computational
Geosciences, 24(1), 17–31. https://doi.org/10.1016/S0098‐3004(97)00082‐4
Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., et al. (2012). Scikit‐learn: Machine learning in Python.
Journal of Machine Learning Research. https://doi.org/10.1007/s13398‐014‐0173‐7.2
Petersen, R. J., Prinds, C., Iversen, B. V., Engesgaard, P., Jessen, S., & Kjaergaard, C. Nitrogen reduction along variable flow pathways in
riparian lowland transects. Submitt. to Water Resour. Res., 2018.
Ransom, K. M., Nolan, B. T., Traum, A. J., Faunt, C. C., Bell, A. M., Gronberg, J. A. M., et al. (2017). A hybrid machine learning model to
predict and visualize nitrate concentration throughout the Central Valley aquifer, California, USA. Science of the Total Environment,
601‐602, 1160–1172. https://doi.org/10.1016/j.scitotenv.2017.05.192
Refsgaard, J. C., Auken, E., Bamberg, C. A., Christensen, B. S. B., Clausen, T., Dalgaard, E., et al. (2014). Nitrate reduction in geologically
heterogeneous catchments—A framework for assessing the scale of predictive capability of hydrological models. Science of the Total
Environment, 468‐469, 1278–1288. https://doi.org/10.1016/j.scitotenv.2013.07.042
Refsgaard, J. C., Christensen, T. H., & Ammentorp, H. C. (1991). A model for oxygen transport and consumption in the unsaturated zone.
Journal of Hydrology, 129(1), 349–369. https://doi.org/10.1016/0022‐1694(91)90058‐P
Refsgaard, J. C., van der Sluijs, J. P., Højberg, A. L., & Vanrolleghem, P. A. (2007). Uncertainty in the environmental modelling process—A
framework and guidance. Environmental Modelling and Software, 22(11), 1543–1556. https://doi.org/10.1016/j.envsoft.2007.02.004
Robertson, W. D., Russell, B. M., & Cherry, J. A. (1996). Attenuation of nitrate in aquitard sediments of southern Ontario. Journal of
Hydrology, 180(1‐4), 267–281. https://doi.org/10.1016/0022‐1694(95)02885‐4
Rodriguez‐Galiano, V., Mendes, M. P., Garcia‐Soldado, M. J., Chica‐Olmo, M., & Ribeiro, L. (2014). Predictive modeling of groundwater
nitrate pollution using random forest and multisource variables related to intrinsic and specific vulnerability: A case study in an agri-
cultural setting (southern Spain). Science of the Total Environment, 476‐477, 189–206. https://doi.org/10.1016/j.scitotenv.2014.01.001
Rodriguez‐Galiano, V., Sanchez‐Castillo, M., Chica‐Olmo, M., & Chica‐Rivas, M. (2015). Machine learning predictive models for mineral
prospectivity: An evaluation of neural networks, random forest, regression trees and support vector machines. Ore Geology Reviews, 71,
804–818. https://doi.org/10.1016/j.oregeorev.2015.01.001
Shen, C. (2018). A trans‐disciplinary review of deep learning research and its relevance for water resources scientists. Water Resources
Research, 54(11), 8558–8593. https://doi.org/10.1029/2018WR022643
Stenger, R., Barkle, G., Burgess, C., Wall, A., & Clague, J. (2008). Low nitrate contamination of shallow groundwater in spite of intensive
dairying: The effect of reducing conditions in the vadose zone‐aquifer continuum. Journal of Hydrology. New Zealand, 47(1), 1–24.
Szatmári, G., & Pásztor, L. (2018). Comparison of various uncertainty modelling approaches based on geostatistics and machine learning
algorithms. Geoderma, 337, 1329–1340. https://doi.org/10.1016/J.GEODERMA.2018.09.008
Tesoriero, A. J., Gronberg, J. A., Juckem, P. F., Miller, M. P., & Austin, B. P. (2017). Predicting redox‐sensitive contaminant concentrations
in groundwater using random forest classification. Water Resources Research, 53, 7316–7331. https://doi.org/10.1002/2016WR020197
Tesoriero, A. J., Terziotti, S., & Abrams, D. B. (2015). Predicting redox conditions in groundwater at a regional scale. Environmental Science
& Technology, 49(16), 9657–9664. https://doi.org/10.1021/acs.est.5b01869
Thayalakumaran, T., Bristow, K. L., Charlesworth, P. B., & Fass, T. (2008). Geochemical conditions in groundwater systems: Implications
for the attenuation of agricultural nitrate. Agricultural Water Management, 95(2), 103–115. https://doi.org/10.1016/j.agwat.2007.09.003
Vaysse, K., & Lagacherie, P. (2017). Using quantile regression forest to estimate uncertainty of digital soil mapping products. Geoderma,
291, 55–64. https://doi.org/10.1016/j.geoderma.2016.12.017
Viscarra Rossel, R. A., Webster, R., & Kidd, D. (2014). Mapping gamma radiation and its uncertainty from weathering products in a
Tasmanian landscape with a proximal sensor and random forest kriging. Earth Surface Processes and Landforms, 39(6), 735–748. https://
doi.org/10.1002/esp.3476
Youssef, A. M., Pourghasemi, H. R., Pourtaghi, Z. S., & Al‐Katheeri, M. M. (2016). Landslide susceptibility mapping using random forest,
boosted regression tree, classification and regression tree, and general linear models and comparison of their performance at Wadi
Tayyah Basin, Asir Region, Saudi Arabia. Landslides, 13(5), 839–856. https://doi.org/10.1007/s10346‐015‐0614‐1