Cultivated Organic soils in Montreal's southwest plain are the most productive soils in the province of Quebec. After their initial drainage to enable farming, Organic soils are susceptible to many forms of degradation and soil loss. In this study, we characterized the physical, chemical, and pedological properties of 114 sites from five peatlands to form soil conservation management zones. We attempted to use the maximum peat thickness (MPT) as a soil degradation proxy. The MPT can be defined as the thickness of the layer of peat until coprogenous or mineral materials are reached. The latter are undesired growing media and are not considered in MPT calculation. A series of multivariate analysis of variance indicated that MPT was moderately related to soil degradation (optimal model's Pillai's trace = 0.495). Three soil degradation groups were defined, separated by two MPT thresholds: 60 and 100 cm. When looking at 17 different depth-property combinations, shallower sites (MPT < 60 cm) showed signs of soil degradation significantly higher than sites with an MPT above 60 cm. The second threshold was proposed for practical purposes. Then, these thresholds were used to separate the study area into spatially distinct management zones. Important spatial contrasts were found. This supports the theory that precision agriculture techniques are needed to target fields to optimize soil conservation interventions. The relationship between the MPT and soil degradation should be further explored to account for other degradation factors, and to better identify degraded soils and soils at risk.
Introduction
In the province of Quebec, cultivated wetland Organic soils account for only 4% of the southern territory but support a strong, prosperous horticultural industry (Groupe AGÉCO 2007; Parent and Gagné 2010). Quebec is the second most important producer of fresh vegetables in Canada with 40% of the total Canadian acreage, of which 40% are cultivated Organic soils, equivalent to 9000 ha (MAPAQ 2018). The Montérégie region, south of Montreal, is home to 56% of the vegetable production acreage of the province. There, organic materials were deposited over many hundreds of years in anaerobic conditions, limiting peat degradation. Once drained for agriculture, Organic soils are in constant transformation, caused by subsidence, oxidation, and erosion (Kroetsch et al. 2011; Vepraskas and Craft 2015), depending on their initial carbon content and land use (Mann 1986). Low density and high volumetric porosity before initial drainage (78%–93%) are favorable conditions that allow important subsidence in the first years post-drainage (Driessen et al. 2000; Ilnicki and Zeitz 2003; van Asselen et al. 2010; Liu et al. 2016). Over time, surface layers will be mineralized and humified (i.e., decomposed) by microorganisms oxidizing peaty materials (Ilnicki and Zeitz 2003; Kroetsch et al. 2011). Combined soil losses, mainly wind erosion, equate to 1–5 cm·yr−1 (lnicki 2003; Esselami et al. 2014). Another threat to productivity comes in the form of an impervious, uncultivable layer made of coprogenous material. Limnic materials can be found between the top peaty layer and the mineral horizon (SCWG 1998). With soil loss, a coprogenous layer, when present, risks being exposed at the surface. Evolution of Organic soils is unidirectional and ultimately leads to degradation if no soil conservation measures are implemented.
Soil degradation has been repeatedly reported as a cause of modifications in physical and chemical properties. Bulk density can be higher with the decomposition of peat soils (Brandyk et al. 2003; Ilnicki and Zeitz 2003; Anshari et al. 2010; Hallema et al. 2015a; Wang et al. 2021). Dessureault-Rompré et al. (2018) reported that the development of a compacted layer is more important with soil degradation. The water retention capacity diminishes (Zeitz and Velty 2002). Moreover, modification of particle size, porosity, and soil structure affect the water retention curve (WRC; Shein et al. 2018). These parameters are modified with peat decomposition (Brandyk et al. 2003; Ilnicki and Zeitz 2003; Okruszko and Ilnicki 2003; Vepraskas and Craft 2015). Saturated hydraulic conductivity (Ksat) is a function of botanical origin as well as peat decomposition level (Ilnicki and Zeitz 2003; Gnatowski et al. 2010; Wang et al. 2021). Other properties are known to be modified by peat decomposition. Increase in pH and total nitrogen content could act as indicators of soil degradation if nitrogen fertilizer inputs are considered (Anshari et al. 2010). Trace elements adsorbed to surface soil could also serve as a degradation indicator, such as inorganic elements — Si, Fe, and Al (Ilnicki and Zeitz 2003). Organic matter and organic carbon contents are reduced with soil degradation (Hallema et al. 2015a; Krüger et al. 2015). To help manage peat soil degradation in wetlands, Sienkiewicz et al. (2019) suggested the use of the ratio between dissolved organic carbon and soil organic carbon as a measure of intensity of peat mineralization, also supporting the use of carbon loss as an indicator. The same authors also mentioned that the ratio could be more sensitive than the C/N ratio and the total nitrogen content as indicators of decomposition. The literature offers a handful of proposed critical thresholds, mainly for physical organic soil properties, that can be useful to classify soils based on their degradation. A meta-analysis showed that the botanical origin and decomposition level (von Post scale) of peat significantly affected hydraulic properties (Liu and Lennartz 2019). A critical value for bulk density was proposed by the authors: 0.2 g·cm−3. Lafond et al. (2015) proposed a limiting value of 10−4 cm·s−1 for the Ksat. Higher values are desired to ensure good agricultural drainage and to avoid root asphyxia. It is also possible to look at artificial growing media, where Caron and Michel (2021) proposed that an air content at −1 kPa inferior to 0.10–0.15 cm3·cm−3 is problematic for crops. Soil compaction can affect water and air flow in the soil matrix. Soil penetration resistance can assess the presence of a compacted layer. Growing roots can exert a maximal force of 0.9–1.3 MPa to colonize the soil (Bengough and Mullins 1990). No clear relationship has been proposed to link the force of growing roots to measured soil penetration resistance in organic soils. Thus, it is safe to assume that 1.3 MPa could serve as a potential soil compaction threshold. However, 0.5 MPa could also serve as an additional potential threshold because this value is associated with perched water table formation, and it is known to lead to low oxygen levels limiting root development.
Since many peat degradation indicators are available, but there are few critical thresholds, there is a need to find a simpler proxy that could simplify soil conservation management operations of Organic soils. In the province of Quebec, three Great Groups can be found with increasing degrees of decomposition: Fibrisols, Mesisols, and Humisols (SCWG 1998). However, a taxonomic distinction in itself is not enough to manage these highly fertile soils. Furthermore, simply relying on rubbed fibers or sodium pyrophosphate index (SCWG 1998; Parent and Caron 2007) to evaluate peat decomposition level is not efficient in managing cultivated peatlands at regional scale. Recently, two large-scale studies were carried out in Montérégie to create groups of soils linked to peat decomposition (Hallema et al. 2015a, 2015b). While the group separation is clearly related to hydraulic properties and to the number of years since land conversion for agriculture, the results imply the use of expensive measurements of soil properties to diagnose and recommend soil conservation measures. Expert knowledge and farmer experience suggest that shallow soils exhibit important degradation signs and limit crop yields. This observation is supported by the fact that freeze–thaw and shrink–swell cycles and microorganism activity are more intense near the surface (Zeitz and Velty 2002; Ilnicki 2003; Okruszko and Ilnicki 2003). Therefore, it was hypothesized that cultivated Organic soils could be divided into at least two groups with differing soil degradation levels related to their cultivation depth, affecting both physical and chemical properties. We are now introducing the maximum peat thickness (MPT) as the thickness of the surface peaty layer down to coprogenous or mineral materials. This definition excludes the coprogenous layer from the cultivable organic soil thickness because of its unsuitable properties for agriculture. The maximum peat thickness appears preferable to organic layer depth as a decision-making tool since it can be converted to a residual number of cultivation years using annual losses. Therefore, it can be translated to arable land value. The specific objectives of this paper were to perform an initial characterization of soils from the study area to allow a temporal monitoring of soil properties and to define management zones based on soil degradation discriminating properties affected by the MPT.
Materials and methods
Site selection
The visited sites are distributed among five drained and cultivated peatlands located in the southwest plain of Montreal, Quebec. These peatlands are in three separate watersheds (Lamontagne et al. 2014). The mean annual rainfall is 1000 mm and the mean annual temperature is 6.0 °C.
The sampling design was centred around contrasted sites belonging to 14 research-partner farms. Thus, expert knowledge was used to cover the study area and to capture the effects of shallow and deep deposits on soil properties. Precise measurements of depth were unknown before sampling. A total of 120 sites were sampled through 2018 and 2019, but 114 were kept for statistical analyses after the removal of sites with missing samples. Due to time and budget constraints, no replicates were collected regarding soil samples for chemical and physical analyses.
Pedological characteristics
At each site, a soil profile was described. Spatial coordinates and elevation were recorded with a differential GPS receiver (GENEQ GNSS F90) and a GPS data collector (SXPad 1000P). Figure 1 shows the location of each site and the study area.
Fig. 1.
Regional view of the final selection of 114 sites visited in Montérégie-Ouest, Quebec, in the summers of 2018 and 2019. [Colour online]
Depth to the mineral layer was measured at every location, even when it exceeded the 1.6 m control section, which is the standard soil profile depth for the Organic soil order (SCWG 1998). Soil pits were dug with a shovel until 70 cm, and then a Macaulay sampler (Eijkelkamp peat sampler) was used to extract 50 cm long soil cores until the underlying mineral layer was reached. At every location, various morphological traits were described, such as thickness of layers, fiber and rubbed fiber content, coarse fragments, structure, and botanical origin of fibers. Afterwards, the MPT was obtained at each site by subtracting the coprogenous layer thickness from the depth to the mineral layer. Soil profile description was based on the Manual for describing soils in the field (Expert Committee on Soil Survey 1982) and the Canadian System of Soil Classification (SCWG 1998).
Physical properties
Saturated hydraulic conductivity and bulk density
An aluminum cylinder with an internal diameter of 82 mm and a height of 55 mm was carefully pushed into the soil matrix to limit soil disturbance. Then, it was extracted with a shovel. A knife was used to trim the excess soil on both sides of the cylinder before a cloth was put in place with a rubber band to contain the soil sample. The cylinder was then placed in a plastic bag and was sealed to keep the sample moisture. Then, the hole was deepened to the next sampling depth. At each site, the soil was sampled at 0–5, 30–35, and 50–55 cm. Once transported to the laboratory, samples were kept at 4 °C before performing lab tests.
Saturated hydraulic conductivity was evaluated with the constant head method (Reynolds 2007). The samples were saturated over 24 hours in a tank. A constant water head was maintained with a Mariotte reservoir. The outflow was measured with calibrated pressure transducers in the Mariotte reservoirs and Ksat was derived using Darcy's law. Data were gathered using Loggernet (version 3.4.1, Campbell Scientific, 2007) software, and the calculation of Ksat was done in Microsoft Excel.
Once Ksat had been measured, the soil was dried and the bulk density was calculated using the equation from Hao et al. (2007). The temperature was adapted to Organic soils, and the soil samples were dried at 65 °C instead of the standard 105 °C (Sheppard and Addison 2007; Vepraskas and Craft 2015).
Water retention curve
The water retention curve can be obtained using one of two general ways: in a laboratory with undisturbed soil samples, or in situ with tensiometers and a time-domain reflectometry probe (Brandyk et al. 2003). In both cases, the water content is studied as a function of the soil matric potential. Since the number of sites and samples were important, the laboratory approach was selected. A second series of three cylinders was extracted side by side at the same depth as that used for Ksat measurements. The soil samples were kept at 4 °C until laboratory manipulations were performed.
For the first WRC points at low tension, a tension table with glass beads was used (CPVQ 1997a). At high matric potential, soil samples were put on a scale with a tensiometer, relying on evaporation losses to complete the curve (Šimůnek et al. 1998). Loggernet software was used to gather and save data. Some samples from 2018 were put in a pressure extractor to obtain the last WRC points instead of using the evaporation method (Reynolds and Topp 2007).
Soil penetration resistance
A numerical penetrometer with a static cone form Eijkelkamp (Eijkelkamp Soil and Water, Giesbeek, Netherlands) was used to collect the soil penetration resistance (PR) in the first 80 cm of soil. The device produces a curve of PR in MPa as a function of depth. Each curve is georeferenced with an integrated GPS receiver. The data can then be exported with the PenetroViewer software (version 6.08, Eijkelkamp, 2011).
Five curves were obtained at each site around the soil profile. Then, they were averaged for each centimeter and the mean PR curve was filtered with Guedessou's (2020) R script to correct for friction noise and experimenter bias. The final product was an average curve with values at 5 cm increments.
Particle size distribution
Particle size distribution is indicative of surface soil erodibility. It was calculated on the first 20 cm of soil (see the chemical analyses section). A sample of at least 100 g of air-dried soil was placed in the upper sieve of a series of nine sieves stacked on a sieve shaker (CPVQ 1997b). Each sieve had a different mesh size ranging from 16 to 0.1 mm. The erodible fraction was assessed based on the percentage of soil mass that went through the 0.84 mm sieve after one minute of shaking (López et al. 2007).
Chemical properties
Concerning chemical analyses, three bags were filled with soil extracted with a Dutch soil sampler probe at different depths: 0–20, 20–40, and 40–60 cm. Soil from four to five locations around the soil profile was sampled and mixed in each bag to obtain a sufficient volume for analysis. The bags were hermetically closed and stocked at 4 °C. Before performing the analyses, the soil samples were dried at 65 °C, ground, and sieved to 2 mm (CPVQ 1988).
One should note that sampling depths were fixed before the beginning of the study, for chemical and physical properties, no matter what the MPT would be. First, this methodology allowed comparison of the properties between sites at the same depths. Second, roots of vegetable crops are mostly found in those depth ranges. Since organic soil degradation and soil loss occur in the topmost layers, it would not have been justified to sample at greater depths (i.e., 70 or 80 cm).
Loss on ignition
The organic matter content was estimated using the loss-on-ignition method. The analysis was performed on 10 g of soil (MDDELCCQ 2017). The soil sample was combusted at 550 °C over 16 hours.
Sodium pyrophosphate index
The sodium pyrophosphate index was obtained with the colorimetric determination procedure (Parent and Caron 2007). The absorbance of a sample of organic soil that reacted with sodium pyrophosphate was read with a UV–VIS Cary 100 Bio spectrophotometer with CarywinUV software (Agilent Technologies Inc.). The degree of peat decomposition can be estimated with this method. With this index, higher values are associated with more degraded peat materials. It is a repeatable laboratory method that can add to the more empirical unrubbed and rubbed fiber field method.
Total nitrogen content and total carbon content
The total nitrogen content and total carbon content were estimated with the LECO macro combustion method with an 828 Series Carbon Nitrogen Protein Determinator (LECO, Michigan, USA) following the manufacturer's protocol.
Soil pH in water and electrical conductivity
The method from CEAEQ (2014) was followed to measure soil pH in water and electrical conductivity. A 1:5 soil–water extraction ratio was used instead of the standard 1:1 and 1:2 extraction ratios. Soils with high organic matter content tend to form a paste instead of a liquid solution, which complicates measurements (Hendershot et al. 2007) A sympHony SB70P pH meter was used (VWR, Pennsylvania, USA) to measure the pH, and a sympHony SB70C conductivity meter (VWR, Pennsylvania, USA) was used to measure electrical conductivity.
Statistical analyses
Spatial manipulations and statistical analyses were performed using QGIS version 3.10.13 (QGIS Development Team 2020) and R version 3.6.3 (R Core Team 2020).
Data preparation and exploration
A total of 53 out of 65 variables (depth-property combinations) were kept in the final matrix. Some pedological properties were removed since they were less relevant and appeared less empirically correlated with the soil degradation phenomenon. Appendix 1 provides a table of descriptive statistics of all final variables. Many variables showed values for skewness and kurtosis higher than 2 (Tabachnick and Fidell 2013). Hence, the R library bestNormalize package (Peterson 2020) was used to automatically find the best transformation for each variable.
The presence of univariate relationships between the variables and the MPT was explored through scatterplots. The objective was to find strong relationships and link them to soil degradation before the multivariate analyses. If multivariate analyses fail, promising univariate relationships could be further explored.
Principal component analysis
A principal component analysis (PCA) was used to select a smaller set of representative and uncorrelated variables among the final selection of 53 variables to obtain a ratio of dependent variables (DVs) to observations (N) lower than 1:10 (Tabachnick and Fidell 2013). The PCA was performed on the correlation matrix of the transformed variables with the psych library (Revelle 2020). A varimax rotation was used. Concerning PCA assumptions, normality and factorability of the correlation matrix were respected.
Following PCA, 10 components were retained. Each of them had an eigenvalue >1, meaning that a component explained at least more variance than a raw variable. For each retained component, the most correlated variable was noted. These 10 chosen variables represented 73.3% of the total variance of the original data set. These variables (and depths) are total nitrogen and total carbon content (0–20 cm), pH (20–40 cm), water content at −3 kPa (0–5 cm), available water content at −10 kPa (0–5 cm), air content at −5 kPa (30–35 and 50–55 cm), total porosity (50–55 cm), depth to the compacted layer, and maximal soil penetration resistance within the first 80 cm. The presence of multivariate outliers in the data set was evaluated. The Mahalanobis test (χ2 for α = 0.001 for 10 df) revealed the absence of outliers for the 114 final observations and the 10 chosen variables.
Multivariate analyses of variance (MANOVAs)
The proposed approach is inspired from Hallema et al. (2015a), who used incremental thresholds of reclamation year to form statistical groups of fields with significantly different properties. To determine the best MPT threshold to define the soil degradation groups, a series of MANOVAs was performed. The 10 variables from the PCA acted as DVs, while dichotomized MPT was the independent variable. The MPT was dichotomized using different thresholds, incremented by 10 cm from 40 to 130 cm to divide the data set into two soil degradation groups. MANOVAs were performed with the manova function (R Core Team 2020). Alpha = 0.05 was used.
The optimal threshold was selected based on two criteria. The first is Pillai's trace. This statistic can have a value between 0 and 1, where a high value reflects a larger effect of the independent variable on the DVs. Pillai's trace is accompanied by its 95% confidence interval (CI) limit derived from 5000 bootstraps. In short, bootstraps generate new data sets of the same sample size, on which MANOVAs will be performed, by randomly sampling with replacement of the original data set. Pillai's trace is then extracted for each MANOVA, and the CI can be derived by aggregating the results from the 5000 iterations and identifying the 2.5th and 97.5th percentiles.
The second criterion is the number of significant analyses of variance (ANOVAs) between the two groups for the 53 variables. The higher the number of properties significantly different between the soil degradation groups, the most relevant an MPT threshold will be considered. Moreover, ANOVAs were used to understand how soil properties differed between soil degradation groups. The Bonferroni adjustment gave control over type I error (Tabachnick and Fidell 2013). In more detail, the 0.05 alpha threshold was divided by the total number of ANOVAs to obtain a new alpha value of 0.000943. Therefore, to consider a property significantly different between soil degradation groups, this severe alpha value was preferred to be over 0.05. Assumptions of normality, independence of observations, and homogeneity of variance were verified for all ANOVAs. Some variables showed signs of heteroscedasticity. In those cases, we used a linear mixed model with heterogeneous variances.
After a first MPT threshold was chosen, a second series of MANOVAs was performed to test if a second MPT threshold could be used to form a third soil degradation group. Increments of 10 cm were tested from 80 to 170 cm. A third group could be needed to facilitate soil conservation. Indeed, the first group would be degraded soils, the second group would be considered at risk of degradation or as priority intervention zones to maintain productivity, whereas a third group would include deposits deep enough not to be considered a priority if conservation resources are limited (i.e., funds or material). Analyses of variance were also performed to find significant differences in soil properties between groups. The Tukey's method was used to adjust multiple comparisons for significant variables.
Results
Relationships between the MPT and soil properties
Plots between the MPT and soil properties (Figs. 2 and 3) were generated as a first univariate approach, before computing MANOVAs, to explore potential relationships. Inverse exponential relationships and inverse relationships can be observed between some properties and the MPT, mostly related to carbon content. At a depth between 15 and 75 cm, a change seems to occur in the scatterplots. It is important to note that only a small number of properties seemed to be related to the MPT and that result might differ in the multivariate feature space.
Fig. 2.
Scatterplots of the maximum peat thickness and (a) the organic matter content between 40 and 60 cm, (b) the total carbon content between 40 and 60 cm, (c) the total nitrogen content between 40 and 60 cm, and (d) the total carbon content between 0 and 20 cm.
Fig. 3.
Scatterplots of the maximum peat thickness and (a) the bulk density at 50–55 cm and (b) the electrical conductivity at 20–40 cm.
Concerning Figs. 2 and 3, some observation points have an MPT smaller than the depth at which a property has been sampled. This is due to the fixed sampling depths. In other words, some samples are from coprogenous or mineral materials. Interestingly, a change in material is not the only factor affecting properties. It is noteworthy that properties measured closer to the surface (0–5 and 0–20 cm) also show signs of degradation linked to depth even if the sampled material is peaty.
Maximum peat thickness threshold based on MANOVAs
A series of MANOVAs were performed with an incremental MPT threshold to form statistically significant groups linked to soil degradation. The results are given in Table 1. Comparison of results indicated that an MPT threshold of 60 cm would be the most appropriate, based on Pillai's trace and the number of statistically significant ANOVAs between groups. At this threshold, the data set was divided into two soil degradation groups of 12 and 102 observations. Combined DVs were significantly different between groups at this MPT threshold, F(10 103) = 10.099, p < 0.001. The results showed a moderate association between soil degradation groups and the combined DVs, Pillai's trace = 0.495 with a 95% CI between 0.36 and 0.68. Homogeneity of variance–covariance matrix assumption was rejected following Box's M test (p < 0.05). Multivariate normality was respected for kurtosis but rejected for skewness. The use of the Pillai's trace metric is suggested when MANOVA assumptions are rejected (Olson 1974). Hence, the conclusions are still valid. Lack of multicollinearity between DVs is assumed.
Table 1.
Series of multivariate analyses of variance (MANOVAs) with group-forming thresholds based on maximum peat thickness incremented between tests to define the first threshold.
Testing a second significant MPT threshold
Once the 60 cm MPT threshold was obtained, a second series of MANOVAs was performed to test the significant addition of a second MPT threshold to form a third soil degradation group. Table 2 provides results from this second series of MANOVAs. Although all MPT thresholds were statistically significant, Pillai's trace ranged from 0.295 to 0.334, remaining lower than two-group MANOVAs (Table 1: 0.260 to 0.498).
Table 2.
Series of multivariate analyses of variance (MANOVAs) with two maximum peat thickness thresholds, a first fixed at 60 cm and a second incremental threshold.
Statistically speaking, there is no evidence to support the need of a second MPT threshold. However, given the fact that consequences of soil degradation are severe at 60 cm, a second MPT threshold could be used to identify some early signs of degradation. Indeed, considering that drain tiles are often placed at a depth of 90 cm, soil between 60 and 100 cm may face drainage limitations due to decomposition and compaction. It would be advised to use a similar depth as an MPT threshold. Hence, we proposed a second threshold at 100 cm. Therefore, as previously stated, this third soil degradation group represents soils that do not need immediate and important soil conservation investments. Their productivity might be less at risk in the short term, unlike soils near 60 cm, but they may need a close follow-up in the near future.
Differences in soil degradation between groups
The MANOVAs revealed a significant difference between the soil properties of the observations above and under the 60 cm MPT threshold. To understand the magnitude of these differences, the results of the 17 statistically significant ANOVAs are given in Table 3.
Table 3.
Mean and standard error (SE) of 17 significantly different variables between Group 1 (<60 cm) and Group 2 (>60 cm).
All variables show signs of degradation in terms of soil quality in group 1. Group 1 comprises sites that are more humified (lower total carbon content, organic matter content, C/N ratio, and fiber and rubbed fiber content), more compact (higher bulk density, lower air content, and available water content), and more saline (higher electrical conductivity). As expected, the presence of mineral soil and coprogenous soil near the surface affected soil properties in the rooting zone.
To add to the previous table of results, Figs. 4–10 provide the density distribution of observations per soil degradation group for the 17 variables through violin diagrams. Violin plots are composed of a boxplot and a symmetrical density histogram. The box contains 50% of observations and is also defined as the interquartile range (IQR). Two whiskers, each usually representing 25% of the observations when no outliers are present, extend on each side of the box. Outliers are defined as the value on each side of the box ± 1.5 IQR. For instance, Fig. 4 shows that the organic matter content distribution within a group is similar except for the last depth (40–60 cm). Indeed, some sites had mineral materials (organic matter content <30%). A similar pattern is observed for the total carbon content (Fig. 10). Concerning bulk density and electrical conductivity, a wider range of values characterizes group 1 at greater depths (Figs. 5 and 6). For Figs. 7 and 9, similar distributions per group are observed at different depths for the C/N ratio and for the air content at −5 kPa. On a general note, one can appreciate the narrower IQR for Group 2, meaning more variability is present in Group 1, except for fiber and rubbed fiber content and the C/N ratio (Figs. 8 and 9). Indeed, the more degraded soils from Group 1 show less variability for these properties.
Fig. 4.
Violin plots of the organic matter content according to the degradation group at three sampled depths: (A) 0–20 cm, (B) 20–40 cm, and (C) 40–60 cm. [Colour online]
Fig. 5.
Violin plots of the bulk density according to the degradation group at two sampled depths: (A) 0–5 cm and (B) 50–55 cm. [Colour online]
Fig. 6.
Violin plots of electrical conductivity according to the degradation group at two sampled depths: (A) 20–40 cm and (B) 40–60 cm. [Colour online]
Fig. 7.
Violin plots (A) of air content at −5 kPa and (B) of available water content at −10 kPa according to the degradation group at a sampling depth of 50–55 cm. [Colour online]
Fig. 8.
Violin plots of the (A) mean fiber content and (B) mean rubbed fiber content of the maximum peat thickness. The fiber and rubbed fiber content (in percentage) of every organic soil layer of the soil profile was multiplied by their corresponding layer thickness. Then, all values were summed up and divided by the maximum peat thickness of that site. [Colour online]
Fig. 9.
Violin plots of the C/N ratio according to the degradation group at three sampled depths: (A) 0–20 cm, (B) 20–40 cm, and (C) 40–60 cm. [Colour online]
Fig. 10.
Violin plots of the total carbon content according to the degradation group at three sampled depths: (A) 0–20 cm, (B) 20–40 cm, and (C) 40–60 cm. [Colour online]
Analyses of variance were also performed on the three-group solution using the 60 and 100 cm MPT thresholds. Eighteen variables were significantly different between groups (Table 4). Sites from Groups 1 and 2 belong to different classes after Tukey's adjustment for multiple comparisons of the means. Group 3 is mostly of the same class as Group 2. Nonetheless, Group 3 is statistically different from Groups 1 and 2 for three soil properties: the organic matter content at 40–60 cm, the C/N ratio at 40–60 cm, and the thickness of the coprogenous layer in the first 160 cm, indicating early signs of degradation. The thicker the cultivable peaty layer, the higher the organic matter content and C/N ratio. The organic matter content doubles from Group 1 to Group 2. Concerning the thickness of the coprogenous layer, it seems that sites from Group 2 have a thicker layer in average.
Table 4.
Mean and standard error (SE) of 18 significantly different variables between Group 1 (<60 cm), Group 2 (60–100 cm), and Group 3 (>100 cm).
Spatial distribution of management zones at a regional scale
Predictive maps of the depth to the mineral layer and of the coprogenous layer thickness were produced with this study's data set combined with other data sets gathered within the last 10 years (Fig. 11). By subtracting the latter from the former map, the MPT was estimated at a regional scale. Both final maps were generated using the cubist machine-learning model with digital soil mapping techniques. Environmental covariates gathered for the study area were spatially intersected with the georeferenced soil observations to act as predictors of the depth to the mineral layer and to the coprogenous layer thickness. One should note that prediction errors added up, augmenting the prediction uncertainty in each location of the final map. The complete workflow can be found in Deragon et al. (2022). Both MPT thresholds of 60 and 100 cm were applied to the final map. Each cell of 10 m spatial resolution was classified to illustrate the distribution of soil degradation groups across the study area that form distinct management zones. Although predictive maps were generated for the complete extent of each peatland, only fields of the 14 research-partner farms are shown in Fig. 11.
Fig. 11.
Regional map of maximum peat thickness management zones based on the proposed thresholds: 60 cm and 100 cm. [Colour online]
First, it is possible to see that soils with a negative depth are predicted (black areas, Fig. 11). These areas can be seen where a thick coprogenous layer was predicted and subtracted from a shallower predicted depth to the mineral soil. Nonetheless, this group only represents 0.6% of the study area (Table 5).
Table 5.
Predicted area corresponding to the maximum peat thickness (MPT) management zones at a regional scale.
Second, degraded soils (Group 1) are mostly found near the borders of peatlands, which was expected, but it is not limited to this spatial pattern. A more degraded zone can be found near the centre of the southwest peatland. This zone is associated with deep organic deposits with a thick coprogenous layer. The peatland to its right is subject to a higher concentration of degraded soils. This concentration can be explained by a shallower deposit or the number of years since land conversion for agriculture of the fields at this location. Land use and annual loss can also affect the speed at which soil resources are lost.
To rephrase the importance of these results for soil conservation practices, 63.7% of the study area has an MPT higher than 1 m (Group 3). A fifth of the study area is at risk (Group 2), while 15.6% of the area is already degraded (Group 1). Also, the 14 farms do not share the same distribution of management zones over their land.
Discussion
The proposed approach, based on the MPT, allowed the definition of a 60 cm MPT threshold related to a significant change in organic soil properties and a 100 cm MPT threshold related to drainage considerations. Therefore, three groups related to soil degradation were formed (MPT < 60, 60 < MPT < 100, and MPT > 100 cm). Indeed, if the 60 cm MPT threshold is the answer to the question “At which maximum peat thickness a soil is sufficiently degraded to affect its properties?”, it does not answer the question “At which maximum peat thickness should we start intensive soil conservation practices?”. Since these practices have a related cost, it is imperative to apply them before a given field gets close to or under an MPT of 60 cm. In other words, soils with an MPT < 60 cm are degraded, while soils between 60 and 100 cm would be targeted for the application of intense and costly soil conservation practices to maintain soil productivity, and soils with an MPT > 100 cm would not be considered priority intervention zones. One should note that basic soil conservation practices should be applied no matter what the MPT depth (i.e., cover crops after harvest, controlled traffic farming, irrigation during windy periods, shelterbelts, etc.), but intensive and more costly measures such as the addition of biomass crop amendments (Dessureault-Rompré et al. 2020, 2022) should be targeting soils that are near 60 cm. Since the practical MPT threshold of 100 cm is not supported by statistical evidence, knowledge of conservation practices efficiency and costs should be integrated in the future to support or modify the proposed threshold. Pillai's trace of 0.495 indicated a moderate relationship between degradation groups and DVs. This result suggests that other degradation indicators should be studied to obtain a stronger relationship. Such indicators could include the number of years since land conversion to agriculture, the botanical origin and type of wetland, and agricultural practices. Furthermore, outliers in Figs. 4C, 5A, and 10C support this observation, hinting toward a need for other properties to obtain a better distinction between soil degradation groups.
The properties presented in Tables 3 and 4 indicated that shallower soils (Group 1) are relatively more degraded than deeper soils (Groups 2 and 3). These properties are mostly related to the carbon content, C/N ratio, bulk density, soil salinity, and WRC. The results are coherent with the expected change in Organic soils properties after drainage. The observed loss of porosity concords with expectations (Brandyk et al. 2003; Ilnicki and Zeitz 2003; Okruszko and Ilnicki 2003; Hallema et al. 2015b). Decomposition of surface layers leads to modifications of the C/N ratio, which was also observed, since nitrogen is mineralized and carbon is released by microorganism respiration (Ilnicki and Zeitz 2003; Anshari et al. 2010; Kruger et al. 2015).
Concerning higher salinity in more degraded soils, this observation could be linked to a migration of solutes by capillarity rise from underlying, saltier mineral and coprogenous layers, fertilization, or saline irrigation water. Salinity sensitivity of vegetable crops can vary significantly. Crops commonly found in the study area, such as carrots (1000 µS·cm−1), onions (1200 µS·cm−1), celery (1800 µS·cm−1), lettuce (2000 µS·cm−1), and spinach (2000 µS·cm−1), are sensitive and moderately sensitive to soil salinity (Machado and Serralheiro 2017). Possible yield loss can be expected when the soil salinity exceeds these thresholds. This was the case in most degraded sites.
The lower fiber content in soils from Group 1 can be attributed to a more advanced stage of degradation. Chemical and physical alterations are more important near the surface with the presence of oxygen and with tillage (Ilnicki and Zeitz 2003; Kroetsch et al. 2011; Pawluczuk et al. 2019). Moreover, the mean surface bulk densities from each of the two groups exceeded the 0.2 g cm−3 threshold proposed by Liu and Lennartz (2019). This result can be explained by soil degradation mostly occurring at this depth. With tillage, soil aggregates and peat residues break down into particles that can rearrange themselves and interlock in a more compacted structure (Okruszko and Ilnicki 2003). The mean bulk density at 50–55 cm shows a different result. Indeed, deeper soils have lower values of bulk density, while shallower soils have values around three times the threshold. A higher density of mineral particles contributed to higher bulk densities in group 1. Wang et al. (2021) found a significant difference in Ksat and bulk density values between degrees of degradation of peat soils, whereas in this study, only bulk density was significantly different between soil degradation groups related to the MPT. Other factors, such as botanical origin, can have a major impact on Ksat and were not accounted for in our study. Soil air content is a physical property related to soil compaction and essential to gas diffusion. The observed values are inferior to the threshold proposed by Caron and Michel (2021) for shallower soils. Accordingly, the available water content is also inferior.
A practical application of the results would suggest that soils under 60 cm, in which vegetable crops will undoubtedly have a shallower rooting depth, also have a lower water retention capacity and limiting gas diffusion properties. This results in an agronomically undesired growing medium for crops with major consequences on irrigation practices: applied quantities and optimal timing for irrigation must be adapted. Soil properties significantly different between groups will be used as discriminant properties to help manage fields in a soil conservation context. That being said, the studied soil properties should not be limited to the 17 variables previously mentioned. For instance, Ksat is a key physical property for drainage plans and should also be integrated in a wider decision frame when it comes to soil management, given the magnitude of degradation of a particular field. Furthermore, our studies were already ongoing when Sienkiewicz et al. (2019) suggested the use of the ratio between dissolved organic carbon and soil organic carbon as a better indicator of soil degradation. This ratio could be explored in the future to perhaps witness a bigger effect of the MPT on the DVs. The MPT management zone map showed visible contrasts between peatlands. The presence of spatial patterns of soil degradation can have major implications on conservation costs. The distribution of the 14 research-partner farms is not homogeneous. For instance, four farms would be severely impacted by soil degradation. Located in the southeast peatland where a great concentration of soil under 60 cm was predicted, this could potentially mean that their soil is already severely affected by soil degradation and that soil resources are limited. Therefore, more conservation practices should be put in place compared with farms completely found in MPT > 100 cm areas.
Limitations of this study and future work
There are limitations to the use of the MPT. Little is known about soil volume changes due to seasonal and shrinkage–swelling cycles (Ilnicki and Zeitz 2003) at a field scale. Uncertainty related to the evaluation of the MPT at a given point over the years is thus unknown. The magnitude and duration of these cycles may likely vary with space, necessitating complex predictive models calibrated with samples (Camporese et al. 2006; Morton and Heinemeyer 2019). Grading operations should also be considered since they can affect the MPT. These variations could not only affect some properties but also change the soil management zone classification of a site. Grading could have had a potential effect on this study for sites visited over two years.
Furthermore, available techniques to measure peat depth do not always discriminate coprogenous materials from peat and the mineral layer. They all have uncertainties, whether it is manual probing or proximal sensing-based methods (Parry et al. 2014; Ji et al. 2019). The predictive map of the MPT only allows visualization of trends and gradients at a regional scale. Its associated uncertainty is high, and it would be advised to use it as a complementary tool until further research is performed with a bigger data set. It is also important to note that our study did not differentiate between soil degradation of shallower deposits near a peatland border (i.e., degradation linked to the peatland's genesis) and degradation related to anthropic activity. Furthermore, one should not assume that human activity alone causes the shallowness of soils observed near peatland borders.
The initial sampling design goal was to provide as much information as possible on the depth and diversity of soil properties at a regional scale by covering both problematic and good fields. Another goal was to obtain data on the depth of the organic deposits in three of the peatlands that were not covered in prior studies. Legacy soil data available for these peatlands, in the form of pedological reports, were acquired at different time periods (1950, 2000, and 2001), and no maps of the depth to the mineral layer or the coprogenous layer thickness were produced. Now that such maps exist, they could be used to elaborate a more robust sampling scheme, meaning one could use a stratified sampling design to sample sites with regular MPT increments of 10 cm. This methodology would increase the number of observations in the shallower soil group and more statistical robustness can be expected. More observations would be needed to obtain more robust conclusions, although they remain valid. Due to time and budget constraints, expert knowledge was used to determine the location of the sites. Although not rigorous, interesting results have been found and will further be investigated.
Conclusion
The MPT of cultivated Organic soils appeared to be a property effectively related to soil degradation. This property could be used to help manage soil conservation practices. A rich database of soil properties at three depths was assembled to help monitor soil degradation. Seventeen soil properties differed significantly between groups formed with two proposed MPT thresholds of 60 and 100 cm. Three management zones were derived to assess soil resource distribution and to facilitate the determination of priority intervention zones. To improve the methodology, it is crucial to consider other data to elaborate an efficient and complete soil conservation plan. For instance, soil degradation factors (soil erosion loss quantification and mineralization rates) and the potential of soil conservation practices (effectiveness and impact on peat decay) must be part of the decision-making process. Indeed, certain deeper zones can experience higher soil loss rates than a shallower one. In other words, the approach should further be explored by integrating other decision-making criteria to determine intervention zones.
Funding information
The authors are thankful for the financial support from the Canadian Graduate Scholarship program by the Natural Sciences and Engineering Research Council of Canada (NSERC) and the Master's Scholarship program (B1X) by the Fonds de recherche du Québec — Nature et technologies granted to R.D. We also acknowledge the financial support of the NSERC through an Industrial Research Chair Grant in Conservation and Restoration of Cultivated Organic Soils (IRCPJ 411630-17) in partnership with Delfland Inc., Productions Maraîchères Breizh Inc., La Production Barry Inc., Les Fermes R.R. et Fils Inc., Le Potager Montréalais Ltée, R. Pinsonneault et Fils Ltée, Patate Isabelle Inc., Les Fermes du Soleil Inc., Les Jardins A. Guérin et Fils Inc., Le Potagers Riendeau Inc., Vert Nature Inc., Fermes Hotte et Van Winden Inc., Production Horticole Van Winden, and Maraichers J.P.L. Guerin and Fils.
References
Appendices
Appendix 1
Table A1. Descriptive statistics of the final selection of 53 depth-property combinations.
