Improved parent material map disaggregation methods in the Saskatchewan prairies using historical bare soil composite imagery

Abstract The major drivers of soil variation in Saskatchewan at scales finer than the existing soil maps are parent material variance, slope position, and salinity. There is therefore a need to generate finer-scale parent material maps as part of updating soil maps in Saskatchewan. As spatially referenced soil point data are lacking in Saskatchewan, predictive soil mapping methods that disaggregate existing soil parent material maps are required. This study focused on investigating important environmental covariates to use in parent material disaggregation, particularly bare soil composite imagery (BSCI). Synthetic point observations were generated using an area-proportional approach based on existing soil survey polygons and a random forest model was trained with those synthetic observations to predict parent material classes. Including BSCI as environmental covariates increased model accuracy from 0.38 to 0.52 and the model Kappa score from 0.19 to 0.35 compared with models where it was not included. Models that included training points from all locations, regardless of whether BSCI was available, and included BSCI as environmental covariates had similar results to the BSCI model with an accuracy of 0.48 and a Kappa value of 0.30. Based on these results, BSCI is an important covariate for parent material disaggregation in the Saskatchewan Prairies. Future work to disaggregate soil classes based on slope position and salinity, and to combine those methods with parent material disaggregation is needed to generate detailed soil maps for the Canadian Prairies.


Introduction
The Province of Saskatchewan had extensive surveying and soil mapping activities conducted between 1958 and 1998. The resulting soil survey maps were produced at scales of 1:100 000 or coarser and are available in an easily accessible digital platform (SKSIS Working Group 2018). While these maps were essential for agriculture management and land use planning in Saskatchewan throughout the late 20th century, there is increasing demand for more detailed soil maps to support targeted land use planning, soil carbon management, and precision agriculture. To meet these needs, predictive soil mapping using machine learning tools, in combination with a suite of environmental covariates, has been an increasing focus of research in Saskatchewan and around the world (McBratney et al. 2003).
As there is currently limited public spatially explicit point data for predictive soil mapping in Saskatchewan, there is a need for mapping approaches that do not rely solely on field-collected point data for model training. Given the extensive, coarse-scale, soil maps available, there is the potential to disaggregate those maps using established techniques to generate improved soil maps across the agricultural region of Saskatchewan. Polygon disaggregation by generating synthetic training points based on polygon labels is one such approach (Holmes et al. 2015;Chaney et al. 2016). Initial polygon disaggregation approaches focused on assigning equal points to polygons and randomly assigning class labels based on their relative proportions within a polygon (Odgers et al. 2014). Improvements have since been made to polygon disaggregation by incorporating area-proportional sampling and informing class assignment based on soil-landscape relationships (Møller et al. 2019). Additionally, the use of random forest models with a single sampling procedure was found to be much more computationally efficient than approaches that used multiple sampling procedures and C5.0 decision trees, with only a slight decrease in predictive accuracy (Møller et al. 2019).
Parent material is one of the key soil-forming factors (Jenny 1941) and a major control of important soil properties in the Canadian Prairies. Currently, in Saskatchewan, the finest-scale parent material maps are the classifications included in Saskatchewan's detailed soil survey. More detailed parent material maps are necessary to map soil variation at finer scales. A previous study in Saskatchewan identified that finer-scale mapping of hydropedological parameters can delineate soil variation associated with topography, but this typically requires high-resolution digital elevation models (DEMs), which are not widely available in Saskatchewan (Pennock et al. 2014;Kiss and Bedard-Haughn 2021). Successful disaggregation of parent material maps may not depend on such data requirements, making it a potential approach to improve soil maps across the Saskatchewan prairies. There have been successful studies in disaggregating existing parent material maps in British Columbia (Heung et al. 2014;Bulmer et al. 2016), which did not rely on high-resolution DEMs for model environmental covariates.
One potential approach to improve soil parent material map disaggregation in the Saskatchewan Prairies is to incorporate bare soil composite imagery (BSCI) as an environmental covariate. BSCI is generated by mining archival remote sensing imagery to isolate pixels where bare soil is present and using the resulting data stack to generate a composite image of bare soil conditions. The creation of BSCI has become a possibility with the development of large-scale cloud computing infrastructure (Gorelick et al. 2017). This is particularly necessary in the Canadian Prairies as pixels with vegetation or dead plant residue cover that exceed 20% make direct estimation of soil properties difficult (Bartholomeus et al. 2007(Bartholomeus et al. , 2008. Vegetation or residue-free images are rare in Saskatchewan in recent decades due to the widespread adoption of conservation tillage (Statistics Canada 2017). BSCI has been used in a variety of locations such as Germany (Rogge et al. 2018), Brazil , and globally  for predicting soil properties such as soil organic carbon. Recently, BSCI has been used to successfully generate historical soil carbon maps and clay content maps for Saskatchewan . Both soil organic carbon and clay content vary with parent material, along with other factors, and therefore there is potential that BSCI could support parent material map disaggregation in Saskatchewan.
This study had two main objectives. The first objective was to evaluate the potential of disaggregation approaches to generate finer-scale parent material maps in Saskatchewan, specifically by using a more computationally efficient random forest model approach with a single set of synthetically generated training data derived from existing soil survey maps of soil parent materials. The second objective of this study was to investigate the value of BSCI for improving parent material disaggregation in Saskatchewan.

Soil parent material context
Nearly the entire extent of Saskatchewan's soils is formed on unconsolidated transported parent materials, predominantly the result of glacial processes (Anderson and Cerkowniak 2010). The most common soil parent materials in Saskatchewan are glacial tills (often enriched with limestonerich bedrock), followed by glaciolacustrine and fluvial deposits (Table 1). Due to the overall slope of the land and drainage to the northeast, large glacial lakes formed that caused extensive glaciolacustrine deposits associated with proglacial lakes (Anderson and Cerkowniak 2010). Additionally, high-energy streams from the melting glaciers were responsible for the creation of fluvial deposits along deltas that formed where streams entered glacial lakes. Minor areas of eolian landscapes are also present, which formed due to postglacial wind erosion. Minor areas of bedrock, eolian, peat, and colluvium are also present. The existing parent material map, from Saskatchewan's detailed soil survey polygons, is provided in Fig. 1.

Remote sensing data
Remote sensing data for the agricultural regions of Saskatchewan ( Fig. 1) were acquired using Google Earth Engine (Gorelick et al. 2017) to be used as environmental Table 1. Distribution of parent material types in the detailed soil survey polygons of Saskatchewan (SKSIS Working Group 2018) and the National Pedon Database parent material class counts (Agriculture and Agri-Food Canada 2016 covariates in the parent material disaggregation models. Landsat 7 Tier 1 surface reflectance data were acquired from 1999 to 2021, with several band indices calculated ( Table 2). The median anthocyanin reflectance index (ARI), which is associated with anthocyanin pigments in plant foliage (Gitelson et al. 2001), was calculated for data from the months of July and August. This index was included as it has been useful for discriminating vegetation differences in Alberta for peatland mapping (DeLancey et al. 2019). The median canopy response salinity index (CRSI) was also calculated using data from the months of July and August, as it has been used to distinguish soil differences associated with salinity (Scudiero et al. 2015). The median normalized difference vegetation index (NDVI) was calculated for the months of July and October, separately. These months were selected as July corresponds to peak photosynthetic activity, and only non-arable cropland would have any photosynthetic activity in October. Non-arable cropland would be expected to be preferentially located on sandier parent materials or those with steeper slopes. The standard deviation of NDVI using data from May to October was also determined. Level 1 C data from Sentinel 2 were acquired for the months of May to October from 2015 to 2020 to generate the median red-edge inflection index (REIP). The REIP was included in the study, as it has been an important parameter for discriminating vegetation differences in Alberta (DeLancey et al. 2019). Each parameter was median focal filtered using a 3 × 3 window size and exported from Google Earth Engine at 30 m spatial resolution. The Google Earth Engine scripts for generating these covariates are available on GitHub (Sorenson 2021a). Data from the Sentinel 1 platform were also retrieved from Google Earth Engine. Processing of the Sentinel 1 data was done using the methodology in Hird et al. (2017). The median vertical-vertical (VV) and vertical-horizontal (VH) backscatter was determined using data from May to October from 2015 to 2020. To focus more on broader landscape trends, rather than fine-scale variation, the data were median focal filtered with a 9 × 9 window size and exported at a 30 m spatial resolution, rather than the default 10 m, to match the resolution of the Landsat 7 data. The Google Earth Engine script for generating these covariates is available on GitHub (Sorenson 2021a).
Bare soil composite imagery BSCI was generated using Google Earth Engine and the methodology from Sorenson et al. (2021). Atmospherically corrected Tier 1 surface reflectance data for Landsat 5 were retrieved for the months of July and August from 1984 to 1995. Landsat 5 data from this time period, rather than more recent Sentinel-2 or Landsat 7 and 8 data, were used because the total number of acres in Saskatchewan under conventional tillage peaked in 1970 and declined to only 18% of farmland by (Statistics Canada 2015. Therefore, more recent imagery is heavily affected by crop residues. The months were limited to July and August to ensure that crop residues from previous years had time to decompose, thereby reducing the influence of crop residues on BSCI. Clouds, shadows, and lowquality pixels were filtered using the Landsat quality assessment band (pixel_qa). The Landsat 5 data had seven bands: Band 1 corresponding to blue light (450-520 nm), Band 2 corresponding to green light (520-600 nm), Band 3 corresponding to red light (620-690 nm), Band 5 (1550-1750 nm), Band 6 (10 400-12 500 nm), and Band 7 (2080-2350 nm). Band 6 was not used in the analysis.
Image collections were filtered using NDVI, normalized burn ratio (NBR 2 ), normalized difference water index (NDWI) (Demattê et al. 2018), and normalized difference index 7 (NDI7) ). Pixels were kept as likely representing the bare soil surface if the NDVI value was less than 0.3 , NBR 2 was less than 0.1 , NDWI was less than 0.5 (Du et al. 2016), and NDI7 was less than 0 (Kokaly et al. 2017). Additionally, any region mapped as pasture or grassland in Agriculture and Agri-Food Canada's 2019 Annual Crop Inventory was masked (Agriculture and Agri-Food Canada 2019). While 2019 does not match the years of the bare soil imagery collections, both land uses tend to be managed for long term. Following filtering and masking, the median reflectance values for each pixel were calculated, and then spatially filtered using a circular 10 × 10 median focal filter. The focal filtering was done to account for the 300 m location uncertainty reported in the
Note: Band ratio equations for each band ratio were used as an environmental covariate in the analysis. The listed bands correspond to Landsat 5 and 7 except for the red-edge inflection point which correspond to Sentinel 2 bands.
National Pedon Database and to smooth out fine-scale variation more likely attributable to factors other than parent material. The final filtered composite image was then aggregated to a spatial resolution of 30 m. Additional details and band ratio equations are available in Sorenson et al. (2021) and the Google Earth Engine script is available on GitHub (Sorenson 2021b). Example results for the BSCI are provided in Fig. 2, which includes true colour red-green-blue images of the bare soil surface.

Terrain attributes
Terrain attributes were determined using the 30 m digital surface model from the Advanced Land Observation Satellite (ALOS) (Japan Aerospace Exploration Agency 2015) retrieved using Google Earth Engine. Two terrain attributes, terrain ruggedness index (TRI) and the standard deviation of elevation, were calculated using a range of median focal filtering window sizes and calculation window sizes (Table 2). Median focal filtering was used, rather than mean, as it is less sensitive to the effects of outliers. These attributes were selected, as they represent coarser scale patterns in landscape variability rather than specific fine-scale patterns related to factors such as local landscape morphometry, hydrological characteristics, and landscape exposure. While these finer-scale attributes were documented to improve predictive parent material mapping at 100 m in British Columbia Fig. 2. Imagery for the four finer spatial resolution example areas. The first row corresponds to the northeast example area, the second row is the northwest example area, the third row is the southwest example area, and the fourth row is the southeast example area. The figures in column (a) are true colour red-green-blue median surface reflectance image (Landsat 7 Bands 3-2-1) for Saskatchewan from July and August from 1999 to 2021. The figures in column (b) are true colour red-green-blue median reflectance bare soil pixel composite images (Landsat 5 Bands 3-2-1) for Saskatchewan from 1984 to 1995. This represents a true colour image of the exposed soil surface. The white areas in column (b) correspond to areas either without bare soil pixels present or permanent pasture areas that have been masked. Coordinates are in UTM Zone 13 N NAD83. (Heung et al. 2014;Bulmer et al. 2016), the characteristics of prairie landscapes make these attributes challenging to accurately calculate with currently available DEMs in Saskatchewan. Satellite-derived DEMs with spatial resolutions of 30 m, like the ALOS model used here, lack the horizontal and vertical accuracy to capture fine-scale variation in the relatively low-relief Canadian Prairies. Specifically, the ALOS 30 m DEM has a reported vertical root mean square error of 1.78 m (Caglar et al. 2018), which is greater than the relief differences driving soil property differences in some hummocky landscapes in Saskatchewan (Landi et al. 2004). Additionally, 30 m spatial resolution is not fine enough to consistently characterize slope and topographic variation in many Saskatchewan landscapes (Landi et al. 2004). However, currently these are the best freely available DEM data for Saskatchewan. Therefore, this study focused on using terrain attributes that accounted for general degrees of variability more broadly.
TRI was calculated using the system for automated geoscientific analyses (SAGA) (Conrad et al. 2015). The TRI was calculated with either a window size of 10 or 20 (Table 2), and with the following median focal filters applied to the DEM prior to processing: no filter, 3 × 3, 5 × 5, and 9 × 9 (Table 2). Applying a 3 × 3 median focal filter to the TRI results with no filter on the DEM was also tested. The standard deviation of elevation was calculated using Google Earth Engine with window sizes of 3 × 3, 5 × 5, 9 × 9, and 21 × 21 following 3 × 3 median focal filtering of the DEM ( Table 2). The standard deviation of elevation was also calculated using window sizes of 21 × 21 and 101 × 101 for a 9 × 9 median focal-filtered DEM. Examples of the terrain features determined as important environmental covariates during model development are illustrated in Fig. 3. TRI was not calculated using a 101 × 101 window because of the computational time required compared to standard deviation. Multiple scales were considered as multi-scale digital terrain analysis has been shown to improve predictive soil mapping results ).

Model development
Predictive models were developed using the detailed soil survey polygons as training data, which were mapped at a 1:100 000 scale. The disaggregation approach was similar to Møller et al. (2019) where a single set of area-proportionally sampled synthetic training points were generated and a random forest model was trained. This approach was selected for the computational efficiency (Møller et al. 2019). Synthetic training data were developed by randomly assigning points within each polygon on an area proportional basis. The number of points generated for each polygon was equal to the natural logarithm of the polygon area divided by the area of the smallest polygon in the data set. This number was then multiplied by 5 for the final number of points per polygon. The rationale for this approach was that it ensured multiple training points are generated for each polygon and that larger polygons have more points generated to reduce the likelihood that they are underrepresented in the training data.
Points were assigned to the parent material class based on the first deposition layer described for the polygon that they were generated in, with the exception of the lacustrine over till class, which reflects both the first and second deposition layers. Polygons with parent material values of undifferentiated indicated that there is too much variability in parent materials at the mapping scale of the detailed soil survey polygons to assign a parent material class, so data from these polygons were not used for model development. As each polygon has a single component, the single component was used to assign the training data label. Multi-component polygons in the parent material maps were labelled as undifferentiated and they were excluded during training data creation. In total, 5 786 772 synthetic training points were generated. The scripts for generating the training data from the parent material polygons are available on GitHub (Sorenson 2021c).
The parent material classes for the synthetic training points generated from the detailed soil survey polygons were simplified. Fluvial eolian and fluvial lacustrine parent materials were grouped as fluvial because of the similar characteristics of these materials. Additionally, minor parent materials with very similar characteristics were removed, specifically glaciolacustrine and glaciofluvial. Fluvial eolian and fluvial lacustrine were grouped as fluvial. Fen and sphagnum peat were grouped as a single peat class. Lacustro-till were relabelled as lacustrine due to its infrequent occurrence and high similarity to lacustrine deposits. Residual and undifferentiated bedrock were combined into a bedrock class. The final parent material classes were: bedrock, colluvium, eolian, fluvial, lacustrine, lacustrine over till, peat, and till.
To improve the computational efficiency of the model building, the total points were randomly subset to three sets with 1 000 000 training points from the original 5 786 772. Separate training subsamples were generated for the three following models: (1) the BSCI-only model, which includes only training points from locations where BSCI is available and includes BSCI as environmental covariates, (2) the exhaustive model, which includes training points from all locations, regardless if BSCI is available and includes BSCI as environmental covariates, and (3) the BSCI-excluded model, which also includes training points from all locations but does not include BSCI as environmental covariates. Where BSCI imagery was not available, a null value of 0 was set for model development purposes.
Following creation of the training data, random forest models were built using the ranger package in R (Wright and Ziegler 2017). As class imbalances exist in the training data, the case weights term in the ranger model was used to balance across classes during model building. The case weights term in the ranger package allows for probabilities to be assigned to each training point that determines the likelihood a point will be sampled during individual tree building in the random forest. As a result, each tree of the random forest can be built with balanced classes, and the overall random forest model will therefore be built on class balanced data even if the original data set has imbalanced classes. Case weights were set so that each parent material class had an equal chance of being subsampled for each individual tree Fig. 3. Terrain attribute imagery for the three terrain parameters selected as important predictors by the feature selection process. The first row provides imagery for the portion of Saskatchewan with conventional soil survey maps at 1:100 000. The bottom row is the imagery for the southwest example area. The figures in column (a) are the ALOS DEM with a 3 × 3 median filter applied and the standard deviation of elevation calculated with a 21 × 21 window. Figs. in column (b) are the ALOS DEM with a 9 × 9 median filter applied and the standard deviation of elevation calculated with a 21 × 21 window. The figures in column (c) are the ALOS DEM with a 9 × 9 median filter applied and the standard deviation of elevation calculated with a 101 × 101 window. Coordinates are in UTM Zone 13 N NAD83. built in the random forest. An initial model was built using all environmental covariates listed in Table 2 to determine feature importance based on Gini index. Features were then selected for the final models using a backward feature selection process, where random forest models were built using sequentially less features. For each iteration, the feature importances were recalculated, and the least important feature was removed. The final features selected were determined based on where the out-of-bag classification error was minimized. In total, the BSCI-only model used 16 features, the exhaustive model used 20 features, and the BSCI-excluded model used 11 features (Table 3). Final random forest models were built for each scenario with the following model hyperparameters: probability set to true, extra trees as the split rule, and case weights set so that subsampling probability was equalized across parent material classes to reduce the influence of class imbalances on the model. The scripts for generating the random forest models and resulting parent material maps are available on GitHub (Sorenson 2021c).

Validation data
Data from the National Pedon Database (Table 1) were used as model validation data (Agriculture and Agri-Food Canada 2016). Historically, the samples were collected to represent modal soils from dominant soil types from a region and were collected by horizon with analyses performed on the homogenized horizon samples. The location accuracy for these samples is low, with specific locations estimated relative to land marks by the soil surveyors. The estimated location uncertainty is approximately ±300 m. The 300 m uncertainty is the location uncertainty reported for the data in the National Pedon Database. These soil pedon observations were used to inform the soil survey parent material mapping and so they are not completely independent from the maps used as training data, but the point observations themselves were not used in the model training. The soil parent material classes were simplified for this study based on the original parent material classes in the survey database using the same approach as the training data. Standard deviation of 9 × 9 median filtered elevation with a 21 × 21 focal window 52 061 Median CRSI 50 442 Note: These models include only training data with BSCI present, using training data with and without BSCI, and without using BSCI. The distribution of parent material classes amongst pedon validation data did not match the distribution of parent material classes reflected by the detailed soil survey polygons (Table 1). This was likely because the pedon data were collected to provide example profiles for the range of soil associations present in Saskatchewan, not to sample parent material classes on a proportional basis. The most frequent parent material in the pedon validation data was fluvial, followed by lacustrine, till, and then lacustrine over till (Table 1). Minor amounts of eolian, bedrock, and colluvium parent materials were present. No peat data are present in the validation data. For the detailed soil survey polygon data, the most common parent material class was till, followed by fluvial and then lacustrine (Table 1).

Model validation
Model performance was evaluated based on overall model accuracy, Cohen's Kappa, which measures inter-rater reliability and considers agreement occurring by chance, and specific class producer accuracy. Accuracy and Kappa values in this study are reported on a scale of 0-1. Producer accuracy refers to the number of correctly predicted observations of a soil class per the total number of observations of that class (Malone et al. 2017), and is referred in the results simply as class accuracy.

Uncertainty estimates
Prediction confidence was determined using the confusion index (Burrough et al. 1997;Odgers et al. 2011;Brungard et al. 2015), which is calculated with the following equation: where CI is the confusion index, u max is the probability value of the most likely prediction based on the random forest proportion of votes in the ensemble, and u (max−1) is the probability value for the second most likely prediction based on the proportion votes in the ensemble. Confidence interval values closer to 1 indicate less certainty and values closer to 0 indicate more certainty in the prediction.

Results and discussion
The results of this study are presented in two sections. The first section is a review of the results for the BCSIonly model, the BCSI-excluded model, and the exhaustive model. The second section provides details on the model features

Model performance
Overall, the BSCI-only model was the best performing model. The accuracy for the most probable parent material was 0.52 and Kappa was 0.35 (Table 4). The model performance for the exhaustive model was very similar, but slightly lower, with an overall accuracy of 0.48 and Kappa of 0.30 (Table  5). The model performance for the BSCI-excluded model was lowest, with an accuracy of 0.38 and a Kappa value of 0.19 (Table 6). Predictive results were comparable between the two models that included BSCI environmental covariates (Fig. 4). However, it is important to note that the BSCI-only model could only generate predicted maps for those regions where BSCI imagery exists, which is a significant limitation (Fig. 5).
In comparison, classification accuracy for regional-scale parent material mapping in British Columbia had an overall accuracy value of 0.779 and Kappa values of 0.697 (Heung et al. 2014). That study focused on topographic derivatives that appear to work well in British Columbia. The greater relief in British Columbia improves results with coarser-scale DEMs with higher vertical error, compared to that in the Canadian Prairies where relief changes related to changes in soil properties can be shallower than the DEM vertical error (Landi et al. 2004). The better results in British Columbia compared to this Saskatchewan study might be explained by the less spatially extensive study area mapped in the British Columbia study (5472 km 2 in British Columbia versus 298 000 km 2 in Saskatchewan). A second British Columbia study targeting a more extensive study area (945 000 km 2 ) had comparable accuracies to this study with values of 0.411 or 0.415, depending on whether a balanced or constrained approach was used for predicting parent material (Bulmer et al. 2016) Cross-class confusion for the BSCI-only model showed variance in performance across the parent material classes. Due to the lack of validation points, the mapping performance for colluvium and peat classes cannot be evaluated as part of this study (Table 4). Eolian and bedrock classes had few validation points resulting in less certainty in the overall predictions for those classes. The bedrock class observations were predicted across the parent material classes. Predictive accuracy for the eolian parent material class was the highest with a class accuracy of 0.75 (Table 4). Fluvial and till parent materials had the majority of the validation points classified correctly, with the primary misclassification for fluvial being lacustrine (Table 4).
Interestingly, an area mapped as till in the southwest corner of the southeast example area was mapped as fluvial deposits rather than till even though it was mapped as till in the original detailed soil polygons (Fig. 4). Without validation data in this area, the accuracy of this mapping cannot be confirmed. However, some of this area could include tills reworked by fluvial processes as they are adjacent to a large glacial lake that was present in this region of Saskatchewan (SKSIS Working Group 2018). Till was primarily misclassified as lacustrine over till, which is the parent material with the most similar characteristics, and surface expression influenced by the underlying Till. Predictions for lacustrine parent material observations were less than 50% correct, with a class accuracy of 0.47. There was a fairly even distribution of misclassifications across fluvial, lacustrine over till, and till observations. Lacustrine over till observations were classified correctly 25% of the time, with the class being primarily misclassified as lacustrine. This is likely because of the importance of bands 5 and 7, which are related to clay content , to the model, and the similar surface textures between these two parent materials.
Cross-class confusion for the exhaustive model (Table 5) was similar to the BSCI-only model. Overall, this model had slightly lower class accuracies for the eolian, fluvial, and till parent material classes. There was a slight increase in class accuracies for the lacustrine and lacustrine over till classes. There were major changes in the performances per class for the BSCIexcluded model. The accuracies of all classes dropped (Table 6). This was particularly the case for the eolian class. This is likely because of the lack of bare soil imagery Bands 5 and 7 to identify the lower clay content areas associated with eolian parent materials. Based on these shifts in overall accuracy, Kappa, and individual class accuracies, the utility of the BSCI-excluded model is low.

Feature importance
Overall, the most important feature for the BSCI-only model was the Sentinel-1 VH backscatter (Table 3). A possible explanation for the importance of this feature is that VH backscatter has been shown to be related to aboveground biomass (Laurin et al. 2018). By comparison, the VV backscatter was less important (11th most important feature). VV is sensitive to soil moisture, but does vary with vegetation characteristics as well (Vreugdenhil et al. 2018). Particularly under Fig. 4. Parent material prediction results for each of the four example areas illustrated in Fig. 1. The first row corresponds to the northeast example area, the second row is the northwest example area, the third row is the southwest example area, and the fourth row is the southeast example area. The figures in column (a)  water-limiting conditions, finer textured soils have increased yields and associated biomass (He et al. 2013), and parent materials having higher average clay contents, and therefore biomass, are likely the reason Sentinel-1 VH backscatter was an important feature. Additionally, increased biomass would be associated with areas managed as pasture, which is also associated with coarser parent materials such as fluvial deposits.
The second most important feature was the standard deviation of 9 × 9 median focal filtered elevation with a 101 × 101 focal window. This terrain feature was the most smoothed and coarsest-scale feature of the terrain attributes (Fig. 3), suggesting that the coarser-scale terrain attributes were more useful for parent material mapping in Saskatchewan, at least given the ALOS DEM used in this study. The coarse smoothing resulted in a terrain feature associated with larger-scale landscape heterogeneity. This likely was helpful in differentiating between glacial till with heterogeneous hummocky landforms compared to level or undulating lacustrine landscapes that would have low amounts of terrain variability.
Optical remote sensing features were then the next most important features, with the median ARI being one of the most important of these features (Table 3). The ARI is a function of anthocyanin pigments in foliage (Gitelson et al. 2001) and is likely corresponding to variation in plant composition, which would be closely related to management practices. Coarser-textured hummocky parent materials, such as fluvial parent materials, are more likely to be managed as pasture compared to annually cropped lacustrine deposits. The standard deviation of NDVI was the most important NDVI feature (third-most important feature overall). Areas with higher NDVI standard deviations indicated more changes in NDVI over the growing season. This is expected to be influenced by factors such as management practices, which, in turn, are influenced by soil and parent material type. The standard deviation of NDVI has been shown to be a useful remote sensing feature for land use monitoring and land classification (Becker et al. 2021). The median October NDVI and Median July NDVI were less important as the eleventh and fourteenth most important features, respectively (Table 3).
While not the most important features, BSCI features were reasonably important: the fifth and eighth most important features were the BSCI bands 5 and 7 (Table 3), which correspond to shortwave infrared bands. The importance of these bands could be due to increased absorption in these wavelengths with higher clay contents considering they were found to be the most important bands for mapping clay content in Saskatchewan with BSCI ). Clay has absorption features within Landsat 5's shortwave bands range (Rossel et al. 2010). However, these bands cover a wide portion of the shortwave infrared portion of the electromagnetic spectrum that may reflect other soil properties such as soil organic carbon, so their relationship to clay content is not certain. The sixth, ninth, and tenth most important features were the visible light bands for the BSCI. These bands were important predictors for SOC prediction in Saskatchewan , along with the near-infrared band which was the 12th most important feature. Silt and clay content, which are affected by parent material type, has been documented to be positively related to soil organic carbon concentration in Saskatchewan (Plante et al. 2006).
Standard deviation of a 3 × 3 median-filtered elevation with a 21 × 21 focal window was the 13th most important feature, which would be reflective of finer-scale terrain patterns compared to the 9 × 9 filtered 101 × 101 standard deviation (Fig. 3). Including a feature that characterizes the terrain patterns at different scales likely helped distinguish different parent materials as the contrasts amongst some parent materials were greater at coarser scales (i.e., lacustrine versus lacustrine over till) and others were greater at finer scales (i.e., lacustrine versus fluvial). While inclusion of finer-scale terrain attributes did improve the results, the coarsest-scale terrain attribute was more important (Table  3), at least for mapping parent material in Saskatchewan with a 30 m DEM, and in conjunction with other remote sensing features. While TRI has value quantifying topographic heterogeneity in other prediction soil mapping studies (Brungard et al. 2015), the simple standard deviation of elevation was determined to be more important by the random forest models in this study for the BCSI-only mapping.
The most important features for the exhaustive model were very similar to the BSCI-only model (Table 3). There were, however, differences in terms of the relative importance of each feature. The importance of Sentinel-1 VV backscatter was much greater for the exhaustive model than for the BSCIonly model. The likely reason for this is that because this model included areas where bare soil composite data were not present, the model needed to rely on other feature to distinguish differences amongst parent material classes. There was also an additional terrain feature included in the final selected feature set, particularly different focal window sizes for the standard deviation of elevation and two TRI features. In contrast with the BSCI-only model, BSCI band 3 was not included in the final model. Overall, BSCI was less important and bands 5 and 7 were found to be the most important BSCI bands. The reduction in importance of the BSCI bands could be partially explained by the fact that a significant portion of the training data would not have variance in bare soil values.
The median REIP and CRSI were also important features for the exhaustive model. The REIP is an approximation of the near-infrared red inflection point in vegetation spectra. The location of REIP is highly correlated with foliar chlorophyll content, photosynthetic activity, and a good predictor of leaf area index in wheat (Herrmann et al. 2010). CRSI has been documented to be related to salinity in other contexts (Scudiero et al. 2015). Salinity is affected by many factors with drainage and hydraulic conductivity being two important factors (Eilers et al. 1997). Both factors vary within the different parent material classes in this study, which explains why they were less important for parent material mapping compared to other features.
For the BSCI-excluded model, Sentinel 1 backscatter, standard deviation of elevation, and NDVI attributes were most important, with the standard deviation of NDVI as most important NDVI feature followed by the median October NDVI (Table 3). As bare soil composite data were not present, the model found features related to land use such as the standard deviation of NDVI to be the best as separating parent material classes. Sentinel 1 VH backscatter was the most important feature followed by the standard deviation of 9 × 9 median filtered elevation with a 101 × 101 focal window, and then Sentinel 1 VV backscatter. Other important features were ARI, REIP, the standard deviation of 3 × 3 median filtered elevation with a 21 × 21 focal window, median July NDVI, and then the standard deviation of 9 × 9 median filtered elevation with a 21 × 21 focal window, followed by the CRSI. For all three models, the coarsest scale terrain attribute was the most important terrain attribute, as the greatest differences amongst parent material classes were apparent in the larger-scale terrain variability-associated features. Other studies have found local and landscape morphometric features along with hydrological characteristics useful (Heung et al. 2017). With a higher quality DEM, these features would also likely be useful in the Canadian Prairies as well; however, with the 30 m DEM available, the coarser scale properties were more useful in the context of this study. This could be because the scale of terrain variation in the Canadian prairies is such that a 30 m DEM cannot consistently detect all slope positions in the landscape, and the coarser scale terrain analysis has enough of the upper and lower slope positions in the DEM to better characterize overall variability if not specific slope positions at a given point.

Conclusion
Parent material disaggregation has an important role in predictive soil mapping efforts in the Canadian Prairies. The major drivers of soil variation at scales finer than the existing soil maps are parent material variance, slope position, and salinity. Disaggregating parent material maps to create finer resolution maps is an important step for creating more detailed soil maps in the Canadian Prairies to support a variety of end uses. Based on these results, the inclusion of BSCI is an important covariate for parent material disaggregation in the Canadian Prairies and is likely essential for generating useful maps when high-resolution DEMs are not available. The model that did not include BSCI were less accurate overall. Therefore, parent material mapping in the Canadian prairies through disaggregation methods may currently need to be limited to areas where BSCI is available to ensure output maps with acceptable accuracy. Future work to disaggregate soil classes based on slope position and salinity, and to combine those methods with parent material disaggregation is needed to generate detailed soil maps for the Canadian Prairies.