Mapping the maximum peat thickness of cultivated organic soils in the southwest plain of Montreal

Abstract Large organic deposits in the southwestern plain of Montreal have been converted to agricultural land for vegetable production. In addition to the variable depth of the organic deposits, these soils commonly have an impermeable coprogenous layer between the peat and the underlying mineral substratum. Estimations of the depth and thickness of these materials are critical for soil management. Therefore, five drained and cultivated peatlands were studied to estimate their maximum peat thickness (MPT)—a potential key soil property that can help identify management zones for their conservation. MPT can be defined as the depth to the mineral layer (DML) minus the coprogenous layer thickness (CLT). The objective of this study was to estimate DML, CLT, and MPT at a regional scale using environmental covariates derived from remote sensing. Three machine-learning models (Cubist, Random Forest, and k-Nearest Neighbor) were compared to produce maps of DML and CLT, which were combined to generate MPT at a spatial resolution of 10 m. The Cubist model performed the best for predicting both features of interest, yielding Lin’s concordance correlation coefficients of 0.43 and 0.07 for DML and CLT, respectively, using a spatial cross-validation procedure. Interpretation of the drivers of CLT was limited by the poor predictive power of the final model. More precise data on MPT are needed to support soil conservation practices, and more CLT field observations are required to obtain a higher prediction accuracy. Nonetheless, digital soil mapping using open-access geospatial data shows promise for understanding and managing cultivated peatlands.


Introduction
Canada accounts for about a quarter of the world's peatland extent (Vepraskas and Craft 2015). Cultivated organic soils cover only 4% of the Province of Quebec's southern region; yet, they greatly contribute to food production and to the economy of the province with exports to northeastern USA (Groupe AGÉCO 2007;Parent and Gagné 2010). These soils are prized for their vegetable production, but they are affected by intense soil loss processes that are unique to the evolution of peat materials. Soil loss occurs primarily by subsidence, oxidation, and erosion after land conversion for agriculture, and is enhanced by the drainage of peatlands (Kroetsch et al. 2011;Vepraskas and Craft 2015). Over the past decades, with an annual estimated soil loss of 2.5 cm (Ilnicki 2003;Esselami et al. 2014), degraded, shallower fields have transitioned to less productive mineral soils, as if the mineral boundary had been moving toward the surface.
Furthermore, coprogenous material can also be found between the peat and mineral layers (Lamontagne et al. 2014) -also reducing the thickness of the cultivable peaty layer. When the material is found within the first 160 cm and has a minimum thickness of 5 cm, it is referred to as a "limnic" layer (SCWG 1998). This gelatinous, impervious material is unsuitable for agricultural production and specific to lacustrine organic deposits (Kroetsch et al. 2011). When it dries, this material shrinks and does not rewet (SCWG 1998). In Poland, calcareous limnic deposits have been shown to limit root growth (Ilnicki 2003). Limnic materials can be of different types: coprogenous earth (sedimentary peat), diatomaceous earth, or marl (SCWG 1998). In the southwestern plain of Montreal, coprogenous earth is most common and sometimes found with a small layer of marl. The latter is effervescent due to the presence of shells and precipitated CaCO 3 , while the former can be a mineral (<17% organic carbon), or organic deposit (≥17% organic carbon) enriched by algae or aquatic life plants transformation products (SCWG 1998). Since pedological surveys have typically focused on the spatial extent of peat and descriptions of soil series, little is known about the depth to the mineral layer (DML) and the coprogenous layer thickness (CLT). To mitigate the impacts of long-term soil loss, soil conservation approaches, such as the addition of biomass crop amendments (Dessureault-Rompre´et al. 2020), the plantation of tree windbreaks, or water table management, are needed; however, these approaches have a related cost and material requirements to be considered. Therefore, priority management zones need to be defined to guide the application of regional soil conservation plans. The effective peat thickness that can be used for agriculture, hereinafter referred to as the "maximum peat thickness" (MPT), could be used to define these priority management zones. In other words, the MPT can be defined as the thickness of the peaty layer of an organic soil, therefore, excluding coprogenous and mineral materials. This definition would better reflect the real long-term agricultural potential of a field than only mapping the DML as other studies have done. It is necessary to understand the spatial distribution of the MPT to better manage shallower soils.
While pedological surveys and taxonomic map products are available at a global scale, new mapping efforts are targeting soil properties and soil functions (FAO 2020). Modern digital soil mapping (DSM) techniques can provide such information by leveraging technological advances for predicting a soil attribute or class using georeferenced field measurements and a suite of environmental covariates obtained via remote or proximal sensing (McBratney et al. 2003). In peatland-related DSM studies, the thickness of the peat and organic carbon content are the most frequently predicted features (Minasny et al. 2019). To the best of our knowledge, predictive mapping of coprogenous materials has never been explored and the relevant covariates are unknown. Yet, the information that could be acquired through a predictive map on a regional scale is crucial to estimating the MPT, and, therefore, the delineation of priority management zones. Manual probing over large areas is labor intensive and relatively slow (Parry et al. 2014). A regional approach, relying on a field calibration data set, could provide useful maps when combined with a relevant set of covariates. Many combinations of covariates are commonly used in regression and classification mapping of peatlands and could be investigated as potential candidates for predicting the thickness of the coprogenous layer. Peat thickness and its extent are often evaluated by combining a digital elevation model (DEM) and its derivatives, airborne gamma radiometric data, electromagnetic data, and satellite data (Rudiyanto et al. 2018;Gatis et al. 2019;Minasny et al. 2019;Siemon et al. 2020). Most of the covariates are the product of remote sensing techniques, while ground penetrating radar, gamma radiometric data, and soil electrical resistivity or conductivity can be obtained with proximal sensing techniques (Rosa et al. 2009;Parry et al. 2014;Comas et al. 2015;Beucher et al. 2020).
Hence, given the lack of data concerning the DML and CLT for the study area and the abundance of covariates, the main objective of this study was to integrate open-access, remote sensing covariates and field data to predict the spatial distribution of the MPT, including mineral and coprogenous materials' depth. Here, the specific objectives of the study were (i) to determine the covariates that most contribute to the prediction of DML and CLT and (ii) to derive a map of MPT from the mineral and coprogenous material predictions as a tool to guide soil conservation practices. Figure 1 shows the methodological framework for this study.

Study area
The study area is comprised of five drained and cultivated peatlands in the southwest plain of Montreal, Quebec, covering approximately 90 km 2 (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.
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.

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 Methodological framework used in this study. Two prediction maps were produced, (1) for the depth to the mineral layer (DML) and (2) for the coprogenous layer thickness (CLT). Then, they were subtracted to obtain the maximum peat thickness map (MPT). Note that the number of covariates and observations in each model differ. VIF = Variance inflation factor, RFE = recursive feature elimination, LOBOCV = leave-one-block-out cross-validation, and HOCV = hold-out cross-validation. Bottom left is a schematic of an organic soil profile. Oh, m, f = Humic, mesic, or fibric peaty layer, Oco = Coprogenous earth layer, and Cg = Mineral soil layer, according to the Canadian System of Soil Classification (SCWG 1998). 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).

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.

Digital elevation model (DEM) data
A DEM is a common data set that can be used to derive a suite of geomorphological and hydrological covariates (McBratney et al. 2003;Minasny et al. 2019). Since organic soils are formed in lowlands, depressions, and water saturated conditions, the topographical covariates are expected to be effective predictors of peat accumulation and the delineation of peatlands. The DEMs were produced and distributed by the Ministère des Forêts, de la Faune et des Parcs, where it was downloaded from the Données Québec repository. Here, the multiple DEM tiles were mosaicked. The DEM data were derived from a LiDAR survey and made available at a 1 m spatial resolution.
The DEM was preprocessed in five steps using the whitebox package in R (Wu 2020). Missing data were first filled-in; following this, a feature preserving smoothing tool was used to remove short-scale noise due to the use of the ultrafine 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. resolution DEM data. Afterward, a mean filter with a 5 m × 5 m window was used to further smooth the DEM and remove artifacts. To reduce the computational demands and memory requirements, the DEM was aggregated to a 10 m resolution. This covariate became the reference to which all other covariates were resampled to match its extent and resolution.
Single cell pits were filled to remove local artifacts and to produce a more continuous output. The Breach Depressions Least Cost tool was used to prepare a DEM for deriving the hydrological covariates (Lindsay 2016). This tool breaches and then fills depressions in the DEM to ensure continuous flow paths through depressions using a least cost pathway method based on a breaching algorithm by Lindsay (2016). After preprocessing, both the RSAGA (Brenning et al. 2018) and whitebox packages were used to generate 70 covariates. Some of these covariates included the difference and deviation from mean elevation, topographic wetness index, multiresolution valley bottom flatness (MRVBF), multiresolution ridge top flatness (MRRTF), catchment areas, hillshade, slope, aspect, and curvature (Beven and Kirkby 1979;Gallant and Dowling 2003;McBratney et al. 2003;Lindsay et al. 2015).

Aerial gamma ray spectroscopy data
Six raster layers were downloaded from Natural Resources Canada's GEOSCAN database (Natural Resources Canada 2019). The layers represented surface concentrations of potassium (K, %), equivalent uranium (URA, ppm), equivalent thorium (THO, ppm), and ratios URA/THO (RUT), URA/K (RUK), and THO/K (RTK). Radiometric data have proven to be useful for differentiating mineral soils from organic soils given their difference in parent material, porosity, and water content (Beamish 2013; Keaney et al. 2013). Water-filled, peat soil attenuates bedrock geology radioactivity and can provide an indication of the magnitude of the deposit. The six raster layers were resampled from a 250 m to 10 m spatial resolution using a bilinear method to match the extent and resolution of the other covariates.

Sentinel-2 L2-A data
Sentinel-2 is a constellation of two satellites with a combined repeat cycle of 5 days over the same area. Common indices and raw bands were tested as potential covariates. Gholizadeh et al. (2018) obtained significant correlations (r = −0.74 to −0.36) between bands 4, 5, 11, and 12 and soil organic carbon; hence, raw bands could serve as discriminant covariates to delineate a peatland's extent. L2-A level of preprocessing generates a bottom-of-atmosphere corrected reflectance raster for all 13 bands. This level of preprocessing includes radiometric, geometric (including orthorectification), and atmospheric corrections. The raster layers were preprocessed by the European Space Agency and made available via the Copernicus Open Access Hub. All final raster layers were resampled to a 10 m spatial resolution using the bilinear method since the raw bands had a native resolution of 10, 20, or 60 m. Lastly, the raster layers were reprojected to match the other covariates.
A multitemporal approach was adopted to capture moisture gradients over the study area (Fathololoumi et al. 2021). Therefore, median and standard deviation layers were produced for all raw bands and indices that were gathered for the following dates: 31 March 2020, 20 April 2020, 25 April 2020, and 20 May 2020. These dates were selected because they had neither snow nor crops covering the fields and, hence, were more effective in differentiating between the organic and mineral soils, as well as forests. The dates were also selected based on data availability and cloud cover over the study area.

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).

Distance and directional covariates
The center of a peatland is often the point where peat accumulation is at its maximum, whereas near the border of the peatland, one can expect thinner deposits. Therefore, project-specific distance layers made in QGIS with the Multi-Distance Buffer plugin were included (Tveite 2018). The first layer was produced by locating a peatland's centroid, then generating 10 m buffers around it until the peatland's border was reached. The product was then rasterized to a 10 m spatial resolution. The second layer was produced by generating 10 m inner buffers from the peatland's border, which was then rasterized. One might think the results would be the same; however, the resulting covariates exhibit significantly different patterns due to the shape of each peatland. Aspect reflects the topographic shape and directionality but is a circular measure and cannot be used directly as a covariate; therefore, we generated raster layers that reflected the eastness and northness of the study area. In addition to these covariates, Euclidean distance fields were generated for the study area to provide spatial context (Behrens et al. 2018).

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 2 ) 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.  Nellis and Briggs 1992 Raw bands 2-8, 8A, 9, 11, and 12 x Note: L = 0.5 (medium density of vegetation) was used in the Soil Adjusted Vegetation Index equation.

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.

Cubist
Cubist is a rule-based model that produces a regression tree based on the M5 model (Quinlan 1992), and was subsequently adapted by Kuhn and Johnson (2013) in R. At each node of the tree, the data set is split based on a set of rules using the value of one or many covariates, and form groups that minimize the within-node variability. The terminal nodes (i.e., leaves) of the resulting tree consist of multivariate linear models, which are applied to the covariates to make predictions. The hybridization of piecewise linear models and the hierarchical tree models allow Cubist to capture both linear and nonlinear relationships between predictors and the response variables (Malone et al. 2017). This model has two hyperparameters to be tuned when training the model (Kuhn and Johnson 2013). Overpredictions and underpredictions can be accounted for with a boosting technique using the committees hyperparameter. It specifies the number of similar trees to be sequentially produced and aggregated to optimize the set of rules. This process is explained in greater details in Kuhn and Johnson (2013). The second hyperparameter is the number of neighbors to be considered in a nearest-neighbor search through the training data set to find the closest observation to the predicted value (Quinlan 1993). This final step is useful for adjusting the value at the predicted location based on the average of its nearest neighbors.

Random Forest
Based on ensemble theory (Zhang and Ma 2012), RF is a nonparametric, decision-tree model, where a set of uncorrelated trees are combined (Breiman 2001). Multiple trees are generated with the use of bagging, where each tree is built on a random bootstrap sample of the original data set (with replacement). Using the bootstrap sample, a series of node-splitting rules are generated with respect to the covariates with the objective of maximizing the within-node homogeneity and the between node heterogeneity. In the caret package, m try is the main tuning parameter and is used to control the number of randomly selected covariate at each node, where the values of m try range from 1 to m number of predictors. Higher values of m try results in a higher likelihood of correlation between the trees; however, it lowers the prediction variance at the same time (Hastie et al. 2009). The trees from the forest are then aggregated to obtain an average prediction for each new observation (Genuer and Poggi 2020). Predictions can be regarded as more effective than those from single tree-learners due to RF's ability to mitigate bias (Kuhn and Johnson 2013).

K-Nearest Neighbors
kNN is a distance-based learner with one hyperparameter, k, which represents the number of neighbors of the unobserved location to be used for the prediction (Hastie et al. 2009). The rationale behind this learner is that observations that have similar properties (i.e., covariates values) will tend to have a similar value for the response variable of interest. Therefore, kNN predictions are made using k observations from training data that have similar values to those of the predicted site in the covariates multivariate feature space. To find neighbors of the predicted value, a distance metric must be used (i.e., Euclidean distance) to assess the closeness between observations. In a regression context, the k closest training data will be averaged, but if k = 1, the predicted value will be assigned the value of the nearest training observation. Since the covariates had different value ranges, they were standardized to ensure nonbiased calculation of the distance metrics (Hastie et al. 2009).
2.6. Model tuning and validation 2.6.1. Spatial cross-validation As stated in a review of peat mapping studies (Minasny et al. 2019), 64 % of the 90 reviewed studies did not perform validation, which may result in overfitting and overoptimistic models' performance metrics. Therefore, each of the machine learners were trained to find the optimal hyperparameters for DML and CLT models. Considering recent papers that had overoptimistic performances of machine learners in DSM applications, spatial cross-validation was used to evaluate these models (Pohjankukka et al. 2017;Meyer et al. 2018;Schratz et al. 2019;Ploton et al. 2020). This approach has proven to be more effective at measuring internal error and controlling for spatial autocorrelation in the data set. When the data are clustered, the model is more likely to correctly predict nearby observations during the training and validation process and, therefore, estimates of model performance can be biased.
Both data sets had clustered observations over the study area. To allow a fair assessment of model performance, leaveone-block-out cross-validation (LOBOCV) was used (Roberts et al. 2017). This form of spatial cross-validation is similar to leave-one-out cross-validation, where one observation from a data set is partitioned out during model training and is used as a test observation. During LOBOCV, one spatial block is put aside as a validation set. The model is trained using all observations except for those from the testing block. The model is then applied to this unseen data and error measures are derived. The process is repeated by changing each validation block for a new one at each fold. The number of folds is equal to the number of blocks, respectively, 45 (CLT) and 44 (DML) for this study. To compare LOBOCV with the conventional k-fold cross-validation approach, the model performance of both final models at the testing stage was also evaluated with a 10-fold cross-validation with 10 repeats.
The CLT data set was used without modification in LOBOCV because of its low number of observations (N = 255), whereas the higher number of observations in the DML data set (N = 4488) made partitioning possible. The traditional method to split a data set for training and testing purposes involves a 70%-30% split, respectively, with the use of random sampling. Since the DML data were highly clustered, random sampling could not divide observations evenly in training and testing data sets for peatlands with less data. Thus, conditioned Latin hypercube sampling on the X and Y coordinates (Minasny and McBratney 2006) was preferred over random sampling. Before merging the two data sets (255 + 4286), they were split based on their spatial coordinates to ensure a similar spatial coverage of the study area between the training and test partitions. The LOBOCV was performed to tune the model with 70% of the data, and to compute internal performance assessment, the test data set (i.e., the remaining 30% yet unseen by the model) was used in a hold-out crossvalidation procedure.
To identify the suitable size for each spatial block unit, the range from the theoretical variograms for the DML and CLT were calculated using the approach described in Oliver and Webster (2014). The experimental and theoretical variograms were computed using the gstat package (Pebesma 2004;Gräler et al. 2016). Weighted least residual sum of squares and visual assessment of the variogram fit were used to select the best model (Oliver and Webster 2014). If all models had an adequate fit via visual assessment, the one with the lowest weighted sum of squared errors was selected. Four outliers were identified in the DML data set; however, they corresponded to deep peat deposits and were kept in the model. Four outliers were identified in the CLT and the data were skewed (Table 1). Since outliers were not mistakes nor did they belong to another population, they were kept in the model and a square root transformation was applied (Oliver and Webster 2014). This transformation was more effective at reducing skewness than a logarithmic transformation. No outliers remained after transformation. The grid was then applied to the area and can be seen in Fig. 2. The same grid was applied to both features prediction under the assumption that CLT and DML would have shared similar autocorrelation if CLT data set had more observations.

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): 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: where s op is the covariance between predicted and observed values, s 2 are their corresponding variance,ō is the mean of observed values, andp 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 2 , 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).

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.

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: where σ MPT is the standard deviation of the MPT for a single cell, σ 2 M is the variance of the mineral predictions, and σ 2 C is the variance of the coprogenous predictions for 100 bootstraps. It is assumed that both variables were dependent (i.e., correlated); therefore, a term corresponding to the covariance between the features was added to the generic equation (Ku 1966). Based on Malone et al. (2017), the mean square error estimate from the validation data was added to the bootstrapping variance σ 2 M and σ 2 C to account for systematic and random errors in the models. Here, σ MC is the covariance, as defined by eq. 5: where M i and C i are the values of the i th cell of one bootstrap raster, andM andC are the mean values of the 100 rasters for that cell, respectively, for mineral and coprogenous maps.
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).
where X MPT is the mean MPT of a given cell obtained from the bootstrapped maps.

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).

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.
To account for spatial autocorrelation during the crossvalidation 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 crossvalidation. 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.

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.
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).

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 peat-lands 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.
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).

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 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.  covariates. A sharp decline in relative importance was observed for the other covariates.

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.

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 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. 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.  Table 3 for a description of covariates' abbreviations. 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 distancebased predictors were relatively important.
Nonetheless, our prediction errors were comparable to other studies. Rudiyanto et al. (2018) tested 14 machinelearning 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 (50 000 ha), while Gatis et al. (2019) obtained a RMSE of 0.31 m for a 40 600 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 2 and we used CCC, we cannot make direct comparisons, but both studies clearly outperformed ours. This may be due to the use of spatial crossvalidation 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.

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 scalerelated 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 gammaradiometric 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.

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.

Conclusion
The use of freely available covariates and DSM techniques provided the first maps of the depth to the mineral layer, coprogenous layer thickness, and maximum peat thickness covering five peatlands at a regional scale. Although their prediction error was comparable to other DSM studies, the CLT model did not achieve a sufficient accuracy to produce CLT and MPT maps of similar precision to site sampling for now and further work is needed. Thus, we were not able to elaborate on the contribution of the individual covariates on the CLT model performance. Although the DML map provides more accurate information as a tool to determine priority intervention zones, MPT still remains a key metric to guide soil conservation practices. The produced MPT map should serve as a baseline to be built upon and improved in future research. For instance, proximal sensing tools could be investigated as a more relevant source of covariate data to produce coprogenous layer thickness maps at a local scale instead of relying on remote sensing tools at a regional scale. Such could provide new insight on the CLT spatial variability and would also improve the accuracy of the MPT map.

Data availability
Data generated or analyzed during this study are not available due to the nature of this research. This data are of strategic importance for the research partners and is considered sensitive information.