Soil sampling was carried out in the French Level II forest monitoring plots (RENECOFOR network). The 102 plots are distributed over the French territory excluding overseas territories, and cover a wide range of climates, soil types and tree species (Figure 1; Table 1). Soil sampling was first conducted between 1993 and 1995, and then repeated between 2007 and 2012. Within the 0.5 ha central zone of each plot, the individual forest floor horizons (litter (Ol), fragmented (Of), and humic (Oh) horizons) and three mineral soil layers (0–10 cm, 10–20 cm and 20–40 cm) were sampled separately on 5 subplots of 13.5 m × 13.5 m each, located at the center and in the four corners (Figure 2). For each subplot, one composite sample per horizon or layer was obtained from five individual sampling points selected among the 16 intersections of a systematic grid (4.5 × 4.5 m), to account for local site variability while avoiding disturbed areas (logging residues, skidder tracks) and proximity to living trees. For the forest floor, the three horizons were collected separately, using a 30 cm × 30 cm frame after discarding the living material (other than fine roots), the branches of a diameter greater than 0.5 cm and the fresh fruits. Three underlying mineral soil layers based on soil depth (0–10, 10–20, 20–40 cm) were also sampled separately by digging a small pit over an approximately 50 cm × 50 cm area. For each mineral layer, an undisturbed soil sample was first taken with a Kopecky cylinder (250 cm³, h = 5 cm, Ø = 8 cm) to measure bulk density, and the remaining soil was then sampled for chemical analyses. The same protocol was followed for both soil surveys, except that the subplots were moved by 1.5 m in a fixed direction per plot during the second campaign to avoid any possible soil disturbance resulting from the previous sampling. For a given plot, the soil was sampled during the same season for both surveys. A more comprehensive description of the sampling grid soil is provided in the guidebook of Ulrich et al. (2009). In addition to this nested sampling design, two soil pits were opened per plot during the first survey to characterize the soil horizons down to a depth of 1 meter and to estimate the volumetric coarse fragment content (see below).

All samples for chemical analysis were first oven-dried at 35°C for 24 h. For the forest floor, the dry to fresh mass ratio was determined on subsamples oven-dried at 105°C until constant weight, and the chemical analyses were performed after grinding. For the mineral soil layers, all samples were passed through a 2-mm sieve and those for fine earth bulk density measurements were further oven-dried at 105°C before being weighed.

All chemical analyses were performed according to French and international standards by one single reference laboratory (Laboratory of Soil Analyses, INRAe Arras, France), using the same methods for both campaigns. This laboratory is routinely involved in ring test analyses, and was selected to assess the uncertainty of SOC concentrations by multi-laboratory analysis in the frame of the BioSoil project (De Vos et al., 2015). Soil pH was determined in the mineral layers according to NF ISO 10390 (2004) using a glass electrode in a 1:5 (V/V) suspension of soil in deionized water (pH H₂O), and in a solution of 0.01 mol.L⁻¹ calcium chloride (pH CaCl₂). SOC content of the forest floor and 0–10 cm layers was determined by dry combustion after decarbonation according to the Dumas method (NF ISO 10694, 1995). For the 10–20 cm and 20–40 cm layers, SOC was determined by sulfochromic oxidation at a temperature of 30°C after decarbonation according to the Anne method (variant of the Walkley-Black method; NF X31-109, 1993). Total nitrogen (N_tot) of the forest floor and 0–10 cm layers was determined by dry combustion according to NF ISO 13878 (1998), while N_tot of 10–20 cm and 20–40 cm layers was measured according to the Kjeldahl method (NF ISO 11261, 1995). Total K, Ca, and Mg were extracted from the forest floor layers by hydrofluoric and perchloric acid according to NF ISO 14869-1 (2001), then measured by flame atomic absorption spectroscopy (FAAS) for Ca and Mg, and by flame atomic emission spectroscopy (FAES) for K.

Extractable phosphorus (P) in the mineral soil layers was quantified either by the Dyer method for samples with pH CaCl₂ < 6.5 (NF X31-160, 1993), or by the Joret-Hébert procedure (NF X31-161, 1993) for samples with a pH CaCl₂ > 6.5. The Dyer method involves a solution of citric acid monohydrate 20 g.L⁻¹ and a 1:5 ratio (w/v), while P is extracted with a solution of ammonium oxalate 0.1 mol.L⁻¹ using a 1/25 ratio (w/v) with the Joret-Hébert method. In both cases, P is determined by colorimetry. Exchangeable cations (K, Ca, Mg, Al) were extracted with a solution of barium chloride (0.1 mol.L⁻¹) in a 1:10 mass/volume ratio according to NF ISO 11260 (1994), then measured by FAAS. Exchangeable Al was measured either by titrimetry (1st campaign) or ICP (2nd campaign), on all samples with a pH CaCl₂ < 7.

Values below the Limit of Quantification (LOQ) were set to either LQ or LQ/2, the latter being used when the parameter was included in sums or ratios. As in practice exchangeable Al concentrations in soils with pH CaCl₂ > 6.5 were generally below the LOQ, the corresponding values were set to 0 for all those samples.

Carbonated soil samples were excluded from temporal change analysis of exchangeable Ca and Mg, since BaCl₂ extraction is known to partly dissolve the carbonates (Dohrmann and Kaufhold, 2009; Nel et al., 2023). Similarly, soil samples involving different analytical methods between campaigns regarding extractable P (Dyer vs. Joret-Hébert procedures) were not considered when analyzing changes in P.

Finally, the whole data set was systematically checked for possible outliers. One individual measurement was considered as outlier and the corresponding (subplot × layer) excluded from further analyses if the two following conditions were met simultaneously: (i) z-score of the individual measurement (calculated based on the 10 individual measurements of a plot × layer combination across the two campaigns) beyond ±2.58 standard deviations and (ii) subplot value beyond the ±2.58 residual standard deviations confidence interval associated with the regression line between the two campaigns.

Effective cation exchange capacity (ECEC) and base saturation (BS) were calculated as follows:

As mentioned above for Ca and Mg, carbonated soil samples were not considered when comparing ECEC and BS among both soil campaigns as extraction with BaCl₂ causes partial dissolution of carbonates, thereby giving rise to overestimation of these two quantities.

For the forest floor, the stocks of elements were calculated by multiplying their concentration by the dry mass of the layer, expressed per unit area for each plot at each sampling campaign. For the mineral layers, the stocks of SOC, N_tot, extractable P and exchangeable cations were calculated by multiplying the soil content of the element (C_x) by the soil bulk density (BD, in g.cm⁻³) and by the thickness of the layer (t, in cm), corrected for the fraction of coarse elements (CE) larger than 2 mm, estimated by visual observation of soil profiles, as follows:

As bulk density measurements of the two soil surveys were not strictly comparable due to a change in the sieving technique (mostly manual for the first soil survey vs mechanical for the second), the use of a campaign-specific BD would bias the comparison of stocks. We therefore compared the BD on a subset of 260 soil samples covering the diversity of soil types and the different soil layers that were sieved mechanically at both soil surveys, and found no significant differences between the first (1.031 g.cm⁻³) and second (1.027 g.cm⁻³) survey (p-value of the paired t-test = 0.71). Based on this, we assumed there was no change in BD and we used BD values of the first campaign (the most complete) as a reference for all stock calculations.

Temporal changes of soil chemical parameters between the sampling campaigns were assessed for each depth individually by linear mixed models accounting for the correlation structure of the data set. These models included the sampling campaign as a fixed effect and accounted for the inter-plot variability with a random factor:

where x is the soil property measure, a s is the parameter estimating the mean value of the soil property at each soil survey s, plot is the random factor accounting for the plot effect with mean value zero and variance σ plot 2 and ε is the residual term with mean value zero and variance σ ε 2 .

A change was considered as significant when the p value associated with the fixed effect was lower than 0.05.

Estimation of model parameters was performed using the method of Restricted Maximum Likelihood (REML). The linearity, homoscedasticity (Levene test), normality [skewness and kurtosis coefficients ranging between −2 and 2 to accept the normality of the distribution, Miller (1986)], and independence of residuals have been checked for all models. All chemical parameters (except for exchangeable Al) were not complying with Gaussian distribution. Data were therefore log-transformed before performing statistical tests.

When comparing large datasets, minor changes are sufficient to generate a significant effect (Lourel et al., 2011). To properly interpret the P levels of the fixed effect, we therefore calculated the effect size with the following indices: Eta² (Cohen, 1973), Cohen’s d (Cohen, 1988), Pearson’s correlation coefficient r (Rosenthal et al., 2000). Effect sizes were classified as small (Eta² < 4%, Cohen’s d < 0.40, r < 0.20), medium (Eta² 4–16%, Cohen’s d 0.40–0.87, r 0.20–0.40) or high (Eta² > 16%, Cohen’s d > 0.87, r > 0.40).

To study the relationships between the changes in the various chemical parameters, principal component analysis (PCA) and Spearman rank correlation were performed. PCA were performed using centered-reduced data at plot scale (individual sub-plot values minus plot mean calculated over the two soil surveys) so that inter-plot variability has been removed, allowing to better highlight the temporal changes.

Finally, for illustrative purposes, we separated the plots into three classes of trophic levels, based on the mean pH H₂O of the 0–10 cm mineral soil layer calculated over the two soil surveys: highly acidic soils with pH < 4.5, intermediate soils with pH between 4.5 and 5.5, and low-acid soils with pH > 5.5. A linear mixed model including the interaction between ‘campaign’ and ‘soil trophic level’ was tested first. In case of a significant (p < 0.05) interaction, a linear mixed model similar to Equation (4) was then run separately for each trophic level to test for the campaign effect.

All statistical analyses were performed using the R version 3.0.2 software (ade4 package, Dray and Dufour (2007) and nlme package, Pinheiro et al., 2013; R Core Team, 2013).

The temporal change of the main chemical parameters of the forest floor between the two sampling campaigns is presented in Table 2. We noted a significant 11.2% increase of SOC concentrations, while N_tot concentrations did not change significantly. Consequently, the C/N ratio increased sharply by 8.1% or 2.6 units. Total calcium (Ca_tot) and total magnesium (Mg_tot) concentrations increased significantly between the two campaigns with, respectively, 8.3 and 10.9% increase. Total potassium (K_tot) concentrations did not display any significant change. As total litter mass did not change significantly over time (data not shown), the changes in stocks reflected the changes in concentrations.

Aside the overall temporal evolution of soil parameters presented above for all soil types together, contrasted change patterns were found among the three considered different classes of trophic levels (Figure 4 Supplementary material S2). Highly acidic soils (pH H₂O < 4.5), corresponding to nearly 50% of the plots, exhibited indices of further acidification process with a significant decrease in soil pHs (1.9 and 2.9% decrease for pH H₂O and pH CaCl₂, respectively) and an increase in exchangeable Al concentrations (10.1% increase) especially pronounced in the 0–10 cm soil layer and tending to weaken with increasing depth. Nevertheless, aside these acidification indices and the general decrease of BS, these soils also exhibited increasing exchangeable Ca (2.7% increase), Mg (12.5% increase) and K (12.5% increase) concentrations over the whole profile, though not significantly for the former element. In contrast, low-acid soils (pH H₂O > 5.5) presented unchanged or increasing (pH H₂O in the 20–40 cm layer: 2.4% increase) pHs between the two campaigns, no significant change in BS, stable or decreasing (20–40 cm layer) exchangeable Al stocks and a general increase in exchangeable Ca, especially for the 0–40 cm layers (14.1% increase on average). Soils with initial pH H₂O ranging between 4.5 and 5.5 showed intermediate change patterns for these properties. Another noticeable difference among trophic levels is the particularly sharp decrease of total nitrogen in the 10–40 cm soil layers of the highly acidic soils (21.5% decrease on average) compared with the corresponding lower reduction for low-acidic soils (6.4% decrease on average, significant for the 20–40 cm layer only), while changes of N_tot stocks in the upper mineral layer and the forest floor were weaker and close for all trophic levels.

For organic carbon, SOC stocks increased significantly in the 0–10 cm mineral layer with similar relative amplitude for all soil types (14.3% increase on average). SOC relative change decreased as moving down in the profile to become significantly negative in the 20–40 cm layer of the highly acidic soils (4.2% decrease) and non-significant for the two other trophic levels. Regarding the forest floor, substantial increase of SOC stock (12.3% increase) was observed for highly acidic soils while the amplitude of the differences among campaigns decreased with increasing trophic level and was found to be non-significant for low-acidic soils. As a result of these patterns for N_tot and SOC, the C/N ratio rose significantly among campaigns in all cases, though with increasing differences from the forest floor (4.8% increase) to the lowest mineral soil layer (21.2% increase) for the more acidic soils and inversely for the lowest ones (forest floor: 14.9% increase, 20–40 cm layer: 9.0% increase). Considering the whole profile, the average temporal increase of the C/N ratio raised with soil acidity.

Significant changes were detected for many chemical parameters during t15-year monitoring period, illustrating the potential of such network to highlight even small changes in soil properties while the risk of detecting “false” trends was limited by accounting for effect size (Sullivan and Feinn, 2012). The correlations between temporal changes in soil parameters were in agreement with well established relationships between soil parameters (positive relationship between ECEC and SOC, pH correlated negatively with exchangeable Al and positively with exchangeable Ca and Mg), which supports the observed trends. The observed trends were also in agreement with the reported evolution of key potential drivers (see below).

The relative or absolute change strongly differed according to the parameter or the soil layer. As expected and also reported by other authors (e.g., Lawrence et al., 2015; Berger et al., 2016; Fraser et al., 2019), the change in the soil chemical properties tended to decrease with depth, with no or only limited effect in the 20–40 cm layer except for total N and the C/N ratio. However, when integrated over the whole soil profile (0–40 cm), the changes in some parameters including increase in C stocks, C/N ratio, ECEC as well as exchangeable K and Mg, were strong enough to be significant and associated to at least a medium size effect.

The mixed models together with the PCAs results showed contrasting evolution patterns between the top mineral 0–10 cm layer and the underlying 10–20 cm and 20–40 cm layers. For the topsoil, results from both analyses showed an increase in SOC, exchangeable cations (ECEC, base cations), acidity (lower pH, higher exchangeable Al) and total N contents. For the underlying layers, the dominant temporal change was an acidification trend (increase in exchangeable Al, decrease in pH) when considering the PCA results, while the mixed-models highlighted a decrease in N_tot, a moderate increase in all exchangeable cations (including Al) and no or little change in SOC. As the increase in ECEC could not be ascribed to SOC changes in these deep layers, the additional exchangeable base cations and Al should have replaced exchangeable sites occupied by a cation not accounted for in the calculation of ECEC. We hypothesize that this could be NH₄⁺, which could also explain the observed decrease of N_tot.

The evolution of soils may be ascribed to change in inputs and outputs of elements under solid, liquid and gaseous forms associated with a range of processes. In the following section we will show how the reported recent change in two major element fluxes, namely, increase in C sequestration (Jonard et al., 2017) and reduced inorganic N inputs (Waldner et al., 2014), contributed to the observed changes in soil chemistry within such a short time span.

Contrasting changes of soil properties as a function of the soil trophic level have been observed by several authors (Kirk et al., 2010; Vanguelova et al., 2010; Cools and De Vos, 2011; Jandl et al., 2012; Hazlett et al., 2020). For all these works, the soil resampling interval covered the period of strong reduction of N and S atmospheric deposition. The comparison of soil types mainly focused on the variations of soil pH over that period, showing generally a more pronounced increase of pH for the acidic soils (pH < 5.5) than for soils with intermediate (pH 5.5–7.0) or high (pH > 7.0) values for which smaller changes, corresponding even sometimes to a pH decrease, were found. Such results suggest that recovery from acidification following the decline in acidic deposition would primarily occur for acidic soils in which the dissolution of aluminosilicates is the main process involved in acid neutralization, than for soils in which cation exchange or carbonate dissolution is the dominant acid buffering mechanism. Yet, these findings differ from our observations of decreasing pH values for the highly acidic soils (especially in the upper mineral soil layer) and of pH stagnation or increase for the intermediate and low-acidic soils (Figure 4; Supplementary material S2). The reasons for the divergences between our results and those of the above works are not straightforward. Investigating changes in forest soil properties at the European level, Cools and De Vos (2011) hypothesized that the observed contrasting temporal trend of soil pH among trophic levels could be due to the spatial patterns of soil types and atmospheric deposition. Indeed, the soils with relatively high pH are better represented in the southern part of Europe and mainly rely on atmospheric deposition for the supply of base cations, which also decreased in parallel with N deposition; in contrast, more acidic soils are mainly located in northwest and northern Europe, and mostly rely on weathering as a source of base cations (Draaijers et al., 1997). Kirk et al. (2010) indicated that pH increase following a decline in acidic deposition is more likely for soils with pH lower than 5.5, in which the dissolution of aluminosilicate minerals is the main weathering mechanism for increasing soil pH, than for soils in which carbonate dissolution (pH ≥ 7) or cation exchange (pH between 5.5 and 7) is the dominating buffering process. In our study, the SOC increase in the top mineral soil associated with the pronounced rise of exchangeable Al for the highly acidic soils would favor complexation of this element with organic acids, which is accompanied with the liberation of protons and would lead to a decrease of pH in these acid soils (Nissinen et al., 1999). Moreover, the nitrification of NH₄⁺ postulated above as a potential explanation for N_tot depletion, which is especially pronounced for highly acidic soils, would also give rise to a further decrease of soil pH in these contexts. Besides, the influence of other factors superimposing on that of the soil trophic level on the evolution of the pH, such as notably the effect of the species composition (Ahrends et al., 2022; Clesse et al., 2022), might also at least partly explain the differences between our results and findings of previous studies. Nevertheless, within the RENECOFOR plots, deciduous species are more frequently associated with eutrophic soils while conifers are mainly found on intermediate and acidic soils, with only few exceptions to this general pattern (Table 1). Therefore, the effect of species can hardly be distinguished from that of the soil trophic level. In addition, as already highlighted by Jonard et al. (2017) when analyzing SOC changes, coniferous stands are on average 30 years younger than deciduous ones within the network, so the effect of age also overlaps with that of species. In other respects, the larger amplitude of the increase of SOC stocks observed in the forest floor of the highly acidic soils compared with the two other soil types (Supplementary material S2) would arise from a reduced decomposition of organic matter in acidic contexts, which is corroborated by the much larger SOC pools at the more acid sites.

While many significant changes were detected, some of them were of very limited amplitude with therefore few expected ecological consequences, except in contexts of nutrient limitation or imbalance for which the increase of the availability of the limiting nutrient or the alleviation of the disequilibrium might have substantial effects on ecosystem functioning. In contrast, the large changes observed in SOC and N_tot stocks are likely to have general consequences. The rise of SOC between 0 and 10 cm resulted in an increase in buffering capacity of soil. Yet, in the current context of global change, such an accumulation of SOC in the surface soil raises questions with regard to its biogeochemical stability. In interface position with the atmosphere and vegetation, the surface soil is more exposed to climate changes and land use changes than the deep mineral layers (Lozano-García et al., 2017). In the French forest soils, the accumulated SOC might be more vulnerable to climate and/or anthropogenic perturbations and be more quickly destabilized than the deep C pool (Soucémarianadin et al., 2018, 2019). However to date, the relative vulnerability of the surface vs deep C stocks to environmental changes remains an open question (Schmidt et al., 2011). Besides, the N decrease in the deep soil layers, particularly marked in the highly acidic soils, raises questions as to the consequences of such changes for tree nutrition. Nitrogen, considered until recently as being in excess in areas having been subject to elevated N deposition (Jonard et al., 2015), might now become limiting for tree growth (Etzold et al., 2020). Nevertheless, while the influence of N exceedance on forests has been widely studied, the response of the ecosystems to N deposition decline is still rather poorly documented and further research should be conducted in that direction (Schmitz et al., 2019).