Abstract

Sediment thickness and bedrock topography are vital for the terrestrial hydrosphere. In this study, we estimated sediment thickness by using information from digital elevation models, geological maps, and public databases. We discuss two different approaches: First, the horizontal distances to the nearest bedrock outcrop were used as a secondary function in kriging and cokriging. Second, we applied Poisson's equation to estimate the local trend of the sediment thickness where bedrock outcrops were used as boundary conditions. Differences between point observations and the parabolic surface from Poisson's equation were minimized by inverse modelling. Ordinary kriging was applied to the residuals. These two approaches were evaluated with data from the Øvre Eiker, Norway. Estimates derived from Poisson's equation gave the smallest mean absolute error, and larger soil depths were reproduced better if the local trend was included in the estimation procedure. An independent cross-validation was undertaken. The results showed the best accuracy and precision for kriging on the residuals from Poisson's equation. Solutions of Poisson's equation are sensitive to the boundary conditions, which in this case were locations of the bedrock outcrops. Bedrock outcrops are available for direct observations; hence, the quality of the estimates can be improved by updating input from high-resolution mapping.

HIGHLIGHTS

  • Estimation of bedrock topography.

  • Spatial distribution of sediment thickness.

  • Use of Poisson's equation and kriging.

  • Utilization of the GRANADA well database.

  • Use of digital elevation models and Quaternary maps.

INTRODUCTION

Questions are put forward as to whether current conceptual models in hydrology are able to predict the impact of environmental change on the terrestrial hydrosphere (Beven 2006; Wagner 2007; Gupta et al. 2012). A fundamental part of the terrestrial environment is the interface between the atmosphere and the unweathered bedrock – labelled the Critical Zone (CZ) in Earth Science (NRC 2001; Giardino & Houser 2015; Richardson 2017). Water saturation in the CZ is a key variable for rainfall-runoff processes (Dunne & Black 1970; Dunne 1978). Evaporation and transpiration from this domain are important for the terrestrial water flux of moisture to the atmosphere (Trenberth 1999; Bierkens & van den Hurk 2008). A major part of the CZ is the unconsolidated material above the bedrock – for simplicity, called the sediment thickness in this paper. The spatial distribution of the sediment thickness and the bedrock topography is interesting in many different applications. In this article, we do not focus on processes in the CZ, but we address different approaches to estimate the geometry of sediment thickness and bedrock topography by utilizing publicly available databases. Before looking more closely into the estimation problem, it is interesting to review briefly some studies that focus on the role of bedrock topography and sediment thickness for rainfall-runoff processes.

Graham et al. (2010) and Graham & McDonnell (2010) used field experiments and modelling results to demonstrate the importance of bedrock topography and permeability for threshold response in runoff events. Above a threshold discharge, the runoff increased dramatically with a given rainfall intensity. In such cases, a major fraction of the subsurface flow took place at the interface between the sediments and the less permeable bedrock. Hopp & McDonnell (2009) used a synthetic dataset to show the importance of surface topography, soil depth, and bedrock topography to simulate the rainfall-runoff response. They concluded that the subsurface connectivity of saturated areas explained the runoff characteristics. Gabrielli et al. (2012) studied the relation between runoff and groundwater in sediments and bedrock in two hillslope catchments in New Zealand and Oregon, USA. In both hillslopes, the bedrock influenced the streamflow, but rapid subsurface stormflow was apparently most important if the bedrock had low permeability compared with the soil. Based on a large dataset, Tromp-van Meerveld & McDonnell (2006a) showed that pipe flow, matrix flow, and total flow had a threshold response to precipitation. Total flow increased two orders of magnitude if precipitation increased to a specified threshold. From their analysis, they concluded that bedrock topography determined the upstream contributing area, not the surface area as usually assumed. In the following study, Tromp-van Meerveld & McDonnell (2006b) showed that micro depressions in the bedrock surface were filled before they contributed to subsurface runoff. Freer et al. (2002) studied the relation between the hydrological response of a storm event and the bedrock topography in the Panola Mountain Research Watershed in Georgia, USA. They concluded that in cases where the bedrock had low permeability, the bedrock topography was important for the groundwater level, which in turn explained the large-scale hydraulic connectivity of the subsurface. Connectivities of macrostructures are important for the estimation of peak flow and recession of water in the catchment. Several other studies also elucidate the importance of the bedrock topography to understand hydrological and geochemical response (Genereux et al. 1993a, 1993b; Hinton et al. 1993; Mulholland 1993; Jencso et al. 2009). In summary, these studies indicate relations between runoff response, sediment properties and bedrock topography that deserve some more attention.

Methods for mapping and/or estimation of sediment thickness or bedrock topography are determined by the purpose of the study. Test drilling yields usually the ultimate point precision, but for resource mapping or hazards assessments, geophysical prospecting methods might give the required precision. In the case of water balance, integrated grid estimates might be sufficient. A thorough review of geotechnical and geophysical methods are beyond the scope of this study, but a short summary of the most relevant estimation approaches is appropriate.

Rempe & Dietrich (2014) suggested to model the general trend of sediment thickness for a hillslope by the analytical solution of differential equations where the transport of sediments was expressed as a diffusive process. The assessment of sediment budget has been explored in the fluvial and glacial environment (Walling & Collins 2008; Bogen 2017) and explored for large mountain areas (Hinderer 2001; Otto et al. 2009; Straumann & Korup 2009), but the aims of such works are usually lumped estimates of large-scale processes and not estimates of the sediment thickness as a function in space. Sulebak & Hjelle (2003) applied Spline Algorithms (SAs) to model landscape geometry with different spatial resolutions. The spline methods are robust smoothing algorithms, and the methods work, therefore, best for areas with gentle variations in the bedrock topography. Clarke et al. (2009) applied an Artificial Neural Network (ANN) to estimate glacial bedrock topography. They evaluated the method by using synthetic input data where the (synthetic) ice thickness was calculated by solving the shallow ice equations. The spatial reproduction of the bedrock topography had a maximum misfit of more than 400 m with a root-mean-square error of 70 m. The relative error of ice volume varied from −45 to +220%. Mey et al. (2015) applied ANN in the same way as Clarke et al. (2009). They used the digital elevation model (DEM) to estimate the sediment thickness in the Rhone valley and reported an average estimation error in sediment fill thickness of less than 27%. The ANN method assumes a ‘sameness’ in the training data and the output data. The assumption implies that the geomorphic signature of the bedrock topography is similar above and below the elevation where the bedrock is covered by sediments or ice. The ANN approach is, therefore, sensitive to the training dataset with respect to statistical homogeneity. A practical challenge arises in cases where the (homogeneous) training data is minor compared with the area of estimation. Neither SA nor ANN has conditional uncertainties as part of the methodology. From a practical point of view, this is a drawback, because in many cases, the uncertainty is equally important as the most likely value.

In this project, we used DEMs, digital Quaternary maps (DQMs), and the Norwegian well database GRANADA (NGU 2020a) to estimate the sediment thickness. The purpose was to evaluate ordinary methods to estimate the bedrock topography in areas covered by sediments. This aim was achieved by using standard methods in geostatistics (kriging and cokriging) and inverse modelling of Poisson's equation. Kriging and cokriging are well documented in textbooks (Isaaks & Srivastava 1989; Journel & Huijbregts 1989; Deutsch & Journel 1998), and a review of these methods is, therefore, redundant. The application of Poisson's equation, however, requires some more attention. The equation describes a steady-state solution of a diffusive process, which can be derived from the continuity of mass and Fickian diffusion. The evolution of landscape shaped by fluvial processes can be described as a diffusive process. The idea can be traced back to Gilbert (1877) who stated that ‘… erosion is most rapid where the slope is steepest’. According to Tucker & Hancock (2010), Gilbert's observation is expressed in an equation by Culling (1960):
formula
(1)
where is the sediment flux, is a constant that includes all properties involved in sediment transport (namely physical properties of the material, climate), and is the gradient of the landscape. Combined with the balance of mass, this gives Boussinesq's equation, which in steady-state is equivalent to Poisson's equation. This approach is widely applied in landscape modelling (Kooi & Beaumont 1994; Braun & Sambridge 1997; Tucker & Hancock 2010; Kitterød 2011). In our case, the bedrock topography is mainly governed by glacial erosion, and the sediment thickness is the result of both glacial and fluvial processes. Empirical equations indicate that glacial erosion is proportional to , where u is the glacial flow velocity. The power term p is usually between 1 and 2 (MacGregor et al. 2009; Herman et al. 2015; Cook et al. 2020). Harbor et al. (1988) applied a numerical model to simulated glacial erosion and calculated a time span of years to develop U-shaped valleys. The glaciation of Scandinavia lasted for more than years, and U-shaped valleys are, thus, common landscape features in this region (Mangerud et al. 2011). Thus, there are physical reasons to model fluvial and glacial landscapes as parabolic functions in a geostatistical framework.

On this background, Poisson's equation offers some interesting properties (Kazhdan et al. 2006; Calakli & Taubin 2011; Liu 2018). First, the solution of Poisson's equation yields parabolic functions in space, which is readily adaptable to U-shaped valleys. Second, Poisson's equation is sensitive to the boundary conditions. This sensitivity can be utilized if the altitude is known. Today, DEMs are available for large parts of the world, thus utilizing DEMs for this purpose is an interesting application. Third, in areas covered by sediments, inverse modelling of Poisson's equation can be used to fit a parabolic surface to the point observations of the sediment thickness.

Kazhdan et al. (2006) applied the numerical solutions of Poisson's equation as a filtering method for surface reconstruction. A conventional application of Poisson's equation assumes that the reconstruction domain is minor compared with the surface that is exposed to observations. The opposite is more common in geo-applications: The unknown domain is usually large compared with the exposed bedrock. One way to overcome this problem is to optimize the source term in Poisson's equation by inverse modelling. The residuals between observed sediment thickness and the optimal parabolic function can be modelled as a stochastic function in space.

The bedrock topography (m a.m.s.l.) is given by
formula
(2)
where (m a.m.s.l.) is the terrain elevation, and (m) is the estimate of the sediment thickness. The challenge is to find that minimize the estimation uncertainty . To do that, it was necessary to find a local trend of the requested variable. Two different approaches were applied to model the local trend of the sediment thickness (m): First, we used relations between and the horizontal distance to the nearest bedrock outcrop (m). Second, we applied inverse modelling of Poisson's equation to adapt a parabolic function to the sediment thickness or the bedrock topography. The boundary condition was obtained by defining an indicator function . In location u, where the bedrock was exposed to the atmosphere , if the bedrock was covered by deposits, then . In the case study presented below, was estimated by using point information from a public well database, and locations of bedrock outcrops from geological maps, but all available information can, in principle, be utilized (namely point information from geotechnical test drilling, line information from geophysical mapping, area information from constructions).

Using in the estimation procedure was inspired by an intuitive procedure among drill operators: Based on previous experiences from a given area, it is often possible to judge the sediment thickness D at a location u, by simple visual inspection of the horizontal distance L, to the nearest bedrock outcrop (Figure 1). The statistical relation between and was used in cokriging and ordinary kriging to evaluate the effect on the estimation variance.

Figure 1

Two principles for including a local trend in the estimation of the bedrock topography. First, by using the relation between the sediment thickness , and the horizontal distance to the nearest bedrock outcrop , where is the locations of the observations. Second, by fitting a parabolic curve to the bedrock topography (red curve) by solving Poisson's equation. Boreholes (white rectangles) penetrate the sediments into the bedrock, while sedimentary wells (blue rectangle) do not penetrate into the bedrock (modified from Kitterød (2017)). Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/nh.2021.102.

Figure 1

Two principles for including a local trend in the estimation of the bedrock topography. First, by using the relation between the sediment thickness , and the horizontal distance to the nearest bedrock outcrop , where is the locations of the observations. Second, by fitting a parabolic curve to the bedrock topography (red curve) by solving Poisson's equation. Boreholes (white rectangles) penetrate the sediments into the bedrock, while sedimentary wells (blue rectangle) do not penetrate into the bedrock (modified from Kitterød (2017)). Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/nh.2021.102.

The quality of publicly available data is interesting to explore because the volume of information is increasing for every year. In this study, was taken from the Norwegian well database GRANADA (NGU 2020a). Data from mainland Norway (denoted the national- or large scale in the following) were used to estimate covariance- and cross-covariance functions of and . Local boreholes were used for conditional estimation and inverse modelling of the constant parameter in Poisson's equation.

The geostatistical methods we evaluate are available in open source and commercial software, but it is not standard procedure to include the bedrock outcrop as a secondary function in kriging or as boundary conditions in Poisson's equation. Most of the algorithms we applied were, therefore, tailored to the data material we used.

Two different types of cross-validation were performed: First, cross-validation was undertaken by leaving one observation out. In this procedure, all available point observations were included. The second cross-validation was done by dividing the observation material into two groups. The first group contained boreholes recorded before 2010, and the second group was sedimentary wells and boreholes recorded after 2010. Because the second cross-validation was done on independent data not included in the parameter inference, this procedure was called jackknife cross-validation.

The results were evaluated according to three simple criteria: The mean absolute error , accuracy , and precision (Goovaerts et al. 2005). In addition, we applied a simple probability score of observations and the continuous ranked probability score (CRPS; Gneiting & Raftery 2004). Before presenting the methods and the results, it is necessary to describe the data material we applied to evaluate the estimation procedures.

DATA MATERIAL

As indicated above, we employed three sources of information for the current study: (i) point observations of sediment thickness, (ii) DQM, and (iii) DEM of the surface. In addition, we also had access to geological maps of the bedrock for the study area. Even though the locations of the different bedrock types were not utilized in the current study, a brief description of the most prominent geological structures in the study area is interesting as background information.

Local geology

The study area Øvre Eiker (9 km 7.1 km) is located about 60 km west of Oslo (Figure 2). Øvre Eiker is a part of Vestfold graben, which is one of the three main graben structures in the Oslo Rift. Large parts of the Vestfold graben come from batholithic intrusions and consist mainly of biotite granite. At the highest altitudes in the eastern parts of the study area, the bedrocks belong to the Glitrevann caldera, which was formed as part of the plutonic intrusion in Permian time (299–252 mill. years BP). At the eastern side of the valley, the bedrocks are mainly shales, marbles, and limestones from Ordovicium (488–440 mill. years BP). The bedrocks at the western side of the valley are mainly quartzites (estimated deposition time about 1,475 mill. years BP) and different kinds of gneiss (estimated intrusion time about 1,500 mill. years BP). At large, the topography of the study area mirrors the character of the bedrock with hard resistant bedrock at high altitudes and softer bedrocks in the valleys. More information on the bedrock geology can be found in Andersen et al. (2008) and NGU (2020b).

Figure 2

Location of the local study area Øvre Eiker (7.2 9 km2), 60 km West of Oslo, Norway, to the left. To the right is a contour map of Øvre Eiker surface topography; locations of boreholes registered before 2010 (red dots); boreholes registered after 2010 (black dots); and sedimentary wells (yellow plus sign). Grey colour indicates areas with bedrock outcrops. [UTM-East, UTM-North] coordinates of the Øvre Eiker study area are [210262.50, 6637562.50] for the lower left corner and [219262.50, 6644687.50] for the upper right corner. Coordinates are given in UTM zone 33 references. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/nh.2021.102.

Figure 2

Location of the local study area Øvre Eiker (7.2 9 km2), 60 km West of Oslo, Norway, to the left. To the right is a contour map of Øvre Eiker surface topography; locations of boreholes registered before 2010 (red dots); boreholes registered after 2010 (black dots); and sedimentary wells (yellow plus sign). Grey colour indicates areas with bedrock outcrops. [UTM-East, UTM-North] coordinates of the Øvre Eiker study area are [210262.50, 6637562.50] for the lower left corner and [219262.50, 6644687.50] for the upper right corner. Coordinates are given in UTM zone 33 references. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/nh.2021.102.

Quaternary maps and surface elevation

DQMs (NGU 2020c) were used to identify locations where the bedrock was exposed to the atmosphere or covered with a thin and patchy layer of organic matter (37% of the study area). For simplicity, the area where is referred to as exposed bedrock in this article. The remaining area where (63%) was covered by the unconsolidated material (Figure 3). Marine deposits are the most frequent sediment type in the area, which cover about 27% of the study area. The marine limit was levelled to 194 m above present mean sea level (Bargel 1987). More than 70% of the surface area covered by sediments was below the marine limit, while more than 75% of the exposed bedrock was located above this altitude (Figure 3).

Figure 3

DQMs of Øvre Eiker (left figure) with boreholes registered before 2010 (red dots), sedimentary wells (black plus sign), and boreholes registered after 2010 (black dots). The hypsographic curves (right figure) indicate the distribution of surface cover as a function of altitude (m a.m.s.l.). The area fraction of the different surface categories was (0) bedrock (36.9%), (1) till and moraine sediments (12.9%), (2) glaciofluvial sediments (6.2%), (3) fluvial sediments (6.9%), (4) marine sediments (27.1%), (5) colluvial or weathered avalanches (8.1%), (6) peat and swaps (1.6%), and (7) anthropogenic material (0.14%). Percentages of the total area are given in parentheses. Coordinates are given in UTM zone 33 references, and distances in the East–West and the North–South directions are 9 and 7.1 km, respectively. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/nh.2021.102.

Figure 3

DQMs of Øvre Eiker (left figure) with boreholes registered before 2010 (red dots), sedimentary wells (black plus sign), and boreholes registered after 2010 (black dots). The hypsographic curves (right figure) indicate the distribution of surface cover as a function of altitude (m a.m.s.l.). The area fraction of the different surface categories was (0) bedrock (36.9%), (1) till and moraine sediments (12.9%), (2) glaciofluvial sediments (6.2%), (3) fluvial sediments (6.9%), (4) marine sediments (27.1%), (5) colluvial or weathered avalanches (8.1%), (6) peat and swaps (1.6%), and (7) anthropogenic material (0.14%). Percentages of the total area are given in parentheses. Coordinates are given in UTM zone 33 references, and distances in the East–West and the North–South directions are 9 and 7.1 km, respectively. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/nh.2021.102.

Till deposits are the second most frequent surface sediment, which are exposed to about 13% of the study area. Most of the tills are deposited as a relatively thin layer above the bedrock, and more than 95% of the exposed tills are located above the marine limit. For simplicity, we also included front moraine in the till category for the present study, which implies that the sediment thickness may be larger in some locations. Glaciofluvial sediments are deposited below 200 m a.m.s.l. and fluvial sediments below 120 m a.m.s.l. (Figure 3). About 8% of the study area is covered by colluvial or weathered material. Peat and swamps are most frequent at higher altitudes and cover less than 2% of the area. A minor fraction of the area (0.14%) is anthropogenic material.

Sediment thickness and the horizontal distance to bedrock outcrop

Point observations of sediment thickness for mainland Norway were taken from the public well database GRANADA (NGU 2020a). DQM was used to calculate the horizontal distance to the nearest bedrock outcrop , for all GRANADA recordings . Semivariogram (covariance and cross-covariance) model parameters used in the current study were based on the GRANADA data registered before 2010 from all over Norway. At that time, the total number of reported boreholes with > 0 were = 19,682 (Kitterød 2017).

GRANADA recordings after 2010 located within the Øvre Eiker study area were included in this study, for primary cross-validation by leaving one observation out. A second cross-validation was done by dividing the dataset into two groups. This procedure is called jackknife cross-validation because it was done on an independent dataset. In addition to point observations of , we used to calculate in a regular grid with a spatial resolution of = 25 m 25 m. The purpose was to evaluate whether modelling results were improved or not if local information of was included as a secondary function. The local relation between and with current available data can be seen in Figure 4.

Figure 4

GRANADA recordings of boreholes and sedimentary wells in the Øvre Eiker study area (cf. Figure 2). The figure shows the relation between the horizontal distance to the nearest bedrock outcrop , and sediment thickness D for boreholes that penetrate the deposits and goes into the bedrock. For sedimentary wells (blue circles), the sedimentary thickness is not recorded, only the length of the well, thus for these wells D is equal to the length of the well, which means that the true sediment thicknesses for these wells are either equal to D shown on the figure or larger. Boreholes recorded before 2010 were used for cross-validation by leaving one well out. After this primary cross-validation, percentiles were conditioned on all boreholes recorded before 2010. Boreholes recorded after 2010 (red crosses) and sedimentary wells (blue circles) were used for jackknife cross-validation. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/nh.2021.102.

Figure 4

GRANADA recordings of boreholes and sedimentary wells in the Øvre Eiker study area (cf. Figure 2). The figure shows the relation between the horizontal distance to the nearest bedrock outcrop , and sediment thickness D for boreholes that penetrate the deposits and goes into the bedrock. For sedimentary wells (blue circles), the sedimentary thickness is not recorded, only the length of the well, thus for these wells D is equal to the length of the well, which means that the true sediment thicknesses for these wells are either equal to D shown on the figure or larger. Boreholes recorded before 2010 were used for cross-validation by leaving one well out. After this primary cross-validation, percentiles were conditioned on all boreholes recorded before 2010. Boreholes recorded after 2010 (red crosses) and sedimentary wells (blue circles) were used for jackknife cross-validation. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/nh.2021.102.

METHODS

The statistical methods we applied for this study were based on the Gaussian theory, which is well documented in textbooks (Isaaks & Srivastava 1989; Journel & Huijbregts 1989; Deutsch & Journel 1998). The equations we apply will, therefore, not be reproduced here, but some more details are given in the Supplementary Appendix.

Experimental semivariograms

Ordinary kriging (), cokriging (), and bedrock kriging () of sediment thickness in the study area were based on semivariogram and cross-semivariogram parameters derived from the national GRANADA database. The estimation parameters should, therefore, be considered as average values for mainland Norway (Kitterød 2017). These estimation parameters are not stationary in space, and local estimation parameters should, in principle, give better results. In this case, the derivation of local parameters was not feasible due to the scarcity of point observations in the study area. For this reason, we preferred to perform local kriging with average large-scale semivariogram model parameters instead of less reliable local model parameters. Before the modelling took place, the GRANADA data required some pre-processing.

Boreholes and wells are usually located in populated areas, which means that most of the GRANADA recordings were from urban areas, while fewer recordings were from rural areas. If such data are used to calculate average properties, uneven sampling will affect the statistics because oversampled domains will dominate. To suppress such clustering effects, each observation received a weight calculated as a function of the distance to other boreholes. Boreholes located close to each other were given less weight than boreholes located far from each other. The calculation of declustering weights was done according to the moving grid algorithm suggested by Deutsch (1989) and applied in Kitterød (2017). It should be noted that declustering weights do not change the value of the observations that we used for conditional estimation. The weights were only applied for the calculation of experimental semivariograms and statistical moments.

We applied two methods to approach Gaussian distributions of the involved variables. For ordinary kriging () and cokriging (), we utilized the normal score transform (Supplementary Appendix A1), while bedrock kriging () was done on logarithmic values. For the sake of completeness, was also performed on logarithmic values. After calculations, the results were back-transformed to sediment depth in metric units for cross-validation, either by the inverse normal score transform or the inverse of the lognormal cumulative distribution function.

Ordinary kriging and cokriging

Expected sediment thickness and estimation variance were obtained by solving the ordinary kriging () and the cokriging () equations. The and equations were expressed in matrix notation and solved simultaneously (cf. Supplementary Appendix A). Several methods exist for including secondary information in the estimation procedure (Deutsch & Journel 1998; Goovaerts et al. 2005). Here, we suggest a method called bedrock kriging (), which is less formal than cokriging. The method is a two-step procedure. First, was done on , and then second on the relation:
formula
(3)
where D (m) is the sediment thickness, L (m) is the horizontal distance to bedrock outcrop, and is the point observations. The motivation for using the simple relation in Equation 3 is that information of is given a priori in every point of interest from the DQMs. See Supplementary Appendix A3 and A4 for more details.

Poisson's equation

Inverse modelling of Poisson's equation P was applied to fit a two-dimensional parabolic surface to the observations. A general version of Poisson's equation is written as follows:
formula
(4)
where denotes the Laplace operator and is the source function where spans the total domain of observations (e.g. mainland Norway).

Three types of boundary conditions can be applied to solve Equation 4: (i) a Dirichlet condition, where elevation of the bedrock surface or soil depth is given; (ii) a Neumann condition, where the gradient is given; or (iii) a Cauchy condition where a weighted average of (i) and (ii) is given (Gosses et al. 2018; Liu 2018). Here, in this study, we applied the Dirichlet condition in two slightly different ways. Because in outcrop locations where the bedrock surface is identical to the terrain , we have two options: First, the topography elevation (m a.m.s.l.) can be used directly as boundary conditions: . In that case, is identical to the local trend in the bedrock topography: , where is given by Equation 6. Second, the soil depth (m) can be used as boundary condition because by definition in outcrop locations , thus: . In that case, is identical to the local trend of the sediment thickness , and the local trend of the bedrock topography is given by: .

In the case study presented below, we fitted , , where is the Øvre Eiker study area. Since , we approximated the right-hand side of Equation 4 to a constant:
formula
(5)
In this context, inverse modelling of means to find a that minimizes the difference between point observations and , which is the solution of Equation 5:
formula
(6)
where is given in Equation 5. First, was calculated by leaving one observation out: , where all observations were included ( = 35), and second for jackknife cross-validation, where only boreholes recorded before 2010 were used ( = 10). This procedure indicates the sensitivity of depending on the left-out observations. Inverse modelling of with all observations included yields a parabolic function with minimum deviation from the observations. The advantage is that the bedrock morphology is honoured directly by using the outcrop locations as boundary conditions. In this way, the complete geometry of the bedrock outcrop is included in the estimation procedure, not only the shortest distance to the bedrock outcrop as in cokriging and bedrock kriging.

We applied an explicit finite difference numerical scheme of second order to solve Equation 5. The explicit scheme implemented the boundary conditions very simply, and the matrix multiplication in MATLAB gave fast convergence to stable solutions. The evaluation of the numerical scheme was beyond the scope of this study, but no numerical instability was revealed, and high numerical precisions were achieved by feasible consumption of CPU time.

To obtain conditional mean and estimation variance, we did ordinary kriging of the residuals between the observations and the optimal solution of Poisson's equation:
formula
(7)
where is the point observation, is the results of inverse modelling (Equation 6), and n is the number of local observations. Ordinary kriging of the residuals from optimal solution of Poisson's equation is called Poisson kriging () in the following. In the procedure, we also included the boundary locations as conditional observations: , where is the locations of the bedrock outcrops.

Evaluation criteria

To keep the model evaluation as simple and consistent as possible, we decided to present the uncertainties in terms of probability density functions by calculation of percentiles (Supplementary Appendix A5).

The modelling results can be expressed as follows:
formula
(8)

The subscript m denotes the estimation model , where is ordinary kriging, is cokriging, is bedrock kriging, is Poisson kriging, and where u is the estimation location, is the percentiles, and denotes all available information.

Three criteria were employed for evaluation of the modelled uncertainty: the mean absolute error , the accuracy , and the precision (Goovaerts et al. 2005; Kitterød 2017). The mean absolute error is the difference between the observed variable and the median of the estimated variable :
formula
(9)
where n is the number of observations used in cross-validation.
The accuracy denotes whether the observed estimate is within two predefined percentile intervals around . For the current case study, we used: . If the observed value is within the predefined percentile interval, then , if it is outside, then , thus:
formula
(10)

By this definition, includes the model uncertainty, which means that in locations with less information and thus larger uncertainty, the model can be accurate (i.e. the estimates are between the given percentile intervals), but not very precise because high uncertainties yield large percentile intervals. An additional criterion for precision is, therefore, necessary.

For this purpose, precision was defined by taking the difference between upper and lower percentiles of the accuracy estimates and sum the inverse difference. By this definition, the modelled precision reads:
formula
(11)
where k is the number of estimates within the predefined accuracy criteria. For the case study presented below: . Hence, a precise estimate has large values because the estimated percentiles are close to the mean value (i.e. small estimation uncertainty). If two modelling results have the same and , the model with the highest is superior.

Gneiting & Raftery (2004) introduced CRPS to measure the success of probabilistic forecasts and estimates. Mean CRPS was calculated by linear interpolation in the MATLAB script provided by Shrestha (2020). CRPS values were calculated for each cross-validation point for 10,000 realizations.

In addition to the criteria above, we calculated a more simple probability score of the observations :
formula
(12)
where is the location of the cross-validated observation, and is the cumulative probability density function (cdf) of , where is the expected sediment thickness, is the estimation uncertainty, and n is the number of observations.

Finally, linear regression was included to indicate the deviation between cross-validated estimates and observations. Linear regression cannot be regarded as a formal evaluation method in this case because the assumption of independence of residuals (i.e. the difference between the regression model and the estimates) are not fulfilled. The correlation coefficients R should, therefore, be interpreted with caution. The slope coefficients () of the linear regression model indicate underestimation (1) or overestimation (1). A successful estimation method should have close to 1.

RESULTS

The modelling results were evaluated by comparing the solutions with respect to the evaluation criteria explained above. This was first done in a primary cross-validation procedure by leaving one observation out, and second, by jackknife cross-validation where the GRANADA observations were divided into two groups: Boreholes recorded before 2010, and sedimentary wells and boreholes recorded after 2010. The first group was used to calculate conditional estimates in the locations of the second group. For all cases presented below, we applied a grid resolution of m.

Semivariogram model parameters

The national GRANADA data were used to calculate experimental semivariograms for ordinary kriging (), cokriging (), and bedrock kriging (). Optimal parameters were obtained by using the simulated annealing algorithm (Kitterød 2017; MATLAB 2020). Semivariogram for Poisson kriging () was based on GRANADA boreholes located within the Øvre Eiker study area. The experimental semivariograms from logarithmic values are shown in Figure 5. For and , the semivariograms were derived from normal score transformed data (cf. Tables 3 and 4, case F in Kitterød (2017)). Estimated model parameters are given in Table 1.

Table 1

Covariance and cross-covariance model parametersa used in ordinary kriging (), cokriging (), bedrock kriging (), and kriging on residuals from Poisson's equation (). The observations were either normal score (*) or lognormal transformed (′)

Case
  0.077 0.726 2371 1.00 
  0.003 0.207 2402 1.01 
  0.003 0.207 2402 1.01 
  0.007 0.638 5865 1.02 
  0.35 0.750 2100 1.0 
  0.20 0.520 8000 1.0 
  20 200 4000 1.0 
Case
  0.077 0.726 2371 1.00 
  0.003 0.207 2402 1.01 
  0.003 0.207 2402 1.01 
  0.007 0.638 5865 1.02 
  0.35 0.750 2100 1.0 
  0.20 0.520 8000 1.0 
  20 200 4000 1.0 

a, and the semivariogram model are given by , where indicates either normal score or lognormal transformed sediment thickness , or the horizontal distance to the bedrock outcrop , , and X denotes the residuals between Poisson's equation and observations (Equation 7). The practical range = log(0.05) for all models.

Figure 5

Semivariograms for to the left and to the right, where D (m) is the sediment thickness, and L (m) is the horizontal distance from the location of the borehole to the location of the nearest bedrock outcrop. Data source: GRANADA database (NGU 2020a). The total number of observations with was 19,682. The parameters in the semivariogram model is the nugget , sill , range a (m), exponential model , and the practical range , in this case 95% of the total variance: .

Figure 5

Semivariograms for to the left and to the right, where D (m) is the sediment thickness, and L (m) is the horizontal distance from the location of the borehole to the location of the nearest bedrock outcrop. Data source: GRANADA database (NGU 2020a). The total number of observations with was 19,682. The parameters in the semivariogram model is the nugget , sill , range a (m), exponential model , and the practical range , in this case 95% of the total variance: .

Inverse modelling of Poisson's equation

Poisson's equation was applied for trend analysis of the local bedrock topography. This was done by inverse modelling of the constant parameter in Equation 5, which means that we search for a that minimizes the mean absolute error, (Equation 9). Since only one parameter was involved, this was obtained by simple incremental stepping of the factor. The sensitivity of was evaluated by successively leaving one observation out and, finally, by including all observations. This procedure was repeated twice. First with all local observations included ( = 35), and second for jackknife cross-validation where the observations were limited to boreholes recorded before 2010 ( = 10). The results show that a minimum was possible to obtain independent on the dataset included, but the absolute value varied between 10 and 11.5 m with all data included, and between 2 and 3 m for boreholes recorded before 2010 (Figure 6). The optimal load parameter varied between 4.0 × 10−5 and 4.5 × 10−5 for all data included and between 4.5 × 10−5 and 6 × 10−5 for boreholes recorded before 2010.

Figure 6

Inverse modelling of Poisson equation (Equation 6) by incremental stepping of the constant factor in Equation 5. One observation is left out successively to estimate sensitivity of . Two inverse modelling results are shown. (a) All observation points within the study area were included () for the optimal solution (magenta curve). In addition, inverse modelling results from three wells are shown (GRANADA id.: 31581, 54057, 4928) to visualize extremes and variance of optimal load parameters ([4.0 × 10−5, 4.5 × 10−5]). (b) The results from boreholes recorded before 2010 were shown (10). Blue curve is the inversion results where all boreholes before 2010 were included. The optimal load parameter 4.5 × 10−5 if all observations were included (magenta curve to the left) but was increased to 5.5 × 10−5 if the observations were limited to boreholes recorded before 2010 (blue curve to the right). Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/nh.2021.102.

Figure 6

Inverse modelling of Poisson equation (Equation 6) by incremental stepping of the constant factor in Equation 5. One observation is left out successively to estimate sensitivity of . Two inverse modelling results are shown. (a) All observation points within the study area were included () for the optimal solution (magenta curve). In addition, inverse modelling results from three wells are shown (GRANADA id.: 31581, 54057, 4928) to visualize extremes and variance of optimal load parameters ([4.0 × 10−5, 4.5 × 10−5]). (b) The results from boreholes recorded before 2010 were shown (10). Blue curve is the inversion results where all boreholes before 2010 were included. The optimal load parameter 4.5 × 10−5 if all observations were included (magenta curve to the left) but was increased to 5.5 × 10−5 if the observations were limited to boreholes recorded before 2010 (blue curve to the right). Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/nh.2021.102.

Primary cross-validation

All available point observations of sediment thickness within the study area were used for primary cross-validation analysis. This cross-validation was performed for observations by leaving successively one observation out and using the remaining boreholes for conditional modelling of the cdf at the location of the left-out borehole: (). For ordinary kriging (), only observations of sediment thickness were used for conditional modelling, while for cokriging () and bedrock kriging (), the secondary variable was used also in the location of the left-out borehole. A DQM was used to identify grid locations u where bedrocks were exposed to the atmosphere or covered by deposits (Figure 7).

Figure 7

Gridded map for Øvre Eiker (cf. Figure 2) showing areas with bedrock outcrops and areas where the bedrock is covered by deposits . Grid resolution: m. UTM(33) coordinates for lower left and upper right corner are [210262.50, 6637562.50] and [219262.50, 6644687.50].

Figure 7

Gridded map for Øvre Eiker (cf. Figure 2) showing areas with bedrock outcrops and areas where the bedrock is covered by deposits . Grid resolution: m. UTM(33) coordinates for lower left and upper right corner are [210262.50, 6637562.50] and [219262.50, 6644687.50].

The primary cross-validation results in Figure 8 show that, in general, the results overestimated observations smaller than the lumped mean sediment thickness and underestimated observations larger than the lumped mean. Most of the estimates were around one standard deviation from the lumped mean observation. The and estimates gave slightly better results, but the best result was obtained from estimates (i.e. kriging on residuals from Poisson's equation, Equation 7). The optimal solution of Poisson's equation (Equations 5 and 6) gave significant deviations from observations (black squares in Figure 8).

Figure 8

Cross-validations by leaving one observation out with all point observations in the study area (35). The median estimation values () are shown for ordinary kriging (; blue and black dots), bedrock kriging (; red dots), cokriging (; magenta dots), and kriging on residuals from Poisson's equation (; open red circles). Optimal solutions of Poisson's equation are indicated by black squares (Equations 5 and 6). Results based on lognormal transformed observations were indicated by ′ and normal score transformed observations by *. The solid line is a 1:1 relation between estimates and observations. Dashed lines indicate percentiles : [0.05, 0.25, 0.75, 0.95] from the cross-validation results. Lumped mean and standard deviation of observed sediment thickness were 14.3 and 15.6 m, respectively. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/nh.2021.102.

Figure 8

Cross-validations by leaving one observation out with all point observations in the study area (35). The median estimation values () are shown for ordinary kriging (; blue and black dots), bedrock kriging (; red dots), cokriging (; magenta dots), and kriging on residuals from Poisson's equation (; open red circles). Optimal solutions of Poisson's equation are indicated by black squares (Equations 5 and 6). Results based on lognormal transformed observations were indicated by ′ and normal score transformed observations by *. The solid line is a 1:1 relation between estimates and observations. Dashed lines indicate percentiles : [0.05, 0.25, 0.75, 0.95] from the cross-validation results. Lumped mean and standard deviation of observed sediment thickness were 14.3 and 15.6 m, respectively. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/nh.2021.102.

A summary of the cross-validation results by leaving one observation out is given in Table 2. The results gave the mean absolute error equal to 6.51 m for normal score transformed data, and 7.46 m for lognormal transform. for and were 6.47 m and 7.26 m, respectively, and 11.30 m for inverse modelling of Poisson's equation () and 2.15 m for kriging on the residuals from Poisson's equation (). Least accuracy was obtained for estimates on log-transformed variables. The estimates had the highest accuracy, and 86% of the cross-validated observations were within predefined percentiles. The precision was a bit lower than the estimates, but the precision is calculated only for the accurate estimates, which was 61% for the estimates.

Table 2

Summary of cross-validation by leaving one borehole out. The number of point observations in the study area was 35. Cross-validations were undertaken by ordinary kriging (), cokriging (), bedrock kriging (), and Poisson equation () with an optimal load parameter estimated by leaving one observation out. Kriging on residuals between observations and Poisson equation () were done for the Poisson solution where the optimal parameter was estimated based on all observations. The superscripts indicate either normal score transform (*) or lognormal transform (′) of observations. In addition to the evaluation criteria , , , , and (see references below), the correlation coefficient R and the slope coefficient from linear regression is included.

 6.51 6.47 7.46 7.26 11.30 2.15 
 0.58 0.61 0.48 0.58 0.49 0.86 
 1.73 4.51 1.41 1.99 0.92 3.97 
mean (0.48 0.47 0.51 0.49 0.66 0.57 
std (0.25 0.28 0.30 0.29 0.22 0.17 
mean (12.10 12.15 5.67 5.95 8.64 1.98 
std (16.01 15.99 10.98 9.55 12.47 1.12 
 0.55 0.59 0.48 0.51 0.40 0.99 
 0.31 0.37 0.23 0.25 0.19 0.89 
 6.51 6.47 7.46 7.26 11.30 2.15 
 0.58 0.61 0.48 0.58 0.49 0.86 
 1.73 4.51 1.41 1.99 0.92 3.97 
mean (0.48 0.47 0.51 0.49 0.66 0.57 
std (0.25 0.28 0.30 0.29 0.22 0.17 
mean (12.10 12.15 5.67 5.95 8.64 1.98 
std (16.01 15.99 10.98 9.55 12.47 1.12 
 0.55 0.59 0.48 0.51 0.40 0.99 
 0.31 0.37 0.23 0.25 0.19 0.89 

, Mean absolute error (Equation 9); , Accuracy (Equation 10); , Precision (Equation 11); , Probability score (Equation 12); , Continuous ranked probability score (Gneiting & Raftery 2004); , Correlation coefficient between observations and estimates; , Slope coefficient estimated by linear regression.

The mean probability score for the cross-validated observations was close to 0.5 for all methods, but this criterion must be judged together with standard deviation which was smallest for the estimates. Mean values of the are similar to with smallest values for the estimates. The standard deviation of the values was also smallest for the estimates.

Even though a simple linear regression is not a formal evaluation method, the slope coefficients ( in Table 2) show an average underestimation between 0.23 and 0.31 for and methods, and 0.37 for . For the method, the underestimation was 0.89, i.e. the regression model gave on average estimates that were 11% less than the observations. The correlation coefficients ( in Table 2) cannot be interpreted as a ‘success factor’, but they indicate the deviation between a linear regression line and the estimates. In this case, the method had estimates closest to the regression line ( = 0.99).

Jackknife cross-validation

Conditional cdfs were calculated in all grid points within the study area for jackknife cross-validation. These cdfs were obtained by using the parameters listed in Table 1. The Poisson's equation was solved by using the optimal load parameter = 5.5 × 10−5 were all boreholes recorded before 2010 were used for inverse modelling (Figure 6(b)). All estimates were conditioned on local observations of from boreholes in the study area recorded before 2010 and by using the horizontal distance to the nearest bedrock outcrop () as secondary information (Figure 7). The results of the estimated cdfs were presented in terms of percentiles: given as fractional values.

Before presenting the jackknife cross-validation results, some characteristic features in Figure 9 should be noted. The results of and in Figure 9(a) and 9(c) show the local median () of the two variables. The estimates shown in Figure 9(e) were dominated by the recorded sediment depth close to the boreholes, while outside the range of the boreholes, the estimates are dominated by L. This effect on the estimates can be seen at the southwestern part of the image where L was large and where the distance to the nearest boreholes were larger than the range (Figure 9(e)). The results showed also the spatial structure of the secondary variable L, but less prominent than the estimates in areas that were distant to the borehole observations (Figure 9(c)). This can be seen in the west-east cross-section that goes through the two boreholes located south of the study area (Figure 10). In general, the uncertainty in the estimates was greater than the estimates. The results, however, show greater estimation variances than the estimates for locations that were close to the bedrock outcrops, but distant from the nearest D observation. The V-shape of the estimates was evident in the and the estimates, while estimates, on the other hand, had the characteristic U-shape as expected (Figure 9(f)).

Figure 9

Images of the Øvre Eiker area showing estimated median () values conditioned on boreholes recorded prior to 2010 (the number of boreholes were 10). (a) Sediment thickness (m) obtained by ordinary kriging (), (b) sediment thickness (m) by cokriging (), (c) (–) by ordinary kriging (), where (–) and L (m) is the horizontal distance to the nearest bedrock outcrop, (d) bedrock kriging weights (–) (Supplementary Appendix Equation A10), (e) sediment thickness (m) obtained by bedrock kriging (), and (f) sediment thickness (m) by the optimal solution of Poisson's equation (Equations 5 and 6). The dotted line through the two southernmost boreholes (id. 31581 and 31588) indicates the cross-section shown in Figure 10.

Figure 9

Images of the Øvre Eiker area showing estimated median () values conditioned on boreholes recorded prior to 2010 (the number of boreholes were 10). (a) Sediment thickness (m) obtained by ordinary kriging (), (b) sediment thickness (m) by cokriging (), (c) (–) by ordinary kriging (), where (–) and L (m) is the horizontal distance to the nearest bedrock outcrop, (d) bedrock kriging weights (–) (Supplementary Appendix Equation A10), (e) sediment thickness (m) obtained by bedrock kriging (), and (f) sediment thickness (m) by the optimal solution of Poisson's equation (Equations 5 and 6). The dotted line through the two southernmost boreholes (id. 31581 and 31588) indicates the cross-section shown in Figure 10.

Figure 10

East-West cross-sections of bedrock topography shown as percentile values for (a) ordinary kriging , (b) cokriging , (c) bedrock kriging , and (d) kriging on residuals from the optimal solution of Poisson's equation . The cross-section goes through the two southernmost boreholes indicated by dotted lines in Figure 9. The 50% percentile is identical to the expected mean value (Supplementary Appendix Equation A5).

Figure 10

East-West cross-sections of bedrock topography shown as percentile values for (a) ordinary kriging , (b) cokriging , (c) bedrock kriging , and (d) kriging on residuals from the optimal solution of Poisson's equation . The cross-section goes through the two southernmost boreholes indicated by dotted lines in Figure 9. The 50% percentile is identical to the expected mean value (Supplementary Appendix Equation A5).

Figure 10 shows clearly the differences in the percentile distribution around the expected median value. The , , and estimates had non-symmetric percentiles around the median values due to the non-Gaussian distribution of the sediment thickness. The percentiles, on the other hand, had symmetric distributions around the median value because these estimates were derived as residuals from the local trend. The cross-sections also illustrate clearly that the results had the smallest estimation variance.

Sedimentary wells and boreholes recorded after 2010 were used for independent cross-validation (Figures 2 and 4). The jackknife cross-validation results show, in general, that small sediment thicknesses were overestimated and large sediment thicknesses were underestimated (Figure 11). It should be noted that there are only minor differences between the parabolic surface obtained by Poisson's equation and the estimates in this case.

Figure 11

Jackknife cross-validation results of boreholes recorded after 2010 (cf. Figure 2). Median value estimates () for ordinary kriging (; blue and black dots), bedrock kriging (; red dots), cokriging (; magenta dots), and kriging on residuals from Poisson's equation, Equation 7 (; open red circles). Optimal solutions of Poisson's equation are indicated by black squares (Equations 5 and 6). indicates optimal solutions of Poisson's equation. All estimates were conditioned on observations recorded before 2010 with 10 boreholes. Dashed lines indicate percentiles : [0.05, 0.25, 0.75, 0.95]. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/nh.2021.102.

Figure 11

Jackknife cross-validation results of boreholes recorded after 2010 (cf. Figure 2). Median value estimates () for ordinary kriging (; blue and black dots), bedrock kriging (; red dots), cokriging (; magenta dots), and kriging on residuals from Poisson's equation, Equation 7 (; open red circles). Optimal solutions of Poisson's equation are indicated by black squares (Equations 5 and 6). indicates optimal solutions of Poisson's equation. All estimates were conditioned on observations recorded before 2010 with 10 boreholes. Dashed lines indicate percentiles : [0.05, 0.25, 0.75, 0.95]. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/nh.2021.102.

Cumulative histograms of probability scores (Equation 12) summarize the results of , , and (Figure 12). Perfect estimates yield = 0.5 with estimation variance equal to zero, which correspond to the vertical bold dashed line in Figure 12. Cross-validation results that overestimated the observations ( had , and vice versa). The cross-validation results by leaving one out gave from 20% to 30% of the observations for the , , and estimates. A smaller fraction of the observations had . The estimates had 86% of the estimates between (corresponding to the = 0.86 in Table 2). In jackknife cross-validation, this number was reduced to 80% for and about 25% for (i.e. 35% over-estimations and 40% under-estimations).

Figure 12

Cumulative probability score () histograms of cross-validated observations (Equation 12). (a) Cross-validation results by leaving one observation out where all observations within the Øvre Eiker study area were included (). (b) Jackknife cross-validation results, where estimates were conditioned on boreholes recorded before 2010 (), and cross-validation was done on boreholes recorded after 2010 (). The black solid curves indicate kriging on residuals from the optimal Poisson's equation () (Equations 5 and 6). Blue solid and dotted lines show results from ordinary kriging (). Red lines show bedrock kriging (), and magenta lines show cokriging () results. The superscripts indicate that calculations were done on lognormal (′) or normal score transformed data (*). Black dashed vertical lines indicate theoretical cross-validation results of observations corresponding to [0.25, 0.5, 0.75] probability scores with zero estimation variance. Results where < 0.5 correspond to overestimation of observations, while results where > 0.5 are underestimation of observations. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/nh.2021.102.

Figure 12

Cumulative probability score () histograms of cross-validated observations (Equation 12). (a) Cross-validation results by leaving one observation out where all observations within the Øvre Eiker study area were included (). (b) Jackknife cross-validation results, where estimates were conditioned on boreholes recorded before 2010 (), and cross-validation was done on boreholes recorded after 2010 (). The black solid curves indicate kriging on residuals from the optimal Poisson's equation () (Equations 5 and 6). Blue solid and dotted lines show results from ordinary kriging (). Red lines show bedrock kriging (), and magenta lines show cokriging () results. The superscripts indicate that calculations were done on lognormal (′) or normal score transformed data (*). Black dashed vertical lines indicate theoretical cross-validation results of observations corresponding to [0.25, 0.5, 0.75] probability scores with zero estimation variance. Results where < 0.5 correspond to overestimation of observations, while results where > 0.5 are underestimation of observations. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/nh.2021.102.

DISCUSSION

The cross-validation results demonstrate that the assumption of second-order stationarity of sediment thickness or bedrock topography is not valid. This conclusion was supported by the relatively low performance of the ordinary kriging cases (, , , and in Table 2). We suggested applying Poisson's equation to estimate the trend in sediment thickness. For this purpose, we used bedrock outcrops as boundary conditions. Ordinary kriging on the residuals between observations and the parabolic surface obtained from Poisson's equation (labelled Poisson kriging in Table 2), reduced the estimation uncertainty. All evaluation criteria applied in this study were consistent with this result. Hence, the method was superior to the other methods.

In this study, we explored the use of public data to estimate the bedrock topography as a continuous function in space. Previous research has underlined the importance of the bedrock topography in rainfall-runoff modelling (Genereux et al. 1993a, 1993b; Hinton et al. 1993; Mulholland 1993; Freer et al. 2002; Jencso et al. 2009; Graham & McDonnell 2010; Graham et al. 2010; Gabrielli et al. 2012).

For runoff response, the relation between sediment thickness and bedrock topography is most relevant in areas with sparse sediment cover and low bedrock permeability. Focusing on flow and threshold effects due to undulations in the bedrock topography may cause non-linear hysteresis effects in rainfall-runoff relations (Tromp-van Meerveld & McDonnell 2006a; Graham & McDonnell 2010; Graham et al. 2010; Gabrielli et al. 2012). To improve such non-linear effects in hydrology, estimates of sediment thickness and bedrock topography need to be included as auxiliary functions.

We modelled the estimation uncertainty of the bedrock topography (Equation 2) by using input data from the GRANADA database and digital information from the local study area (Øvre Eiker area, Figure 2). The purpose was to minimize the uncertainty by using DQMs, DEMs, and point observations of sediment thicknesses . The methods we applied relied on multi-Gaussian statistics, which is a simplification of the reality. We also assumed that the information in the DQM could be simplified to a spatial indicator function , where in areas where the bedrock was exposed (0 m), and in areas covered by deposits (0 m). Below, we discuss briefly some implications of the statistical and geological simplifications we made, but before doing so, it is necessary to comment on some general problems of using public data sources like the GRANADA database, for geostatistical modelling.

Clustering and bias of empirical data

Like most geo- and environmental databases, the GRANADA database contains clustered observations. This is due to the fact that spatial data is usually sampled for a specific purpose, which implies a large number of observations in a few limited areas, while most of the sampling domain has less frequent observations. In this case, boreholes and sedimentary wells were drilled in populated locations, which means that there is a high frequency of recordings in some places and less in other places. Uneven sampling known as clustering effects, ought to be controlled before statistical inference are undertaken. In this project, declustering was done by the moving grid method (Deutsch 1989; Deutsch & Journel 1998). Results from declustering of sediment thickness recorded in the GRANADA database show that small is more clustered than observations of large . Small sediment thicknesses should, therefore, receive less weights than large sediment thicknesses, and vice versa.

In addition to simple clustering, the GRANADA database has most likely another bias effect. The primary purpose for the GRANADA database was groundwater management, not sampling of sediment thicknesses. Because boreholes and sedimentary wells are usually drilled for economic reasons, locations with large sediment thickness are usually avoided due to the casing costs. Thus, boreholes may be abandoned for economic reasons if the sedimentary thickness is too large. The opposite is true for sedimentary wells. They are usually drilled for domestic water supply, which imply a desire for high water extraction capacity. The water flux into the well is governed by the hydraulic conductivity and the thickness of the water conductive layer. Preferred locations for such wells are, therefore, areas with large sediment thicknesses. Thus, there are different preferences with respect to sediment thickness D, for boreholes and sedimentary wells. The recordings of D are, therefore, prone to bias, with an over-representation of small D for boreholes and vice versa for sedimentary wells. These preferences have also impact on the location of boreholes and sedimentary wells with respect to horizontal distance to the nearest bedrock outcrop L, and the preferred Quaternary deposit. These effects are illustrated in the dataset from Øvre Eiker (Figure 2). Only one of the boreholes was drilled at a location where km, and 11 of the 25 boreholes had 100 m (Figure 4). Most of the sedimentary wells were also located quite close to the bedrock outcrop (100 m), but the sedimentary wells were located where D was expected to be large. All sedimentary wells within the study area were drilled in fluvial deposits. The nine clustered wells belonged to the same waterworks. The information attached to the well recordings indicate that these wells were located after a seismic survey, which means that a larger area was scanned and the most useful locations were selected. Since the water pumping capacity is proportional to the product of the permeability and the thickness of the sediment, the preferred location had large sediment thicknesses. Thus, if the relation is going to be utilized for modelling of the bedrock topography, data sources in addition to the GRANADA data should be considered. The GRANADA recordings are still useful for local conditioning, and because more and more boreholes and wells are registered the value of the database will increase.

Cross-validation results

Cross-validation of observations from the study area confirmed that kriging on the residuals from Poisson's equation gave the smallest mean absolute error, and best accuracy and precision. Estimation uncertainty was reduced by including locations of the bedrock outcrops, either as the boundary condition in Poisson's equation or as a secondary function in cokriging or bedrock kriging. The uncertainties were still significant as illustrated in Figures 8 and 11. This is not unexpected because the magnitude of the small-scale variance () was large compared with the total variance (). This fact explains the minor differences between a simple (regression) average and the kriging results. The large also explains the large spread in the percentile values at the points of observations. If the small-scale variance was zero (), all percentile values would have been equal to the observed value (Figure 10).

Bedrock kriging () and cokriging () rely on the relation between sediment thickness D and horizontal distance to the nearest outcrop L. The observations from the study area, however, show a large variance in the relation (Figure 4). This explains the large uncertainties in the and estimation results. These uncertainties were confirmed by the jackknife cross-validation exercise. It should also be noted that the relation in the national GRANADA database had a similar high variability. The percentile values for the national declustered relation (i.e. corrected for clustering effects) had the median value of = 4.7 [1.2, 18.0] m/km (0.25 and 0.75 percentile values in square brackets, Figure 13). In comparison, data from the study area Øvre Eiker had the median of 30 m/km, which corresponds to the 0.8 percentile value of the national dataset. This indicates that the local estimates should be utilized instead of the large-scale national percentiles. Histogram of from the national GRANADA dataset shows that for 50% of the boreholes with horizontal distance to bedrock outcrop 1 km had sediment thickness 4.7 m (Figure 13). Non-declustered median value of the same data material gave 6.7 m, which shows that there is a clustering of minor values in the GRANADA database. Declustered data shows that 25% of the wells with 1 km had sediment thickness 18.0 m and 10% had 60.8 m. Corresponding values for raw data (not declustered) was 20.2 and 54.7 m. For Øvre Eiker, the median value was 30.7 m/km. This shows that local conditioning should be used for modelling of bedrock topography if L is used as a secondary variable, which was the motivation for using cokriging and bedrock kriging in the present study. Even though the variance in the relation is high, the observations shown in Figure 4 indicate a physical relation: If L is large, there is no small D, and vice versa.

Figure 13

Histogram and cumulative probability density function (cdf) of for the national GRANADA database. The histogram (blue) and the cdf-obs (yellow) are derived from declustered observations. The brown curve (pdf-fit) and the purple curve (cdf-fit) were fitted to an optimal Gaussian pdf. The table was derived from inverse lognormal percentiles. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/nh.2021.102.

Figure 13

Histogram and cumulative probability density function (cdf) of for the national GRANADA database. The histogram (blue) and the cdf-obs (yellow) are derived from declustered observations. The brown curve (pdf-fit) and the purple curve (cdf-fit) were fitted to an optimal Gaussian pdf. The table was derived from inverse lognormal percentiles. Please refer to the online version of this paper to see this figure in colour: http://dx.doi.org/10.2166/nh.2021.102.

The cumulative histograms of (Equation 12) indicate that the observed sediment thickness was not estimated well by ordinary kriging , bedrock kriging , or cokriging (Figure 12). Poisson kriging , however, gave better results. The main reason is that Poisson's equation utilizes more of the available information than the other kriging methods. In and , information from a single point only was utilized for estimation (i.e. the distance from the estimation point to the nearest outcrop), while in , the total information of the outcrop geometry was included. This ability of Poisson's equation is due to the fact that the boundary condition is incorporated in the solution. Another reason for using Poisson's equation is that only one parameter – the load factor in Equation 4 is necessary to optimize to find a parabolic function, which can be adapted to the variable of interest. In this context, we regard the load factor as a stationary stochastic function , which is equivalent to the assumption of the ‘sameness’ referred to in the ANN application (Clarke et al. 2009; Mey et al. 2015). In this case, the variable of interest is the result of glacial and fluvial processes, which have been modelled by others as diffusive processes. Thus, the parabolic functions from Poisson's equation are intuitively convenient to apply.

A final warning might be added on the application of Poisson's equation: If is given as a deterministic function, and the point observations are used as boundary conditions in the same way as the outcrop locations, it is possible to fit a parabolic surface that goes perfectly through all the point observations. This is bad modelling practice because the uncertainty in is ignored. Thus, we reemphasize that instead of striving for a ‘perfect’ fit, it is better modelling practice to appreciate the uncertainty because, in most cases, the uncertainty is equally (or more) important than the best fit.

Public data are always affected by uncertainties. The areal information of DQMs is not always consistent with independent point information. At least in one case, the bedrock surface was registered as 7 m below the ground in the GRANADA database, but from the DQM, the borehole was located on exposed bedrock. Thus, geographical uncertainties should be included as part of the estimation procedure, and not hidden in the results.

Statistical simplifications

Modelling of the bedrock topography was first done by conditioning on normal score transformed data, and second, by using the lognormal transform. The normal score transform ranks the observations from minimum to maximum and finds the corresponding entries in a standard normal pdf (Deutsch & Journel 1998). If the empirical dataset does not cover extreme values, the normal score transform relies on extrapolations. Hence, the normal score transform is sensitive to the extreme values in the empirical data. The lognormal transform is by definition independent of the empirical data and should, therefore, yield more robust results with regard to extreme values. The estimation variance, however, usually becomes too large because the expected value is part of the variance (Supplementary Appendix Equations A13 and A14). This drawback explains why the lognormal modelling results had low precision. In the present case study, the modelling results based on the normal score transform did not differ very much from the lognormal transformed results. One reason is that even though the univariate pdf by definition is perfect Gaussian, the same is not automatically true for the bivariate pdf. This problem is visualized in Figure 14, where the bivariate normal score transforms of D and L were compared with a synthetic bivariate Gaussian variable with the covariance structure given in Table 1 for separation distance = 0 m. As can be seen, the Gaussian character is roughly preserved, but the bivariate pdf is noisy.

Figure 14

Synthetic multi-Gaussian probability density function to the left (a) normal score transformed sediment thickness D and the horizontal distance to the nearest bedrock outcrop L in the middle (b), and the difference between the synthetic pdf and the empirical to the right (c). The synthetic pdf was generated with zero mean and covariance function equal to the empirical data for separation distance (cf. Table 1).

Figure 14

Synthetic multi-Gaussian probability density function to the left (a) normal score transformed sediment thickness D and the horizontal distance to the nearest bedrock outcrop L in the middle (b), and the difference between the synthetic pdf and the empirical to the right (c). The synthetic pdf was generated with zero mean and covariance function equal to the empirical data for separation distance (cf. Table 1).

It should be underlined that the percentiles of sediment thickness and the bedrock topography yield a smoothed version of the reality (Figure 10). Thus, the percentiles do not give a realistic picture of the bedrock topography itself. The percentiles, however, is used to visualize the estimation uncertainty of the bedrock topography. If the physical structure of the bedrock topography or the sediment thickness plays an important role, then equally probable realizations should be generated. This subject was beyond the scope of the present paper, and it is, therefore, left for further studies.

Multi-Gaussian methods maximize entropy, which implies that the connectivity of extreme values is underestimated (Gómez-Hernández & Wen 1998). Subsurface canyons for example, which are extremes with respect to topography, will not be captured well by multi-Gaussian modelling. This is a drawback in areas with glacial erosion. Glacial flow may give rise to significant differences in topography over short distances perpendicular to the flow direction, and the longitudinal extent of the valleys may, therefore, be significant. In geological terms, it means that steep and narrow valleys cannot be simulated unconditionally by multi-Gaussian methods.

Geological simplifications

Geological maps hold vast amounts of information that might be utilized to estimate the bedrock topography. In this study, however, the only information we used was the indicator function , where in locations where the bedrock was exposed to the atmosphere, and else . Another relevant question is the relation between the geological setting (i.e. the sediment type) and the bedrock topography. It is common knowledge that the glacial history of the area explains the main character of the landscape. The weight of the glacier suppressed the continental crust, and the isostatic rebound after the melting explains why marine deposits cover the bedrocks in the lower part of the study area. Glacial transgression and retreat explain the deposits of the area, but surface mapping of Quaternary deposits does not usually include thickness of sediments. Some sediment types (in specific areas) usually have significant thickness (like marine sediments), while other kinds of deposits usually are patchy and sparse (like colluvial or weathered material). Thus, a relevant project for further research is to match the GRANADA recordings with DQMs to find statistical relations between sediment category and sediment thickness. This kind of information might be utilized for future (soft) conditioning of the bedrock topography. For sediment category and sediment thickness D within the study area, it can be noted that boreholes located in fluvial deposits have the largest sediment thickness D. If there exists a statistically significant relation between sediment category and D, then modelling of D might be improved by using sediment category as an auxiliary function.

CONCLUSIONS

With this paper, we wanted to draw attention toward modelling of sediment thickness and the associated bedrock topography (Equation 2). We explored the opportunities of combining information from the national well database (GRANADA), DQMs, and DEMs to reduce the estimation uncertainty.

Several MATLAB scripts were tailored to process the dataset for this study. This included declustering of the national data, normal score transformation, semivariogram calculations, and routines for inverse modelling of Poisson's equation. In this case, the small-scale variance (the nugget) was large compared with the total variance (the sill), which explains why the mean absolute error and the estimation variance were large. The large-scale national semivariogram model for the sediment thickness showed a correlation length of about 2 km, and the cross-semivariogram between the sediment thickness and the horizontal distance to the bedrock outcrop had a correlation length of more than 5.5 km (Table 1). These correlation lengths indicate that the estimation of the bedrock topography or the sediment thickness can capitalize on observations in the GRANADA database if the distances to the nearest GRANADA recording are less than 2–5.5 km depending on the estimation procedure.

Poisson's equation was used to estimate the local trend in the sediment thickness. This was done by inverse modelling of the load parameter in Poisson's equation (equivalent to the sink/source term in groundwater modelling). By this approach, we fitted a parabolic function to the bedrock topography, which minimized the differences between the point observations and the numerical solutions of Poisson's equation. At this stage, we solved Poisson's equation by using the Dirichlet boundary conditions and the sediment thickness as the primary variable. Identical results can be obtained by using terrain elevation as the Dirichlet boundary instead of the sediment thickness. Alternatively, the Neuman condition might be applied by using the gradient of the exposed bedrock topography, but this approach was not evaluated in the current study.

Local trend analysis improves the mean absolute error, but the independent (jackknife) cross-validation shows that the estimation uncertainty is large. Because the number of recorded boreholes, wells, and geotechnical test drillings increases every year, the conditional estimation uncertainties will be reduced. Point databases (as GRANADA and/or NADAG (NGU 2020d)) in combination with line and/or areal estimates (e.g. from geophysical surveys) will provide plentiful information in coming years. Thus, an important job for the future is to develop simple and robust estimation procedures that can benefit from the increasing amount of information continuously coming in.

AUTHOR CONTRIBUTIONS

The main part of the work was done by N.O.K. with significant contributions from É.L.

COMPETING INTERESTS

The authors declare no conflict of interests.

ACKNOWLEDGEMENTS

Thanks to the Geological Survey of Norway (NGU) for providing the dataset and to the Norwegian Institute of Bioeconomy Research (NIBIO) for supporting the data harvest. Thanks to anonymous and non-anonymous reviewers for critical and constructive comments, and finally big thanks to Camille Jouin who prepared the data.

DATA AVAILABILITY STATEMENT

All relevant data are available from the Geological Survey of Norway (geo.ngu.no/).

REFERENCES

REFERENCES
Andersen
T.
,
Trønnes
R. G.
,
Nilsen
O.
&
Larsen
A. O.
2008
Alkaline rocks of the Oslo Rift, SE Norway: a field trip with emphasis on felsic to intermediate intrusive rocks and their associated mineralizations. August 1–5, 2008. Eurogranites 2008/IGCP 510 field trip 2008. Publication of the Department of Geosciences, University of Oslo, 54 p. Available from: https://www.nhm.uio.no/english/about/organization/research-collections/people/rtronnes/4/oslorift-excursion-articl/eurogranites-2008-oslorift.pdf (accessed 13 March 2021)
.
Bargel
T. H.
1987
Hokksund 1714 I. Kvartærgeologisk kart M 1:50 000, Geological Survey of Norway. Available from: http://www.ngu.no/FileArchive/198/K17141.pdf (accessed 13 March 2021)
.
Beven
K.
2006
A manifesto for the equifinality thesis
.
Journal of Hydrology
320
.
https://doi.org/10.1016/j.jhydrol.2005.07.007
.
Bierkens
M.
,
van den Hurk
B.
2008
Feedback mechanisms: precipitation and soil moisture
. In:
Climate and the Hydrological Cycle
(
Bierkens
M. F. P.
,
Dolman
A. J.
&
Troch
P. A.
, eds).
IAHS Special Publication
,
Wallingford, UK
.
Bogen
J.
2017
Erosion rates and sediment yields of glaciers
.
Annals of Glaciology
22
,
48
52
.
doi:10.1017/s0260305500015202
.
Calakli
F.
&
Taubin
G.
2011
SSD: smooth signed distance surface reconstruction
.
Computer Graphics Forum
30
(
7
),
1993
2002
.
Clarke
G. K. C.
,
Berthier
E.
,
Schoof
C. G.
&
Jarosch
A. H.
2009
Neural networks applied to estimating subglacial topography and glacier volume
.
Journal of Climate
22
,
2146
2160
.
doi:10.1175/2008JCLI2572.1
.
Cook
S. J.
,
Swift
D
,
Kirkbride
A.
,
Knight
M. P.
&
Waller
P. G. R. I.
2020
The empirical basis for modelling glacial erosion rates
.
Nature Communications
11
,
759
.
doi:10.1038/s41467-020-14583-8
.
Culling
W. E. H.
1960
Analytical theory of erosion
.
Journal of Geology
68
,
336
344
.
Deutsch
C. V.
&
Journel
A. G.
1998
GSLIB. Geostatistical Software Library and User's Guide
, 2nd edn.
Oxford University Press
,
Oxford, New York
.
ISBN 0 19 510015 8
.
Dunne
T.
1978
Field studies of hillslope flow processes, Chapter 7
. In:
Hillslope Hydrology
(
Kirkby
M. J.
, ed.).
Wiley
,
London
, pp.
227
293
.
Dunne
T.
&
Black
R. D.
1970
An experimental investigation of runoff production in permeable soils
.
Water Resources Research
6
,
170
191
.
Freer
J.
,
McDonnell
J. J.
,
Beven
K. J.
,
Peters
N. E.
,
Burns
D. A.
,
Hooper
R. P.
,
Aulenbach
B.
&
Kendall
C.
2002
The role of bedrock topography on subsurface storm flow
.
Water Resources Research
38
(
12
).
doi:10.1029/2001WR000872
.
Gabrielli
C. P.
,
McDonnell
J. J.
&
Jarvis
W. T.
2012
The role of bedrock groundwater in rainfall-runoff response at hillslope and catchment scales
.
Journal of Hydrology
450–451
,
117
133
.
doi:10.1016/j.jhydrol.2012.05.023
.
Genereux
D. P.
,
Hemond
H. F.
&
Mulholland
P. J.
1993a
Spatial and temporal variability in streamflow generation on the west fork of Walker Branch Watershed
.
Journal of Hydrology
142
,
137
166
.
Giardino
J. R.
&
Houser
C.
2015
Introduction to the critical zone
. In:
Principles and Dynamics of the Critical Zone
(J. R. Giardino & C. Houser, eds). Vol.
19
.
Elsevier
,
Amsterdam
,
Netherlands
, pp.
1
14
,
doi:10.1016/B978-0-444-63369-9.00001-X
.
Gilbert
G. K.
1877
Geology of the Henry Mountains, USGS Unnumbered Series
.
Government Printing Office
,
Washington, DC
.
Available from: http://pubs.er.usgs.gov/publication/70038096 (accessed 13 March 2021)
.
Gneiting
T.
&
Raftery
A. E.
2004
Strictly Proper Scoring Rules, Prediction, and Estimation. Technical Report No. 463
.
Department of Statistics, University of Washington
,
Seattle, Washington
.
doi:10.1198/016214506000001437
.
Gómez-Hernández
J. J.
&
Wen
X.-H.
1998
To be or not to be multi-Gaussian? A reflection on stochastic hydrogeology
.
Advances in Water Resources
21
(
1
),
47
61
.
Goovaerts
P.
,
AvRuskin
G.
,
Meliker
J.
,
Slotnick
M.
,
Jacquez
G.
&
Nriagu
J.
2005
Geostatistical modeling of the spatial variability of arsenic in groundwater of southeast Michigan
.
Water Resources Research
41
,
W07013
.
doi:10.1029/2004WR003705
.
Gosses
M.
,
Nowak
W.
&
Wohling
T.
2018
Explicit treatment for Dirichlet, Neumann and Cauchy boundary conditions in POD-based reduction of groundwater models
.
Advances in Water Resources
115
,
160
171
.
doi.org/10.1016/j.advwatres.2018.03.011
.
Graham
C. B.
&
McDonnell
J. J.
2010
Hillslope threshold response to rainfall: (2) development and use of a macroscale model
.
Journal of Hydrology
393
,
77
93
.
doi:10.1016/j.jhydrol.2010.03.008
.
Graham
C. B.
,
Woods
R. A.
&
McDonnell
J. J.
2010
Hillslope threshold response to rainfall: (1) a field based forensic approach
.
Journal of Hydrology
393
,
65
76
.
doi:10.1016/j.jhydrol.2009.12.015
.
Gupta
H. V.
,
Clark
M. P.
,
Vrugt
J. A.
,
Abramowitz
G.
&
Ye
M.
2012
Towards a comprehensive assessment of model structural adequacy
.
Water Resources Research
48
,
W08301
.
doi:10.1029/2011WR011044
.
Harbor
J. M.
,
Hallet
B.
&
Raymond
C. F.
1988
A numerical model of landform development by glacial erosion
.
Nature
333
,
347
.
Herman
F.
,
Beyssac
O.
,
Brughelli
M.
,
Lane
S. N.
,
Leprince
S.
,
Adatte
T.
,
Lin
J. Y. Y.
,
Avouac
J.-P.
&
Cox
S. C.
2015
Erosion by an Alpine glacier
.
Science
350
(
6257
),
193
195
.
doi:10.1126/science.aab2386
.
Hinderer
M.
2001
Late Quaternary denudation of the Alps, valley and lake fillings and modern river loads
.
Geodinamica Acta
14
,
231
263
.
Hinton
M. J.
,
Schiff
S. L.
&
English
M. C.
1993
Physical-properties governing groundwater-flow in a glacial till catchment
.
Journal of Hydrology
142
,
229
249
.
Hopp
L.
&
McDonnell
J. J.
2009
Connectivity at the hillslope scale: identifying interactions between storm size, bedrock permeability, slope angle and soil depth
.
Journal of Hydrology
376
,
378
391
.
doi:10.1016/j.jhydrol.2009.07.047
.
Isaaks
H.
&
Srivastava
R. M.
1989
An Introduction to Applied Geostatistics. Oxford University Press, New York, pp. 561
.
Jencso
K. G.
,
McGlynn
B. L.
,
Gooseff
M. N.
,
Steven
M.
,
Wondzell
S. M.
,
Bencala
K. E.
&
Marshall
L. A.
2009
Hydrologic connectivity between landscapes and streams: transferring reach- and plot-scale understanding to the catchment scale
.
Water Resources Research
(
45
),
W04428
.
doi:10.1029/2008WR007225
.
Journel
A. G.
&
Huijbregts
C. J.
1989
Mining Geostatistics
, 4th edn.
Academic Press
,
London
.
Kazhdan
M.
,
Bolitho
M.
,
Hoppe
H.
2006
Poisson surface reconstruction
. In:
Eurographics Symposium on Geometry Processing
(
Polthier
K.
&
Sheffer
A.
, eds). The Eurographics Association.
Kitterød
I. H.
2011
Coupled Tectonic Surface Process Modelling of Large Scale Tectonic Geomorphology of Rifted Margins
.
MSc Thesis
,
Dep. Geosci. Earth Science, University of Bergen
,
Bergen, Norway
.
Kitterød
N.-O.
2017
Estimating unconsolidated sediment cover thickness by using the horizontal distance to a bedrock outcrop as secondary information
.
Hydrology and Earth System Sciences
21
,
4195
4211
.
https://doi.org/10.5194/hess-21-4195-2017
.
Kooi
H.
&
Beaumont
C.
1994
Escarpment evolution on high-elevation rifted margins: insights derived from a surface processes model that combines diffusion, advection, and reaction
.
Journal of Geophysical Research
99
(
12
),
191
209
.
Liu
C.-S.
2018
Fast solving the Cauchy problems of Poisson equation in an arbitrary three-dimensional domain
.
Computer Modeling in Engineering and Sciences
114
(
3
).
doi:10.3970/cmes.2018.114.351
.
MacGregor
K. R.
,
Anderson
R. S.
&
Waddington
E. D.
2009
Numerical modeling of glacial erosion and headwall processes in alpine valleys
.
Geomorphology
103
,
189
204
.
doi:10.1016/j.geomorph.2008.04.022
.
Mangerud
J.
,
Gyllencreutz
R.
,
Lohne
Ø.
&
Svendsen
J. I.
2011
Glacial history of Norway. In: Developments in Quaternary Sciences, Vol. 15, pp. 279–298. Quaternary Glaciations – Extent and Chronology: A Closer Look (J. Jürgen Ehlers, P. L. Gibbard, P. D. Hughes, eds). https://doi.org/10.1016/B978-0-444-53447-7.00022-2
.
MATLAB, Release
2020
The MathWorks Inc
.
Mey
J.
,
Scherler
D.
,
Zeilinger
G.
&
Strecker
M. R.
2015
Estimating the fill thickness and bedrock topography in intermontane valleys using artificial neural networks
.
Journal of Geophysical. Research: Earth Surface
120
,
1301
1320
.
https://doi.org/10.1002/2014JF003270
.
NGU
2020a
Geological Survey of Norway, GRANADA – National Groundwater Database
.
Available from: http://geo.ngu.no/kart/granada/ (accessed 13 March 2021)
.
NGU
2020b
Bedrock – National Bedrock Database
.
Available from: http://geo.ngu.no/kart/berggrunn (accessed 13 March 2021)
.
NGU
2020c
Superficial Deposits – National Database
.
Available from: http://geo.ngu.no/kart/losmasse/ (accessed 13 March 2021)
.
NGU
2020d
Geological Survey of Norway, NADAG – National Database for Geotechnical Surveys
.
Available from: http://geo.ngu.no/kart/nadag/ (accessed 07 April 2020)
.
NRC (National Research Council, USA)
2001
Basic research opportunities in earth science. In: The Critical Zone: Earth's Near-Surface Environment, pp. 35–45. The National Academies Press. Available from: https://www.nap.edu/read/9981/chapter/4 (accessed 13 March 2021)
.
Otto
J. C.
,
Schrott
L.
,
Jaboyedoff
M.
&
Dikau
R.
2009
Quantifying sediment storage in a high Alpine valley (Turtmanntal, Switzerland)
.
Earth Surface Processes and Landforms
34
,
1726
1742
.
doi:10.1002/esp.1856
.
Rempe
D. M.
&
Dietrich
W. E.
2014
A bottom-up control on fresh-bedrock topography under landscapes
.
PNAS
111
(
18
),
6576
6581
.
Richardson
J. B.
2017
Critical zone
. In:
Encyclopedia of Geochemistry. Encyclopedia of Earth Sciences Series
(
White
W.
, ed.).
Springer
,
Cham
.
doi:10.1007/978-3-319-39193-9_355-1
.
Shrestha
D. L.
2020
Continuous Rank Probability Score. MATLAB Central File Exchange. Available from: https://www.mathworks.com/matlabcentral/fileexchange/47807-continuous-rank-probability-score (retrieved 17 November 2020)
.
Straumann
R. K.
&
Korup
O.
2009
Quantifying postglacial sediment storage at the mountain-belt scale
.
Geology
37
,
1079
1082
.
doi:10.1130/G30113A.1
.
Sulebak
J. R.
&
Hjelle
Ø.
2003
Multiresolution spline models and their applications in geomorphology
. In:
Concepts and Modelling in Geomorphology: International Perspectives
(I. S. Evans, R. Dikau, E. Tokunaga, H. Ohmori & M. Hirano, eds).
TERRAPUB
,
Tokyo
, pp.
221
237
.