# Soil texture explains soil sensitivity to C and N losses from whole-tree harvesting in the boreal forest

Rock Ouimet, Nathalie Korboulewsky, Isabelle Bilger

<figure><img src="image_1.png" /></figure>

## To cite this version:

Rock Ouimet, Nathalie Korboulewsky, Isabelle Bilger. Soil texture explains soil sensitivity to C and N losses from whole-tree harvesting in the boreal forest. Soil Systems, 2023, 7 (2), pp.39. 10.3390/soil-systems7020039. [https://doi.org/10.3390/soil-systems7020039](https://doi.org/10.3390/soil-systems7020039). [hal-04211043](https://HAL.inrae.fr/hal-04211043)

HAL Id: hal-04211043

https://hal.inrae.fr/hal-04211043

Submitted on 19 Sep 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

Article

# Soil Texture Explains Soil Sensitivity to C and N Losses from Whole-Tree Harvesting in the Boreal Forest

Rock Ouimet<fnref n="1" />*, Nathalie Korboulewsky<fnref n="2" /> and Isabelle Bilger<fnref n="2" />

**Abstract:** The use of forest biomass to produce energy is increasingly viewed as a means to reduce fossil fuel consumption and mitigate global warming. However, the impact of such practices on soils in the long term is not well known. We revisited forest plots that were subjected to either whole-tree (WTH, n = 86) or stem-only (SOH, n = 110) harvesting 30 years ago in the boreal forest in Quebec, Canada. The objective of the present study was to find soil properties that could explain the lower soil C and N stocks at the sites subjected to WTH compared to SOH after 30 years. Compared to SOH, lower soil C and N stocks attributable to WTH occurred when soil particle content <20 µm was below 30%. The theoretical separation of soil organic matter into two fractions according to soil particle content <20 µm—a recalcitrant and a labile fraction—could explain the observed pattern of soil C and N differences between WTH and SOH. Imperfect or poor soil drainage conditions were also associated with lower soil C and N in WTH compared to SOH. Limiting additional biomass harvesting from these sites would help to preserve soil C and N from potential losses.

**Keywords:** boreal forest soils; carbon; drainage; mineral-associated organic matter; nitrogen; soil fine particles; whole-tree harvesting

<figure><img src="image_4.png" /></figure>

check for
updates

**Citation:** Ouimet, R.; Korboulewsky, N.; Bilger, I. Soil Texture Explains Soil Sensitivity to C and N Losses from Whole-Tree Harvesting in the Boreal Forest. *Soil Syst.* **2023**, *7*, 39. https://doi.org/10.3390/soilsystems7020039

Academic Editor: Klaus Von Wilpert

Received: 18 February 2023

Revised: 16 March 2023

Accepted: 5 April 2023

Published: 18 April 2023

<figure><img src="image_5.png" /></figure>

**Copyright:** © 2023 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article
distributed under the terms and
conditions of the Creative Commons
Attribution (CC BY) license (https://
creativecommons.org/licenses/by/
4.0/).

## 1. Introduction

According to the IPCC [1], the provision of biomass for bioenergy and other bio-
based products represents an important share of the total mitigation potential for global
warming associated with agriculture, forestry, and other land uses. However, the use
of bioenergy can lead to either reduced or increased emissions, depending on the scale
of deployment, conversion technology, the amount of fossil fuel displaced, and on how
and where the biomass is produced [2]. This is particularly true regarding forest biomass.
Literature reviews and associated meta-analyses show that whole-tree harvesting (WTH)
or forest residue harvesting can reduce soil carbon (C) and nitrogen (N) reserves and
other mineral nutrients compared to stem-only harvesting (SOH) [3–5]. Since forest soils
are among the largest terrestrial C reservoirs [6,7], C (and N) losses therein caused by
certain biomass harvesting practices may become an unforeseen source of atmospheric
CO₂, thereby impeding the ultimate objective of residual biomass harvesting for energy
production, which is to reduce global CO₂ emissions.

A striking example of the soil C and N losses associated with biomass harvesting has been observed for WTH compared to SOH in boreal forest stands located in distinct soil landscapes in Quebec, Canada [8]. In this study, we sampled 196 forest plots in the boreal forest in four so-called ‘soil provinces’ with contrasting soil types. The plots had been harvested 30 years before with either WTH or SOH. Mineral soil C and N stocks (0–60 cm soil depth) were on average 34% lower with WTH than with SOH in the Laurentian and the Mistassini Highlands after 30 years, while there was virtually no difference between the two harvesting intensities in the Appalachians or the Abitibi and James Bay Lowlands. In

<footnote n="1">¹ Direction de la Recherche Forestière, Ministère des Ressources Naturelles et des Forêts du Québec, Complexe Scientifique, 2700 Einstein, Quebec, QC G1P 3W8, Canada</footnote>

<footnote n="2">² INRAE—UR Écosystèmes Forestiers, Domaine des Barres, 45290 Nogent-sur-Vernisson, France</footnote>

<footnote>* Correspondence: rock.ouimet@mffp.gouv.qc.ca; Tel.: +1-418-643-7994</footnote>
the two former soil provinces, the estimated extra biomass harvested with WTH compared to SOH amounted to 21 Mg ha⁻¹ on average, which is approximately 10.5 Mg C ha⁻¹ and 0.15 ± 0.01 Mg N ha⁻¹, but this increase in biomass retrieval increased soil C loss to approximately 40 Mg C ha⁻¹ and 2 Mg N ha⁻¹ [8].

Soil C and N can take several decades to recover following harvest; according to James et al. [9], some Podzol C can only be recovered after at least 75 years. Therefore, research must define the characteristics that distinguish soils "sensitive" to extra forest biomass harvesting from those able to preserve their C and N stocks. Carbon management is the primary impetus for justifying the use of forest biomass as energy feedstock and soil carbon losses should be addressed in forest biomass harvesting policies and guidelines. Most existing guidelines have been established to encourage, or try to ensure, that biomass removal is environmentally sustainable [10]. Most of these guidelines address many issues, such as biodiversity, social values, water and aquatic ecosystems, and soil chemical fertility, acidity, erosion, structure, moisture retention, and compaction. However, the indicators they mention for the long-term storage of C and N derived from soil organic matter often remain vague and general due to the indicators' unknown reactivity or stability [11].

Soil C and N accumulation are influenced by soil texture and other ecological factors, such as vegetation type [7], soil drainage class [12], pedogenesis per se [13,14], and time since major disturbance [15,16]. In addition, soil organic matter (SOM) encompasses a vast array of highly diverse organic residues depending on their size and composition, whether they are in a free state, associated with mineral particle surfaces to varying degrees, or buried within aggregations of mineral particles [17]. These different forms of SOM lead to different properties, which especially influence their stability and lability [18]. In order to improve our understanding of SOM persistence, Lavallée et al. [19], and many others, have proposed separating SOM into two operational physical fractions: "particulate organic matter" (POM) and "mineral-associated organic matter" (MAOM). The two SOM fractions have very distinct compositions and properties [19,20]. For example, mean residence time for POM is from <10 years to several decades, while for MAOM, it ranges from decades to centuries. Therefore, POM is much more labile than MAOM, and is logically more at risk of loss after disturbances such as forest biomass harvesting. In addition, there is a strong relationship between the maximum SOM fraction as MAOM can accumulate in the soil and mineral soil particle size; MAOM is associated with fine particles <20 µm [21,22]. As total SOM is the sum of MAOM and POM, the estimation of potential maximum MAOM allowed us to estimate the corresponding minimum POM fractions from soil particle analyses.

The objective of this study was to find soil properties that could explain the soil C and N differences between WTH and SOH we had found in our earlier study [8], and therefore, to define indicators to identify soils sensitive to C and N loss caused by forest biomass harvesting. In particular, two hypotheses were designed to explain the differences in soil C and N between WTH and SOH:

1. The difference could be related to the existing C and N stocks in the soils;

2. The difference could be related to fine particle content.

In addition, the decomposition process of organic matter is hampered in poorly drained soils [23,24]. Therefore, we also tested if differences in total soil C and N between WTH and SOH were associated with the soil drainage class.

## 2. Materials and Methods

The experimental set-up and relevant soil sampling protocol are briefly described below. More details on forest stand composition and sampling can be found in [8].

### 2.1. *Study Area*

The study area was mainly located in the balsam fir–white birch bioclimatic domain in eastern Quebec and the southern part of the black spruce–feather moss bioclimatic domain in western Quebec, Canada, between 47° and 50° N latitudes. The study plots included four of the five soil provinces recognized in southern Quebec [25]: the Appalachians (B);
the Laurentians (C); the Abitibi and James Bay Lowlands (D); and the Mistassini Highlands (E) (Supplementary Materials, Figure S1). These soil provinces are distinguished according to the main parent material (geology and geomorphology), altitude (e.g., areas invaded by post-glacial seas), topography (slopes and landforms), soil texture and climate (temperature, rainfall, degree days, etc.; Supplementary Materials, Table S1). Black spruce (*Picea mariana* (Mill.) B.S.P.), balsam fir (*Abies balsamea* (L.) Mill.), and Jack pine (*Pinus banksiana* Lamb.) are the main tree species growing in this boreal forest.

## 2.2. *Experimental Design*

To determine the effect of harvesting intensity on natural regeneration, from 1982 to 1988, 562 circular plots (0.04 ha each, 11.28 m radius) in a completely randomized design were established just before clear-cut harvesting took place throughout the boreal forest of Quebec. The felling system was either WTH (feller buncher, mechanical harvester and cable skidder, or manual felling and cable skidder) or SOH (manual or mechanical felling and delimbing on-site, then stem extraction by cable skidder). No differences in tree species composition or abundance after natural regeneration were observed between WTH and SOH plots, even up to 10 years after harvesting [26].

In 2011, we selected 196 of these plots among the four soil provinces (n-SOH = 110; n-WTH = 86) based on their accessibility and their similarities regarding pre-harvest stand composition, surface deposit types and depth, and drainage class. Within each soil province, we considered that the WTH and SOH plots had a same pedogenesis and vegetation type, and that there was no difference between pre-harvest and post-regeneration stand composition. Soil characteristics were, therefore, likely to be the main factor explaining the difference in C and N contents between harvesting treatments. The selected plots were revisited for tree measurements and soil sampling from 2011 to 2013, a period representing a median of 30 years after treatment (interquartile range: 2 years).

## 2.3. *Field Sampling*

The mineral soil was sampled quantitatively at depth intervals of approximately 15 cm, down to 100 cm when possible, with a bipartite 8 cm volumetric root auger at each of the four cardinal points (subsamples every 100 cm, or approximately 16 m). The precise sampling depth was recorded for each sample. When it was not possible to quantitatively sample a soil layer, we used a standard Edelman auger. We also attributed a soil drainage class according to provincial standards [27]. In all, 16 plots were classified as having excessive drainage (Class 1), 96 were well drained (Class 2), 57 were moderately well drained (Class 3), and 27 were imperfectly to poorly drained (Classes 4 and 5).

## 2.4. *Laboratory Analyses*

We collected a total of 3482 soil samples. All the samples were kept frozen (–18 °C) until preparation for laboratory analysis. The samples were air-dried for at least four days, after which the quantitative samples were weighed individually. The air-dried mineral soil subsamples were passed through a 2 mm mesh sieve. The fine fraction was analyzed for remaining humidity and the concentration of organic matter was assessed with the loss on ignition method. Then, the subsamples of the fine fraction were further ground to 250 µm to determine total C and N with the dry combustion method (LECO CR-412, LECO Corporation, St. Joseph, MI, USA). Since soils did not contain any carbonates, total C was considered as organic C. Particle size analysis was carried out on all 0–15 cm depth soil samples (hydrometer method [28]) and size classes from the Canadian Soil Classification Working Group [29] were used: sand (50–2000 µm); coarse silt (20–50 µm); fine silt (5–20 µm); clay (0–5 µm). The >2 mm fraction of the quantitative mineral soil samples was weighed to determine the proportion of coarser fragments (f) when present.
## 2.5. Computations

Total C and N contents in a given mineral soil layer were calculated following Equation (1)

$$Q_m = 0.1 \times [x] \times D_b \times E_e, \qquad (1)$$

where

$Q_m$ is the element content in the mineral soil (kg·ha⁻¹);

[x] is the element concentration (mg·kg⁻¹) in oven-dried samples;

$D_b$ is bulk density (g·cm⁻³); and

$E_e$ is the effective thickness of the layer (cm), i.e., the thickness of the soil layer without fragments (f).

Effective horizon thickness was calculated following Equation (2)

$$E_e = E \times (1 - f) \qquad (2)$$

where

$E_e$ is the effective thickness of the layer (cm);

E is the measured thickness of the layer (cm); f is the coarse fraction of the volumetric sample (>2 mm) (% / 100).

We used a quantitative relationship obtained from the measured Db and SOM concentrations of the quantitative mineral soil samples (n = 1368) to estimate the bulk density of the soil layers that had not been sampled quantitatively (n = 2114) for each great soil texture class (sand, loam, and clay), following Federer et al.'s methodology [30] (Supplementary Materials, Table S2 and Figure S2). We used the aqp R package v. 1.17 [31] to compile soil C and N concentrations and stocks, and C/N ratios at fixed depths of 0–15, 15–30, 30–45, and 45–60 cm (measured thickness). C and N mean concentrations and ratios were calculated for the whole depth (0–60 cm), and soil C and N pools in each layer were summed. All soil concentrations and pools are reported for oven-dried mass.

Since no contrasting changes or discontinuities in soil texture were observed down to 60 cm depth at any of our soil sampling points, we considered the soil texture data from the upper sampled layer (0–15 cm) to be representative of the whole 0–60 cm soil sample. The clay fraction of the top mineral layer of forest soils, including the upper B and C horizons, is quite similar throughout Quebec (r = 0.86, p < 0.001) [32].

We used the relationship obtained with the boundary-line approach by Feng et al. [22], that is, 0.78 times the percentage of fine soil particles (0–20 µm) to estimate the potential maximum MAOM–C concentrations in our soil samples. The factor 0.78 represents the upper envelope of the relationship between MAOM–C and the percentage of fine soil particles in various soils, so the calculated MAOM–C and –N represent the potential maximal C and N fraction that can be associated with fine particles. To estimate the potential MAOM–N concentrations, we simply divided the MAOM–C concentrations by the average C/N ratio in the 0–60 cm soil layer. Soil POM–C and POM–N were then estimated as follows:

$$[POM-x] = [x] - \min([x], [MAOM-x]) \qquad (3)$$

where [x] is observed soil C or N concentrations.

C and N fractions as percentage of total C and N concentrations were also computed. It should be noted that this approach calculates MAOM–C and –N as the potential maximal C and N fractions that can be associated with fine particles, while POM–C and –N represent the minimum C and N fractions that are associated with particulate organic matter.

We calculated the relative loss in C and N soil stocks with WTH compared to SOH as the ratio of the stocks in each individual WHT plot over the average C and N stock values in the SOH plots within a given soil province, as follows:

$$\text{Relative } \times \text{ stock}=100 \times \text{ stock in WTH} / \text{ Average } \times \text{ stock in SOH} \qquad (4)$$

where x = soil C or N.
## 2.6. Statistical Analysis

We used a generalized linear mixed model to analyze soil C and N fractions and pools for the whole 0–60 cm mineral soil depth range. Harvesting treatments (SOH and WTH), soil provinces, proportion of soil fine particles, and their third-order interactions were considered fixed effects. A second set of analyses was performed using soil C and N pools with harvesting treatment and drainage class, and their interactions, as fixed effects. As some combinations of harvesting treatment and drainage class were missing for this second analysis, we combined them and selected a priori contrasts to test treatment combination effects. The homoscedasticity and normality of sample distributions were verified through residual plot analysis. We then plotted the standardized residuals of the models against all the independent variables to detect possible heterogeneity in their variances. If present, variance heterogeneity was corrected by allowing for different variances per strata. We selected the models for which the variance function structures had the lowest Akaike information criterion (AIC) scores. Adjusted (predicted) means were computed with the *emmeans* R package v. 1.4.2 [33], and statistical comparisons were made through specific a priori contrasts between harvesting treatments within soil provinces. The analyses were performed with the *nlme* R package v. 3.1-141 [34].

We estimated the proportion of fine particles needed and its 95% confidence interval (CI) so that relative soil C and N stocks with WTH attained 90% of the average stock for SOH within a given soil province. We used the *mod_alcc* function of the *soiltestcorr* R package [35] to assess the relationship between relative stocks and soil fine particle content. This function builds an arcsine-log calibration curve based on a standardized major axis; it is a bivariate regression model that assumes that both axes are random variables. To compare the C and N proportions for the sites in the two harvesting treatment groups, we applied chi-square tests for equality of proportions with Yate's continuity correction (*prop.test* function of the *stats* R package) [36].

# 3. Results and Discussion

## 3.1. Soil C and N Stocks with WTH Compared to SOH

Soil stocks ranged from 18 to 268 Mg ha⁻¹ for C, and from 0.4 to 13.4 Mg ha¹ for N for all the studied plots. For both C and N, some provinces only showed differences between SOH and WTH and when observed, losses were greater for WTH. Figure 1 shows the mean differences in C and N soil stocks between SOH and WTH, and average C and N stocks in SOH by soil province. Average differences in soil C stocks between SOH and WTH ranged from around zero for Soil Provinces B and D (*p* ≥ 0.407) to as much as 42 ± 11 Mg ha⁻¹ in favor of SOH in Soil Provinces C and E (*p* ≤ 0.004, Table S3). For Soil Provinces C and E, the degree of soil C loss due to WTH was similar, but the soil C stocks in SOH were quite different (Soil Province C: 148 ± 10 Mg C ha⁻¹ vs. E: 100 ± 15 C Mg ha⁻¹, *p* ≤ 0.007). A similar pattern was observed for soil N stocks; there was an average difference of 1.92 ± 0.46 Mg ha⁻¹ in favor of SOH in Soil Provinces C and E (*p* ≤ 0.006), although soil N stocks in SOH were distinct (Soil Province C: 6.23 ± 0.42 Mg ha⁻¹ vs. E: 4.81 ± 0.57 Mg ha⁻¹; *p* = 0.029). Meanwhile, in Soil Provinces B and D, no difference in soil N stocks was found between SOH and WTH (*p* ≥ 0.582). These results question the hypothesis that the soil C and N loss attributable to WTH is directly proportional to the initial soil C and N stocks.

A limitation of this study was the lack of pre-harvest soil data. However, the high replication of sites compensated for this issue, at least partially. In general, SOH leads to an increase in SOM pools as compared to non-harvested sites, while WTH does the opposite [37]. Therefore, the differences in C and N stocks between WTH and SOH sites that we observed in soils of the two Soil Provinces C and E can be attributed, at least in part, to an increase in these pools in SOH sites and a loss in WTH sites.
<figure><img src="image_6.png" /><figcaption>**Figure 1.** Differences in soil (A) C and (B) N stocks (0–60 cm soil depth) between stem-only (SOH) and whole-tree (WTH) harvesting as a function of soil C and N stocks for SOH by soil province. Adjusted means ± 95% confidence intervals. Soil provinces: B: Appalachians; C: Laurentians; D: Abitibi and James Bay Lowlands; E: Mistassini Highlands.</figcaption></figure>

3.2. *Relationship between Soil C and N Stocks with WTH and Fine Particle Content*

The four soil provinces were distinct in terms of soil particle content in the fine earth (Figure 2). Soil Provinces C and E were the richest in sand (mean ± SD: 69 ± 12% and 60 ± 11%, respectively) and the poorest in clay (9 ± 3% and 8 ± 5%, respectively). Conversely, Soil Provinces B and D were the poorest in sand (38 ± 11% and 43 ± 28%, respectively) and the richest in clay (25 ± 6% and 33 ± 26%, respectively). In terms of fine particles <20 µm, Soil Provinces B and D were the richest (49 ± 10% and 48 ± 29%, respectively) and Soil Provinces C and E the poorest (18 ± 6% and 25 ± 9%, respectively) (*p* < 0.001).

For Soil Provinces B and D, our analysis of the relationship of soil C and N stocks in the 0–60 cm mineral soil depth as a function of soil province, harvesting treatment, and soil fine particle content revealed no significant effect of harvesting treatment or of soil fine particle content, nor did their interaction have any effect (p ≥ 0.432). However, a different pattern was observed for Soil Provinces C and E, where a significant interaction between fine particle content and harvesting treatment was apparent for soil C, and to a lesser extent, for soil N (0.006 ≤ *p* ≤ 0.093). The differences in total soil C and N stocks between WTH and SOH, expressed as relative stocks for WTH over average stocks for SOH, reached 90% of the average stocks for SOH within each soil province when the soil fine particle content reached ≥ 30% (Figure 3). Its confidence intervals (CI95%) range from 26.1% to 34.4% for 90% soil C stocks and from 28.0% to 37.4% for 90% soil N stocks.
<figure><img src="image_7.png" /><figcaption>**Figure 2.** Soil texture diagram (Canadian soil texture triangle) with forest sampling points according to soil province. Soil provinces: B: Appalachians; C: Laurentians; D: Abitibi and James Bay Lowlands; E: Mistassini Highlands.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 3.** Relationship between soil fine particle content and (A) relative soil C stocks and (B) relative soil N stocks (0–60 cm soil depth) after whole-tree harvesting (WTH) (relative to average stem-only harvesting, SOH) by soil province. The vertical dashed line represents estimated fine particle content (and its 95% confidence intervals (CI)) required to reach 90% of the relative soil C and N stocks after SOH with the arcsine–log calibration curve; r is the square root of the proportion of variation explained by the curve. Soil provinces: B: Appalachians; C: Laurentians; D: Abitibi and James Bay Lowlands; E: Mistassini Highlands. Curves: relative soil C stocks = 100(sin(−0.776 + 0.596 ln(x))²); relative soil N stocks = 100(sin(−0.637 + 0.543 ln(x))²).</figcaption></figure>
Our results are in agreement with the meta-analysis by Wan et al. [38], who reported a general increase in mineral soil C losses with WTH as compared with SOH only for forest soils with a clay content of < 20%. This result is very similar to our threshold of 30% for fine particles <20 µm, which corresponds to a clay content of 12% to 17%. In similar boreal forests, Morris et al. [39] found that 20 years after various harvesting intensities, mineral soil C and N stock losses were greater only in coarser-textured sandy soils.

Although deadwood in SOH sites does not appear to be a significant source of SOM in coarse-textured soils [40], the slash left on the ground after SOH can be abundant, averaging 34 Mg ha⁻¹, compared to 5 Mg ha⁻¹ for WTH, for woody debris ≤7 cm in diameter [41]. Abundant slash may play a significant role in preserving POM in soils partly through its insulating and cooling properties [42–46]. The forest soil temperature regime and C balance can change more drastically after WTH compared to SOH [47–49], with consequences lasting several years. This is also true when the soil surface is exposed [50,51]. In addition, the natural forest regeneration established after harvesting may have played a role in SOM accumulation. Although the natural tree regeneration growth and abundance were not affected for up to 10 years after the harvesting treatments [26], lower tree growth, in the order of 15% in stem basal area increment, was observed for WTH compared to SOH in Soil Provinces C and E during the first 30 years [8]. This lower tree growth and the resulting reduced below and aboveground SOM input may have contributed to the differences in soil C and N stocks between WTH and SOH after 30 years.

### 3.3. *The Role of C and N Fractions*

The effect of soil fine particles on soil C and N stocks can be explained by the status of SOM, which can either be rather stable in the MAOM fraction, or more labile in the POM fraction [19,20]. The role of MAOM in soil C and N preservation is shown in Figure 4. In this figure, we plotted potential MAOM-C and MAOM-N concentrations estimated by the proportion of fine particles <20 µm in the soil against actual soil C and N concentrations according to soil provinces. The lines represent a 1:1 relationship; data above that line indicate that the maximal capacity to sequester C and N in the MAOM fraction is greater than the actual C and N measured. The figure indicates that almost all the sites in the two soil provinces with no loss of soil C and N attributable to WTH (i.e., Soil Provinces B and D) are above the 1:1 line.

This means that most of the soil C and N could potentially be locked in the MAOM fraction. In fact, Soil Provinces B and D have the capacity to sequester two to four times more C and N in their fine particles than their actual current concentrations.

That is not the case for Soil Provinces C and E, where losses of soil C and N attributable to WTH compared to SOH were observed. In these soils, a large part of the C and N could not be sequestered by the MAOM fraction, and therefore, remaining as soil POM, which, theoretically, is more prone to decomposition. In Soil Province C (n = 69), 51% of the sites did not have the potential to sequester all the C and N in their fine particle fraction (i.e., points falling below the 1:1 line). In Soil Province E, 27% of the sites did not have sufficient sequestration potential.

Obviously, it may be unrealistic to believe that the potential maximal MAOM fraction was attained in these soils because the 0.78 factor converting fine soil particle content to maximum MAOM-C fraction may not apply to all of them. However, applying the exact factor would not change the general pattern observed among soil provinces; that is, Soil Provinces B and D had a higher capacity to sequester SOM in a stable form than Soil Provinces C and E because of their higher soil fine particle content.

Since there is no theoretical limit to the accumulation of (more labile) POM in soils [52,53], the potential for soil C and N loss should be influenced more directly by the amount of POM rather than by total SOM. We examine this hypothesis in Figure 5, which shows the relationships between the estimated minimum POM–C and –N fractions and their corresponding actual C and N stocks according to soil provinces and harvesting intensity. Minimum POM–C and –N fractions increased with their corresponding total stocks only
in Soil Provinces C and E, where there were C and N losses with WTH compared to SOH. Meanwhile, there was no increase in minimum POM–C and –N fractions with their corresponding total stocks in Soil Provinces B and D, where no such losses in C and N occurred when WTH was compared to SOH.

<figure><img src="image_9.png" /><figcaption>**Figure 4.** Potential mineral-associated organic matter (MAOM) (A) C concentrations and (B) N concentrations as a function of measured C and N concentrations in the soils, according to forest harvesting treatment and soil province (n = 196). Lines represent a 1:1 relationship. Harvesting treatments: whole-tree harvesting (WTH); stem-only harvesting (SOH). Soil provinces: B: Appalachians; C: Laurentians; D: Abitibi and James Bay Lowlands; E: Mistassini Highlands.</figcaption></figure>

Therefore, it appears that actual soil C and N losses attributable to WTH are not entirely associated with their initial total stocks in the soil, but rather that they are more related to the abundance of soil POM, at least in these boreal forest soils. This conclusion is in agreement with what Achat et al. [3] report in their meta-analysis; they state that soil C losses in the upper mineral layers after WTH increased with increasing initial soil C stocks; this is probably because a higher fraction of the C stocks is held as POM in soils with higher soil C stocks, as we found for Soil Provinces C and E.

### 3.4. Relationship between Soil C and N Stocks with WTH and Drainage Conditions

Soil C and N stocks differed significantly between SOH and WTH for imperfectly or poorly drained soils: –50% for both C and N stocks (drainage class 4–5, *p* = 0.043, Figure 6). Soil N stocks also differed between SOH and WTH in moderately well-drained soils (–30%, drainage class 3, *p* = 0.001). We first assumed that the poorer the drainage, the higher the soil fine particle content. However, our initial analyses revealed the same trends even when the data were split by soil fine particle content (<30% and ≥30% of particles <20 µm). In addition, we found no relationship between drainage class and soil fine particle <20 µm content within the soil provinces (Chi² ≤ 2.92, *p* ≥ 0.232), excluding a possible confounding factor. Nevertheless, our results are based on relatively few observations (n = 10 for the WTH/drainage class 4–5 combination, n = 17 for the SOH/drainage class 4–5 combination), so these findings must be taken with caution.
<figure><img src="image_10.png" /><figcaption>**Figure 5.** Estimated minimum soil (A) C fractions and (B) N fractions in the form of particulate organic matter (POM) as a function of measured soil C and N stocks for whole-tree harvesting (WTH) and stem-only harvesting (SOH) in Soil Provinces B and D (left panels) and C and E (right panels). The shaded area around the regression lines represents ± 1 SE. Soil Provinces B and D: Appalachians, and Abitibi and James Bay Lowlands; C and E: Laurentians and Mistassini Highlands.</figcaption></figure>

The lower C and N stocks we observed in the WTH plots associated with poor soil drainage conditions suggest that SOM might have been higher in the POM fraction in those soils due to reduced SOM decomposition [54]. Indeed, slower-draining soils usually accumulate more C than well-drained soils for similar soil textures [12,55]. In well-drained soils, fine particles are a major factor of SOM stabilization. In poorly drained soils, on the other hand, cold temperatures, the lack of oxygen and substrate quality interact to stabilize SOM [23,24]. Therefore, as this SOM is not protected by mineral particles and other ligands, poorly drained soils may be prone to C and N losses in case of a major soil surface disturbance, such as WTH. The persistence of SOM in soils is ultimately controlled by interacting microscale processes shaped by the moisture regime [56]. More research is needed in order to better understand the impact of additional biomass harvesting on soil C and N stability in poorly drained soils.

### 3.5. *Management and Policy Implications*

Our analyses indicate that the differences in total soil C and N between SOH and WTH were not directly proportional to their initial stocks, but were more related to soil fine particle content and associated SOM. These two factors made a clear differentiation between the soil provinces having preserved, or lost, their soil C and N after WTH. Therefore, soil C and N sensitivity to loss after additional biomass harvesting appears to be closely related to the amount of fine mineral particles the soil contains, as stated in our second hypothesis.
Finally, differences in soil C and N stocks between SOH and WTH also appear to be associated with a particular soil drainage class, namely imperfectly or poorly drained soils. However, more observations are needed on the impact of WTH on soils in poor drainage conditions before making more general recommendations concerning additional biomass harvesting.

<figure><img src="image_11.png" /><figcaption>**Figure 6.** Relationships between (A) soil C and (B) soil N stocks (0–60 cm) and drainage class according to forest harvesting treatment (WTH: whole-tree harvesting; SOH: stem-only harvesting) after 30 years. * indicates a significant difference between harvesting treatments at probability *p* ≤ 0.05 and ***p* ≤ 0.001.</figcaption></figure>

Currently the forest sustainability policy for upland forests in Quebec relies only on soil acidification criteria to classify soil sensitivity to additional biomass harvesting. Interestingly, the soil provinces where critical loads of acidity are exceeded [57] are the same soil provinces we found to be the most sensitive to C and N losses after WTH compared to SOH. Existing policy targets forest ecological types with certain basic site attributes (e.g., coarse-textured, thin or organic soils). No specific guidelines for soil C and N preservation exist. The analyses presented in this paper allow us to propose some guidelines that may help fill this void and to preserve soil C and N stocks from the impacts of additional biomass harvesting in boreal forests in the long term. We found that a minimum threshold of 30% fine particles <20 µm ensures that WTH will have a similar impact on soil C and N to SOH in the long term. Thirty percent of fine particles <20 µm corresponds to a threshold of 45% silt + clay (Supplementary Materials, Figure S3). Limiting additional biomass harvesting to sites where forest soils contain ≥45% silt + clay and have a drainage class <4 (i.e., upland soils) would help the forest sector preserve soil C and N from potential losses. These guidelines could be applied readily as soil drainage and texture maps are now available for the whole forest territory in Quebec [32,58,59]. It is, therefore, possible to determine soil sensitivity to C and N losses after additional biomass harvesting at the scale of the forest stand in Quebec's boreal forest.
The extrapolation of these results to temperate broad-leaved forests is not straightforward. In general, coniferous forests store more SOM than do broad-leaved forests, and have a higher proportion of soil C in the POM fraction [7]. In coniferous forests, the MAOM-C fraction, on the other hand, is comparable to that of broad-leaved forests, so in broad-leaved forests, soil C and N should be less impacted by forest biomass harvesting [60]. This hypothesis remains largely undocumented, however.

**4. Conclusions**

In this study, we hoped to find soil properties that could explain the higher losses of soil C and N at sites subjected to WTH compared to sites subjected to SOH 30 years before. We found that lower soil C and N stocks with WTH than with SOH were not directly associated with corresponding initial stocks, but to the abundance of soil fine particles. C and N stocks were maintained after WTH in soils with at least 30% of fine particles <20 µm. A low fine particle <20 µm content leads to the saturation of the MAOM fraction and a greater abundance of POM in the soil, and this could explain why soil C and N were lower in WTH soils compared to SOH soils in the Laurentians and Mistassini Highlands soil provinces. Conversely, a higher fine particle <20 µm content increases the ability of the MAOM fraction to store C and N; this explains why there was no difference in soil C and N between WTH and SOH in the Appalachians and Abitibi and James Bay Lowlands provinces. Furthermore, imperfectly and poorly drained soils appeared to be sensitive to additional biomass harvesting regardless of soil texture, but this finding warrants further investigation of the mechanisms at play. In the context of the fight against climate change, it is important to stress that soil C and N reserves should be conserved, and if possible even enhanced, by forest operations in order to reduce greenhouse gases in the atmosphere.

**Supplementary Materials:** The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/soilsystems7020039/s1, Table S1. Characteristics of the four soil provinces in Quebec included in the study. Table S2. Parameter estimates for the relationship between soil bulk density (Db) and organic matter (OM) concentrations according to broad soil texture groups. Table S3. Total carbon (C) and total nitrogen (N) stocks in the first 60 cm of soil (forest floor excluded) as a function of harvesting treatment (stem-only harvesting (SOH) and whole-tree harvesting (WTH)) in four soil provinces in Quebec: B = Appalachians; C = Laurentians; D = Abitibi Lowlands; E = Mistassini Highlands. Data presented are model-adjusted means ± SE. Figure S1. Map of sampling plot numbers and locations. Colored areas show the different soil provinces in Quebec (Canada). Green lines delimit vegetation zones according to bioclimatic domains: SM-NH: sugar maple northern hardwoods; BF-YB: balsam fir-yellow birch domain; BF-WB: balsam fir-white birch domain; BS-FM: black spruce-feather moss domain. Figure S2. Relationship between soil bulk density (Db) and organic matter (OM) concentrations according to broad soil texture groups. Lines show values predicted by the model with a 95%CI. Figure S3. Relationship between soil fine particle content (silt + clay) and the relative average C and N stocks after whole-tree harvesting (WTH) (relative to stem-only harvesting, SOH) by soil province. The dashed line in the arcsine-log calibration curve represents the silt + clay content (95% confidence intervals (CI)) required to reach 90% of the soil C and N stocks after SOH [35]. Curves: relative C stocks = 100(sin(-1.384 + 0.698 ln(x))²); relative N stocks = 100(sin(-1.191 + 0.636 ln(x))²).

**Author Contributions:** Conceptualization, R.O., N.K. and I.B.; methodology, R.O.; validation, N.K. and I.B.; formal analysis, R.O.; writing—original draft preparation, R.O.; writing—review and editing, N.K. and I.B.; funding acquisition, R.O., N.K. and I.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Ministère des Ressources naturelles et des Forêts du Québec (project No. 142332122), by the Quebec 2030 Plan for a Green Economy, and by Le STUDIUM, Loire Valley Institute for Advanced Studies, Orléans & Tours, France, for a research fellowship to R.O. at INRAE.
**Data Availability Statement:** The soil samples for further analyses and data presented in this study are available on request from the corresponding author; data are also available at https://doi.org/10.5281/zenodo.7654866.

**Acknowledgments:** The authors gratefully acknowledge the research fellowship support by Le STUDIUM, Loire Valley Institute for Advanced Studies. We also thank Vicki Moore for editing the English in the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; and in the writing of the manuscript.

## References

1. IPCC. *Climate Change 2022: Impacts, Adaptation and Vulnerability*. Contribution of Working Group II to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change; IPCC: Cambridge, UK; New York, NY, USA, 2022; p. 3056.

2. Calvin, K.; Cowie, A.; Berndes, G.; Arneth, A.; Cherubini, F.; Portugal-Pereira, J.; Grassi, G.; House, J.; Johnson, F.X.; Popp, A.; et al. Bioenergy for climate change mitigation: Scale and sustainability. *GCB Bioenergy* **2021**, *13*, 1346–1371. [CrossRef]

3. Achat, D.L.; Deleuze, C.; Landmann, G.; Pousse, N.; Ranger, J.; Augusto, L. Quantifying consequences of removing harvesting residues on forest soils and tree growth—A meta-analysis. *For. Ecol. Manag.* **2015**, *348*, 124–141. [CrossRef]

4. Clarke, N.; Kjær, L.P.; Janne Kjønaas, O.; Bárcena, T.G.; Vesterdal, L.; Stupak, I.; Finér, L.; Jacobson, S.; Armolaitis, K.; Lazdina, D.; et al. Effects of intensive biomass harvesting on forest soils in the Nordic countries and the UK: A meta-analysis. *For. Ecol. Manag.* **2021**, *482*, 118877. [CrossRef]

5. James, J.; Page-Dumroese, D.; Busse, M.; Palik, B.; Zhang, J.; Eaton, B.; Slesak, R.; Tirocke, J.; Kwon, H. Effects of forest harvesting and biomass removal on soil carbon and nitrogen: Two complementary meta-analyses. *For. Ecol. Manag.* **2021**, *485*, 118935. [CrossRef]

6. Lal, R. Forest soils and carbon sequestration. *For. Ecol. Manag.* **2005**, *220*, 242–258. [CrossRef]

7. Lugato, E.; Lavallee, J.M.; Haddix, M.L.; Panagos, P.; Cotrufo, M.F. Different climate sensitivity of particulate and mineral-associated soil organic matter. *Nat. Geosci.* **2021**, *14*, 295–300. [CrossRef]

8. Ouimet, R.; Duchesne, L.; Tremblay, S. Long-term soil fertility and site productivity in stem-only and whole-tree harvested stands in boreal forest of Quebec (Canada). *Forests* **2021**, *12*, 583. [CrossRef]

9. James, J.; Harrison, R. The effect of harvest on forest soil carbon: A meta-analysis. *Forests* **2016**, *7*, 308. [CrossRef]

10. Stupak, I.; Lattimore, B.; Titus, B.D.; Tattersall Smith, C. Criteria and indicators for sustainable forest fuel production and harvesting: A review of current standards for sustainable forest management. *Biomass Bioenergy* **2011**, *35*, 3287–3308. [CrossRef]

11. Titus, B.D.; Brown, K.; Helmisaari, H.-S.; Vanguelova, E.; Stupak, I.; Evans, A.M.; Clarke, N.; Guidi, C.; Bruckman, V.J.; Varnagiryte-Kabanskiene, I.; et al. Sustainable forest biomass: A review of current residue harvesting guidelines. *Energy Sustain. Soc.* **2021**, *11*, 10. [CrossRef]

12. Paré, D.; Banville, J.L.; Garneau, M.; Bergeron, Y. Soil Carbon Stocks and Soil Carbon Quality in the Upland Portion of a Boreal Landscape, James Bay, Quebec. *Ecosystems* **2011**, *14*, 533–546. [CrossRef]

13. De Vos, B.; Cools, N.; Ilvesniemi, H.; Vesterdal, L.; Vanguelova, E.; Carnicelli, S. Benchmark values for forest soil carbon stocks in Europe: Results from a large scale forest soil survey. *Geoderma* **2015**, *251*, 33–46. [CrossRef]

14. Kögel Knabner, I.; Amelung, W. Soil organic matter in major pedogenic soil groups. *Geoderma* **2021**, *384*, 114785. [CrossRef]

15. Jandl, R.; Lindner, M.; Vesterdal, L.; Bauwens, B.; Baritz, R.; Hagedorn, F.; Johnson, D.W.; Minkkinen, K.; Byrne, K.A. How strongly can forest management influence soil carbon sequestration? *Geoderma* **2007**, *137*, 253–268. [CrossRef]

16. Mayer, M.; Prescott, C.E.; Abaker, W.E.A.; Augusto, L.; Cécillon, L.; Ferreira, G.W.D.; James, J.; Jandl, R.; Katzensteiner, K.; Laclau, J.; et al. Tamm Review: Influence of forest management activities on soil organic carbon stocks: A knowledge synthesis. *For. Ecol. Manag.* **2020**, *466*, 118127. [CrossRef]

17. Baldock, J.A.; Broos, K. Soil organic matter. In *Handbook of Soil Sciences*; Huang, P.M., Li, Y., Sumner, M.E., Eds.; CRC Press: Boca Raton, FL, USA, 2011; Volume 1, pp. 1–37.

18. Cotrufo, M.F.; Lavallee, J.M. Chapter One—Soil organic matter formation, persistence, and functioning: A synthesis of current understanding to inform its conservation and regeneration. In *Advances in Agronomy*; Sparks, D.L., Ed.; Academic Press: Cambridge, MA, USA, 2022; Volume 172, pp. 1–66.

19. Lavallee, J.M.; Soong, J.L.; Cotrufo, M.F. Conceptualizing soil organic matter into particulate and mineral-associated forms to address global change in the 21st century. *Glob. Change Biol.* **2020**, *26*, 261–273. [CrossRef]

20. Sokol, N.W.; Sanderman, J.; Bradford, M.A. Pathways of mineral-associated soil organic matter formation: Integrating the role of plant carbon source, chemistry, and point of entry. *Glob. Change Biol.* **2019**, *25*, 12–24. [CrossRef]

21. Hassink, J. The capacity of soils to preserve organic C and N by their association with clay and silt particles. *Plant Soil* **1997**, *191*, 77–87. [CrossRef]

22. Feng, W.; Plante, A.F.; Six, J. Improving estimates of maximal organic carbon stabilization by fine soil particles. *Biogeochemistry* **2013**, *112*, 81–93. [CrossRef]

23. Hobbie, S.E.; Schimel, J.P.; Trumbore, S.E.; Randerson, J.R. Controls over carbon storage and turnover in high-latitude soils. *Glob. Change Biol.* **2000**, *6*, 196–210. [CrossRef]
24. Freeman, C.; Ostle, N.; Kang, H. An enzymic 'latch' on a global carbon store. *Nature* **2001**, *409*, 149. [CrossRef] [PubMed]

25. Lamontagne, L.; Nolin, M.C. *Cadre Pédologique de Référence pour la Corrélation des Sols*; Équipe Pédologique du Québec, Centre de Recherche et de Développement sur les Sols et les Grandes Cultures, Agriculture et Agroalimentaire Canada: Sainte-Foy, QC, Canada, 1997; p. 69.

26. Pothier, D. Évolution de la régénération après la coupe de peuplements récoltés selon différents procédés d'exploitation. *For. Chron.* **1996**, *72*, 519–527. [CrossRef]

27. Saucier, J.P.; Berger, J.P.; D'Avignon, H.; Racine, P. *Le point D'observation Écologique: Normes Techniques*; Service des Inventaires Forestiers, Direction de la Gestion des Stocks Forestiers, Ministère des Ressources Naturelles du Québec: Québec, QC, Canada, 1994; p. 116.

28. Bouyoucos, G.J. Hydrometer method improved for making particle size analyses of soils. *Agron. J.* **1962**, *54*, 464–465. [CrossRef]

29. Soil Classification Working Group. *The Canadian System of Soil Classification*, 3rd ed.; Agriculture and Agri-Food Canada: Ottawa, ON, Canada, 1998; p. 187.

30. Federer, C.A.; Turcotte, D.E.; Smith, C.T. The organic fraction-bulk density relationship and the expression of nutrient content in forest soils. *Can. J. For. Res.* **1993**, *23*, 1026–1032. [CrossRef]

31. Beaudette, D.E.; Roudier, P.; O'Geen, A.T. Algorithms for quantitative pedology: A toolkit for soil scientists. *Comput. Geosci.* **2013**, *52*, 258–268. [CrossRef]

32. Duchesne, L.; Ouimet, R. Digital mapping of soil texture in ecoforest polygons in Quebec, Canada. *PeerJ* **2021**, *9*, e11685. [CrossRef]

33. Lenth, R. Emmeans: Estimated Marginal Means, Aka Least-Squares Means; R Package Version 1.4.2; 2019. Available online: https://github.com/rvlenth/emmeans (accessed on 1 September 2018).

34. Pinheiro, J.; Bates, D.; DebRoy, S.; Sarkar, D.; R Core team. nlme: Linear and Nonlinear Mixed Effects Models. R Package Version 3.1-137. 2018. Available online: http://cran.r-project.org/package=nlme (accessed on 12 June 2018).

35. Correndo, A.A.; Salvagiotti, F.; García, F.O.; Gutiérrez, Boem, F.H. A modification of the arcsine-log calibration curve for analysing soil test value–relative yield relationships. *Crop Pasture Sci.* **2017**, *68*, 297–304. [CrossRef]

36. R Core Team. *R: A Language and Environment for Statistical Computing*; R Foundation for Statistical Computing: Vienna, Austria, 2020.

37. Achat, D.L.; Fortin, M.; Landmann, G.; Ringeval, B.; Augusto, L. Forest soil carbon is threatened by intensive biomass harvesting. *Sci. Rep.* **2015**, *5*, 15991. [CrossRef]

38. Wan, X.; Xiao, L.; Vadeboncoeur, M.A.; Johnson, C.E.; Huang, Z. Response of mineral soil carbon storage to harvest residue retention depends on soil texture: A meta-analysis. *For. Ecol. Manag.* **2018**, *408*, 9–15. [CrossRef]

39. Morris, D.M.; Hazlett, P.W.; Fleming, R.L.; Kwiaton, M.M.; Hawdon, L.A.; Leblanc, J.D.; Primavera, M.J.; Weldon, T.P. Effects of biomass removal levels on soil carbon and nutrient reserves in conifer-dominated, coarse-textured sites in Northern Ontario: 20-year results. *Soil Sci. Soc. Am. J.* **2019**, *83*, S116–S132. [CrossRef]

40. Senez Gagnon, F.; Thiffault, E.; Paré, D.; Achim, A.; Bergeron, Y. Dynamics of detrital carbon pools following harvesting of a humid eastern Canadian balsam fir boreal forest. *For. Ecol. Manag.* **2018**, *430*, 33–42. [CrossRef]

41. Bessaad, A.; Bilger, I.; Korboulewsky, N. Assessing Biomass Removal and Woody Debris in Whole-Tree Harvesting System: Are the Recommended Levels of Residues Ensured? *Forests* **2021**, *12*, 807. [CrossRef]

42. Yole, D. Effect of Alternative site Preparation Treatments on Soil Chemistry, Physical Properties, Climate and Seedling Growth on a Mesic Site in the Northern Interior of British Columbia. Ph.D. Thesis, University of British Columbia, Vancouver, BC, Canada, 1996.

43. Smith, C.T.; Lowe, A.T.; Skinner, M.F.; Beets, P.N.; Schoenholtz, S.H.; Fang, S. Response of radiata pine forests to residue management and fertilisation across a fertility gradient in New Zealand. *For. Ecol. Manag.* **2000**, *138*, 203–223. [CrossRef]

44. Waldrop, M.P.; McColl, J.G.; Powers, R.F. Effects of Forest Postharvest Management Practices on Enzyme Activities in Decomposing Litter. *Soil Sci. Soc. Am. J.* **2003**, *67*, 1250–1256. [CrossRef]

45. Moroni, M.T.; Carter, P.Q.; Ryan, D.A.J. Harvesting and slash piling affects soil respiration, soil temperature, and soil moisture regimes in Newfoundland boreal forests. *Can. J. Soil Sci.* **2009**, *89*, 343–355. [CrossRef]

46. Augusto, L.; Achat, D.L.; Bakker, M.R.; Bernier, F.; Bert, D.; Danjon, F.; Khlifa, R.; Meredieu, C.; Trichet, P. Biomass and nutrients in tree root systems–sustainable harvesting of an intensively managed Pinus pinaster (Ait.) planted forest. *GCB Bioenergy* **2015**, *7*, 231–243. [CrossRef]

47. Mahendrappa, M.K.; Pitt, C.M.; Kingston, D.G.O.; Morehouse, T. Environmental impacts of harvesting white spruce on Prince Edward Island. *Biomass Bioenergy* **2006**, *30*, 363–369. [CrossRef]

48. Mahendrappa, M.K.; Kingston, D.G.O. Intensive harvesting impacts on soil temperature and solution chemistry in the maritimes region of Canada. *N. Z. J. For. Sci.* **1994**, *24*, 402–414.

49. Lavoie, M.; Kellman, L.; Risk, D. The effects of clear-cutting on soil CO2, CH4, and N2O flux, storage and concentration in two Atlantic temperate forests in Nova Scotia, Canada. *For. Ecol. Manag.* **2013**, *304*, 355–369. [CrossRef]

50. Kubin, E.; Kemppainen, L. Effect of soil preparation of boreal spruce forest on air and soil temperature conditions in forest regeneration areas. *Acta For. Fenn.* **1994**, *244*, 56. [CrossRef]

51. Jiménez Esquilín, A.E.; Stromberger, M.E.; Shepperd, W.D. Soil Scarification and Wildfire Interactions and Effects on Microbial Communities and Carbon. *Soil Sci. Soc. Am. J.* **2008**, *72*, 111–118. [CrossRef]
52. Barré, P.; Angers, D.A.; Basile-Doelsch, I.; Bispo, A.; Cécillon, L.; Chenu, C.; Chevallier, T.; Derrien, D.; Eglin, T.K.; Pellerin, S. *Ideas and perspectives: Can we use the soil carbon saturation deficit to quantitatively assess the soil carbon storage potential, or should we explore other strategies?* *Biogeosci.* **2017**, *2017*, 1–12. [CrossRef]

53. Cotrufo, M.F.; Ranalli, M.G.; Haddix, M.L.; Six, J.; Lugato, E. Soil carbon storage informed by particulate and mineral-associated organic matter. *Nat. Geosci.* **2019**, *12*, 989–994. [CrossRef]

54. Keiluweit, M.; Wanzek, T.; Kleber, M.; Nico, P.; Fendorf, S. Anaerobic microsites have an unaccounted role in soil carbon stabilization. *Nat. Commun.* **2017**, *8*, 1771. [CrossRef]

55. Tan, Z.X.; Lal, R.; Smeck, N.E.; Calhoun, F.G. Relationships between surface soil organic carbon pool and site variables. *Geoderma* **2004**, *121*, 187–195. [CrossRef]

56. Schlüter, S.; Leuther, F.; Albrecht, L.; Hoeschen, C.; Kilian, R.; Surey, R.; Mikutta, R.; Kaiser, K.; Mueller, C.W.; Vogel, H.; et al. Microscale carbon distribution around pores and particulate organic matter varies with soil moisture regime. *Nat. Commun.* **2022**, *13*, 2098. [CrossRef]

57. Ouimet, R.; Arp, P.A.; Watmough, S.A.; Aherne, J.; DeMerchant, I. Determination and mapping critical loads of acidity and exceedances for upland forest soils in Eastern Canada. *Water Air Soil Pollut.* **2006**, *172*, 57–66. [CrossRef]

58. Gouvernement du Québec. Carte écoforestière à Jour. Données Québec. 2022. Available online: https://www.donneesquebec.ca/recherche/dataset/carte-ecoforestiere-avec-perturbations (accessed on 10 January 2023).

59. Sylvain, J.; Daniel; Anctil, F.; Thiffault, É. Using bias correction and ensemble modelling for predictive mapping and related uncertainty: A case study in digital soil mapping. *Geoderma* **2021**, *403*, 115153. [CrossRef]

60. Hume, A.M.; Chen, H.Y.H.; Taylor, A.R. Intensive forest harvesting increases susceptibility of northern forest soils to carbon, nitrogen and phosphorus loss. *J. Appl. Ecol.* **2018**, *55*, 246–255. [CrossRef]

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.