[go: up one dir, main page]
More Web Proxy on the site http://driver.im/
Next Article in Journal
Towards Better Visualisation of Alpine Quaternary Landform Features on High-Resolution Digital Elevation Models
Previous Article in Journal
BiFA-YOLO: A Novel YOLO-Based Method for Arbitrary-Oriented Ship Detection in High-Resolution SAR Images
You seem to have javascript disabled. Please note that many of the page functionalities won't work as expected without javascript enabled.
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Remote Sensing
Volume 13
Issue 21
10.3390/rs13214210
remotesensing-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A High-Resolution, Random Forest Approach to Mapping Depth-to-Bedrock across Shallow Overburden and Post-Glacial Terrain

by
Shane Furze
1,2,
Antóin M. O’Sullivan
1,2,*,
Serge Allard
3,
Toon Pronk
3 and
R. Allen Curry
1,2,4
1
Faculty of Forestry and Environmental Management, University of New Brunswick, Fredericton, NB E3B 5A3, Canada
2
Canadian Rivers Institute, University of New Brunswick, Fredericton, NB E3B 5A3, Canada
3
New Brunswick Department of Natural Resources and Energy Development, Fredericton, NB E3B 5H1, Canada
4
Department of Biology, Faculty of Science, University of New Brunswick, Fredericton, NB E3B 5A3, Canada
*
Author to whom correspondence should be addressed.
Remote Sens. 2021, 13(21), 4210; https://doi.org/10.3390/rs13214210
Submission received: 14 June 2021 / Revised: 20 September 2021 / Accepted: 21 September 2021 / Published: 20 October 2021
Browse Figures
Figure 1
<p>Visualization of the geographic location of New Brunswick, Canada (46.5653° N, 66.4619° W) outlining elevation change and physiographic regions (<b>A</b>) and an overview of bedrock lithology (<b>B</b>), overlain on World Ocean Basemap [<a href="#B30-remotesensing-13-04210" class="html-bibr">30</a>] in ArcGIS 10.3 software [<a href="#B31-remotesensing-13-04210" class="html-bibr">31</a>]. Bedrock geology was acquired from [<a href="#B28-remotesensing-13-04210" class="html-bibr">28</a>] (Scale 1:2,500,000).</p> "> Figure 2
<p>Overview of sample locations for each DTB source: boreholes (<b>A</b>); drillholes (<b>B</b>); pedons (<b>C</b>); site cards (<b>D</b>); well logs (<b>E</b>); and bedrock outcrops (<b>F</b>). Total sample size = 170,719 (see also <a href="#remotesensing-13-04210-t002" class="html-table">Table 2</a>).</p> "> Figure 3
<p>Histogram and continuous density distribution (red line) of bedrock depth for boreholes (<b>A</b>), drillholes (<b>B</b>), pedons (<b>C</b>), site cards (<b>D</b>), and well logs (<b>E</b>), and all data combined (<b>F</b>), including sample size, minimum, maximum, and mean DTB values for each source. Note: Depth intervals are (i) 0.25 m for pedons and site cards, and 5 m for the remaining sources, and (ii) the total sample size with minimum, maximum, and mean exclude the rock outcrop samples since all of these samples have a DTB of 0.</p> "> Figure 4
<p>The framework for developing the DTB model covariates. The multi-scale generation of covariates beginning with the 10 m DEM is shown as an example (blue frame) and the correlation analyses that followed to reduce the number of covariates (red frame).</p> "> Figure 5
<p>Scree plot used to select the number of PCs to capture a minimum of 90% variability (<span class="html-italic">n</span> = 11).</p> "> Figure 6
<p>The RF variable importance plot highlighting the explanatory strength of the 20 most significant covariates with importance measured as decreasing from 100% (most important) to 40%. Covariate abbreviations are explained in <a href="#remotesensing-13-04210-t003" class="html-table">Table 3</a>.</p> "> Figure 7
<p>Actual vs. fitted model results for the training and validation subsets ((<b>A</b>,<b>B</b>), respectively) including 95% confidence limits, and histograms displaying densities of residual errors for DTB model on training and validation subsets ((<b>C</b>,<b>D</b>), respectively).</p> "> Figure 8
<p>Predicted depth to bedrock (DTB) at 10 m<sup>2</sup> resolution for the Province of New Brunswick, Canada (<b>A</b>) with an example of the fine-scale resolution showing sediment accumulation across an upland–hillslope–valley bottom transition (<b>B</b>).</p> ">
Review Reports Versions Notes

Abstract

:
Regolith, or unconsolidated materials overlying bedrock, exists as an active zone for many geological, geomorphological, hydrological and ecological processes. This zone and its processes are foundational to wide-ranging human needs and activities such as water supply, mineral exploration, forest harvesting, agriculture, and engineered structures. Regolith thickness, or depth-to-bedrock (DTB), is typically unavailable or restricted to finer scale assessments because of the technical and cost limitations of traditional drilling, seismic, and ground-penetrating radar surveys. The objective of this study was to derive a high-resolution (10 m2) DTB model for the province of New Brunswick, Canada as a case study. This was accomplished by developing a DTB database from publicly available soil profiles, boreholes, drill holes, well logs, and outcrop transects (n = 203,238). A Random Forest model was produced by modeling the relationships between DTB measurements in the database to gridded datasets derived from both a LiDAR-derived digital elevation model and photo-interpreted surficial geology delineations. In developing the Random Forest model, DTB measurements were split 70:30 for model development and validation, respectively. The DTB model produced an R2 = 92.8%, MAE = 0.18 m, and RMSE = 0.61 m for the training, and an R2 = 80.3%, MAE = 0.18 m, and RMSE = 0.66 m for the validation data. This model provides an unprecedented resolution of DTB variance at a landscape scale. Additionally, the presented framework provides a fundamental understanding of regolith thickness across a post-glacial terrain, with potential application at the global scale.

1. Introduction

Bedrock is the consolidated rock that comprises the Earth’s surface. It is overlain by an amalgamation of unconsolidated materials, or regolith, including fractured bedrock, deposited (glacial, water, wind) and weathering minerals, and soil or organic residues at various stages of decomposition [1,2]. The composition of regolith can have varying heterogeneity depending on a combination of geology (mineralogy and age), topography, and past depositional processes [3,4]. The depth of the regolith, i.e., depth to bedrock (DTB) similarly varies significantly from fully exposed (outcrops) to depths exceeding 1 km beneath the surface in some locations [3,5]. Regolith composition and depth, in combination with bedrock fracturing and transmissivity, are fundamental regulators of surface (on-ground) and near-surface (below-ground) landscape forming processes [6,7]. Weathering combined with past restructuring in post-glacial landscapes and ongoing hydrological processes dictate the types of landforms including slopes and surface/sub-surface water pathways, and the types of soils [8,9]. These physical features in turn regulate the biota which occur within the regolith, on its surface, and in the surface waters that follow topography, and feedback within [10,11]. As a consequence, the regolith composition and depth are critical characteristics that influence many human needs and activities, including water supplies, agriculture, mineral extraction, forest harvesting, and most physical constructions on which humans depend, e.g., roads and buildings [12]. Despite its importance, DTB is often misunderstood and ignored because describing and modelling it at broad scales across the landscape is technologically challenging and costly [3,4,5].
Estimating the DTB shows promise with emerging geospatial technologies and the inclusion of various kinds of remote sensing information, such as digital elevation models (DEMs), associated topo-hydrological derivatives (e.g., slope, aspect, flow accumulation) and orthophotography. Table 1 outlines previous studies that predicted DTB for varying study extents, spatial resolutions, and statistical methods. Of the global DTB models, Shangguan et al. [5] explained 59% of DTB variability at 250 m2 using Random Forest (RF) and gradient boosting (GB), whereas Pelletier et al. [13] explained 34% on upland hillslopes at 1 km2 using algebraic equations applied to topographic derivatives based on conceptual models. Of the regional models, Yan et al. [2] followed a similar statistical approach to [5] and explained 57% of the DTB variability at 100 m2 resolution while [4] compared linear regression and ordinary kriging (OK) for predicting DTB and captured 68% of the variability at 20 m2. Metelka et al. [14] used artificial neural networks (ANN) capturing 94.2% of DTB variability at 30 m2. Additionally, Wilford et al. [15] and Wilford et al. [16] applied the cubist model to predict 38% and 64% of DTB variability at 90 m2 and 10 m2, respectively. Gomes et al. [3] achieved an RMSE of 1.4 m via Bayesian modeling for their predictions at 4 m2 whilst also utilizing linear regression, predicting DTB at 2 m2 with an average RMSE of 2 m. Finally, Shafique et al. [17] applied stepwise regression to account for 91.7% of the variability at 30 m2, while Kuriakose [18] tested regression kriging (RK) for their DTB model at 20 m2, explaining 59% of the variability.
Varying interpolation techniques have also been applied for DTB predictions. Karlsson et al. [19] utilized Inverse Distance Weighting (IDW), while Kuriakose [18] applied RK for predicting DTB. Alternatively, other studies have developed models for predicting soil depth, but these typically do not specify whether it is soil depth to parent material or bedrock [20]. Typically, bedrock has been represented in the form of broad-scale, tessellated vector datasets outlining the bedrock formation, age, permeability, etc., but lacks information pertaining to DTB [21,22]. Alternatively, DTB has been measured at specific locations via soil profiles (capturing shallow DTB), drill holes and boreholes for mineral exploration, well logs for water-well development, and transects for locating areas of rock outcrops, with these data represented in point vector format. These data types are commonly used in spatial modeling for predicting DTB across the landscape [2,19]. In some locations, spatial predictions of soil depth can be, and have been, used as surrogates for DTB [4,16].
During the last glacial maximum (18,000 years (18 kyr)) ~30% of the Earths’ land surface was covered in ice [23]. As the climate warmed, this extent receded to ~10% of the Earths’ land surface, limiting ice to high latitudes and altitudes [24,25,26]. Post-glacial landscapes comprise a cacophony of features such as eskers, kames, kettles, moraines, outwash plains, terraces, deltas, etc., and the depth-to-bedrock in post-glacial landscapes can vary from >100 m, e.g., Northern Alberta, Canada, to <5 m, e.g., New Brunswick, Canada. Understanding the DTB at a landscape scale could significantly advance geological and hydrological sciences and applications in natural resource management and development, e.g., water supplies, mineral exploration, road building, and forest management. The objective of this study was to develop a high-resolution spatial map of DTB in a post-glacial landscape characterized by shallow overburden, covering 72,000 km2 region of New Brunswick, Eastern Canada as a case study. This was attempted by building an RF model utilizing field-acquired DTB measurements, high-resolution—LiDAR-derived—topographic and hydrographic derivatives, and geomorphological parameters. New Brunswick was selected for this analysis due to the variability in landforms, elevation, and geological features.

2. Methods

2.1. Study Region

The province of New Brunswick (NB), situated in the Atlantic zone of Eastern Canada, spans an area of ~72,000 km2. NB houses an array of diverse geomorphologies resulting from a total glacial coverage during the Pleistocene 40,000 years before the present (ybp), as well as ongoing erosional and depositional processes, resulting in numerous, continuously varying physiographic regions, topographies, landforms, and landscapes which are faced with ongoing erosional and depositional processes (Figure 1A) [12]. Deposits across NB are represented by a combination of sedimentary, igneous, and metamorphic lithological bedrock types ranging from Early Neoprotozoic (850 Ma—1 Ga) to Early Cretaceous (99–145 Ma) with an average DTB of 1.3 m (range of 0–184 m) [27] (Figure 1B). Physiographic regions within the province range from highlands in the southeast (Caledonia Highlands) and northwest (Chaleur Uplands, Miramichi Highlands, and Edmundston Highlands) to lowlands in the east (New Brunswick lowlands), and lowlands/foothills in the southwest (Saint John River Valley/Highland Foothills) [28]. The variability in bedrock type and geomorphological features across the landscape has resulted in a heterogeneous DTB with exposed bedrock occurring within the highlands, uplands, and foothills, to deep post-glacial deposits and glacial landforms [29].

2.2. Regolith and DTB Data Sources

Spatially explicit DTB measurements were acquired and amalgamated from publicly available datasets retrieved from governmental sources as outlined in Table 2 and Figure 2 and Figure 3, including:
  • Boreholes: A borehole database was accessed from NB Department of Energy and Resource Development (NBERD) via a national topographic system (NTS) index [32]. These data represent the results of petroleum, potash, and coal production/exploration (n = 332).
  • Drillholes: This database is an amalgamation of drill holes collected as part of diamond drilling mineral exploration within the province. DTB was recorded in exploration reports. These data were accessed from Service New Brunswick’s GeoNB data catalogue [21] (n = 9,828).
  • Pedons: Soil samples collected as part of multiple research projects, including the development of county-based soil surveys by the Canadian Soil Information Service [33] and a permanent sample plot database from NBERD [34]. From the databases, samples were queried and selected for those with DTB recorded (n = 199).
  • Site Cards: A set of till geochemistry site cards exist that represent till geochemistry surveys conducted throughout NB. Site cards were retrieved from the geoscience publication search query [35] and limiting search results to open file reports. These data were amalgamated and only samples with DTB measurements were selected for analyses (n = 324).
  • Well Logs: This dataset represents the recorded data for all new and deepened drinking water wells in NB since 1994, amalgamated by the New Brunswick Department of Environment and Local Government (NBELG). The depth at which bedrock was contacted was also recorded in the dataset. This information was retrieved from GeoNB data catalogue [21] (n = 33,187).
  • Rock Outcrops: the rock outcrops represent the results of on-ground transects conducted throughout the province to locate areas in which bedrock was exposed at the surface (NBERD) (n = 126,849).

2.3. Data Quality Assessment

Prior to statistical analyses, a procedure was developed to locate potential outliers within each DTB dataset. First, each dataset was assessed to determine if any point locations had multiple DTB values (two measurements for the same geographic location). This was completed by searching the database for rows that shared identical latitude and longitude. If identical, both samples were removed from the dataset, resulting in the removal of 2 samples from rock outcrops, 12 samples from site cards, 4 samples from boreholes, 1099 samples from drill holes, and 5346 samples from well logs. Second, points that occurred within open quarries and mines were excluded from the analysis. These features were located using the GeoNB non-forest inventory dataset (SLU = Mi and SLU = QU) [36]. This resulted in the removal of 877 drill holes, 4 boreholes, 347 rock outcrops, 1 site card, and 5568 well logs. Next, the aggregate vector data set, acquired from NBERD granular aggregate database [37] was used to identify DTB samples occurring on deep alluvial and glaciofluvial deposits. From this, bedrock outcrop points occurring within the delineated granular aggregates were removed due to the contradictory nature of each source (n = 347). At this time, the well logs were removed from the spatial database due to inconsistencies arising throughout the samples. These data were collected as a by-product, with the focus being the depth-to-ground water for each well. As a result, what was believed to be the surface of bedrock may have been a clay belt, hardpan, large coarse fragment deposit, or fractured bedrock. Thus, the actual DTB may not be adequately captured within this dataset. In addition, some well log recordings had multiple DTB values for any single location, raising caution when incorporating this dataset.

2.4. Modelling Framework

The digital elevation model (DEM) utilized in this study is a seamless dataset at 10 m resolution derived from the amalgamation of light detection and ranging (LiDAR) acquisitions at 1 m resolution. The LiDAR digital elevation models (DEM) were combined, ensuring correct projections, datums, and spatial extents, and re-sampled to 10 m using the mean statistic [31]. The original airborne LiDAR acquisitions for NB were collected between 2015 to 2018 at a point density of 6 points∙m−2. The individual LiDAR datasets are provided in .laz format with a tile size of 1 km2. This open-sourced data was retrieved from the GeoNB LiDAR Download web app (https://www.arcgis.com/apps/webappviewer/index.html?id=a6d24f0cc5f149dfb8b569c4ab8aa9f4 (accessed on 4 June 2017)). This form of elevation data collection results in DEMs at higher quality and resolution than counterparts derived from other forms of remote sensing (i.e., satellites, photogrammetry, or contours).
A framework was developed using several DEMs in combinations of increasing cell size (resolution) and localized averaging of neighboring cells (focal smoothing) of the original LiDAR-derived 10 m DEM to produce myriad topographic datasets via algorithms within ArcGIS and SAGA GIS software (Figure 4, Appendix A). In addition, each DEM was hydro-conditioned to fill depressions and ensure correct hydrological flow patterns and used to develop the hydrological datasets [38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73]. All gridded datasets (topo-hydrological and geological), hereafter referred to as the ‘covariates’ (n = 398), were resampled to 10 m resolution so that they shared: (i) spatial projection (NAD 1983 CSRS New Brunswick Stereographic); (ii) resolution (10 m); and (iii) spatial extent (number of columns and rows). Similarly, the point vector dataset (field-based DTB measurements) shared the same spatial projection as covariates. Within the modeling framework, the covariates created were annotated with their DEM resolution, e.g., 10 m resolution DEM abbreviated as ‘dem10′, whereas the 10 m DEM averaged using a 5-cell search window was abbreviated as ‘dem105′. The result was eight variant models for each of the DEM scales, e.g., 10 m (‘dem10′), 10 m smoothed (‘dem105′), 30 m (‘dem30′), 30 m smoothed (‘dem305′), 60 m (‘dem60′), 60 m smoothed (‘dem605′), 100 m (‘dem100′), and 100 m smoothed (‘dem1005′).

2.5. Model Parameters

2.5.1. Geological Data

The available tessellated vector dataset representing bedrock geology [28] at a scale of 1:500,000 was converted to a gridded format to represent bedrock age at 10 m resolution. In addition, available tessellated surficial geology delineations [29], also at a 1:500,000 scale, were utilized to develop gridded datasets (at 10 m resolution), representing the mode of deposition and primary lithology, and from these, mineral hardness (derived from Moh’s hardness scale), dominant rock type, and dominant grain size of the surficial geology. Additionally, the vector dataset representing aggregate deposits, at a 1:50,000 scale, was converted to a 10 m gridded dataset in binary form, where 1 was assigned if within an aggregate and 0 assigned elsewhere. In addition, aggregate locations were used as reference for calculating Euclidean distances from the aggregate deposits as an attempt to highlight possible un-delineated aggregate deposits. This technique was also applied to the rock outcrop locations, delineated waterbodies, and wetlands [21].
The surficial geology datasets were included since the study extent was heavily influenced by past glaciation. The surficial aggregate data sets (“agg_dist” and “agg_10”) were included because they represent the deep deposits with highest DTBs. The lithological data sets (“geo_class”, “grain_class”, “lith_class”, and “hardness”) were selected because they represent the parent materials’ resistance to weathering. Finally, the landform datasets (“land_class” and “wetland”) were included since they represent the different parent material modes of deposition, with these also influencing DTB. Definitions for abbreviations are provide in the supplementary data (see Appendix A).

2.5.2. Topo-Hydrological Data

Bedrock plays a large role in influencing current landform patterns and complex changes in topography. As such, topographic datasets are typically used in predicting DTB [3]. These, in combination with delineations of surficial geology, provide the necessary information required to make predictions on DTB across a landscape. With that said, topographic change is not adequately captured via a single spatial resolution, e.g., 10 m DEM. At this resolution, there are trade-offs between improved vertical accuracy of elevation and inability to visualize large-scale changes in topography. At 10 m resolution, it is difficult to interpret the entire hillslope from ridge to valley; instead, small-scale variability in topography is captured with such high-resolution DEMs. Thus, a multi-scale approach was used on a LiDAR-derived 1 m DEM resampled to 10 m in which the same topographic derivatives are produced at multiple DEM resolutions and scales (Figure 4). Recreating the same topographic derivatives at larger resolutions, while applying focal smoothing, better captures the change in topography from ridge to valley. This study created the topo-hydrological derivatives at 10, 30, 60, and 100 m2 resolution with each set smoothed using a 5-cell focal search window. The topo-hydrological datasets utilized in this study were chosen due to their ability to represent landscape evolution.

2.5.3. Covariate Reduction

The output covariates from each DEM-scaled model and all geological covariates were matched to the field-measured DTB observations. Null values, locations where the corresponding covariate value for any observation = −9999, were excluded from the analyses. The result was 398 covariates for 176,658 DTB observations. Principal Components Analysis (PCA) was conducted on the covariates using the ‘prcomp’ function in the ‘stats’ package in R statistical software [74]. Predictor variables were standardized (value/1 standard deviation) and centered (value—average). A Scree plot was then used to determine the suitable number of principal components (PCs) which capture 90% of the variability (n = 11) (Figure 5). Variance of 90% was selected for two reasons: (i) it captures adequate variability while keeping the number of PCs low and (ii) adding additional PCs has very few improvements in variability captured. Adversely, adding additional PCs slows down computational efficiency.
From here, the loading values for each covariate were retained for the first 11 PCs, while the others were excluded. Pearson correlation analysis was conducted on the remaining loading values to assess the degree of multicollinearity (absolute correlation values ≥ 70% were considered significant). Counts were conducted to determine how many covariates met this condition. The covariate with the highest count was retained while the correlated covariates were removed. This counting and removal procedure was implemented iteratively until no correlated covariates remained. This covariate reduction technique was chosen for two reasons: (i) the inherent multicollinearity between topo-hydrological datasets produced from the same DEM (although derived at different resolutions and scales), (ii) potential model overfitting, and (iii) the computational limitations arising from the management of 398 covariates, each at 10 m resolution and each covering ~72,000 km2.
Of the topographic datasets, slope and curvature datasets help represent zones of accumulation and loss with steep convex slopes representing zones of loss, while relatively flat concave slopes represent zones of accumulation. In addition, longer and flatter slopes will have fewer erosional forces over time in comparison to shorter, steeper slopes. Datasets representing surface exposure to the atmosphere (how exposed a position is) (i.e., ‘diu’ [45], ‘aspect’ [38,39], ‘sky’, ‘trn’, ‘view’, ‘vis’ [58,59,60], ‘tpo’, and ‘tno’ [67,68]) represent areas which may be exposed to prolonged periods of weathering or increased biological activity, which will increase DTB. For example, south-facing slopes receive more direct sunlight and warmer temperatures then north-facing slopes, while areas of exposed ridges may receive more direct temperature fluctuations and precipitation increasing the rate of rock and mineral weathering. Datasets representing overall change in topography (‘mrrtf’, ‘mrvbf’ [57], ‘tpi’ [39,66], ‘rsa’ [57], ‘tri’ [64], ‘tst’ [65], and ‘vrm’ [70]) represent the change in slope, curvature, and landforms to broader extents, which helps quantify zones of additions and losses. The hydrological datasets are important for understanding the power of water in carving landscapes. For example, ‘ca’ [40,41], ‘cs’ [42,43], ‘sca’ [40], and ‘scs’ [42,43] represent how much water is flowing into any given location and how fast that water is flowing through the landscape. These, in combination with ‘st_hdist’, ‘st_vdist’, ‘swi’ [71], ‘val’ [69], and ‘mrvbf’, show where steeply incised stream channels have developed, as well as where the eroded sediments have likely deposited, hence representing shallow and thicker DTB, respectively.

2.6. Statistical Modeling

Random Forest (RF) is a machine learning algorithm with which multiple decision trees are produced with each decision tree developed using a random subset (bootstrap sampling) of the original dataset [75,76]. Prior to model development, the dataset was split so that 30% of the data was set aside for validation. Of this 30%, an equal number of samples from each source were selected. The remaining 70% of the dataset was used for model development. The RF model was developed and applied utilizing the ‘caret’ and ‘raster’ packages in R statistical software [75,76,77]. This model was set up using a tuning grid, k-fold cross validation, and controlling the number of trees to be developed. The tuning grid tested a matrix of covariates at each node within each regression tree, ranging from 3 to 15 at iterations of three. The maximum of 15 variables was selected to ensure that it remained under the ‘p/3′ rule [75]. K-fold cross validation was set at 10 folds with 20 repeats [78]. The 10 folds were selected due to the large population size of the training subset (n = 94,385) and as an attempt to limit bias during model development. In addition, 50 trees were selected for RF model development to avoid computational delays. The resulting RF model was applied to the spatial covariates using the ‘write raster’ function, within the ‘raster’ package in R statistical software [76], yielding a spatial prediction for DTB at 10 m resolution for NB. Model performance was assessed using a coefficient of determination (R2), mean absolute error (MAE), and root mean squared error (RMSE).

3. Results

3.1. Covariate Reduction

The PCA-PCoA covariate reduction technique was applied to the database resulting in the removal of 354 predictor covariates, from 398 to 44 covariates (see Appendix A). Table 3 highlights the remaining covariates used for model development and the number of correlated covariates of each.

3.2. Model Results

Variable importance for the 20 most significant covariates in the RF model is highlighted in Figure 6, with the full list provided in Appendix A. This shows the explanatory strength of the covariates which explain the highest amount of variability, with scores standardized to equal 100. The DEM at 10 m resolution was the strongest predictor, followed by valley depth (10 m), terrain surface texture (‘tst’, 100 m smoothed with 5 m search window), flow line curvature (‘flw’) at 100 m smoothed and at 100 m, then multi-resolution index of ridgetop flatness (‘mrrtf’) at 100 m smoothed. The RF for the training data produced an R2 = 92. 8%, MAE = 0.18 m, and RMSE = 0.61 m, and an R2 = 80.3%, MAE = 0.18 m, and RMSE = 0.66 m for the validation data. The RF model results are visualized in Figure 7 with the resulting spatial DTB predictions for NB, at 10 m2 resolution, displayed in Figure 8. As illustrated in Figure 7, model performance decreases with increasing depth, with all samples having a mean DTB of 1.3 m (Figure 3), and 98.9% of samples (both training and validation) with an absolute error of ±2 m (Table 4).

4. Discussion

4.1. Depth-to-Bedrock Model

Post-glacial landscapes account for ~20% of the Earth’s surface landscapes. These landscapes, which occur mostly at high latitudes, are also the most at risk from a changing climate [23,24]. This presents an increasing necessity to develop an understanding of the regional-scale mechanisms driving abiotic and biotic processes across, and in, these landscapes. A foundational control on these processes is DTB [79,80]. This study demonstrates that DTB in a post-glacial landscape with shallow overburden can be accurately mapped using a combination of field measurements, multi-scale, high resolution topo-hydrological derivatives, geologic parameters, and a Random Forest model (R2~0.93). The data underlying the DTB produced here are typically available from regional and local resources managers, and LiDAR data are becoming more ubiquitous. In contrast to most literature on predicting DTB, this study introduces the concept of incorporating delineations of surficial geology into the modelling procedure due to the extensive impact past glacial events have had on shaping the contemporary bedrock composition and geomorphology. The framework presented here provides a pathway for future studies attempting to better understand how DTB varies across post-glacial and other landscapes with shallow DTB.
Geomorphic features that characterize post-glacial landscapes can dilute the relationship between DTB and surface topographic derivatives [81,82]. Terminal moraines, for example, are formed by glaciers, pushing eroded material downslope and depositing upon glacier retreat [78,83]. Sub-glacial tunnel melt waters that deposit glaciofluvial material become eskers as glaciers retreat [9,84], while a low gradient valley bottom set in a paleo-valley may comprise of deposits >40 m [8,85]. Such geomorphic features are driven by erosion and deposition processes, and undoubtedly confound the relationship between topographic controls on DTB [86]. We sought to capture some of this variance by utilizing multi-resolution topographic derivatives, geological parameters such as bedrock, surficial geology, and aggregate deposition in valleys by applying a non-linear machine learning model (RF). Unsurprisingly, the model’s capacity to capture DTB diminished with increasing depth (Figure 7); however, the model illustrated some ability to delineate zones of deposition, and their deeper DTB, in valley bottoms (Figure 8). The results outlined in this study are similar to those of [14] and [17] in terms of model performance, while outperforming the others. Of the models which utilized RF, the model produced in this study outperformed those ([5] and [2]) by 21 and 23%, respectively. The performance of this model is likely due to the high resolution of the input predictor data sets as well as the distribution of the DTB, with most measurements being <2 m.
Some insights into controls on DTB in this post-glacial landscape are also evident. In NB, DTB is shallowest (0–1 m) in the Highlands and deepens moving eastward towards the NB Lowlands (100+ m) (Figure 1 and Figure 8). The variable importance results from the RF model show that elevation was the strongest covariate for predicting DTB. Areas of higher elevation in NB represent dense igneous bedrock within the upland and highland physiographic regions (Figure 1). Additionally, highland areas in NB also have tors, which have been attributed to the presence of cold-based glaciers [87]. Cold-based glaciers are characterized by basal temperatures below 0 °C, and can be fused to the bedrock, limiting deposition of eroded materials relative to warm-based glaciers, and potentially leading to shallow DTB in these areas [88,89]. Contrastingly, areas with prolonged, subtle changes in elevation (such as the New Brunswick Lowlands—Figure 1) may represent regions which have been more eroded and/or more influenced by glacial activity. These areas also represent locations with sedimentary bedrock, which is more susceptible to weathering [90,91]. Valley depth was the second most important predictor, followed by the real surface area (‘rsa’) and terrain surface texture (‘tst’). Areas with steeply incised stream channels tend to have shallower DTB (Figure 8). The geomechanical characteristics of the glacial deposits in NB can be broadly classified as having low angles of repose, excluding basal tills and clays [28,92]. Weathering, and fluxes in effective stresses underpinned by pore water pressure in this northern temperate region likely lead to recurrent slope failures; thus, shallow DTB across steep slopes in NB.
Flowline curvature (‘flw’, at three resolutions and scales) and topographic negative openness (‘tno’) represented four of the strongest covariates. Flowline curvature is a function of elevation, slope, and rate of change of the slope. We suggest that these areas result in shallow DTB due the eroding capacity of water. Conversely, areas defined as ‘open’ have lower slopes, are mostly delimited by wetlands and glaciofluvial outwash areas, and result in deep DTB [28,93]. This is clearly illustrated in the depositional zone in Figure 8, where a large open valley, comprising outwash, is demarcated by deep DTB. This runs congruent to the open, upland areas which are also characterized by deeper DTB (see Figure 1). The final hydrological covariate identified as a strong predictor is modified catchment area (‘mca’), representing the upriver catchment contributing area. The largest upriver catchment-contributing areas occur at the lowest elevations. A decrease in potential energy in these lowland areas leads to deposition in these river reaches [93]. The New Brunswick Lowlands are populated with paleo-terraces that defined the hydrography of the ancient Northumbrian Valleys (see [94]). Hence, deeper DTB is to be expected in this physiographic setting (Figure 1 and Figure 8).
The exposure of ridge tops and upper portions of hills is represented by a multi-resolution index of valley bottom flatness (‘mrrtf’), sky view factor (‘sky’, at three resolutions and scales), real surface area (‘rsa’), and diurnal anistropic heating (‘diu’). Exposed ridge tops and upper hillslopes tend to have shallower DTB than lower lying areas (e.g., [95,96]). This is likely caused by glaciers removing material from exposed areas and compacting the remaining material as basal tills [97,98]. The curvature covariates: tangential curvature (‘tan’), longitudinal curvature (‘lng’), local downslope curvature (‘locdn’), maximum curvature (‘max’), and upslope curvature (‘up’), all represent the rate of change in elevation. Concave landscapes accumulate material, whereas their upslope, convex counterparts tend to lose material as weathering, gravity, glaciers, and runoff erode materials [99,100]. Finally, the grain size (‘grain_size’) of the lithology of the parent material was a dominant predictor. While this may be pointing towards weathering rates associated with differing lithological units, it is possible that this is also a function of geographic position. Most of the sedimentary parent material occurs in the NB Lowlands (Figure 1). This area is characterized by low relief and wetlands [101,102]. Low potential energy in this environment facilitates deposition, hence the identification of deeper DTB along with grain size (Figure 1 and Figure 8).

4.2. Limitations & Uncertainties

Using multi-sourced data always requires an emphasis on the quality assessment of each dataset, i.e., (i) the source of each dataset, (ii) methods of measurement, (iii) units of measurement, and (iv) assessing the spatial location of point observations. This study has highlighted techniques to develop a homogenous database from different sources of data, but inconsistencies still occur which cannot be addressed. For example, drill hole information does not provide the angle in which the drill pierced the regolith to reach the bedrock. It is assumed that this occurred at 90° through the regolith, and if not, then the DTB was corrected for the drilling angle. Additionally, when drilling, it is difficult to distinguish the contact zone between consolidated bedrock and overlying regolith. The regolith can be composed of fractured bedrock, dense basal till, or even physically or chemically induced compacted or cemented mineral soil horizons, such as hardpans and orstein horizons, respectively. Additionally, contradictions arose when comparing delineations of surficial geology to the spatial locations for field-based measurements, e.g., rock outcrops occurring within or adjacent to aggregate delineations. These points were excluded, but emphasize the critical importance of quality the assessment for each data set.
Although high-resolution datasets pertaining to topo-hydrology from a LiDAR-derived DEM were utilized in this analysis, it is difficult to fully capture the DTB variability across the landscape. This is due to the influence of past glacial events and the ongoing erosional and depositional processes, which have been acting on the soil and parent material for the past centuries since the Pleistocene glaciation. Varying surficial geological deposits have resulted in DTB values ranging from 0 m (exposed bedrock) on ridges and steep banks, to 80+ m in sediment-filled valleys. As such, existing delineations of surficial geology are unable to capture the full variability in changing surficial deposits and more information is required in the form of specific glacial deposit types (e.g., drumlins, kames, and eskers) and lithology of different deposits (the type of rock and grain size), since this will influence how easily weathered the material will be. These are reflected in the model results in Figure 8 with the maximum DTB prediction error reaching 30 m. Although lacking information, existing surficial geology datasets have importance in the DTB predictions since water-deposited materials (i.e., alluvial floodplains and glaciofluvial outwash plains) have more accumulated sediments on top of bedrock than some glacial and residual deposits (e.g., basal till and residual material). In any case, the ability to predict DTB at increasing depths will become more difficult since there are no single attributes explaining these depths. For example, the aggregate deposits do not specify the depth of the aggregates; thus, the bedrock could be 10 m from the surface or 80+ m. This is reflected in Figure 6 in which model performance decreases with DTB exceeding 5 to 10 m.
Applying such a model to other geographic locations poses additional difficulties. This model was developed for a post-glacial landscape in which topography is largely a replicate of shallow DTB. The transferability of this model is likely limited to regions that have witnessed similar landscape and climatic evolution, e.g., the eastern United States and Scotland. Inversely, this modelling exercise would most likely perform poorly in areas with deep overburden, where topography is not a replicate of topography, e.g., north-central United States and Alberta, Canada [5]. Additionally, the database used for model development is biased to shallow DTB measurements (Figure 3), with 94.6% of the samples recording a DTB of ≤1 m. Hence, model results and spatial prediction are biased toward shallow DTB. A more homogenous sample distribution would likely lead to poorer performance results due to the inherent difficulty in predicting deeper DTB.

5. Conclusions

This study aimed to model DTB at 10 m2 resolution for the Province of New Brunswick, Canada, by way of RF modeling on open-sourced data. New Brunswick was selected due to its variability in both bedrock and geomorphological features, as well as the availability of open-sourced data. Using this information, the model captured 80.3% of the variability on a randomly selected validation subset of the DTB measurements with an MAE and RMSE of 0.18 and 0.66 m, respectively. At 10 m resolution, the DTB model derived in this study can directly aid in spatial analyses which attempt to quantify natural processes across the landscape. These predictions were derived by comparing LiDAR-derived topo-hydrological, and geology-derived gridded datasets (from tessellated vector delineations) to an amalgamated, spatially explicit dataset of DTB measurements. The final model was then applied to the covariates, yielding a spatial prediction of DTB for NB at 10 m resolution (Figure 8).
The DTB model developed in this study provides two main benefits to further spatial analyses: (i) a 10m DTB spatial dataset and (ii) a framework for similar analyses to be conducted in other geographic settings, or with updated data and statistical techniques as these come to light. This analysis highlights the applicability and limitations of utilizing multiple open-sourced datasets in spatial statistical modeling studies. We continue to work on improving datasets, defining covariates, and testing alternative statistical approaches, and we look forward to testing the framework on other glaciated and non-glaciated landscapes. The framework we present for building landscape scale maps of regolith depth is a promising step to advance science and resource management in post-glacial terrains with shallow overburden.
Understanding DTB at a landscape scale is certain to advance earth science studies, such as those investigating regolith structure, soils, hydrogeology, and hydrology. Mapping DTB and defining its characteristics at landscape scales will also advance applications in mineral exploration, e.g., aggregates, hydrogeology of water supplies, and engineering human structures. These DTB maps are also consequential to ecological studies, for instance, examining the spatial variability of habitats and biota, e.g., forest type and productivity; thus, they are also useful for regional natural resource management and conservation which rely heavily on landscape scale data. Unfortunately, it is often too difficult and costly to obtain detailed depictions of DTB. This study addressed this limitation by developing an RF model to derive a high-resolution depiction of DTB utilizing publicly available data.

Author Contributions

S.F.: conceptualized, collated data, built models, wrote the first draft; A.M.O.: conceptualized, wrote, reviewed and edited; S.A.: provided data, reviewed and edited; T.P.: provided data, reviewed and edited; R.A.C.: wrote; reviewed, edited, and attained funding. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by NSERC (RGPIN-2018-06015—RAC); Atlantic Salmon Conservation Foundation (128044); New Brunswick Innovation Foundation (128133).

Institutional Review Board Statement

Not Applicable.

Informed Consent Statement

Not Applicable.

Data Availability Statement

The data used for this study can be downloaded from the following sources: LiDAR: https://www.arcgis.com/apps/webappviewer/idex.html?id=a6d24f0cc5f149dfb8b569c4ab8aa9f4 (accessed on 4 June 2017); Drillholes: http://www.snb.ca/geonb1/e/DC/Drill.asp (accessed on 21 May 2019); Boreholes: https://www1.gnb.ca/0078/GeoscienceDatabase/Borehole/Search.asp?_ga=2.157373114.1653925072.1605800013-910544185.1605800013 (accessed on 21 May 2019); Pedons: extracted from individual soil survey reports retrieved from: http://sis.agr.gc.ca/cansis/publications/surveys/nb/index.html (accessed on 21 May 2019); Site Cards: http://dnr-mrn.gnb.ca/ParisWeb/PublicationSearch.aspx (accessed on 21 May 2019); Well Logs: https://www.elgegl.gnb.ca/0375-0001/pointRadiusSearch.aspx?queryType=3&userType=1&provinceWide=1 (accessed on 21 May 2019).

Acknowledgments

The authors would like to acknowledge the dedication and hard work of the government agencies who have spent countless hours collecting, recording, and editing the data used in this study. A special thank you to the N.B. Department of Natural Resources and Energy Development, Service New Brunswick, N.B. Department of Environment and Local Government, and the Canadian Soil Information service for making DTB measurements publicly available. We specifically thank W. Gilmore of the N.B. Department of Energy and Mines. Finally, the authors would like to thank P. Arp, UNB for providing the computation requirements necessary for conducting this study. AMO’S is funded by grants from the Atlantic Salmon Conservation Foundation, New Brunswick Innovation Foundation, and with RAC, NSERC RGPIN-2018-06015. AMO’S would also like to thank Lord Pisces.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

Table
Table A1. Covariates assessed for use in bedrock modeling with description, reference, and software utilized for development.
Table A1. Covariates assessed for use in bedrock modeling with description, reference, and software utilized for development.
CovariateAbbreviationDefinitionReferenceSoftware Used
AggregateAgg1_0Binary representing if a cell falls within a delineated aggregate deposit or not.This studyArcGIS
Distance to AggregateAgg_DistEuclidean (horizontal) distance to aggregate depositThis StudyArcGIS
AspectAspCardinal direction of slope (measured in degrees)[38,39]SAGA GIS
Bedrock AgeBed_ageAverage age of bedrock geology types across New BrunswickThis StudyArcGIS
Catchment AreaCaUpslope area contributing discharge into any given cell, calculated from flow accumulation[40,41]ArcGIS
Catchment SlopeCsSlope of catchment, derived as intermediate with SAGA Wetness Index[42,43]SAGA GIS
Convergence Index (gradient)CvgrDetermines the convergence or divergence based on neighboring cells (calculated using gradient)[44]SAGA GIS
Convergence Index (aspect)CvasDetermines the convergence or divergence based on neighboring cells (calculated using aspect)[44]SAGA GIS
Cross-sectional CurvatureCrsLike planar curvature[38]SAGA GIS
Distance to WaterbodyWb_DistEuclidean (horizontal) distance to delineated waterbodyThis studyArcGIS
Distance to WetlandWetlandEuclidean (horizontal) distance to delineated wetlandThis study
Diurnal Anistropic HeatingDiuInfluence of topography on diurnal heat balance[45]SAGA GIS
Curvature ClassificationCurRepresentation of surface curvature/9 geometric forms of hillslopes[46,47]SAGA GIS
Downslope CurvatureDwnLike profile curvature[40]SAGA GIS
Downward Distance GradientDwdMeasures impact of local slope on hydraulic gradient[49]SAGA GIS
Digital Elevation ModelDemRepresentation of elevation change across landscape, height above sea level[50,51]ArcGIS
Filled DEMFiDepression-filled version of the DEM[52]SAGA GIS
Flowline CurvatureFlwLike profile curvature[40]SAGA GIS
General CurvatureGenCombination of planar and profile curvature[39]SAGA GIS
Dominant Geology TypeGeo_classGridded representation of changing geological types (sedimentary, igneous, metamorphic) of parent materialThis studyArcGIS
Dominant Grain SizeGrain_classGridded representation of changing mineral sizes from very fine to coarse based on lithology of parent materialThis studyArcGIS
Horizontal Distance to StreamSt_hdistEuclidean distance to nearest streamThis studyArcGIS
LandformLand_ClassGridded representation of parent material (surficial geology) mode of depositionThis studyArcGIS
LatitudeLatGridded representation of changing latitude from north to south (decimal degree format)This studyArcGIS
Local CurvatureLocCalculates the total gradient to neighboring cells[40]SAGA GIS
Local Downslope CurvatureLocdnSame as local curvature but only looks at downslope cells[40]SAGA GIS
Local Upslope CurvatureLocupSame as local curvature but only looks at upslope contributing cells[40]SAGA GIS
LongitudeLongGridded representation of changing longitude from west to east (decimal degree format)This studyArcGIS
Longitudinal CurvatureLngLike profile curvature[38]SAGA GIS
LithologyLith_classGridded representation of parent material (surficial geology lithology)This studyArcGIS
Maximum CurvatureMaxCalculates the maximum slope of the slope on a defined search radius (secondary derivative)[38,53,54,55]SAGA GIS
Mineral HardnessHardnessGridded representation of the hardness of lithologic types based on Moh’s hardness scaleThis studyArcGIS
Minimum CurvatureMinCalculates the minimum slope of the slope on a defined search radius (secondary derivative)[38,53,54,55]SAGA GIS
Modified Catchment AreaMcaAdjustment to catchment area calculation to correct for flow in low-lying flat areas.[42,43]SAGA GIS
Modified Specific Catchment AreaMscaModified version of Specific Catchment Area[40]SAGA GIS
Multi-resolution Index of Valley Bottom FlatnessMrvbfDelineation of valley bottom flatness from surrounding hillslopes[56]SAGA GIS
Multi-resolution Index of Ridge Top FlatnessMrrtfDelineation of ridgetops flatness from surrounding hillslopes, intermediate of MRVBF[56]SAGA GIS
Outcrop DistanceOut_DistEuclidean (horizontal) distance to rock outcropsThis studyArcGIS
Planar CurvaturePlnCurvature perpendicular to slope[38]SAGA GIS
Profile CurvaturePrfCurvature parallel to slope[38]SAGA GIS
Real Surface AreaRsaCalculation of ‘real’ cell area[57]SAGA GIS
Dominant Rock TypeRock_classGridded representation of changing rock types based on lithology of parent materialThis studyArcGIS
SAGA Wetness IndexSwiComputes a modified topographic wetness index[42,43]SAGA GIS
Sky View FactorSkyAmount of sky hemisphere visible from the ground[58,59,60]SAGA GIS
Slope (%)SlpRate of elevation change between adjacent cells (as %)[38,39]SAGA GIS
Slope LengthSlDetermine the length of slope[61]SAGA GIS
Slope Length and SteepnessLsCalculates slope length factor, typically used in Revised Universal Soil Loss Equation (RUSLE)[43,62,63]SAGA GIS
Specific Catchment AreaScaCatchment area divided by cell width[40]SAGA GIS
Specific Catchment SlopeScsCatchment area divided by cell width[42,43]SAGA GIS
Tangential CurvatureTanDetermines areas of convex and concave flows, calculated from planar and profile curvature[39]SAGA GIS
Terrain Ruggedness IndexTriTotal change in elevation within a defined radius of any given cell[64]SAGA GIS
Terrain Surface ConvexityTscvDetermines percentage of upward and convex cells within a defined radius[65]SAGA GIS
Terrain Surface ConcavityTsccDetermines percentage of upward and concave cells within a defined radius, same algorithm as terrain surface convexity[65]SAGA GIS
Terrain Surface TextureTstDetermines the variability in frequency and intensity of pits and peaks within a defined radius[65]SAGA GIS
Terrain View FactorTrnOutput covariate from Sky View Factor[58,59,60]SAGA GIS
Multi-scale Topographic Position IndexTpiPosition on hillslope in relation to adjacent cells[39,66]SAGA GIS
27Topographic Positive OpennessTpoRepresents landscape exposure to atmosphere (positive)[67]SAGA GIS
Topographic Negative OpennessTnoRepresents landscape exposure to atmosphere (negative)[67,68]SAGA GIS
Total CurvatureTotCurvature of the surface[39]SAGA GIS
Upslope CurvatureUpAverage local curvature of upslope contributing area for any given cell[40]SAGA GIS
Valley DepthValLike inverse of vertical distance to channel network[69]SAGA GIS
Vector Ruggedness MeasureVrmQuantifies landscape ruggedness via slope and aspect[70]SAGA GIS
Vertical Distance to StreamSt_vdistVertical distance from nearest stream based on increasing elevation[71,72] ArcGIS
View DistanceViewAverage distance to horizon, output from Sky View Factor[58,59,60]SAGA GIS
Visible SkyVisThe percentage of unobstructed hemisphere, output from Sky View factor[58,59,60]SAGA GIS

References

  1. Weil, R.R.; Brady, N.C. The Nature and Properties of Soils, 15th ed.; Fox, D., Gilfillan, A., Dimmick, L., Eds.; Pearson Education: Upper Saddle River, NJ, USA, 2017; ISBN 9780133254488. [Google Scholar]
  2. Yan, F.; Shangguan, W.; Zhang, J.; Hu, B. Depth-to-bedrock map of China at a spatial resolution of 100 meters. Science 2020, 7, 1–13. [Google Scholar]
  3. Gomes, G.J.C.; Vrugt, J.A.; Vargas, E.A., Jr. Toward improved prediction of the bedrock depth underneath hillslopes: Bayesian inference of the bottom-up control hypothesis using high-resolution topographic data. Water Resour. Res. 2016, 52, 3085–3112. [Google Scholar] [CrossRef] [Green Version]
  4. Devkota, S.; Shakya, N.M.; Sudmeier, K.; McAdoo, B.G.; Jaboyedoff, M. Predicting soil depth to bedrock in an anthropogenic landscape: A case study of Phewa Watershed in Panchase region of Central-Western Hills, Nepal. J. Nepal Geol. Soc. 2018, 55, 173–182. [Google Scholar] [CrossRef]
  5. Shangguan, W.; Hengl, T.; de Jesus, J.M.; Yuan, H.; Dai, Y. Mapping the global depth to bedrock for land surface modeling. J. Adv. Model. Earth Syst. 2017, 9, 65–88. [Google Scholar] [CrossRef]
  6. Perron, J.T.; Dietrich, W.E.; Kirchner, J.W. Controls on the spacing of first-order valleys. J. Geophys. Res. Earth Surf. 2008, 113, 1–21. [Google Scholar] [CrossRef] [Green Version]
  7. Chen, A.; Darbon, J.Ô.; Morel, J.M. Landscape evolution models: A review of their fundamental equations. Geomorphology 2014, 219, 68–86. [Google Scholar] [CrossRef]
  8. Rampton, V.N.; Gauthier, R.C.; Thibault, J.; Seaman, A.A. Quaternary Geology of New Brunswick; Geological Survey of Canada: Ottawa, ON, Canada, 1984; ISBN 0660115921.
  9. Burke, M.J.; Brennand, T.A.; Sjogren, D.B. The role of sediment supply in esker formation and ice tunnel evolution. Quat. Sci. Rev. 2015, 115, 50–77. [Google Scholar] [CrossRef]
  10. Corenbilt, D.; Steiger, J. Vegetation as a major conductor of geomorphic changes on the Earth surface: Toward evolutionary geomorphology Vegetation. Earth Surf. Process. Landf. 2009, 34, 891–896. [Google Scholar] [CrossRef]
  11. Phillips, J.D. Biological energy in landscape evolution. Am. J. Sci. 2009, 309, 271–289. [Google Scholar] [CrossRef]
  12. Pronk, A.G.; Ruitenberg, A.A. Applications of surficial mapping to forest management in New Brunswick. Atl. Geol. 1991, 27, 209–212. [Google Scholar] [CrossRef] [Green Version]
  13. Pelletier, J.D.; Broxton, P.D.; Hazenberg, P.; Zeng, X.; Troch, P.A.; Niu, G.-Y.; Williams, Z.; Brunke, M.A.; Gochis, D. A gridded global data set of soil, intact regolith, and sedimentary deposit thicknesses for regional and global land surface modeling. J. Adv. Model. Earth Syst. 2016, 8, 41–65. [Google Scholar] [CrossRef]
  14. Metelka, V.; Baratoux, L.; Jessell, M.W.; Barth, A.; Jezek, J. Automated regolith landform mapping using airborne geophysics and remote sensing data, Burkina Faso, West Africa. Remote Sens. Environ. 2018, 204, 964–978. [Google Scholar] [CrossRef]
  15. Wilford, J.R.; Searle, R.; Thomas, M.; Pagendam, D.; Grundy, M.J. A regolith depth map of the Australian continent. Geoderma 2016, 266, 1–13. [Google Scholar] [CrossRef]
  16. Wilford, J.; Thomas, M. Predicting regolith thickness in the complex weathering setting of the central Mt Lofty Ranges, South Australia. Geoderma 2013, 206, 1–13. [Google Scholar] [CrossRef]
  17. Shafique, M.; Van Der Meijde, M.; Rossiter, D.G. Geophysical and remote sensing-based approach to model regolith thickness in a data-sparse environment. Catena 2011, 87, 11–19. [Google Scholar] [CrossRef]
  18. Kuriakose, S.L.; Devkota, S.; Rossiter, D.G.; Jetten, V.G. Prediction of soil depth using environmental variables in an anthropogenic landscape, a case study in the Western Ghats of Kerala, India. Catena 2009, 79, 27–38. [Google Scholar] [CrossRef]
  19. Karlsson, C.S.J.; Jamali, I.A.; Earon, R.; Olofsson, B.; Mörtberg, U. Comparison of methods for predicting regolith thickness in previously glaciated terrain, Stockholm, Sweden. Geoderma 2014, 226–227, 116–129. [Google Scholar] [CrossRef]
  20. Hengl, T.; de Jesus, J.M.; MacMillan, R.A.; Batjes, N.H.; Heuvelink, G.B.M.; Ribeiro, E.; Samuel-Rosa, A.; Kempen, B.; Leenaars, J.G.B.; Walsh, M.G.; et al. SoilGrids1km—Global soil information based on automated mapping. PLoS ONE 2014, 9, e105992. [Google Scholar] [CrossRef] [Green Version]
  21. Service New Brunswick GeoNB Data Catalogue. Available online: http://www.snb.ca/geonb1/e/DC/catalogue-E.asp (accessed on 21 May 2019).
  22. Gleeson, T.; Smith, L.; Moosdorf, N.; Hartmann, J.; Dürr, H.H.; Manning, A.H.; Van Beek, L.P.H.; Jellinek, A.M. Mapping permeability over the surface of the Earth. Geophys. Res. Lett. 2011, 38, 1–6. [Google Scholar] [CrossRef] [Green Version]
  23. Hughes, P.D.; Gibbard, P.L.; Ehlers, J. Timing of glaciation during the last glacial cycle: Evaluating the concept of a global “Last Glacial Maximum” (LGM). Earth-Sci. Rev. 2013, 125, 171–198. [Google Scholar] [CrossRef]
  24. Haeberli, W.; Hoelzlem, M.; Paul, F.; Zemp, M. Integrated monitoring of mountain glacier as key indicators of global climate change: The European Alp’s by Haeberli and others. Ann. Glaciol. 2007, 46, 150–160. [Google Scholar] [CrossRef] [Green Version]
  25. Comiso, J.C.; Vaughan, D.G.; Allison, I.; Carrasco, J.; Kaser, G.; Kwok, R.; Mote, P.; Murray, T.; Paul, F.; Ren, J.; et al. Observations: Cryosphere. In Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change; Stocker, T.F., Qin, D., Plattner, G.K., Tignor, M., Allen, S.K., Boschung, J., Nauels, A., Xia, Y., Bex, V., Midgley, P.M., Eds.; Cambridge University Press: Cambridge, UK; New York, NY, USA, 2013; p. 84. [Google Scholar]
  26. Dalton, A.S.; Margold, M.; Stokes, C.R.; Tarasov, L.; Dyke, A.S.; Adams, R.S.; Allard, S.; Arends, H.E.; Atkinson, N.; Attig, J.W.; et al. An updated radiocarbon-based ice margin chronology for the last deglaciation of the North American ice sheet complex. Quat. Sci. Rev. 2020, 234, 106223. [Google Scholar] [CrossRef]
  27. Pronk, A.G.; Allard, S. Landscape Map of New Brunswick; Map NR–9. Scale 1:770,000; New Brunswick Department of Natural Resources and Energy, Minerals, Policy and Planning Division: Fredericton, NB, Canada, 2003.
  28. Bedrock Geology of New Brunswick; Map NR–1, (2008 Edition) Scale 1:500,000; New Brunswick Department of Natural Resources Policy and Planning Division: Fredericton, NB, Canada, 2008.
  29. Rampton, V.N. Generalized Surficial Geology Map of New Brunswick; Map NR–8. Scale 1:500,000; New Brunswick Department of Natural Resources and Energy, Minerals, Policy and Planning Division: Fredericton, NB, Canada, 1984.
  30. Environmental Systems Research Institute World Ocean Basemap. 2018. Available online: https://www.arcgis.com/home/item.html?id=6348e67824504fc9a62976434bf0d8d5 (accessed on 7 July 2018).
  31. Environmental Systems Research Institute ArcGIS Deskop 2018. Release: 10.3; Environmnetal Systems Research Institute: Redlands, CA, USA, 2018.
  32. New Brunswick Department of Energy and Mines New Brunswick Borehole Database. Available online: https://www1.gnb.ca/0078/GeoscienceDatabase/Borehole/Search.asp?_ga=2.246412036.1432594313.1632484669-910544185.1605800013 (accessed on 21 May 2019).
  33. Government of Canada; Agriculture and Agri-Food Canada. Canadian Soil Information Service. Available online: https://sis.agr.gc.ca/cansis/ (accessed on 21 May 2019).
  34. Porter, K.B.; Maclean, D.A.; Beaton, K.P.; Upshall, J. New Brunswick Permanent Sample Plot Database (PSPDB v1.0): User’s Guide and Analysis; Natural Resources Canada, Canadian Forest Service, Atlantic Forestry Centre: Fredericton, NB, Canada, 2001; ISBN 0662299884.
  35. New Brunswick Department of Energy and Mines Geoscience Publication Search Query. Available online: http://dnr-mrn.gnb.ca/ParisWeb/PublicationSearch.aspx (accessed on 21 May 2019).
  36. New Brunswick Department of Energy and Resource Development New Brunswick Non-Forest Inventory 2019. Available online: http://www.snb.ca/geonb1/e/dc/non-forest.asp (accessed on 21 May 2019).
  37. New Brunswick Department of Energy and Resource Development New Brunswick Granular Aggregate Database. Available online: http://www1.gnb.ca/0078/GeoscienceDatabase/GranularAgg/GranAgg-e.asp (accessed on 21 May 2019).
  38. Zevenbergen, L.W.; Thorne, C.R. Quantitative Analysis of Land Surface Topography. Earth Surf. Process. Landf. 1987, 12, 47–56. [Google Scholar] [CrossRef]
  39. Wilson,, J.P.; Gallant,, J.C. (Eds.) Terrain Analysis—Principles and Applications; John Wiley & Sons Inc.: Hoboken, NJ, USA, 2000; ISBN 978-0-471-32188-0. [Google Scholar]
  40. Freeman, G.T. Calculating Catchment Area with Divergent Flow Based on a Regular Grid. Comput. Geosci. 1991, 17, 413–422. [Google Scholar] [CrossRef]
  41. Tarboton, D.G. A New Method for the Determination of Flow Directions and Upslope Areas in Grid Digital Elevation Models. Water Resour. Res. 1997, 33, 309–319. [Google Scholar] [CrossRef] [Green Version]
  42. Boehner, J.; Koethe, R.; Conrad, O.; Gross, J.; Ringeler, A.; Selige, T. Soil Regionalisation by Means of Terrain Analysis and Process Parameterisation; Soil Classification, 2001, Micheli, E., Nachtergaele, F., Montanarella, L., Eds.; European Soil Bureau: Luxembourg, 2002. [Google Scholar]
  43. Boehner, J.; Selige, T. Spatial prediction of soil attributes using terrain analysis and climate regionalisation. In SAGA—Analyses and Modelling Applications; Boehner, J., McCloy, K.R., Strobl, J., Eds.; Goettinger Geographische Abhandlungen: Goettingen, Germany, 2006. [Google Scholar]
  44. Koethe, R.; Lehmeier, F. SARA—System zur Automatischen Relief-Analyse, 2nd ed.; User Manual; Department of Geography, University of Goettigen: Goettingen, Germany, 1996. [Google Scholar]
  45. Hengl, T.; Hannes, I.R. (Eds.) Geomorphometry: Concepts, Software, Applications; Elsevier: Amsterdam, The Netherlands, 2009. [Google Scholar]
  46. Dikau, R. Entwurf Einer Geomorphographisch-Analytischen Systematik von Reliefeinheiten, 5th ed.; Heidelberger Geographische Bausteine: Heidelberg, Germany, 1998. [Google Scholar]
  47. Birkeland, P. Soils and Geomorphology, 3rd ed.; Oxford University Press: Oxford, UK, 1999. [Google Scholar]
  48. Hugget, R.J. Fundamentals of Geomorphology, 2nd ed.; Routledge Fundamentals of Physical Geography Series; Routledge: New York, NY, USA, 2007. [Google Scholar]
  49. Hjerdt, K.N.; McDonnel, J.J.; Seibert, J.; Rodhe, A. A New Topographic Index to Quantify Downslope Controls on Local Drainage. Water Resour. Res. 2004, 40, 1–6. [Google Scholar] [CrossRef] [Green Version]
  50. Moore, I.D.; Grayson, R.B.; Ladson, A.R. Digital Terrain Modelling: A Review of Hydrological, Geomorphological, and Biological Applications. Hydrol. Process. 1991, 5, 3–30. [Google Scholar] [CrossRef]
  51. Yue, T.-X.; Du, Z.-P.; Song, D.-J.; Gong, Y. A New Method of Surface Modeling and its Application to DEM Construction. Geomorphology 2007, 91, 161–172. [Google Scholar] [CrossRef]
  52. Wang, L.; Liu, H. An efficient method for identifying and filling surface depressions in digital elevation models for hydrologic analysis and modeling. Int. J. Geogr. Inf. Sci. 2007, 20, 199–213. [Google Scholar] [CrossRef]
  53. Irvin, B.J.; Ventura, S.J.; Slater, B.S. Fuzzy and Isodata Classification of Landform Elements from Digital Terrain Data in Pleasant Valley, Wisconsin. Geoderma 1997, 77, 137–154. [Google Scholar] [CrossRef]
  54. Petry, F.E.; Robinson, V.B.; Cobb, M.A. (Eds.) Fuzzy Modeling with Spatial Information for Geographic Problems; Springer: Berlin/Heidelberg, Germany, 2005. [Google Scholar]
  55. Olaya, V. Basic land-surface parameters. In Geomorphometry: Concepts, Software, Applications; Developments in Soil, Science; Hengl, T., Reuter, H.I., Eds.; Elsevier: Amsterdam, The Netherlands, 2006. [Google Scholar]
  56. Gallant, J.C.; Dowling, T.I. A Multiresolution Index of Valley Bottom Flatness for Mapping Depositional Areas. Water Resour. Res. 2003, 39, 1–6. [Google Scholar] [CrossRef]
  57. Grohmann, C.H.; Smith, M.J.; Riccomini, C. Multiscale Analysis of Topographic Surface Roughness in the Midland Valley, Scotland. IEEE Trans. Geosci. Remote Sens. 2011, 49, 1200–1213. [Google Scholar] [CrossRef] [Green Version]
  58. Oke, T.R. Boundary Layer Climates; Taylor and Francis: New York, NY, USA, 2000. [Google Scholar]
  59. Häntzschel, J.; Goldberg, V.; Bernhofer, C. GIS-based regionalisation of radiation, temperature and coupling measures in complex terrain for low mountain ranges. Meteorol. Appl. 2005, 12, 33–42. [Google Scholar] [CrossRef]
  60. Bohner, J.; Antonie, O. Land-Surface Parameters Specific to Topo-Climatology. In Geomorphometry: Concepts, Software, Applications; Developments in Soil Science; Hengl, T., Reuter, H.I., Eds.; Elsevier: Amsterdam, The Netherlands, 2009. [Google Scholar]
  61. Olaya, V. Slope Length. SAGA GIS Slope Length Module in Terrain Analysis—Hydrology Toolset. 2004. Available online: http://www.saga-gis.org/saga_tool_doc/2.2.6/ta_hydrology_7.html (accessed on 8 March 2018).
  62. Desmet, P.J.J.; Govers, G. A GIS Procedure for Automatically Calculating the USLE LS Factor on Topographically Complex Landscape Units. J. Soil Water Conserv. 1996, 51, 427–433. [Google Scholar]
  63. Kinnell, P.I. Alternative Approaches for Determining the USLE-M Slope Length Factor for Grid Cells. Soil Sci. Soc. Am. J. 2005, 69, 674–680. [Google Scholar] [CrossRef]
  64. Riley, S.J.; DeGloria, S.D.; Elliot, R. A Terrain Ruggedness Index that Qauntifies Topographic Heterogeneity. Int. J. Sci. 1999, 5, 23–27. [Google Scholar]
  65. Iwahashi, J.; Pike, R.J. Automated Classifications of Topography from DEMs by an Unsupervised Nested-means Algorithm and a Three-part Geometric Signature. Geomorphology 2007, 86, 409–440. [Google Scholar] [CrossRef]
  66. Guisan, A.; Weiss, S.B.; Weiss, A.D. GLM versus CCA Spatial Modeling of Plant Species Distribution. Plant Ecol. 1999, 143, 107–122. [Google Scholar] [CrossRef]
  67. Yokoyama, R.; Shirasawa, M. Visualizing Topography by Openness: A New Application of Image Processing to Digital Elevation Models. Photogramm. Eng. Remote Sens. 2002, 68, 251–266. [Google Scholar]
  68. Prima, O.D.A.; Echigo, A.; Yokoyama, R.; Yoshida, T. Supervised Landform Classification of Northeast Honshu from DEM-derived Thematic Maps. Geomorphology 2006, 78, 373–386. [Google Scholar] [CrossRef]
  69. Conrad, O. Valley Depth. SAGA GIS Valley Depth Module in Terrain Analysis—Channels Toolset. 2012. Available online: http://www.saga-gis.org/saga_tool_doc/2.1.3/ta_channels_7.html (accessed on 8 March 2018).
  70. Sappington, J.M.; Longshore, K.M.; Thompson, D.B. Quantifying Landscape Ruggedness for Animal Habitat Analysis: A Case Study Using Bighorn Sheep in the Mojave Desert. J. Wildl. Manag. 2007, 71, 1419–1426. [Google Scholar] [CrossRef]
  71. O’Callaghan, J.F.; Mark, D.M. The Extraction of Drainage Networks from Digital Elevation Data. Comput. Vision Graph. Image Process. 1984, 28, 323–344. [Google Scholar] [CrossRef]
  72. Nobre, A.D.; Cuartas, L.A.; Hodnett, M.; Renno, C.D.; Rodrigues, G.; Silveira, A.; Waterloo, M.; Saleska, S. Height Above the Nearest Drainage—A Hydrologically Relevant New Terrain Model. J. Hydrol. 2011, 404, 13–29. [Google Scholar] [CrossRef] [Green Version]
  73. Saunders, W.K.; Maidment, D.R. A GIS Assessment of Nonpoint Source Pollution in the San Antonio-Nueces Coastal Basin; Center for Research in Water Resources, University of Texas at Austin: Austin, TX, USA, 1996. [Google Scholar]
  74. Florinsky, I. Digital Terrain Analysis in Soil Science and Geology, 1st ed.; Elsevier: Amsterdam, The Netherlands, 2012; ISBN 9780874216561. [Google Scholar]
  75. Kuhn, M. Building Predictive Models in R using the Caret Package. Caret: Classification and Regression Training 2017. J. Stat. Softw. 2008, 28, 1–26. [Google Scholar] [CrossRef] [Green Version]
  76. Hijmans, R.J.; van Etten, J. Raster: Geographic Data Analysis and Modeling 2019. R package version 2.0-12. Available online: https://cran.r-project.org/web/packages/raster/index.html (accessed on 16 April 2018).
  77. Kuhn, M.; Johnson, K. Applied Predictive Modeling, 1st ed.; Springer Science + Business Media: New York, NY, USA, 2013; ISBN 978-1-4614-6848-6. [Google Scholar]
  78. Benediktsson, Í.Ö.; Schomacker, A.; Johnson, M.D.; Geiger, A.J.; Ingõlfsson, Ó.; Guomundsdõttir, E.R. Architecture and structural evolution of an early Little Ice Age terminal moraine at the surge-type glacier Múlajökull, Iceland. J. Geophys. Res. Earth Surf. 2015, 120, 1895–1910. [Google Scholar] [CrossRef]
  79. Phillips, J.D.; Marion, D.A. Biomechanical effects of trees on soil and regolith: Beyond treethrow. Ann. Assoc. Am. Geogr. 2006, 96, 233–247. [Google Scholar] [CrossRef] [Green Version]
  80. Hahm, W.J.; Rempe, D.M.; Dralle, D.N.; Dawson, T.E.; Lovill, S.M.; Bryk, A.B.; Bish, D.L.; Schieber, J.; Dietrich, W.E. Lithologically controlled subsurface critical zone thickness and water storage capacity determine regional plant community composition. Water Resour. Res. 2019, 55, 3028–3055. [Google Scholar] [CrossRef]
  81. Winter, T.C.; Harvey, J.W.; Franke, O.L.; Alley, W.M. Ground Water and Surface Water: A Single Resource; USGS Publications North Dakota Water Science Center: Grand Forks, ND, USA, 1998. [CrossRef] [Green Version]
  82. Ferone, J.M.; Devito, K.J. Shallow groundwater–surface water interactions in pond–peatland complexes along a Boreal Plains topographic gradient. J. Hydrol. 2004, 292, 75–95. [Google Scholar] [CrossRef]
  83. Winkler, S.; Matthews, J.A. Observations on terminal moraine-ridge formation during recent advances of southern Norwegian glaciers. Geomorphology 2010, 116, 87–106. [Google Scholar] [CrossRef]
  84. Livingstone, S.J.; Lewington, E.L.M.; Clark, C.D.; Storrar, R.D.; Sole, A.J.; Mcmartin, I.; Dewald, N.; Ng, F. A quasi-annual record of time-transgressive esker formation: Implications for ice-sheet reconstruction and subglacial hydrology. Cryosphere 2020, 14, 1989–2004. [Google Scholar] [CrossRef]
  85. Pronk, A.G. Surficial Mapping in the Caledonia Highlands of Southern New Brunswick: Mineral Exploration and Land Use Applications of a Till Sampling Program. New Brunswick Department of Natural Resources and Energy, Minerals and Energy Division: Fredericton, NB, Canada, 1996. [Google Scholar]
  86. Phillips, J.D.; Marion, D.A.; Luckow, K.; Adams, K.R. Nonequilibrium regolith thickness in the Ouachita Mountains. J. Geol. 2005, 113, 325–340. [Google Scholar] [CrossRef]
  87. Wilson, R.A.; Parkhill, M.A.; Carroll, J.I. New Brunswick Appalachian Transect: Bedrock and Quaternary Geology of the Mount Carleton—Restigouche River Area; New Brunswick Department of Natural Resources: Bathurst, NB, Canada, 2005.
  88. Hodson, A.J.; Ferguson, R.I. Fluvial suspended sediment transport from cold and warm-based glaciers in Svalbard. Earth Surf. Process. Landf. 1999, 24, 957–974. [Google Scholar] [CrossRef]
  89. Cuffey, K.M.; Conway, H.; Gades, A.M.; Hallet, B.; Lorrain, R.; Severinghaus, J.P.; Steig, E.J.; Vaughn, B.; White, J.W.C. Entrainment at cold glacier beds. Geology 2000, 28, 351–354. [Google Scholar] [CrossRef]
  90. Small, E.E.; Anderson, R.S.; Repka, J.L.; Finkel, R. Erosion rates of alpine bedrock summit surfaces deduced from in situ 10BE and 26Al. Earth Planet. Sci. Lett. 1997, 150, 413–425. [Google Scholar] [CrossRef]
  91. Wan, J.; Tokunaga, T.K.; Williams, K.H.; Dong, W.; Brown, W.; Henderson, A.N.; Newman, A.W.; Hubbard, S.S. Predicting sedimentary bedrock subsurface weathering fronts and weathering rates. Sci. Rep. 2019, 9, 1–10. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  92. Broster, B.E.; Munn, M.D.; Pronk, A.G. Inferences on glacial flow from till clast dispersal, waterford area, New Brunswick. Geogr. Phys. Quat. 1997, 51, 29–39. [Google Scholar] [CrossRef] [Green Version]
  93. Leopold, L.B.; Wolman, M.G.; Miller, J.P.; Wohl, E. Fluvial Processes in Geomorphology; Dover Publications: Mineola, NY, USA, 1964; ISBN 0486685888. [Google Scholar]
  94. Ganong, W.F. Notes on the Natural History and Physiography of New Brunswick; Barnes &, Co.: Saint John, NB, Canada, 1906. [Google Scholar]
  95. Rempe, D.M.; Dietrich, W.E. A bottom-up control on fresh-bedrock topography under landscapes. Proc. Natl. Acad. Sci. USA 2014, 111, 6576–6581. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  96. Patton, N.R.; Lohse, K.A.; Godsey, S.E.; Crosby, B.T.; Seyfried, M.S. Predicting soil thickness on soil mantled hillslopes. Nat. Commun. 2018, 9, 1–10. [Google Scholar] [CrossRef]
  97. Truffer, M.; Harrison, W.D.; Echelmeyer, K.A. Glacier motion dominated by processes deep in underlying till. J. Glaciol. 2000, 46, 213–221. [Google Scholar] [CrossRef] [Green Version]
  98. Christoffersen, P.; Piotrowski, J.A.; Larsen, N.K. Basal processes beneath an Arctic glacier and their geomorphic imprint after a surge, Elisebreen, Svalbard. Quat. Res. 2005, 64, 125–137. [Google Scholar] [CrossRef]
  99. Montgomery, D.R.; Dietrich, W.E.; Torres, R.; Anderson, S.P.; Heffner, J.T.; Loague, K. Hydrologic response of a steep, unchanneled valley to natural and applied rainfall. Water Resour. Res. 1997, 33, 91–109. [Google Scholar] [CrossRef]
  100. Dietrich, W.E.; Perron, J.T. The search for a topographi signature of life. Nature 2006, 439, 411–418. [Google Scholar] [CrossRef] [PubMed]
  101. Monk, W.; Wilbur, N.M.; Allen Curry, R.; Gagnon, R.; Faux, R.N. Linking landscape variables to cold water refugia in rivers. J. Environ. Manag. 2013, 118, 170–176. [Google Scholar] [CrossRef]
  102. O’Sullivan, A.M.; Devito, K.J.; Ogilvie, J.; Linnansaari, T.; Pronk, T.; Allard, S.; Curry, R.A. Effects of topographic resolution and geologic setting on spatial statistical river temperature models. Water Resour. Res. 2020, 56, 1–23. [Google Scholar] [CrossRef]
Remotesensing 13 04210 g001 550
Figure 1. Visualization of the geographic location of New Brunswick, Canada (46.5653° N, 66.4619° W) outlining elevation change and physiographic regions (A) and an overview of bedrock lithology (B), overlain on World Ocean Basemap [30] in ArcGIS 10.3 software [31]. Bedrock geology was acquired from [28] (Scale 1:2,500,000).
Figure 1. Visualization of the geographic location of New Brunswick, Canada (46.5653° N, 66.4619° W) outlining elevation change and physiographic regions (A) and an overview of bedrock lithology (B), overlain on World Ocean Basemap [30] in ArcGIS 10.3 software [31]. Bedrock geology was acquired from [28] (Scale 1:2,500,000).
Remotesensing 13 04210 g001
Remotesensing 13 04210 g002 550
Figure 2. Overview of sample locations for each DTB source: boreholes (A); drillholes (B); pedons (C); site cards (D); well logs (E); and bedrock outcrops (F). Total sample size = 170,719 (see also Table 2).
Figure 2. Overview of sample locations for each DTB source: boreholes (A); drillholes (B); pedons (C); site cards (D); well logs (E); and bedrock outcrops (F). Total sample size = 170,719 (see also Table 2).
Remotesensing 13 04210 g002
Remotesensing 13 04210 g003 550
Figure 3. Histogram and continuous density distribution (red line) of bedrock depth for boreholes (A), drillholes (B), pedons (C), site cards (D), and well logs (E), and all data combined (F), including sample size, minimum, maximum, and mean DTB values for each source. Note: Depth intervals are (i) 0.25 m for pedons and site cards, and 5 m for the remaining sources, and (ii) the total sample size with minimum, maximum, and mean exclude the rock outcrop samples since all of these samples have a DTB of 0.
Figure 3. Histogram and continuous density distribution (red line) of bedrock depth for boreholes (A), drillholes (B), pedons (C), site cards (D), and well logs (E), and all data combined (F), including sample size, minimum, maximum, and mean DTB values for each source. Note: Depth intervals are (i) 0.25 m for pedons and site cards, and 5 m for the remaining sources, and (ii) the total sample size with minimum, maximum, and mean exclude the rock outcrop samples since all of these samples have a DTB of 0.
Remotesensing 13 04210 g003
Remotesensing 13 04210 g004 550
Figure 4. The framework for developing the DTB model covariates. The multi-scale generation of covariates beginning with the 10 m DEM is shown as an example (blue frame) and the correlation analyses that followed to reduce the number of covariates (red frame).
Figure 4. The framework for developing the DTB model covariates. The multi-scale generation of covariates beginning with the 10 m DEM is shown as an example (blue frame) and the correlation analyses that followed to reduce the number of covariates (red frame).
Remotesensing 13 04210 g004
Remotesensing 13 04210 g005 550
Figure 5. Scree plot used to select the number of PCs to capture a minimum of 90% variability (n = 11).
Figure 5. Scree plot used to select the number of PCs to capture a minimum of 90% variability (n = 11).
Remotesensing 13 04210 g005
Remotesensing 13 04210 g006 550
Figure 6. The RF variable importance plot highlighting the explanatory strength of the 20 most significant covariates with importance measured as decreasing from 100% (most important) to 40%. Covariate abbreviations are explained in Table 3.
Figure 6. The RF variable importance plot highlighting the explanatory strength of the 20 most significant covariates with importance measured as decreasing from 100% (most important) to 40%. Covariate abbreviations are explained in Table 3.
Remotesensing 13 04210 g006
Remotesensing 13 04210 g007 550
Figure 7. Actual vs. fitted model results for the training and validation subsets ((A,B), respectively) including 95% confidence limits, and histograms displaying densities of residual errors for DTB model on training and validation subsets ((C,D), respectively).
Figure 7. Actual vs. fitted model results for the training and validation subsets ((A,B), respectively) including 95% confidence limits, and histograms displaying densities of residual errors for DTB model on training and validation subsets ((C,D), respectively).
Remotesensing 13 04210 g007
Remotesensing 13 04210 g008 550
Figure 8. Predicted depth to bedrock (DTB) at 10 m2 resolution for the Province of New Brunswick, Canada (A) with an example of the fine-scale resolution showing sediment accumulation across an upland–hillslope–valley bottom transition (B).
Figure 8. Predicted depth to bedrock (DTB) at 10 m2 resolution for the Province of New Brunswick, Canada (A) with an example of the fine-scale resolution showing sediment accumulation across an upland–hillslope–valley bottom transition (B).
Remotesensing 13 04210 g008
Table
Table 1. Comparison of reviewed literature for predicting overburden depth-to-bedrock, DTB, where RF = Random Forest, GB = gradient boosting, ANN = artificial neural network, OK = ordinary kriging, and RK = regression kriging.
Table 1. Comparison of reviewed literature for predicting overburden depth-to-bedrock, DTB, where RF = Random Forest, GB = gradient boosting, ANN = artificial neural network, OK = ordinary kriging, and RK = regression kriging.
ReferencesLocationSpatial Coverage
(km2)		Resolution (m2)Ratio of Resolution/Area CoveredRecently Glaciated
(Y/N)		Statistical Procedure
Shangguan et al. (2017) [5]Global510,100,0002501.2 × 10−10PartialRF and GB
Pelletier et al. (2016) [13]Global510,100,00010002.0 × 10−3PartialAlgebraic Equations
Yan et al. (2018) [2]China9,579,0001001.0 × 10−9YRF and GB
Wilford et al. (2016) [15]Australian continent7,692,000902.3 × 10−10NCubist
Wilford and Thomas (2013) [16]Mt. Lofty range, South Australia1280108.0 × 10−8NCubist
Karlsson et al. (2014) [19]Three municipalities, Stockholm, Sweden98622.0 × 10−9YLinear Regression
Metelka et al. (2018) [14]Burkina Faso, West Africa686301.3 × 10−6NANN
Shafique et al. (2011) [17]Three cities, Kashmir, India470301.9 × 10−6NStepwise Regression
Devkota et al. (2018) [4]Phewa watershed, Central-Western Hills, Nepal111203.6 × 10−6NLinear Regression and OK
Kuriakose et al. (2009) [18]Western Ghats of Kerala, India9.5204.2 × 10−5NRK
Gomes et al. (2016) [3]Papagaio river basin, Rio de Janeiro, Brazil543.2 × 10−6NBayesian Modeling
Table
Table 2. Overview of DTB datasets acquired for analysis including source, sample size, and overall representation of each dataset including the overall dataset used for modelling. Note: Sample size represents all samples in the original datasets prior to data filtering.
Table 2. Overview of DTB datasets acquired for analysis including source, sample size, and overall representation of each dataset including the overall dataset used for modelling. Note: Sample size represents all samples in the original datasets prior to data filtering.
Data SetSourceSample Size% of Total
BoreholesNBERD3320.17
DrillholesNBERD98285.15
PedonsNBERD, CANSIS1990.10
Site CardsNBERD3240.17
Well LogsNBELG33,18717.4
OutcropsNBERD126,84976.99
Table
Table 3. The final set of covariates including the number of correlated covariates. Refer to Appendix A for full name.
Table 3. The final set of covariates including the number of correlated covariates. Refer to Appendix A for full name.
Covariate# of Correlated CovariatesCovariate# of Correlated CovariatesCovariate# of Correlated Covariates
Ca10518Max60515tan1056
Ca10051Max10052Tan3050
Crs6051Mca30521Tan6022
Cs101Min100523Tno100
Cur1059Min3020Tno3050
Cvgr31005Mrrtf10056Tot102
Dem1017Mrvbf10015Tscv100511
Diu6016Pln1020Tst10050
Flw308Prf3054Val1017
Flw1000Rsa102Val605
Flw10050Sky1014View3015
Grain_size6Sky1006Vis100512
Lng602Sky306Wb_dist0
Locdn101Swi10052Up300
Ls603Swi3013
Table
Table 4. Frequency distribution comparison of residual errors between training and validation subsets for the DTB model.
Table 4. Frequency distribution comparison of residual errors between training and validation subsets for the DTB model.
Range (m)TrainingValidation
Count%Count%
−12–−820.00200.000
−8–−420.00220.005
−4–−2500.053220.049
−2–069,73973.88830,11274.439
0–223,60425.008989824.469
2–46920.7332880.712
4–82280.2421040.257
8–12440.047210.052
12–24180.01940.010
24–3650.00520.005
36–7210.00110.001
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Furze, S.; O’Sullivan, A.M.; Allard, S.; Pronk, T.; Curry, R.A. A High-Resolution, Random Forest Approach to Mapping Depth-to-Bedrock across Shallow Overburden and Post-Glacial Terrain. Remote Sens. 2021, 13, 4210. https://doi.org/10.3390/rs13214210

AMA Style

Furze S, O’Sullivan AM, Allard S, Pronk T, Curry RA. A High-Resolution, Random Forest Approach to Mapping Depth-to-Bedrock across Shallow Overburden and Post-Glacial Terrain. Remote Sensing. 2021; 13(21):4210. https://doi.org/10.3390/rs13214210

Chicago/Turabian Style

Furze, Shane, Antóin M. O’Sullivan, Serge Allard, Toon Pronk, and R. Allen Curry. 2021. "A High-Resolution, Random Forest Approach to Mapping Depth-to-Bedrock across Shallow Overburden and Post-Glacial Terrain" Remote Sensing 13, no. 21: 4210. https://doi.org/10.3390/rs13214210

APA Style

Furze, S., O’Sullivan, A. M., Allard, S., Pronk, T., & Curry, R. A. (2021). A High-Resolution, Random Forest Approach to Mapping Depth-to-Bedrock across Shallow Overburden and Post-Glacial Terrain. Remote Sensing, 13(21), 4210. https://doi.org/10.3390/rs13214210

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop