Study areas

Soil profiles were described and sampled from three Saskatchewan study areas near Swift Current, St. Denis, and Smith Creek (Fig. 2). General characteristics for the study areas are summarized in Kiss and Bedard-Haughn (2021). The climate gradient spanned by the Swift Current, St. Denis, and Smith Creek study areas corresponds to the Brown, Dark Brown, and Black soil zones, respectively, which reflect general soil organic carbon storage trends influenced by the climatic gradient (Pennock et al. 2011). Similar soil catenas are found at the three sites, where Regosols and thinner Chernozems are commonly found in the eroded hillslope shoulder positions, Chernozems with thicker solums develop in mid-slope positions, and eluviated and gleyed Chernozems and Gleysols form in the depressional positions (Pennock et al. 2011). Soils in the Smith Creek area are commonly strongly enriched with CaCO₃ due to the limestone-rich parent materials (Saskatchewan Soil Survey Staff 1991). LiDAR-derived DEMs were available for all three sites (meta-data for these are described in the supplemental material of Kiss and Bedard-Haughn (2021)). Data collection was conducted in fall 2015 and 2016 for the Swift Current study area; fall 2015 for the St. Denis study area; and fall 2016 for the Smith Creek study area. By fall, water levels within the wetlands were at their lowest. Data collection and predictive mapping were restricted to areas mapped in the soil survey (Saskatchewan Land Resource Unit 2009) as containing glacial till parent material as it is most representative of the PPR landscape.

Fig. 2.

Study area locations within approximate extents of the Prairie Pothole Region (Ducks Unlimited 2005) and the Lake Winnipeg Watershed Basin (Government of Canada 2008). Training and testing sample point locations shown within each study area. Satellite RGB imagery for Swift Current from Agriculture and Agri-Food Canada (2009), satellite RGB imagery for St. Denis from Saskatchewan Geospatial Imagery Collaborative (2016), and satellite RGB imagery for Smith Creek from Copernicus (2017). [Colour online.]

Sampling design

Wetland discharge rings represent a small portion of the overall landscape; randomly placed sampling points within a PPR landscape are likely to miss these features, so alternative sampling approaches were required. Wetland soil characteristics only develop within depressions, but not all depressions develop wetland soil characteristics. Individual wetlands could be interpreted from water and vegetative extents present in aerial imagery to inform the sampling design, but those can change drastically over short periods of time as a function of climate and management. Instead, the focus of this study was on the topographic characteristics that influence the development of wetland soil characteristics. By basing the sampling design on depressions rather than only on confirmed wetlands, the data could provide information on the differences between depressions that form wetland soil characteristics and those that do not.

Individual depression boundaries were based on closed topographic depressions determined from the DEM, which are groups of raster grid cells sharing internal drainage (Lindsay 2016). Depressions were selected for sampling using a stratified random sampling design. Depressions were stratified based on size (area and depth) and relative elevation in the watershed (upper and lower) (Table 1). Size classes were based on Millar’s (1976) size classifications for wetlands in the Canadian PPR. The smallest depressions (≤4047 m²) were further stratified based on depth of the closed topographic depression to capture very small depressions in the sampling design. Relative elevation was defined by two categories, upper and lower halves of the watershed, which were determined using Jenks natural breaks classification method (Jenks 1967).

Table 1.

Criteria for depression size classes and number of depressions sampled per study area per relative position in watershed and size stratifications.

Sampled depressions were located in sites with a range of land uses including cultivated, pasture, and native prairie. Some of the smaller depressions within the cultivated sites were tilled. The wetland soil characteristics of gleying and discharge rings are expected to persist regardless of land use. Tillage erosion can result in redistribution of CaCO₃ to wetland positions, but this effect can be differentiated from natural, hydrologically influenced CaCO₃ distributions with field inspections by determining if the CaCO₃ distributions have abrupt horizontal boundaries that indicate they were redistributed through tillage. Wetland vegetation has been found to contribute to the formation of the discharge ring (Hayashi et al. 1998). However, the discharge rings developed over long periods of time (Knuteson et al. 1989), so the current vegetative characteristics do not necessarily reflect the long-term vegetative regimes for a depression. Wetlands that had undergone obvious mechanical drainage were avoided because the soils can be substantially disturbed through such processes. There has been extensive drainage of wetlands throughout the lower portion of the Smith Creek watershed and so sampling was done only in the upper half of that watershed (Table 1).

Each depression was sampled with a single transect with points spaced according to specified changes in elevation (dependent on the depth of the depression) with an average of five points per depression. Refer to Appendix A for a detailed description of the transect sampling design. Transect points were sampled with a truck-mounted hydraulic corer (Giddings Machine Company Ltd., Windsor, CO) where possible. Otherwise, transect points were sampled by hand auger. Sample point locations were collected with a Trimble GeoExplorer 2005 Series GeoXT GPS (Trimble, California, USA). At each point, sample cores were taken to at least 50 cm to facilitate classification and to 90 cm when possible. The profiles at each point were described and classified according to the Canadian System of Soil Classification (CSSC) (Soil Classification Working Group 1998) and were additionally classified as either calcareous wetland, noncalcareous wetland, upland, or buried/depositional soil types according to the criteria defined in Table 2, as the CSSC does not include classifications that reflect these soil types. Appendix A includes a detailed description and justification of the criteria that define the classifications used in this study.

Table 2.

Soil type classification criteria.

Model training datasets

The data collection took place over a period of relatively high precipitation for these areas (Brown et al. 2017a), which meant that the larger, more permanent ponds were often inundated with water. This limited sampling access to the wetland basin centers. An additional dataset from the St. Denis study area was used as supplementary training data to better reflect the soil class distributions in these landscape positions. This dataset was collected between 2000 and 2002 during a period of relative drought when many of the larger wetlands had dried out. A description of this dataset and the criteria used to classify its soil point observations according to the soil types of this study is provided in Appendix A.

Table 3 shows the number of observations per wetland soil type classification from each study area. These observations were used to train the predictive models. Buried/depositional soils were not considered in the models because it would not be possible to accurately predict their distribution based on the current elevation surface (Bedard-Haughn and Pennock 2002). There were fewer calcareous wetland soil type observations as compared to noncalcareous wetland soil type observations, which is representative as they are less abundant within a typical PPR landscape.

Table 3.

Number of observations per soil type in the training dataset.

Models were trained using two training data approaches—a generalized approach and a site-specific approach. For the generalized approach, training data from all three sites were used to train a predictive model and the same model was applied to predictively map the soil types at the three sites. For the site-specific approach, models were trained only using the data from the site to be mapped.

Independent validation datasets

Two datasets of soil point observations from previous studies were used as independent validation datasets to test model predictive accuracy. Independent validation soil point observations were sampled from wetland depressions that were not sampled as part of the model training datasets of this study. These included data from the BIOCAP study at St. Denis (Pennock et al. 2014) and the Brown et al. (2017a) study at Smith Creek. The BIOCAP dataset consisted of three wetlands sampled with two transects, each with 10–12 points per transect. The three sampled wetlands consisted of two recharge wetlands and a discharge wetland. This sample set was considered a rigorous test for the predictive models due to the high density of samples per wetland and the inclusion of a fully discharge wetland, which contained calcareous wetland soils throughout its basin.

The Brown dataset consisted of samples from 42 wetlands; however, 32 of those wetlands had been mechanically drained at some point. The data from the undrained and drained wetlands were separated into two validation sets. Each of the 42 wetlands was sampled with two points, one at the wetland toe-slope and one in a mid-slope position just above the discharge ring. The Undrained dataset did not include any calcareous wetland soil observations and the Drained dataset included only four. These datasets effectively tested the boundary surrounding the discharge ring. The number of observations per independent validation dataset are shown in Table 4. The point observations of the independent validation datasets included CSSC classifications (Soil Classification Working Group 1998) and profile descriptions that enabled their classifications into the target soil types of this study, except there were six observations that were left unclassified as the data required to distinguish them as either calcareous wetland or upland soil types were not available.

Table 4.

Number of soil type observations per independent validation set.

Spatial resolution

The PPR landscape is characterized by topographic variation at scales where hilltops and depressions can occur within 5–10 m of each other and very different soils develop in these positions (Pennock et al. 1987). These features would be smoothed out and lost in low-resolution DEMs and so the DEMs were kept at resolutions that correspond to the scale of variation for terrain derivative calculation. Separate models were tested using terrain derivatives generated at 2, 5, and 10 m to assess differences in predictive accuracy at varying spatial resolutions. Mapping was tested at 5 m spatial resolution to stay consistent with Bedard-Haughn and Pennock's (2002) PSM study on the distribution of PPR wetland soil types. Due to the specific sampling design used in our study, at a resolution of 5 m, many adjacent sampled point pairs shared the same pixel. Based on recommendations in Hengl (2006), a spatial resolution of 2 m was also tested because at this resolution, less than 5% of closest point pairs shared pixels. The 2 m DEMs provide a more detailed representation of the land surface and may capture features that would be smoothed out at 5 m resolution. Ten metre DEMs were also tested to assess if similar predictive performance could be achieved at coarser spatial resolution, which is typically easier to obtain and reduces computational workloads.

DEM processing

LiDAR collection for the three sites occurred either late summer or fall when water levels are lowest; however, at St. Denis and Smith Creek sites, standing water was still present in the largest wetlands. This caused the DEM surface in these positions to reflect the water surface and not the underlying sediment. The water surface causes a noisy DEM in these positions. A method for identifying and reinterpolating water surfaces in the DEM is described in Appendix A, as are the methods used to smooth and resample the DEM to the target spatial resolutions following reinterpolation.

Predictor variables

As the focus of this study was to determine the effectiveness of using terrain derivatives derived from high-resolution DEMs as predictor variables for predictively mapping wetland soil types and to determine which terrain derivatives were most useful, predictor variables were limited primarily to terrain derivatives. It is likely that the predictive mapping of wetland soil types would be improved by including predictor variables derived from other remote and proximal sensing data, like synthetic aperture radar (SAR), optical remote sensing data, band indices, electromagnetic surveys, etc. but they were not included to enable more straightforward evaluation of terrain derivatives as predictor variables.

A total of 35 predictor variables were calculated at each DEM resolution (2, 5, and 10 m) (Table 5). Most of these were calculated using the System for Automated Geoscientific Analysis (SAGA) (Conrad et al. 2015) or Whitebox Geospatial tools (Lindsay 2014). Measures of Euclidean distance were included based on Behrens et al. (2018) and were calculated using R code published by Saurette (2021). Cartographic depth-to-water (Murphy et al. 2009) was calculated with GRASS GIS (GRASS Development Team 2017) using R code published by Schönauer and Maack (2021).

Table 5.

List of predictor variables and their source.

Several terrain derivatives were specifically developed to better represent the unique morphological characteristics of the PPR that were expected to relate to wetland soil type distributions (Fig. 3). Three of these derivatives were based on a grid cell's position relative to a local depression. These were developed to capture relationships identified by Pennock et al. (2014) between soil type distribution and topographic position relative to depression bottoms and maximum water levels.

Fig. 3.

3D scene for an example area in the Swift Current study site (3× vertical exaggeration based on DEM) showing an air photo and the terrain derivatives developed for this study. [Colour online.]

The contributing watersheds were determined for each closed topographic depression within the DEM. Grid cells were related to their local depression based on the contributing basin they fell within. Values were calculated for a grid cell’s elevation above the associated depression’s basin minimum (elevation from depression minimum derivative) and elevation above or below the depression's spillover elevation (elevation from depression spillover derivative). The spillover elevation reflects the maximum elevation of a closed depression before water would spill externally from it. For the depression depth derivative, the depression depth value was ascribed to all pixels within the depression's contributing basin. The measure of elevation from depression minimum per depth was included to potentially capture differences in relationships of soil distribution for wetlands of varying sizes. The watershed max catchment area derivative ascribes the maximum catchment area for a depression’s contributing area to all pixels within that contributing area.

Wetland Strahler order minimum and maximum were calculated for each depression following the methodologies described in Kiss and Bedard-Haughn (2021). For these variables, the characteristics were quantified for each closed depression; all pixels within the depression were ascribed the same value and all pixels not included in closed topographic depressions were given values of zero. R scripts for calculating the terrain derivatives developed in this study and parameters used to calculate existing derivatives are available on GitHub (Kiss 2022). Both the newly developed terrain derivatives and existing derivatives were calculated using a suite of Whitebox Geospatial (Lindsay 2014) and SAGA GIS (Conrad et al. 2015) tools.

Prior to model training, recursive feature elimination was conducted to reduce the number of predictor variables used in the final models to reduce the risk of model overfitting. Feature reduction has been found to improve model performance in PSM studies focused on soil class mapping (Brungard et al. 2015; Hounkpatin et al. 2018). The recursive feature elimination was implemented using the caret package (Kuhn 2008) in R (R Core Team 2020), using a random forest model as the classifier and 5-fold cross validation with five iterations using Kappa score as the evaluation metric to determine the optimal predictor variables to include in the final models. Importance of each predictor variable in the final models was determined using the Gini index (Wright and Ziegler 2017). Overall variable importance was determined by summing the Gini index of each variable within the final random forest models.

Model training, tuning, and validation

The objective of this study was to assess the performance of various terrain derivatives as predictor variables for mapping wetland soil types and not to determine the optimal machine learning model type for this problem. Random forest is a tree-based machine learning model capable of modelling complex nonlinear relationships between numerous predictor variables and categories (Strobl et al. 2009). In addition to its frequently cited success in the PSM literature (Brungard et al. 2015; Heung et al. 2016; Jeune et al. 2018; Assami and Hamdi-Aїssa 2019), random forest models provide an intuitive evaluation of predictor variable importance, which makes it an ideal candidate for this study.

Following the recursive feature elimination step, the caret package (Kuhn 2008) was used again to optimize model parameters and train the final random forest models through K-folds cross-validation. Each model was trained with the training dataset using 5-fold cross-validation with 20 iterations using Kappa score as the evaluation metric. From this, cross-validation accuracy and Kappa scores were determined as well as the optimal mtry parameter setting per model, which is the number of predictor variables available for each split in a random forest model. Each random forest was grown with 1000 trees (ntree).

Model performance was assessed based on overall accuracy, Kappa score, producer’s accuracy, and user’s accuracy as determined from the cross-validation and independent validation tests. Cross-validation performance was assessed in addition to independent validation tests because independent validation data were not available for the Swift Current study area, so cross-validation performance allowed for a fairer assessment of model performance across all three sites. Overall accuracy refers to the total number of correct classification predictions per the total number of observations. Kappa score corrects the overall accuracy by accounting for the possibility of a correct prediction based on chance alone (Landis and Koch 1977). Producer’s accuracy, which is reported per soil type class (Table 8), refers to the number of correctly predicted soil classifications of a particular soil type per the total number of observations of that soil type (Malone et al. 2017). User’s accuracy, also reported per soil type, refers to the number of correctly predicted soil classifications of a soil type per the total number of predictions of that soil type.

Visual assessment of model performance

The predictive accuracies based on the internal and independent validations provide an excellent performance summary of the predictive models. However, it is crucial to ensure that in addition to achieving acceptable accuracy scores, model predictions of soil distributions make conceptual sense. Maps were generated for each site from each of the models and were visually inspected using expert knowledge to determine if the distribution of soils matched the concepts described in the introduction. This involved inspecting areas at each site where wetland types were known either from previous studies or from the sampling conducted in this study.

For most wetlands, inspection was conducted to determine if noncalcareous wetland soil types were predicted to be distributed throughout the wetland basin floor with calcareous wetland soils formed in a discharge ring pattern in the foot to mid slope positions surrounding the wetland basin. At the St. Denis study area, there were numerous wetlands known to have calcareous wetland soils distributed not only in discharge rings surrounding the basin, but throughout their wetland basins; these were inspected to determine if calcareous wetland soils were predictively mapped throughout these wetland basins. Distributions of upland soils mapped within known wetlands or wetland soils mapped in locations well beyond the boundaries of discernable wetlands were deemed erroneous. Soil types are expected to be distributed continuously and therefore single pixel or small clusters of soil types were typically considered erroneous artefacts, except for the models generated based on predictor variables of 10 m resolution, as wetland depressions may be small enough to be represented by a single 10 m pixel.

Tuned model parameters and internal cross-validation results

Table 6 shows the optimized model parameters and cross-validation accuracy per model. Based on the cross-validation results, the models developed on predictor variables of 2 m resolution had the highest predictive accuracies, but cross-validation accuracy was generally similar between the models based on different resolutions (ranging from 75% to 82% accurate) (Table 6). According to the cross-validation results, the site-specific models performed slightly better than the generalized models; the average Kappa scores across all three sites (All (mean)) for the site-specific models were slightly higher than the generalized models at all resolutions (Table 6). Overall, there were only minor differences in performance across all models shown in Table 6 based on the cross-validation evaluation. The models vary widely in the total number of predictor variables used (from 7 to all 35). There were no noticeable trends in the difference of performance of models that incorporated many or few predictor variables.

Table 6.

Model parameters, cross-validation overall accuracy, and Kappa score per model per site.

Independent validation results

The models were tested for predicting the soil types of the independent datasets from the St. Denis and Smith Creek study areas. Table 7 shows the overall accuracy and Kappa scores based on the independent validation tests for the models. Table 8 shows the producer’s and user’s accuracy for each soil type classification for the BIOCAP–St. Denis and Brown–Smith Creek–Undrained independent validation tests for the models.

Table 7.

Overall accuracy and Kappa score per model per independent validation test sets.

Table 8.

Producer and user accuracy per soil type classification per model per independent validation test set (Brown–Smith Creek–Drained not included).

Independent validation results—BIOCAP–St. Denis

The best-performing model for the BIOCAP–St. Denis dataset was the 2 m site-specific model, which was built using predictor variables with a spatial resolution of 2 m and used only training data from the St. Denis study site. This model achieved an overall accuracy of 76% and a Kappa score of 0.63. The site-specific models consistently performed better than the generalized models for the St. Denis independent validation test.

The models consistently had low producer's accuracy and high user’s accuracy for calcareous wetland soil types. This communicates that this soil type is underestimated but where it is predicted, the map user can have confidence they will find that soil type there. The 2 and 10 m site-specific models achieved the most balanced results for the calcareous wetland soil types with producer’s accuracies of 54% and 65% and use’s accuracy of 94% and 90%, respectively. Conversely, the models had very high producer’s accuracy and low user’s accuracy for noncalcareous wetland soil types. This means that all noncalcareous wetland observations were correctly predicted as such, but that they were overpredicted elsewhere.

Independent validation results—Brown–Smith Creek

The models generally had high accuracies and Kappa scores for predicting the soil types of the Undrained independent validation set; some achieving perfect (100% accurate, Kappa score = 1.0) and near perfect (90% accurate, Kappa score = 0.80) results (Table 7). The 10 m resolution models had the lowest accuracies for this validation test (ranging from 84% to 85% accurate). Site-specific and generalized models performed relatively similar for the Undrained dataset, but for the Drained test set, the generalized models consistently performed better, with some Kappa scores 0.30 greater. Predictions for the Drained validation set were generally less successful than for the Undrained validation set but still typically had high levels of agreement. This indicates that despite the impact drainage may have on the landscape (and the derived DEM), the distribution of wetland soils may still be predictable in landscapes that have undergone mechanical drainage, meaning these methods could potentially be used in restoration projects in drained landscapes to identify historic wetland extents.

The Brown–Smith Creek–Undrained validation dataset did not contain any calcareous wetland soil type observations due to the sampling design. Therefore, models that underpredicted the extent of calcareous wetland soil types were likely to perform well for this dataset. The 2 and 5 m resolution-based models performed similarly, sometimes overpredicting noncalcareous wetland soil types and sometimes overpredicting upland soil types (as indicated by high producer’s accuracy and lower user’s accuracy (Table 8)).

Predicted soil type maps—visual assessments

Predicted soil maps for the full study areas from each model for St. Denis and Smith Creek are shown in Figs. 6 and 8. The full map of the Swift Current study area is not included as differences in soil distributions between the models are not obvious at coarse scale. Full-sized maps for each study area from each model are provided in the   Supplementary file (cjss-2022-0034suppla.pdf). Figures 4, 5, and 7 show predicted maps for example areas from each site.

Fig. 4.

3D scene (3× vertical exaggeration based on DEM) of predicted soil type maps and air photo for an example area at the Swift Current study site based on the generalized and site-specific models built at each spatial resolution. Satellite RGB imagery from Agriculture and Agri-Food Canada (2009). [Colour online.]

Fig. 5.

3D scene (3× vertical exaggeration based on DEM) of predicted soil type maps and air photo for an example area at the St. Denis study site based on the generalized and site-specific models built at each spatial resolution. Satellite RGB imagery from Saskatchewan Geospatial Imagery Collaborative (2010). [Colour online.]

Generally, the predicted maps generated by the models were acceptable; wetland soils (calcareous and noncalcareous soil types) were typically mapped within or directly adjacent to depressional positions and upland soils were mapped in upslope positions. Wetland basin centers were often mapped with noncalcareous wetland soils with rings of calcareous wetland soils mapped surrounding the wetlands reflecting the discharge ring phenomenon. Wetland soils were most frequent at the Smith Creek study area and least frequent at the Swift Current study area, which coincides with the general climatic moisture trends for these sites. The St. Denis study area had the most wetlands predicted to have calcareous wetland soils throughout their basins, which is consistent with historic site observations (Pennock et al. 2014). Calcareous wetland soil types were much less frequently predicted within the Swift Current study area, which is consistent with the small amount of calcareous wetland soils encountered in the field data collection (Table 3).

At the Swift Current study area, with exception to the 10 m based models, all models generated relatively consistent maps (Fig. 4). The largest wetlands were mapped to have calcareous wetland soils surrounding the wetland basin, but most other wetlands were mapped with only noncalcareous wetland soil types. The generalized models predicted more occurrences of calcareous wetland soil types than the site-specific models. The maps generated by the 10 m based models typically had more single-cell and small cluster distributions of soil types compared to those of the finer-scale resolutions. These scattered soil type predictions were often wetland soil types erroneously mapped in upland positions.

The predicted maps generated for the St. Denis study area were relatively consistent across all models, except that the site-specific models predicted more frequent occurrences of wetland soils generally and more frequent occurrences of calcareous wetland soils (Fig. 6). This contributed to its better performance when predicting for the BIOCAP validation test set. When mapping at 10 m, the rings of calcareous wetland soils surrounding the wetlands become blocky and less detailed as a result of the pixel size. As discussed, larger depressions were more frequently predicted to have calcareous wetland soils within their basin floors at the St. Denis study area. All models predicted Pond 90 (the largest wetland in the top right of the example area scene (Fig. 5)) to be predominantly covered with calcareous wetland soils, which is appropriate as the largest low-lying wetlands in this portion of the St. Denis watershed are known to have calcareous wetland soils throughout their basins (Pennock et al. 2014). However, resulting maps often include irregular distributions of noncalcareous soils within these basins (as seen in ponds adjacent to Pond 90 in Fig. 5), which is not expected to occur. Additionally, not all wetlands at St. Denis known to have calcareous wetland soils throughout their basins were consistently predicted as such.

Fig. 6.

Predicted soil type maps for the full St. Denis study site based on the generalized and site-specific models built at each spatial resolution. Full-sized versions of each map are provided in supp1a. Satellite RGB imagery from Saskatchewan Geospatial Imagery Collaborative (2016). [Colour online.]

Compared to the other sites, the resulting maps for the Smith Creek study area differed substantially between some models, especially in the eastern portion of the study area where training data were lacking (Fig. 8). The site-specific models predicted greater occurrences of calcareous wetland soil types and wetland soils generally (Figs. 7 and 8). Noncalcareous wetland soils appear to extend beyond the wetland boundaries observable in the air photo in the maps generated by the site-specific models (Fig. 7) and may be overestimated. It is possible that this is appropriate in an area as dense with wetlands as shown in the example area in Fig. 7; however, in the eastern portion of the study area where wetland drainage is more common, it is likely that wetland soils are being overestimated. This explains the poorer performance of the site-specific models for predicting the classes for the Drained validation test set.

Fig. 7.

3D scene (3× vertical exaggeration based on DEM) of predicted soil type maps and air photo for an example area at the Smith Creek study site based on the generalized and site-specific models built at each spatial resolution. Satellite RGB imagery from Saskatchewan Geospatial Imagery Collaborative (2014). [Colour online.]

Fig. 8.

Predicted soil type maps for the full Smith Creek study site based on the generalized and site-specific models built at each spatial resolution. Full-sized versions of each map are provided in supp1a. Satellite RGB imagery from Copernicus (2017). [Colour online.]

The general high accuracies of each model based on the independent validation tests despite the apparent differences in resulting maps is likely due to the nature of the validation test sets from the area. Point observations in these datasets were collected from distinctly different slope positions and were likely easier to predict than if the points spanned the gradient between these positions. The questionable results generated by the site-specific models for the eastern portion of the study area became apparent only through the visual assessment, as they achieved acceptable predictive accuracies based on the validation datasets.

The map generated by the 10 m based generalized model for Smith Creek contains many instances of single pixel and small clusters of the soil types which results in a speckled pixelated-looking output as seen in Fig. 7e. The map generated by the 10 m based site-specific model follows generally the same trends mapped by the 2 m site-specific model but because of its coarseness, the mapped soil type distributions do not follow the contours of the wetlands and depressions with the same level of precision (Figs. 7b and 7f).

Variable importance

Figure 9 shows a plot of the overall importance of every predictor variable within the random forest models. Of the top 10 most important variables, three were developed specifically for this study (elevation from depression minimum per depth, elevation from depression spillover, and elevation from depression minimum) and two were developed for the related study focused on predicting wetland solute-richness classes in PPR landscapes (Wetland Strahler order minimum and maximum) (Kiss and Bedard-Haughn 2021). Elevation from depression minimum per depth and elevation from depression spillover were the most important predictor variables overall. Four of the top 10 predictor variables were calculated using Whitebox Geospatial tools, including stochastic depression analysis (Lindsay and Creed 2005) and maximum elevation deviation (local, meso, and broad scales) (Lindsay et al. 2015). SAGA-calculated terrain derivatives slope height and normalized height (t = 1000) (Böhner and Selige 2006) were the seventh and 11th most important predictor variables, respectively. Euclidean distance from NE site corner, Watershed max catchment area, and Aspect (northness) had the lowest Gini index scores, indicating that they were the least important predictor variables.

Fig. 9.

Predictor variable importance based on Gini index summed across random forest models.

Overall model performance

The relatively high overall accuracies and Kappa scores according to the cross-validation and independent validation tests, and additionally the acceptable map outputs as determined through visual assessment for many of the models indicate their success in predictively mapping wetland soil types at the PPR study sites. Based on all model evaluation criteria, the best-performing model overall was the 2 m site-specific model, which achieved overall accuracies and Kappa scores ranging from 80%–82% and 0.62–0.72, respectively, based on cross-validation (Table 6) and 76%–100% and 0.63–1.0, respectively, based on independent validation tests (Table 7). The related study focused on predictively mapping wetland solute-richness classes at the same PPR study sites achieved similar accuracies: independent validation accuracy of 69%–80% and Kappa scores of 0.40–0.41 (Kiss and Bedard-Haughn 2021). The model developed for that study was simpler and based on only a few terrain-derived predictor variables and its objective was only to predict a single solute-richness class for each wetland, whereas in this study, detailed mapping of soil types within and surrounding individual wetlands was accomplished.

Recent studies focused on predictively mapping wetland classes achieved similar or higher overall accuracies based on independent validation tests: Mahdianpari et al. (2017) achieved 94% overall accuracy predictively mapping land cover types (primarily wetland classes) in the Avalon Peninsula of Newfoundland using remote sensed multispectral imagery; Jahncke et al. (2018) achieved overall accuracy of 89.2% mapping wetland classes in an area south of Halifax, Nova Scotia by using various remote sensed data, including LiDAR data; and Millard and Richardson (2014) achieved overall accuracy of 72.8% predictively mapping land cover classes in a peatland region near Ottawa, Ontario using SAR and LiDAR data. These studies focused on mapping wetland classes that were typically defined by their vegetative characteristics, which can be closely correlated to spectral reflectance bands and vegetative indices, making them more easily predicted. Achieving similar overall accuracies when predictively mapping soil classes, particularly in agricultural landscapes, is often more difficult.

Other studies focused on soil class mapping in North America achieved similar or poorer overall accuracies based on independent validation tests to those achieved in this study: Goldman et al. (2020) achieved overall accuracies of 77.1% and 70.6% mapping soil drainage groups and texture groups, respectively, in a low-relief wetland-rich landscape in the Choptank River Watershed in Maryland and Delaware, USA; Heung et al. (2016) achieved overall accuracies of 70% and 58% mapping soil taxonomic great groups and soil orders, respectively, in the Lower Fraser Valley, British Colombia; and Brungard et al. (2015) achieved Kappa scores of 0.53, 0.32, and 0.19 mapping soil taxonomic subgroups at study sites in Wyoming, New Mexico, and Utah, USA, respectively.

The success of the predictive models of this study indicates that it is possible to accurately predict the spatial distribution of the targeted wetland soil types using primarily terrain derivatives. This corresponds to the findings of Pennock et al. (2014) who were able to relate key hydropedological features to terrain and elevation-based characteristics due to the substantial influences of topography and hydrology on wetland soil formation in this region. This is not to discount the importance of the SCORPAN model (soil as a function of soil, climate, organisms, relief, parent material, and spatial position) (McBratney et al. 2003), which has proven to be a guiding principle for predictive digital soil mapping throughout the world. The predictive models built in this study would likely improve by incorporating predictor variables that more directly reflect the other elements of SCORPAN like remote sensing multispectral imagery, climate data, and parent material information. The success of the discussed wetland class studies that included remote sensed data beyond terrain data highlights their importance for mapping wetland resources; studying the effects of incorporating these types of features as predictor variables is a crucial next step for research on mapping wetland soils and wetlands more generally in the Canadian PPR.

Predictor variable importance

Terrain derivatives generated for this study were some of the most important predictor variables for accurately mapping wetland soil types, specifically the elevation from depression minimum per depth and elevation from depression spillover derivatives, which were the first and second most important predictor variables, respectively. Pennock et al. (2014) found that the boundary between the distribution of recharge (noncalcareous) wetland soils within wetland centers and the surrounding discharge ring of calcareous wetland soils closely matched the extent of maximum recorded water levels. Their finding was determined at St. Denis where pond water levels have been recorded annually since 1968. This information does not exist for most locations but elevation from depression spillover works as a proxy for this information because the depression spillover elevation reflects the maximum possible water level for a wetland. The elevation from depression minimum per depth derivative modified the elevation from depression minimum derivative to reflect different relationships between depression size and soil distributions. A gleysolic soil may form on the edge of a large wetland at an elevation of 1 m above the depression bottom because a large wetland is likely to store enough water to inundate that position, whereas at a position of 1 m along a hillslope above a small depression, it is less likely that a gleysolic or wetland-related soil will form (Fig. 10).

Fig. 10.

Schematic diagram showing varying values for the elevation from depression minimum per depth given a point that is 1 m above the depression minimum for a shallow depression (left) and a deep depression (right). Different soils are expected to form in these two positions. The elevation from depression minimum per depth derivative modifies the elevation from depression minimum derivative by dividing by the depression's overall depth, which helps to quantify the relationships between these two characteristics for depressions of varying sizes. [Colour online.]

Wetland Strahler order calculations were developed and successfully used to map the distributions of solute-rich wetlands in the same regions of the Canadian PPR (Kiss and Bedard-Haughn 2021). Solute-rich wetlands typically correspond to those with calcareous wetland soils throughout their basin, so these derivatives were likely important in the random forest models for differentiating between wetlands that would and would not form such soil distributions.

The maximum elevation deviation derivatives and the normalized height and slope height derivatives aim to quantify a grid cell’s positions in relation to the local and (or) broader landscape. These can reflect similar characteristics to the depression-focused terrain derivatives developed for this study like the elevation from depression minimum derivative but do not limit the calculating window to the watershed of each depression, which helps to better reflect nuanced pedological relationships at the watershed edges. Additionally, hydropedological forces occur over a range of scales (Pennock et al. 2014) and features that capture characteristics reflective of varying scales, like these, are clearly important for hydropedological mapping.

Other studies focused on mapping wetland soils and wetland classes incorporated similar terrain derivatives that quantify a grid cell’s position in relation to its surrounding landscape and found them to be important predictor variables. Goldman et al. (2020) incorporated the deviation from mean elevation derivative, which calculates the difference in a grid cell’s elevation compared to the mean elevation from its 200 m radius divided by standard deviation. This was consistently one of the top two most important predictor variables in their study on mapping soil drainage groups and texture classes in a wetland-rich landscape in the Atlantic Coastal Plain of the USA. This derivative is similar to the maximum elevation deviation derivatives used in this study, except that it only considers the deviation compared to the single spatial scale of 200 m rather than across multiple spatial scales. Jahncke et al. (2018) found that, in addition to slope, topographic position index was an important predictor variable for mapping wetland classes in Nova Scotia. Topographic position index assigns an index to a grid cell according to its relative position compared to local valleys and hilltops. The normalized height and slope height derivatives similarly quantify a grid cell’s position compared to local valleys and hilltops but can account for that relationship across various neighbourhood scales.

Whitebox Geospatial’s stochastic depression analysis (Lindsay and Creed 2005) derivative (fourth most important variable) estimates the likelihood that a position belongs to a closed topographic depression. Depressions were accurately captured in the outputs generated by this tool, which were likely useful for in distinguishing the distributions of wetland and upland soils in the random forest models. This terrain derivative is recommended for depression or wetland landscape studies in this region.

Spatial resolution

Predicted maps generated based on predictor variables derived at the highest tested resolution (2 m) generally had the highest predictive accuracies based on cross and independent validation. The resulting soil maps generated at this resolution were acceptable and matched the established understanding of the distribution of these soil types without excessive errors or artefacts. Calculating terrain derivatives and applying predictive models to stacks of predictor variable rasters at this resolution is computationally intensive, which is an important consideration. Predictive mapping at 5 and 10 m resolution is exponentially less computationally intensive. Model performance at these resolutions was only slightly poorer than mapping at 2 m resolution. Predictive accuracy for cross and independent validation tests were similar between the 5 and 10 m based models, but visual inspection determined that the maps generated at 10 m resolution typically lacked the detail achieved at higher resolutions and soil types were often distributed in blocky and discontinuous patterns. Maps generated at 10 m resolution were still considered successful but working at higher resolutions produced better outputs generally.

Generalized versus site-specific models

The site-specific models performed slightly better overall compared to the generalized models. They had higher cross-validation Kappa scores on average at all resolutions and generally achieved higher or similar accuracy based on independent validation tests, with the exception for the Brown–Smith Creek–Drained independent validation test, in which the generalized models consistently had higher accuracy.

A notable outcome of the generalized models is that the resulting maps were more consistent across sites than the site-specific models. In this case, this effect may be not desirable because the distribution and frequency of wetland soils at each study area varies substantially. Wetland soils and calcareous wetland soils occur less frequently at the Swift Current study site, which causes the generalized models to underpredict these soil types at the other study sites. The generalized models did, however, achieve higher predictive accuracy according to the Smith Creek–Drained validation test and produced output maps considered more appropriate in the eastern portion of the Smith Creek study area because they predicted fewer occurrences of wetland soils there. The poorer performance of the site-specific models in this instance is likely due to the lack of training data in that region of the study area, as training data were only collected in the western and northern portions (Fig. 2). With better distribution of training data across the entire study area, it is expected that the site-specific models would perform better. This highlights an advantage of the generalized models: because they were built based on more generalizable patterns, they were more generally applicable and are less variable based on local training data.

Most of the terrain derivatives used as predictor variables are relative measures and cannot be used to differentiate between sites in the generalized models. The exception to this is the elevation variable, as each site falls within nearly entirely separate bands of elevation values (Swift Current: 706–832 m; St. Denis: 540–600 m; Smith Creek: 520–545 m). Elevation was a relatively useful predictor variable (Fig. 9), which was likely due to its ability to split the training data based on study area within the random forest models. Incorporating more predictor variables that can be used to distinguish sites may improve generalized model performance. Incorporating climatic data would provide a variable that directly relates to a key difference between these sites: their moisture regimes.

In the case where predictive mapping models can be trained using either only local or local and global training data, it is recommended to test both and determine on a case-by-case basis which performs better. Extrapolating predictive models to areas where no training data were collected should be done so cautiously. In the case of this study, a model trained using only data from Swift Current is very likely to underpredict the extent of wetland soils at Smith Creek and vice versa. Although St. Denis represents something of a midpoint in terms of frequency of wetland soils, models trained at the other two sites would not be capable of predicting the fully calcareous wetlands at this site.