2.1..

Study area

The study area is comprised of five drained and cultivated peatlands in the southwest plain of Montreal, Quebec, covering approximately 90 km² (45.1°N to 45.3°N latitude and −73.3°W to −73.7°W longitude). Figure 2 shows the study area and the spatial distribution of sampling sites. Natural forests are still present but are likely affected by agricultural drainage. Most of the soils are used for horticulture, while a small proportion is used for producing bags of gardening soil. The five peatlands are part of three separate watersheds (Lamontagne et al. 2014). The elevation ranges between 45 and 70 m above mean sea level; however, the agricultural fields are leveled across the peatlands.

Fig. 2.

Study area and sampling sites categorized by the nature of the information collected. Overlaid is a square grid of 1.5 km size used to spatially partition the data sets to allow spatial cross-validation. Numbers beside the basin limits are peatlands identification numbers used in this study.

Fieldwork was carried out on 14 partnering research farms, which were converted to agricultural use between 1950 and 2010. Most fields were classified as Humisols and Mesisols, while a small area of recently converted fields was classified as Fibrisols (SCWG 1998). The organic materials found in the southwest plain of Montreal were deposited in channels or in depressions. Basin bogs or shore swamps, the latter having shallower deposits, were the result of hundreds of years of accumulation (Lamontagne et al. 2014). Peatlands in this region are mainly composed of forest peat or herbaceous materials (LaSalle 1963; Grenon 1988). Limnic materials are found throughout the study area but not in every peatland: the northwestern-most peatland appeared to be devoid of limnic material, whereas the other peatlands had regions with limnic layers.

2.2..

Soil data

Two data sets were used to build the predictive models. The first data set consisted of 255 sites that were sampled between 2019 and 2021 (shown as white dots in Fig. 2). The CLT and the DML were obtained by manually extracting soil cores using a Macauley corer (Eijkelkamp peat sampler) until the mineral horizon was reached.

The second data set consisted of 4488 observations. Here, 4286 locations were manually probed (shown as grey dots in Fig. 2) between 2010 and 2014, and were combined with the first data set. Manual probing involved inserting a thin metal rod in the soil and when a change of soil resistance to the rod penetration was felt, the depth was recorded as the DML. It should be noted that coprogenous materials did not offer notable resistance when a rod was inserted into the soil in comparison to the apparent change in density when the probe reached the underlying mineral soil. Manual probing was susceptible to measurement errors due to buried wood pieces (Parry et al. 2014). Nonetheless, manual probing could be done rapidly and required less energy to do manually than the Macauley sampling technique. Furthermore, the sites were clustered in three of the five peatlands. It was believed that the addition of this data set to this study could contribute to the predictive power of the models when generating DML predictions even though manual probing might have been more prone to measurement bias or error. Since some of the sampling locations were close to each other, observations that fell into the same raster cell were averaged. The summary statistics of the two data sets are presented in Table 1 and were also computed for both data sets for each individual peatland (Table A1).

Table 1.

Summary statistics of the coprogenous layer thickness (CLT) and of the depth to the mineral layer (DML).

2.3..

Environmental covariates

When generating a training data set to be fitted with a predictive model, the georeferenced soil observations are spatially intersected with a suite of geospatial environmental layers (i.e., covariates or predictors). The predictive models are used to establish the relationships between the covariate and soil properties to predict the spatial phenomenon of interest. Here, we present the five categories of covariates used, which were selected based on the SCORPAN model (i.e., soil, climate, organisms, relief, parent material, time, and spatial location) in McBratney et al. (2003) and a review of the digital mapping of peatlands in Minasny et al. (2019). The description and preprocessing details are provided for the final selection of covariates.

2.3.4..

Landsat 8 OLI_TIRS C2 Level-2 data

Landsat 8 is a satellite with a 16-day repeat cycle over the same area. C2 Level-2 data are referred to as “analysis ready” data—being preprocessed by the USGS Earth Resources Observation and Science (EROS) Center and then made freely available. After Level-2 processing, surface reflectance (bottom of atmosphere) values were obtained. More information on the Landsat missions and preprocessing can be found in Young et al. (2017).

This satellite provides information on the Earth’s temperature and land surface. The data acquired from two thermal infrared sensors (100 m spatial resolution) and nine spectral bands (30 m spatial resolution with one at 15 m resolution) are sensitive to different wavelength ranges. Similar to the Sentinel-2 data, the Landsat 8 raw bands and indices were resampled to a 10 m spatial resolution. A multitemporal approach was also used to capture seasonal and moisture gradients. Three series of satellite images were chosen based on cloud percentage, the absence of snow and crops. Images from 6 May 2015, 14 May 2018, and 30 May 2018 were used. The median and standard deviation over the three dates were used as covariates. The remote sensing indices were generated from both satellites, using their respective bands and were summarized in Table 2. All indices were computed in QGIS (QGIS.org, version 3.10).

Table 2.

Bands and indices from Landsat 8 and Sentinel-2 satellites explored as potential covariates.

2.4..

Variance inflation factor analysis

The 142 covariates produced were reduced using a stepwise variance inflation factor (VIF) procedure implemented with the vifstep function (Naimi et al. 2014) to reduce multicollinearity (O'Brien 2007). Since it required the use of standardized covariates, the covariates were scaled using the raster package. The VIF approach takes each of the standardized covariate and uses the remaining ones as independent predictors in a multiple linear regression via ordinary least squares. Then, the coefficient of determination of the covariate acting as the dependent variable (R_i²) is calculated. This process is repeated for every covariate. In eq. 1, the VIF is computed for every i^th predictor, where a high proportion of variance explained by a given combination of predictors will result in a higher VIF score. The model then evaluates if one or more covariates has a VIF score that exceeds a predetermined threshold value. In our case, a score of 5 was selected although 10 is also a common threshold (O'Brien 2007; James et al. 2014; Bian et al. 2020). The covariate with the highest VIF is removed, and the process is repeated until all remaining covariates are below the threshold.

(1)

2.5..

Modeling approaches

Three commonly used machine learning algorithms in DSM were compared for predicting CLT and DML (Minasny et al. 2019; Rudiyanto et al. 2018): Cubist, Random Forest (RF), and k-Nearest Neighbor (kNN). Although this section will briefly summarize each model, details concerning the learners are provided in Heung et al. (2016), where the models were reviewed and compared. All the modeling was carried out using the R statistical software (R Core Team 2020) and the caret package, which included all the tested models (Kuhn 2020). Final maps were generated in QGIS at a 10 m spatial resolution.

2.6..

Model tuning and validation

2.6.2..

Model performance

Two accuracy metrics were used to select the best hyperparameters based on accuracy and precision: root mean square of error (RMSE, eq. 2) and Lin’s concordance correlation coefficient (CCC, eq. 3):

(2)

where O is the observed value and P the corresponding predicted value and N is the total number of observations. CCC was calculated as follows:

(3)

where s_op is the covariance between predicted and observed values, s² are their corresponding variance, ō is the mean of observed values, and is the mean of predicted values.

CCC can vary between −1 and 1, while a value near 0 indicates a lack of concordance between two variables. This index is similar to R², because it measures agreement between two variables but with the notable difference that it corrects for systematic bias if the relationship departs from the 1:1 line (Lawrence and Lin 1989). CCC can never be higher than the absolute value of the Pearson correlation coefficient in the presence of a bias. The CCC and RMSE statistics were obtained using the onsoilsurvey package in R (Saurette 2021).

2.7..

Model simplification and variable importance analysis

Recursive feature elimination (RFE) was performed with the rfe function in the caret package (Kuhn 2020) to further reduce the number of covariates for training the final models. This model simplification procedure was needed to aid in interpreting soil–environmental relationships and improving generalizability. Briefly, RFE incorporates resampling to measure model performance with a reduced selection of covariates. It is a backward feature selection method (i.e., the model initially includes all predictors), whereby the predictors are iteratively removed to simplify the model, and in some cases, improve its performance. To fit the RFE models, the caretFuncs function was used with 5-fold cross-validation repeated five times. This allowed the computation of RMSE (see eq. 2) that was used as the metric to determine the best model. The most notable difference between VIF and RFE is that VIF is carried out independently of the response variable to address multicollinearity, whereas the RFE reduces covariates as a function of the model performance in predicting the response variable to remove irrelevant predictors. In other words, the VIF was performed once, but RFE was performed six times (i.e., two response variables with three machine learners). To assess the contribution of the remaining selection of covariates, relative variable importance plots were made for both final models using the VarImp function in the caret package.

2.8..

Prediction uncertainty

Many peat mapping studies do not include local estimates of the prediction uncertainty (Minasny et al. 2019). To produce a prediction uncertainty map, a bootstrapping approach was adapted from Malone et al. (2017). Here, 100 realizations were produced for both CLT and DML final models. Each raster was generated by randomly sampling with replacement the original data set with an equal number of observations. Afterward, the 0.05 and 0.95 quantiles for each cell were obtained to generate the lower and upper bound of the 90% confidence interval of the predicted property.

Since the MPT map was the result of the difference between the DML and CLT maps, propagation of error had to be considered to adequately capture uncertainty. As such, the mean, variance, and covariance of bootstrap predictions were calculated, assuming a normal distribution of the statistics for each cell. Equation 4 describes how the standard deviation value was calculated for each cell of the final map, using the bootstrap-produced maps:

(4)

(5)

Finally, to obtain the upper bound (95^th percentile) and lower bound (5^th percentile) of the 90% prediction interval, the corresponding z-value of 1.645 was multiplied to the standard deviation of the MPT. It was then subtracted and added to the predicted MPT value of a given cell (eq. 6).

(6)

3.1..

Variance inflation factor and recursive feature elimination analyses

The VIF procedure reduced the number of covariates from 142 to 59, while the RFE procedure further reduced the number of covariates for the DML and CLT models to numbers between 5 and 25 covariates. Table 3 lists the final selection of covariates for all models. Only two Landsat 8 covariates and four Sentinel-2 covariates were retained in the final selection, while a larger number of DEM derivatives were retained. The distance and gamma radiometric layers were also retained in the final selection. Both selected kNN models show the lowest number of covariates after the RFE procedure.potassium (K, %), equivalent uranium (URA, ppm), equivalent thorium (THO, ppm), and ratios URA/THO (RUT), URA/K (RUK), and THO/K (RTK).

Table 3.

Final selection of covariates for the six models after stepwise variance inflation factor and recursive feature selection. DEM = digital elevation model; kNN = k-Nearest Neighbor; RF = Random Forest.

3.2..

Leave-one-block-out cross-validation

A spherical model was used to fit the CLT data, while an exponential model was used to fit the DML data to evaluate spatial autocorrelation (Fig. 3). The range for the CLT was 4733 m and for the DML was 957 m. The range provides the maximum distance at which two measurements of the same property are related. The results indicated the presence of autocorrelation among both features that needed to be captured during model tuning. One should note that for the exponential model, it was the effective range (i.e., lag corresponding to 95 % of the sill) that was defined.

Fig. 3.

Experimental variograms (dots) of the two predicted features and their best fit model (line). (A) Variogram of the coprogenous layer thickness using 255 square-root-transformed observations. (B) Variogram of the depth to the mineral layer using 4488 observations. The lag was 100 m for both variograms.

To account for spatial autocorrelation during the cross-validation procedure, the ranges obtained above were used to generate spatial clusters of sites. The study area was divided in squares using a grid of 1000, 1500, and 2000 m for comparison. For each model, LOBOCV was performed. This approach was compared to a non-spatial 10-fold cross-validation. The latter returned overoptimistic values of concordance and RMSE compared to LOBOCV (data not shown), reinforcing the need of a spatial cross-validation approach. It was also determined that a distance of 1500 m produced better results according to the accuracy metrics from the cross-validation procedure and was used to form spatial blocks. The spatial blocks are shown in Fig. 2.

3.3..

Accuracy assessment and model selection

Internal validation results of the three machine learners are summarized in Table 4. Models were compared based on their CCC and RMSE obtained after LOBOCV for the CLT and after hold-out cross-validation for DML. It was clear that CLT models underperformed compared to DML models; CLT models had lower RMSE, but far lower CCC compared to DML models. The Cubist model had the best performance for CLT predictions, with a CCC = 0.07 and RMSE = 30 cm. For the DML predictions, model performance of the Cubist model (CCC = 0.43 and RMSE 48 cm) was similar to the RF model (CCC = 0.40 and RMSE = 44 cm). The difference between the observed and predicted range also led to the selection of the Cubist model as the final DML model. The range of predicted values by Cubist were closer than RFs to the actual range of observed values, although Cubist overpredicted slightly.

Table 4.

Tuning and leave-one-block-out cross validation results for the coprogenous layer thickness (CLT) and hold-out cross validation for the depth to the mineral layer (DML) models. Standard deviation of the concordance and of the root mean square of error (RMSE) are in parentheses. RMSE, observed, and predicted ranges are in centimeter. kNN = k-Nearest Neighbor; RF = Random Forest.

Furthermore, a plot of observed and predicted values was generated for both final models (Fig. 4). DML predictions seemed to follow the 1:1 line; however, a shift in the trend was observed for values above 200 cm (Fig. 4A). Concerning the CLT, the datapoints did not follow the 1:1 trendline. Moreover, many sites without a coprogenous layer (observed = 0 cm) had a predicted CLT up to 120 cm (Fig. 4B).

Fig. 4.

Plots of the correlation between observed and predicted values for the Cubist model (A) after hold-out cross-validation for the depth to the mineral layer (N = 1346), and (B) after leave-one-block-out cross-validation for the coprogenous layer thickness (N = 255). A perfect fit is represented by the red 1:1 line.

3.4..

CLT and the DML predictions and uncertainty estimates

The final CLT map and its 90% confidence interval bounds were produced (Fig. 5). According to Fig. 5B, thicker layers of coprogenous soil were predicted in the two lower peatlands and tended to be thicker towards the center of those peatlands. This trend was consistent with field observations. Nevertheless, no uniform gradient was found, indicating localized accumulation sites across the peatlands. Furthermore, while the peatland in the northwest was reportedly exempt of coprogenous material according to our field sampling survey, the Cubist model predicted a CLT between 0 and 55 cm. The prediction uncertainty maps derived from the bootstrapping technique showed a 90% confidence interval range of 11–304 cm.

Fig. 5.

Maps of the predicted coprogenous layer thickness (cm) across the study area using the Cubist model and 90% prediction limits derived using bootstrapping (N = 100). (A) Lower prediction limit (i.e., 5^th percentile), (B) prediction, (C) upper prediction limit (i.e., 95^th percentile), and (D) prediction interval width.

Concerning the predicted DML final map (Fig. 6B), the spatial pattern of the thickness of organic material deposits shared similarities with the final CLT map. Yet, high DML values did not always correspond to high CLT values. Furthermore, higher DML values were not exclusively found near the center of each peatland (i.e., the southwest and the center peatlands). Spatial artifacts can be seen, mainly around forested areas, and are related to the covariates used in the models. The southeast peatland was the shallowest on average, with a higher concentration of low DML predictions. Compared to CLT predictions, the DML 90% confidence interval width was narrower (i.e., 2–190 cm).

Fig. 6.

Maps of the predicted depth to the mineral layer (cm) across the study area using the Cubist model and 90% prediction limits derived using bootstrapping (N = 100). (A) Lower prediction limit (i.e., 5^th percentile), (B) prediction, (C) upper prediction limit (i.e., 95^th percentile), and (D) prediction interval width.

3.5..

Variable importance analysis

The relative variable importance plots for both final Cubist models are shown in Fig. 7. These plots show covariate importance in all condition and (or) linear model included in the Cubist tree. They provide information concerning the number of times that a covariate was used in the final model's tree (Kuhn 2020). The contribution of each variable was evaluated with model specific metrics since the final models were both Cubist. For the CLT, MRRTF (100 %) contributed the most to the final model, followed by DIST_MID (86 %), MRVBF (81 %), MSP (69%), and MaxElevDevScale2187 (64 %) for the top 5 covariates. Following that, a gradual decline in variable importance can be seen. Dist_MID (100 %) was the most important covariate in the final DML model, followed by THO (86 %), MRRTF (65 %), DIST_Y (62 %), and Center (49%) for the top 5 covariates. A sharp decline in relative importance was observed for the other covariates.

Fig. 7.

Relative importance of each covariate (%) on the final Cubist model (A) of the depth to the mineral layer and (B) of the coprogenous layer thickness. Refer to Table 3 for a description of covariates’ abbreviations.

3.6..

Predictions of maximum peat thickness and uncertainty estimates

The MPT map at a regional scale was derived by subtracting the CLT map from the DML map to reflect the real thickness of the potentially arable peaty layer (Fig. 8B). MPT ranged from −79 to 367 cm and approximately 0.27% (or 32.05 ha) of the MPT map are cells with a negative MPT. The lower MPT prediction limit ranged between −167 and 271 cm (Fig. 8A), while the upper prediction limit ranged between 9 and 464 cm (Fig. 8C). With the bootstrapping technique and the propagation of error, uncertainty tended to be relatively high and variable across the study area (Fig. 8D). Indeed, the 90% prediction interval varied between 130 and 231 cm.

Fig. 8.

Maximum peat thickness map obtained by subtracting the coprogenous layer thickness to the depth to the mineral layer and 90% prediction limits derived using bootstrapping (N = 100). (A) Lower prediction limit (i.e., 5^th percentile), (B) prediction, (C) upper prediction limit (i.e., 95^th percentile), and (D) prediction interval width.

4.1..

Interpreting model performance

The use of regional, open access covariates has proven useful for predicting DML; however, an underperformance of the CLT models was observed in this study. This was likely due to size and spatial distribution differences between the data sets for the two predicted soil properties. Table A1 further highlights differences found across each peatland with regards to the descriptive statistics of both predicted features. Here, the number of observations, the range, and the coefficient of variation of each peatland may have affected the model performance. For the sake of simplicity, a single model was trained for the entire study area, although the CLT and DML spatial gradients differed between peatlands. The DML model was able to moderately generalize across all peatlands (CCC = 0.43 and RMSE = 48 cm), which was not the case for the CLT model (CCC = 0.07 and RMSE = 30 cm). The high coefficient of variation of the CLT (109%, Table 1) indicated high heterogeneity in the data set, compared to that of the DML (52%). This can be explained by the fact that CLT could be considered as a zero-inflated variable with a wide range of values. A large proportion of the data set was composed of sites without a coprogenous layer, as seen in Fig. 4b. Furthermore, data clustering and poor spatial coverage of the study area's extent could have limited the model performance in areas with a lower sampling density, especially when considering variable importance which showed that spatial distance-based predictors were relatively important.

Nonetheless, our prediction errors were comparable to other studies. Rudiyanto et al. (2018) tested 14 machine-learning models to map the peat thickness in Indonesian peatlands, where DML was predicted with a RMSE ranging from 1.8 to 2.8 m at a regional scale (50000 ha), while Gatis et al. (2019) obtained a RMSE of 0.31 m for a 40600 ha study area. These RMSE values were of the same magnitude as the RMSE values observed in this study (Table 4). Moreover, Rudiyanto et al. (2018) had better results with Cubist and RF models compared to the other tested models. Considering the number of observations (1779 and 159) and the peat depth range (0–7 and 0–12 m) in Gatis et al. (2019) and Rudiyanto et al. (2018), respectively, our CLT and DML models showed similar model performances when treated separately. Since both studies reported R² and we used CCC, we cannot make direct comparisons, but both studies clearly outperformed ours. This may be due to the use of spatial cross-validation instead of standard cross-validation. As recently described in Wadoux et al. (2021), spatial cross-validation of autocorrelated data are not always the best solution to obtaining a representative assessment of the bias of a model. Standard cross-validation might yield better results without being overoptimistic. This could partially explain why the CCC obtained with LOBOCV for the CLT model was as low as 0.07, but standard cross-validation yielded a CCC = 0.65, and RMSE = 51.6 cm. The DML model had a higher CCC = 0.95 and a lower RMSE = 21.6 cm using standard cross-validation than when using LOBOCV. The latter cross-validation technique was preferred due to its more conservative estimates of model performance and to account for spatial autocorrelation of the data.

4.2..

Variable importance analysis

The most important DML predictors aligned with other studies, while the interpretation of CLT predictors was limited by the model performance. Results from Rudiyanto et al. (2018) showed that MRVBF was ranked as the fourth most important covariate among those tested. This supports the importance of this DEM-derived covariate on the prediction of deposited materials. Some soil-forming factors affect limnic and organic deposits differently depending on the scale at which phenomenon are studied (Behrens et al. 2018). As stated by Gallant and Dowling (2003), MRRTF and MRVBF can be important predictors for hydrologic and geomorphic processes that are related to valley, depressions, and slopes, like sedimentary soil deposits (i.e., coprogenous material and subsequently peat deposits). Multiresolution covariates can provide this information to the model to allow a better understanding of the magnitude of the deposits. This was likely why so many topographic and scale-related covariates were retained in the final models of this study.

The Euclidean distance to the center of the study area (DIST_MID), the Euclidean distance to the north of the study area (DIST_Y) and the distance to the center of a given peatland (Center) were also important predictors, meaning higher deposits are generally found near the center and thinner deposits near the peatland border, but also that there is a general gradient from a peatland to another. Furthermore, we expected that gamma radiometric covariates would contribute to delineating the peat extent and predicting peat thickness (Minasny et al. 2019). Indeed, the most important gamma-radiometric covariate was ranked as the second most important variable (Fig. 7A) in the DML model. The same authors also suggested that multi-temporal satellite covariate might provide useful information on peat thickness. In this study, most of the Landsat 8 and Sentinel-2 predictors were highly correlated, possibly due to a short period from when the imagery was acquired. Moreover, many Landsat and Sentinel covariates were not retained from RFE, possibly due to their low effectiveness as predictors. While the multi-temporal median of the soil moisture index from Sentinel-2 had a variable importance of 35%, little can be interpreted from this due to the low accuracy of the CLT model. A weak relationship between CLT and the suite of covariates could suggest that the model failed to capture the variability of CLT at a regional scale. Other sampling surveys should be made at a smaller spatial extent (i.e., individual farm) and the potential benefits of proximal soil sensing tools for producing more relevant covariates should be evaluated.

The study area encompassed four pedological surveys for which digital maps were freely available online. Legacy soil data were considered at first due to their popular use in DSM but were not included as covariates due to temporal incoherence between the surveys (1950, 2000, 2001, and 2014). The evolution of organic soils over time can lead to substantial biases and render maps less useful, even though they can be updated with modern techniques (Kempen et al. 2009). Otherwise, this covariate could have been useful to predict areas with coprogenous soil and to delineate peat extent.

4.3..

Predictions of maximum peat thickness and uncertainty estimates

The results suggested a potential limited use of the MPT map produced at a regional scale for soil management purposes at the field-scale given the high uncertainty of the combined CLT and DML predictions. Artifacts in the DML and CLT maps created sharp transitions between neighbouring cells that were not observed in nature and may require further field verifications. These issues may affect interpretation of the final MPT map for fields near the artifacts and the propagation of error may occur when combining maps with such artifacts. These errors may limit the use of such mapping products for precision agriculture and conservation projects. Moreover, the accuracy was not sufficient to predict at which depth tile drainage could be installed or if soil conservation MPT thresholds for management are reached (for details, see Deragon et al. 2022).

Negative values were observed in a small portion of the final map. This was due to a predicted coprogenous layer being thicker than the predicted depth to the mineral layer. As previously stated, higher imprecision in the CML model may have been linked to a poor coverage of the study area; furthermore, a weak relationship between covariates and the predicted feature could be responsible for model uncertainty. Peatlands differed quite significantly on their average MPT. For instance, the southeast peatland showed the lowest average MPT compared to the other ones. This has major implications from a soil conservation perspective, because degraded and shallow soils are not uniformly distributed nor confined to a peatlands’ border. Therefore, farmers from different peatlands are not affected equally by the shallowness of their soil and will have to use different soil conservation measures accordingly to their situation.

Many sources of uncertainty were present in this study (Heuvelink 2017). Two methods of sampling with a different precision were combined. In addition, the data set from 2010 was mainly clustered in two peatlands. Since model tuning inevitably fits the model on a majority of points coming from those peatlands, the final model might not be as effective in predicting the features found in the other three peatlands if their formation and evolution differs significantly (i.e., botanical origin, soil management, groundwater level variations based on differences related to the watershed and elevation, average number of years since conversion to agriculture, etc.). Preprocessing of the original data to produce the covariates may have introduced error in the model as well (i.e., resampling, projecting, and smoothing raster layers). As stated in Samuel-Rosa et al. (2015), the number and the location of calibration points matter; furthermore, the covariates are only approximations of the real-world soil-forming conditions and thereby inherently prone to errors. These errors are complicated if not impossible to quantify and are propagated at each new step of the workflow. A better sampling design would have been needed from the start to cover the full feature space of DML and CLT of the study area, and a better spatial distribution of observation points. Unfortunately, little was known about the DML and CLT in three of the peatlands prior to the start of the study. Therefore, the MPT map should be used with caution, as one decision tool among others to manage organic soil conservation. Yet, not considering coprogenous materials would greatly overestimate the remaining soil resource that can be used to produce crops. Moreover, despite important errors in depth estimates, digital mapping results were consistent with the observed distribution of the DML and CLT.