Subsoils, but not toeslopes, store millennia-old PyC in a gently sloping catchment under temperate climate after centuries of cultivation

Pyrogenic carbon (PyC) is the carbonaceous solid residue of incomplete combustion of biomass. It is a continuum of mostly condensed and aromatic molecules. PyC persists for longer in soils relative to non-PyC organic carbon. However, estimates of PyC residence time vary greatly. The time and spatial scales investigated are not always adapted to the long-residence time and vertical and lateral mobility of PyC in the soil profile and the landscape. In addition, agricultural land-use and shallow slopes are under-represented in the PyC literature.

We measured the concentrations and stocks of PyC down to 60 cm along three toposequences in a small agricultural catchment with shallow slopes and homogeneous soil parent material in the west of France. We used two methods (chemo-thermal oxidation – CTO and hydropyrolysis – HyPy) of PyC quantification that cover the intermediate to highly condensed part of the PyC continuum, and also measured the radiocarbon values in both total soil organic carbon (SOC) and the PyC fraction. There was likely little persistent PyC inputs to the catchment in the last 150 years which gave us access to the resultant, long term PyC distribution in the landscape. In particular, we aimed to investigate whether the vertical and horizontal distribution of PyC were similar or differed from SOC and whether they were affected by the soil types along the slope.

Topographic position was not the main driver of PyC stocks in this landscape. The stock of PyC_CTO to 60 cm depth averaged 2.5 ± 0.22 t ha⁻¹ across topographic positions. PyC stocks were the highest in a Solimovic Cambisol at the toeslope (3.3 ± 0.26 t ha⁻¹), likely formed following changes in erosion dynamics with land-use. Contrary to previous reports, erosion redistributed already aged PyC without enrichment or depletion. PyC_HyPy concentrations in the topsoil decreased from upslope (median = 1.6, IQR = 0.22 g C kg⁻¹ soil) to downslope positions (median = 1.10, IQR = 0.40 g C kg⁻¹ soil), which we tentatively attribute to PyC_HyPy leaching following the destabilisation of mineral associations with iron oxides in the water-table affected portion of the transects. The subsoil (30–60 cm) represented between 37 % and 51 % of the PyC_CTO stock. PyC_HyPy proportion in SOC increased with depth and reached an average of 11 ± 3.3 % at 50–60 cm depth. PyC_HyPy had an uncalibrated radiocarbon age of 2520 to 9600 years BP at this depth, significantly older than bulk SOC at the same depth and than PyC_HyPy at 0–10 cm (1530 to 2630 years BP). These results confirm the long persistence of PyC in soils and point to a slow advection of PyC towards the soil depth under the pedoclimatic conditions of our study area.

Future studies should assess whether erosion modalities and age and quality of PyC affect its fate during erosion events. Identifying the proportion of PyC produced which is quickly transported away from the watershed and that which remains and is stabilised in soils for millennia after a fire is an important knowledge gap that still needs to be investigated to close the terrestrial PyC budget.

1 Introduction

Pyrogenic carbon (PyC) is the carbonaceous solid residue of incomplete combustion or pyrolysis of biomass and fossil fuels. It is present in all compartments of the earth system (Santín et al., 2016; Jones et al., 2019; Wagner et al., 2018). PyC is a continuum of molecules in which the degree of aromaticity and condensation is positively correlated with production temperature (Keiluweit et al., 2010; Wiedemeier et al., 2015). Different PyC detection and quantification methods cover a different part of this continuum (Hammes et al., 2007; Hammes and Abiven, 2013). Across methods, PyC represents on average 15 % of soil organic carbon (SOC) (Reisser et al., 2016).

The stable chemical structure of PyC allows most of it to persist for longer in soils relative to non-PyC, however the range of PyC residence time estimates vary widely (Bird et al., 2015). Large scale modelling of PyC fluxes (Bowring et al., 2022) indicates a persistence in the order of 5000 years. Methods based on radiocarbon dating of charcoal fragments (Carcaillet, 2001; Hajdas et al., 2007; Gavin et al., 2023; Liang et al., 2008) or of the PyC fraction of soil organic carbon (Butnor et al., 2017; Schiedung et al., 2024) have shown that PyC can remain in soils for millennia and up to more than 10 000 years. However, much younger ¹⁴C ages have also been reported (Kane et al., 2010; Ohlson et al., 2009; Krull et al., 2006; Czimczik et al., 2005; Bellè, 2023). These could be due to the rejuvenation of the PyC pool by recent fires or high fire frequency (Kane et al., 2010; Bellè, 2023), or to methods of PyC isolation and detection that may exclude older PyC (Ohlson et al., 2009; Krull et al., 2006), but a centennial PyC turnover rate cannot be excluded. Extrapolations based on laboratory or field incubations (see Azzi et al., 2024, for an extensive list) give estimates in the range of decades to millennia, with significant methodological drawbacks related to model choice and non-representative conditions and time scales (Leng et al., 2019; Azzi et al., 2024). Long term field studies with re-sampling (Hammes et al., 2008b; Lutfalla et al., 2017) or based on chronosequences (Alexis et al., 2012; Lehndorff et al., 2014; Selvalakshmi et al., 2018; Nguyen et al., 2009; Cheng et al., 2008) tend to yield shorter estimates, from decades to centuries. However, PyC residence times estimated from the measurement of PyC losses at the plot scale do not provide true estimates of residence time at the landscape scale. In addition to mineralisation, which represents a net loss of PyC at all scales, PyC may also disappear from the monitored plot via fragmentation to smaller particles and decomposition to less condensed, more oxidised molecules that can evade PyC quantification (Krull et al., 2006), downward vertical transport out of the soil profile, and lateral transport (erosion by surface run-off and leaching by subsurface flow) to nearby or further away locations (Abney and Berhe, 2018). PyC that has been displaced or transformed to another form of soil carbon still resides in the landscape, although its fate may be altered (i.e. it may become more or less subject to mineralisation depending on its chemistry and the environmental conditions at its new location).

Downward vertical transport of PyC has been observed in controlled experiments (Schiedung et al., 2020; Santos et al., 2022; Hilscher and Knicker, 2011) and in the field under different pedo-climatic contexts (Maestrini et al., 2014; Singh et al., 2014, 2015; Schiedung et al., 2023; Dai et al., 2005; Alexis et al., 2012; Bellè, 2023; Bonhage et al., 2022; Vasilyeva et al., 2011). In addition, many studies found that the proportion of PyC in SOC is greater in the subsoil relative to the topsoil, which is attributed in part to transport of PyC to depth (Velasco-Molina et al., 2016; Rodionov et al., 2006; Soucémarianadin et al., 2019). The processes by which PyC can be redistributed vertically in the soil profile have been reviewed by Hobley (2019). Leaching in dissolved form is usually limited (Abiven et al., 2011; Maestrini et al., 2014; Major et al., 2010; Hilscher and Knicker, 2011; Schiedung et al., 2020; Abney et al., 2024) but transport as colloids and larger particles may be significant, in particular when soil texture is coarse, soil structure is loose, sufficient water is present and/or mixing processes (bioturbation, peloturbation, cryoturbation) are intense (Bellè, 2023; Lehndorff et al., 2016; Schiedung et al., 2020; Rodionov et al., 2006; Leifeld et al., 2007). Whether PyC will be retained in the subsoil or transported further down to the altered bedrock and groundwaters depends on its interaction with the minerals and its stabilisation in the subsoil (Brodowski et al., 2007; Guggenberger et al., 2008; Vasilyeva et al., 2011; Santos et al., 2017).

Recent reviews have highlighted the role of erosion in redistributing carbon (C) at the slope, watershed, and global scale (Berhe et al., 2018; Van Oost and Six, 2023; Zheng et al., 2025). These processes are likely exacerbated for PyC (Abney and Berhe, 2018). Erosion rates often increase after fires due to removal of the vegetation and litter layer, and in some cases, increased soil hydrophobicity (Belcher, 2013; Girona-García et al., 2024). PyC has a low density which favours its detachment and transport, even in apparently flat landscapes (Pyle et al., 2017). Manipulative rainfall experiments showed that more than 50 % of deposited PyC can be eroded after a single rain event (Bellè et al., 2021; Rumpel et al., 2009). PyC was enriched in eroded sediments relative to the mineral soil at the plot scale (about 1 m², Rumpel et al., 2006b) and at the micro catchment scale (Chaplot et al., 2005; Cotrufo et al., 2016b; Abney et al., 2019) independent of the slope, but depleted at the larger catchment scale (Chaplot et al., 2005) showing the existence of deposition and storage areas in the landscape. These areas may include alluvial deposits at riverbanks (Cotrufo et al., 2016b), flat areas along the slope or at the toeslope (Abney et al., 2017), patches presenting high surface roughness favouring infiltration of runoff (Boot et al., 2015), where unburned vegetation may retain the charcoal particles (Galanter et al., 2018), or where moderate to low burn severity did not result in increased erosion rates (Abney et al., 2019). Indeed, slope and burn severity seem to have an interactive effect in controlling PyC redistribution by erosion post-fire (McGuire et al., 2021; Galanter et al., 2018). Most of the above-cited studies took place within 4 years post fire. On the decadal time scale, Güereña et al. (2015) showed erosion-related redistribution of mineral soil PyC 10 years after deforestation associated with on site charcoal production in equatorial humid region, but the absence of a flat toeslope (convex toposequence) may have prevented the accumulation of PyC near the stream 62 years after deforestation. Abney et al. (2017) showed that PyC eroded from the slope and deposited on the surface of a riparian area 1 year post-fire had been either decomposed, buried, or remobilised (laterally or vertically, by erosion and/or leaching) 10 years post-fire. Only a few chronosequences of mineral soil PyC have investigated geomorphic factors and landscape positions on longer time scales (Czimczik et al., 2005; Sass and Kloss, 2015). Centennial to millennial time scales are relevant as they are commensurate with the estimates of PyC residence time in soils. Studies on large time- and spatial-scales are needed to capture the balance of erosion processes on the C cycle (Van Oost and Six, 2023). Since PyC dynamics is often related to its interactions with the hydrological cycle (Masiello and Berhe, 2020), the watershed seems to be a relevant scale to trace the fate of PyC. PyC dynamics in the landscape have been studied in temperate to subtropical forests, often with monsoonal precipitation regime and/or in mountainous regions (Galanter et al., 2018; Boot et al., 2015; Matosziuk et al., 2020; McGuire et al., 2021; Abney et al., 2017), in grasslands (Dai et al., 2005; Rodionov et al., 2006), in boreal forests (Kane et al., 2007; Guggenberger et al., 2008) or in ecosystems actively managed by fire (Selvalakshmi et al., 2018; Rumpel et al., 2006a; Alexis et al., 2012; Nicolay et al., 2024, see Fig. S1 in the Supplement for a map of existing studies). Data on other climate zones, geomorphic settings, land-use and vegetation, and pyromes are needed to better understand the drivers of PyC redistribution.

In this work, we investigated whether the vertical and horizontal distributions of soil PyC differed from that of SOC in a shallow sloping, agricultural watershed under oceanic temperate climate. Cropland has dominated the catchment for at least 150 years, implying no or little recent PyC inputs and allowing us to study the resultant, long-term PyC distribution in the landscape. We hypothesised that, even on shallow slopes, PyC produced during past fires was transported downslope by erosion due to its light nature and enhanced post-fire erosion, and deposited at the toeslope where it was protected from further decomposition by burial and/or unfavourable conditions for microbial activity due to frequent water-logging. Owing to the long term stability of PyC, these processes would result, to this day, in higher stocks and higher radiocarbon ages at this position compared to upslope (1). We also made the hypothesis that vertical downward transport of PyC in the soil profile, combined with its high stability, would result in increasing proportions of PyC in SOC with increasing soil depth (2a), and that PyC would be older than SOC at all depths (2b). We measured the concentrations and stocks of PyC down to 60 cm along three toposequences through the convex–concave of the slope. We used two methods of PyC quantification (chemo-thermal oxidation – CTO and hydropyrolysis – HyPy) that cover the intermediate to highly condensed part of the PyC continuum, and also measured the radiocarbon (¹⁴C) values in both total SOC and the PyC fraction.

2 Material and methods

2.1 Study site

Our study site is located in a catchment of 120 ha in the region of Brittany in the west of France (Fig. 1). Over the period 1994–2020, the mean annual temperature was 12 °C (January 6.1 °C–July 18.5 °C) and the mean annual precipitation was 1043.7 mm, well distributed over the year with only slightly wetter autumn-winter and drier spring-summer. The Kervidy-Naizin catchment is part of the ORE (Observatoire de Recherche en Environnement) AgrHys and the critical zone observatories network OZCAR (Observatoires de la Zone Critique: Application et Recherche). It has been studied since the 1970's and was chosen to represent shallow sloping catchments on sedimentary basement typical in the Armorican Massif, that have experienced agricultural intensification over the last 70 years. It feeds the Coët-Dan, an intermittent stream that often dries out in summer. Slopes are shallow (median 3.0 %, first quartile (Q1) 1.7 % – third quartile (Q3) 4.1 %) and present a convex–concave profile. Soils are developed on loess derived locally during the Last Glacial Period (∼115 000 to 11 700 BP, maximal loess deposition in the Northern European Loess Belt around 21 800 BP, Bosq et al., 2023, see also Appendix A1), overlaying weathered siltstone and sandstone, except in the immediate vicinity of the stream where alluvium is the parent material (Walter and Curmi, 1998). The soil texture is silt loam (FAO soil triangle) and soils are acidic (pH_H2O 3.7–6.9). Soils are organised along the convex–concave slope in a gradient of increasing soil evolution towards the stream, typically following the sequence (from summit to toeslope): Cambisol – Haplic Luvisol – Stagnic (solimovic) Luvisol – Stagnosol – Stagnic Fluvisols (transposed to the World Reference Base IUSS Working Group WRB (2022) from Walter and Curmi, 1998). This spatial distribution of soils (sequence and extent of each soil type) is representative of catchments with similar topography and geology in Brittany (Chaplot et al., 2003). The uphill domain is well drained with mainly vertical water flow through the soil volume. Going down slope, the winter water table comes closer to the soil surface, the water flow develops a small horizontal component, and in the soils with stagnic properties, a perched water table above a dense albic horizon can sometimes be observed (Walter and Curmi, 1998). Surface flow and erosion are present in the current agricultural land-use but of limited extent and intensity. There is little export of material out of the catchment (Cros-Cayot, 1996).

Figure 1 (a) Location of the three transects (T1 to T3) and individual sampling sites within the study area. Shading represents topography. Sites are colored by soil type and labelled according to their position along the toposequence (Up = summit, flat and shoulder, Mid = midslope, Down = footslope and toeslope). (b) Organisation of the soil horizons along the lower part of a convex–concave toposequence. Modified from Cheverry (1998). Inset: location of the study area in France.

Because PyC has a long residence time in soil, we also gathered information on land-use change over several millennia (see also Appendix A1). Little is known about the vegetation in Brittany in the periglacial landscape at the end of the last glacial period (11 700 BP) and how it transitioned to the Atlantic oak forest that was likely predominant in Brittany around 7500 BP, when it started to be opened by early human populations for pasture. It is estimated that most of the primary Atlantic forest had been cleared by the end of the Iron Age (Gaudin et al., 2014). Heath and cropland increased globally from the Bronze Age to the end of the Middle Ages (Gaudin et al., 2014). Fire regimes in these successive landscapes are poorly known (Appendix A1). There are no nowadays equivalent to the Atlantic forest, and although a little more is known about managed fire in heathland (Hobbs and Gimingham, 1987) there are only two values of PyC production for heath in the extensive database from Jones et al. (2019). Fire could also have been used to re-open land in a context where arable plots shifted often (Tonnerre, 1992). The transition from a landscape dominated by heath and secondary growth woodland to one where cropland are predominant occurred towards the end of the 19th century (Astill and Davies, 1997). Unintentional fires are virtually inexistent in European croplands under temperate, oceanic climate (Forrest et al., 2024). Agricultural practices such as spreading of plant ashes (Antoine, 2000) or on-site stubble burning (banned from most of Western Europe in the 1980's) could have contributed relatively labile (low temperature) PyC to croplands (see also Appendix A2). The study site is not located downwind of any major road nor, to the best of our knowledge, any current or historical industry that could have been a source of a significant amount of fossil fuel derived PyC (soot). We assume that the current PyC distribution in the landscape integrates the effects of the changes in land use and fire regime since the deposition of the loess that forms the soil parent material (i.e. in the last ∼ 20 000 years, see Appendix A1) and interpret our data accordingly.

Most of the catchment area has been under cultivation or pasture since before 1833, except for a small wood at the south-west border. Plot limits did not vary between 1833 and 1952 and many were delimited by hedgerows and/or paths (see Appendix A2 for details on how we derived this information). In the 1970's, plots were pooled together to form larger fields, some paths were abandoned and hedgerows uprooted. This affected especially the area where transect T1 and T3 are located whereas the vicinity of T2 seemed already devoid of hedgerows by 1952 (see Fig. 1 for transects location).

Table 1 Summary statistics of soil characteristics by slope position and layer. CEC_Co – cation exchange capacity determined using a cobaltihexamine chloride solution, C_fine – concentration of fine soil, Fe_DCB – Iron extracted by the dithionate-citrate-bicarbonate solution, Fe_oxalate – Iron extracted by the oxalate/oxalic acid buffer solution, Subsoil – > 30 cm depth (50–60 cm for all soils except for the Solimovic Cambisol where the subsoil was sampled at 70–80 cm, and an additional sample at 30–40 cm to account for the strongly differentiated horizons in the Stagnic Luvisols and Stagnosol), min. – minimum, max. – maximum.

2.2 Study design and sampling

We sampled along three topographic transects (toposequences) on north-east (T1), south to south-east (T2) and south-west (T3) facing slopes, at upslope (6 sites – 2 per transect), mid-slope (4 sites – 2 at T1, 1 at T2 and T3) and downslope (6 sites – 2 per transect) positions (Fig. 1). Upslope sites were located at the summit (T1) or at the end of a long, shallow (< 2 %) slope running from the summit, just before or at the shoulder (convex portion, T2 and T3). Mid-slope sites were located along the steepest sloping part of each transect (backslope), whereas downslope sites were located at the foot- and toeslope where the slope becomes shallow again (concave portion, see the elevation profiles along the three transects in Appendix Fig. A1). All sites were located in croplands except for the uppermost and lowermost sites of transect 1 which were in a woodland and a cultivated grassland, respectively. Existing sampling sites from previous field campaigns in the study area were relocated based on their GPS coordinates. The GPS coordinates of new sampling sites were recorded. At each sampling site, one lined core was taken using a hydraulic corer and three hand-held auger cores were made about 1.5 m away from each other around the hydraulic core. The hand auger cores were sampled every 10 cm and composited together. Samples and cores were refrigerated within 48 h and kept at 4 °C until analysis. A field soil description was established based on the hand auger core. Selected soil characteristics are summarised in Table 1.

2.3 Soil characteristics

2.3.1 Sample preparation

The soils were dried at 40 °C for a minimum of 48 h and sieved at 2 mm. A 20 g sub-sample of the fine soil was finely ground (≤ 200 µm) using a ball mill (1 min at 25 s⁻¹). The finely ground samples were used for all carbon analyses.

2.3.2 Physical characteristics

Particle size distribution was determined by the Laboratoire d'analyse des sols (LAS, INRAE) according to NF X31 107. In short, the sample was exposed to heated H₂O₂ to break down organic matter prior to particle size evaluation. Coarse sand (> 200 µm) was separated by dry sieving. The remaining material was dispersed in solution and fine fractions (< 50 µm) were determined using Robinson's pipetting method. Fine sand (between 50 and 200 µm) was separated by wet sieving under water flow.

The concentration of fine soil (< 2 mm) was measured for a representative subset of sites using the hydraulic cores. Compaction during coring was evaluated by comparing the depth of soil in the cores to the recorded coring depth. Compaction was considered homogeneous when < 5 % and increasing with depth otherwise. Depth increments equivalent to 10 cm of uncompacted soil were delimited on the outside of the transparent core liner before cutting it open longitudinally and slices were cut out with a knife following the delimitations. The soil was dried at 105 °C for 24 h and weighed (m_tot). Coarse elements > 2 mm were separated from the fine soil by wet sieving, dried again and their weight was recorded (m_rocks). The concentration of fine soil (C_fine, in g cm⁻³) was calculated according to Eq. (1). Note that this is not the bulk density of fine soil as we did not correct for the volume of coarse fragments, but it is enough to calculate the correct carbon stocks as per Eqs. (8) and (9) of Poeplau et al. (2017).

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{C}_{\text{fine}}& ={\displaystyle \frac{(m_{\text{tot}}-m_{\text{rocks}})}{V_{\text{eq\_slice}}}}\\ \text{with }{V}_{\text{eq\_slice}}& =\mathit{\pi }\cdot (\mathrm{7.54}{)}^{\mathrm{2}}\cdot \mathrm{10}\end{array}\end{array}$$

C_fine was estimated for missing sites using a linear model fitted on the measured values with depth, horizon and transect as predictors ($R_{\text{adj}}^{\mathrm{2}}$= 0.47, residual standard error = 0.15 g cm⁻³, p< 0.001).

2.3.3 pH and cation exchange capacity (CEC)

pH and CEC were determined by the LAS (INRAE) according to NF EN ISO 10390 and ISO 23470, respectively. pH was measured in water and in KCl (1 mol L⁻¹) at a 1:5 ($v/v$) soil to solution ratio. CEC was determined based on the exchange of cobaltihexamine ions at soil pH. 50 mL of cobaltihexamine solution at 50 $\mathrm{mmol}+{\mathrm{L}}^{-\mathrm{1}}$ was mixed with 2.5 g of soil. The concentration of cobaltihexamine remaining in solution after exchange was measured by spectrophotometry at 475 nm and CEC was calculated by difference with the initial concentration.

2.3.4 Iron and aluminium oxides and total elemental composition

Iron (Fe) and aluminium (Al) were extracted at the LAS (INRAE) by the Tamm method in the dark (also called oxalate extraction, described in  Courchesne and Turmel, 2008), the Mehra–Jackson method (also called dithionate-citrate-bicarbonate (DCB) extraction Mehra and Jackson, 1960) and hydrofluoric acid (HF) digestion (NF ISO 14869-1) and measured by ICP-AES. The Tamm method extracts iron in organo-metallic complexes and amorphous or poorly crystalline oxides (Fe_oxalate). It is conducted at pH 3, 20 °C and in the dark to avoid iron reduction by UV radiations. 50 mL of oxalate buffer solution was mixed with 1.25 g of finely ground soil during 4 h. The Mehra–Jackson method extracts iron and aluminium in oxides and oxyhydroxides (amorphous and crystalline, Fe_DCB). 0.5 g of finely ground soil was exposed to 25 mL of sodium tri-citrate (0.267 mol L⁻¹) and sodium bicarbonate (0.111 mol L⁻¹) and 1.5 mL of sodium dithionite (200 g L⁻¹) for 30 min at 80 °C. The amount of iron in crystalline oxides was calculated by difference: Fe_crystalline = Fe_DCB – Fe_oxalate. The HF digestion breaks down all silicate minerals and gives access to the total iron (Fe_tot) and aluminium (Al_tot) content of the soil. Other elements measured by ICP-AES on the HF extract include calcium (Ca_tot), potassium (K_tot), magnesium (Mg_tot), manganese (Mn_tot), sodium (Na_tot) and phosphorus (P). From these data we calculated the weathering index (WI, Eq. 2) and the chemical index of alteration (CIA, Eq. 3). The WI is the ratio of Fe_crystalline to total iron. Iron forms in soils are strongly related to pedogenetic processes such as alteration of crystalline iron oxides present in the parent rock, clay formation and breakdown and secondary precipitation of more or less crystalline iron oxides. The CIA was proposed by Nesbitt and Young (1982) as an indicator of the alteration of feldspars from the parent rock (low CIA) to clay minerals (high CIA), during which the mobile elements (Na, Ca, K) are lost whereas Al is retained. Note that K fertilisation and liming will modify CIA in the topsoil independent of weathering.

$$\begin{array}{}\text{(2)}& {\displaystyle }\text{WI}={\displaystyle \frac{{\mathrm{Fe}}_{\mathrm{DCB}}-{\mathrm{Fe}}_{\mathrm{oxalate}}}{{\mathrm{Fe}}_{\text{tot}}}}\text{(3)}& {\displaystyle }\begin{array}{rl}\text{CIA}& ={\displaystyle \frac{{\mathrm{Al}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{3}}}{{\mathrm{Al}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{3}}+{\mathrm{Na}}_{\mathrm{2}}\mathrm{O}+{\mathrm{K}}_{\mathrm{2}}\mathrm{O}+\mathrm{CaO}\ast }}\\ & ={\displaystyle \frac{\frac{{\mathrm{Al}}_{\mathrm{tot}}}{\mathrm{2}M\left(\mathrm{Al}\right)}}{\begin{array}{c}\frac{{\mathrm{Al}}_{\mathrm{tot}}}{\mathrm{2}M\left(\mathrm{Al}\right)}+\frac{{\mathrm{Na}}_{\mathrm{t}}\mathrm{tot}}{\mathrm{2}M\left(\mathrm{Na}\right)}+\frac{{\mathrm{K}}_{\mathrm{tot}}}{\mathrm{2}M\left(\mathrm{K}\right)}\\ +\left(\frac{{\mathrm{Ca}}_{\mathrm{tot}}}{M\left(\mathrm{Ca}\right)}-\frac{{\mathrm{CaCO}}_{\mathrm{3}}}{M\left({\mathrm{CaCO}}_{\mathrm{3}}\right)}\cdot \mathrm{10}\right)\end{array}}}\end{array}\end{array}$$

where X_tot are the concentrations of the elements in weight % (or g X/100 g soil), M(X) is the atomic mass of the respective elements in g mol⁻¹, the factors 2 are to convert from elemental to oxide molar content and CaCO₃ is the carbonate content in g kg⁻¹ soil.

2.3.5 Soil organic carbon and pyrogenic carbon

SOC content was measured by dry combustion – cavity ring down spectroscopy (Picarro G2101-i Isotopic CO₂ and combustion module Costech) calibrated against an in-house standard. Each sample was measured in triplicate. The soils contained no carbonates so total carbon was considered as organic carbon.

PyC content was determined using the chemo-thermal oxidation (CTO-375) method (Gustafsson et al., 1997; Caria et al., 2011) and the hydrogen pyrolysis (or hydropyrolysis, HyPy) method (Meredith et al., 2012). The CTO-375 method removes non-pyrogenic organic matter and the least condensed fraction of the PyC continuum and targets soot-like PyC (Hammes et al., 2007; Gustafsson et al., 2001; Elmquist et al., 2004). HyPy isolates PyC of > 7 fused aromatic rings (Meredith et al., 2012; Wurster et al., 2013), which covers a wider range of the PyC continuum than CTO-375 (Hammes and Abiven, 2013).

Each sample was prepared and measured in duplicate by CTO-375. 40 mg of finely ground sample were weighed in two silver capsules (three if SOC < 0.2 %). The samples were placed in a stainless steel plate and heated in a muffle furnace to 350 °C at 10 °C min⁻¹ then to 375 °C at 1 °C min⁻¹ and kept at this temperature for 24 h. This heating program has been found to prevent in-situ formation of PyC during the procedure (M. Schiedung, personal communication, 2022). The ventilation holes in the oven door were kept open throughout to ensure sufficient oxygen supply for complete combustion of thermally labile organic matter. After cooling to room temperature, the samples were acidified with 20 µL of 1 % ($v/v$) hydrochloric acid (HCl) and exposed to concentrated HCl vapours in a bell jar for 8 h. Samples were then dried successively in a desiccator containing silica beads and under vacuum for 24 h and in a vacuum oven at 40 °C for 3 d. The two or three silver capsules per sample were then composited in a tin capsule and analysed for C content (noted PyC_CTO) as described above for SOC. An in-house standard was submitted to the same procedure with every batch of samples to check for analytical deviation and calculate the coefficient of variation of the PyC_CTO content (CV = 11 %).

Additionally, a subset of samples were analysed by HyPy. We selected the 0–10 cm (top of the A horizon) and 50–60 cm (base of the lowermost B horizon, transition towards the C horizon) depth increments for all profiles except for the Solimovic Cambisol where we selected the base of the solimovic material at 40–50 cm and the base of the structural B horizon at 70–80 cm. To represent the eluvial E horizon an additional intermediate depth increment (20–30 or 30–40 cm) was analysed for mid- and downslope sites presenting this horizon. The finely ground samples were mixed with ammonium dioxydithiomolybdate (5 % $m/m$) and suspended in a 20 % ($v/v$) water-methanol solution to dissolve the catalyst. The solution was then evaporated to dryness at 80 °C. This allowed coating of the catalyst on the soil particles for higher contact surface. Samples were then weighed into pre-combusted quartz inserts and were pyrolysed at 550 °C under a hydrogen pressure of 150 bar and sweep gas flow of 5 L min⁻¹ (ATP). Full details of the pyrolysis programme, including temperature ramp and hold rates, can be found in Meredith et al. (2012). The carbon content of samples before and after hydropyrolysis was measured by elemental analysis (EA). The PyC measured via HyPy is expressed as the proportion of SOC that is made up of PyC (noted PyC_HyPy (%SOC)). PyC_HyPy content was obtained by multiplicating the PyC_HyPy proportion by SOC content.

2.3.6 PyC and SOC radiocarbon dating

Among the samples selected for HyPy, we further restrained the radiocarbon analysis to upslope and downslope sites of each transect.

For ¹⁴C measurement of the PyC fraction, samples were prepared following the HyPy procedure as described above. After hydropyrolysis, the quartz crucible containing PyC was converted to CO₂ by combustion in a sealed quartz tube, as described in Ascough et al. (2024). The sample CO₂ was cryogenically purified and then reduced to graphite by sequential reaction with Zn and Fe. The sample preparation was carried out at the Environmental Radiocarbon Laboratory (Scottish Universities Environmental Research Center, East Kilbride, UK). The ${}^{\mathrm{14}}\mathrm{C}/{}^{\mathrm{12}}\mathrm{C}$ ratio of the graphitised samples was measured by accelerator mass spectrometry (AMS) at the Keck Carbon Cycle AMS Facility (University of California, Irvine, USA). Stable carbon isotope ratios (${}^{\mathrm{12}}\mathrm{C}/{}^{\mathrm{13}}\mathrm{C}$) were measured online on the AMS to correct the measured radiocarbon values for kinetic fractionation effects.

For the bulk ¹⁴C analysis, the homogenised soils were fumigated with 37 % HCl for 72 h at 60 °C to remove all inorganic carbon. After fumigation, the excess acid was neutralised by placing the samples in a closed bell jar with NaOH pellets at 60 °C for a minimum of 48 h. The samples were subsequently wrapped in tin boats and analysed using a coupled elemental analyser-accelerator mass spectrometer (EA-AMS) system (vario MICRO cube, Elementar; Mini Carbon Dating System MICADAS, Switzerland) described in Synal et al. (2007). The ¹⁴C analysis was carried out at the Laboratory of Ion Beam Physics at the Swiss Federal Institute of Technology (ETH), Zürich, Switzerland.

Measured ${}^{\mathrm{14}}\mathrm{C}/{}^{\mathrm{12}}\mathrm{C}$ ratios are reported as F¹⁴C, as described in Reimer et al. (2004). We did not calculate calibrated ages for PyC ¹⁴C data as PyC_HyPy results from the accumulation of fires over time and not a single event that can be dated. We also did not attempt to calculate turnover times for either carbon pool. SOC is not an homogeneous pool with a single turnover time, and even assuming PyC_HyPy is a homogeneous pool, the turnover time calculation requires either a steady state assumption (inputs = outputs, unlikely for PyC produced from sparse fire events), known input rates, or several time points (Torn et al., 2009), neither of which is met.

2.4 Data processing and statistical analysis

Slope position and soil type were determined based on cartographic information from SAS laboratory visualised in QGIS v3.34.8 and confirmed by field observations and laboratory analysis. This led us to modify the grouping of sites from the original design. In particular, we singled out the summit site of transect 1 (Hyperdystric Cambisol, Fig. 1), located in a small woodland and strongly acidic (pH 3.7–4.6), from the other upslope sites (Haplic Cambisols, n=5) under cropland since at least the first half of the 19th century and located either at intermediate flat position along a longer slope or at the slope shoulder. We also singled out the lowest downslope site of transect 3, the only soil showing strong colluvial features (Solimovic Cambisol), and we further subdivided the remaining downslope positions into footslope sites (n= 3), at the transition from Cambisol to Luvisol and toeslope sites (n=2) in the Stagnic Luvisol and Stognosol domain. Mid-slope Cambisols were not subdivided (n= 4). Elevation above sea-level was determined in QGIS based on the GPS coordinates of the sites and the 1 m resolution digital elevation model (DEM) RGE ALTI^® (IGN, 2013). For each site, we calculated its altitude relative to the site closest to the stream in the respective transect, as a proxy for the slope length from that site to the toeslope. Slope was calculated in QGIS from the 25 m resolution DEM BD ALTI^® (IGN, 2015). We used this lower resolution DEM to calculate slope to smooth out the effect of small, localised topographic features. The overall slope were similar between the 3 transects and local slope at the sampling sites ranged from 0 % to 3 %.

Carbon stocks were calculated for each 10 cm depth increment and summed to 60 cm (total stocks), the deepest common depth to all sites, which corresponds to the apparition of the C horizon (altered parent material) at the shallowest profiles (the deepest profiles extended deeper than 90 cm). When available, analytical error was propagated to the final results. In particular, for total stocks, the error was propagated from each depth increment to the sum then to the mean. The degrees of freedom of the variance obtained from error propagation seemed to be erroneously estimated by the Welch-Satterthwaite equation for our small sample size, as observed by O'Hagan et al. (2021), which prevented us from constructing satisfactory confidence intervals for these quantities. Thus we chose to represent the standard deviation, unless specified otherwise. The effect of slope position (original design) was assessed using the Kruskal–Wallis test, a non parametric alternative to ANOVA, after excluding the only solimovic site. Slope positions were compared by the Wilcoxon–Mann–Whitney test and p-values were adjusted for multiple comparisons with the Bonferonni–Holm correction. Some soil groups from the revised design did not contain enough observations to allow meaningful statistical analysis. Trends were assessed visually, taking into account the propagated error. Linear regression was fitted to the concentration and stock data as a function of altitude relative to the toeslope. Diagnostic plots were examined for strong deviations to the normality and homoscedasticity assumptions and identification of high leverage outliers. Correlations between PyC content, F¹⁴C and other soil properties were assessed by Kendall's rank correlation coefficient (τ). Significance was assessed at the 95 % confidence level. Analytical and field replicates were averaged as necessary. All data processing, statistical analysis and figure production was done with R version 4.4.0 (R Core Team, 2021) in Rstudio version 2024.04.1+748 (Rstudio, 2021), using the packages tidyverse (Wickham et al., 2019) and corrplot (Wei and Simko, 2021). All data used to produce the figures presented below can be found in Lebrun Thauront et al. (2025).

3 Results

3.1 PyC in the soil profile

The evolution of PyC_CTO and PyC_HyPy proportion in SOC with depth for the six soil groups is represented in Fig. 2. For all soil groups but the Hyperdystric Cambisol under woodland at the summit, the proportion of PyC_CTO in SOC remains the same in the first 20 to 30 cm due to homogenisation of the soil material by ploughing. Below 30 cm, an increasing proportion of SOC is made up by PyC_CTO as depth increases – between 3.4 ± 0.23 % and 3.83 ± 0.26 % at 30–40 cm and between 6.7 ± 0.40 % and 15.0 ± 2.6 % at 70–80 cm – except for the Solimovic Cambisol at the hill foot which shows a nearly flat profile down to 60 cm depth. The Stagnic Albic Luvisol and Albic Fluvic Stagnosol at the toeslope presented a peak at 20–30 cm (albic horizon) due to a sharper drop in SOC than PyC. Although PyC proportion in SOC increases, PyC concentration decreases with depth in all profiles (Fig. S2).

Figure 2 Profiles of PyC_CTO (dots) and PyC_HyPy (triangles) proportion in SOC with depth by soil group, according to the gradient of soil evolution with slope position. CM – Cambisol, LV – Luvisol, ST – Stagnosol, jd – hyperdystric, ha – haplic, lv – luvic, rw – relectistagnic, st – stagnic, ab – albic, fv – fluvic, sv – solimovic. Error bars represent standard error (typical analytical error when n=1)

The PyC_HyPy profiles show the same trends as PyC_CTO profiles, with the proportion of PyC_HyPy in SOC higher than the proportion of PyC_CTO in SOC by a median of 4.3 %, consistent with the commonly held assumption that the HyPy method measures a wider range of the PyC continuum than the CTO-375 method (Hammes and Abiven, 2013).

3.2 PyC along the toposequence

The amount of PyC in the first 10 cm of the soil, from the toeslope upward, is shown in Fig. 3. PyC_CTO (panel a) tended to be higher at upslope and downslope positions (median = 0.5 and 0.6 g C kg⁻¹ soil, interquartile range (IQR) = 0.18 and 0.15 g C kg⁻¹ soil, respectively) and lower at midslope (median = 0.4 and IQR = 0.04 g C kg⁻¹ soil) but these differences were not significant (χ² = 5.3, p= 0.072). PyC_CTO content in the topsoil decreased significantly with increasing slope (R² = 0.29, p= 0.02, Fig. S3). These tendencies remained visible in the subsoil (50–60 cm depth) although the variability was greater (not shown). The PyC_HyPy content (panel b) and proportion in SOC (not shown) linearly decreased (R² = 0.46, p= 0.01 and R² = 0.51, p= 0.01, respectively) from the summit to the toeslope in the topsoil, but not in the subsoil (not shown). Slope position was a significant factor for PyC_HyPy content (χ² = 7.2, p= 0.027), with significantly more PyC_HyPy in the topsoil at upslope positions (median = 1.6, IQR = 0.22 g C kg⁻¹ soil) than at downslope positions (median = 1.10, IQR = 0.40 g C kg⁻¹ soil). Contrary to PyC_CTO, PyC_HyPy content in the topsoil did not vary systematically with slope (Fig. S3).

Figure 3 (a) PyC_CTO and (b) PyC_HyPy content in the topsoil (0–10 cm) with altitude relative to the toeslope. Points are colored by soil group, according to the gradient of soil evolution with slope position. CM – Cambisol, LV – Luvisol, ST – Stagnosol, jd – hyperdystric, ha – haplic, lv – luvic, rw – relectistagnic, st – stagnic, ab – albic, fv – fluvic, sv – solimovic. In panel a points joined by a line of the same type belong to the same transect, solid line – transect 1, dashed line – transect 2, dotted line – transect 3. Error bars represent the usual analytical error.

3.3 PyC stocks

The stocks of PyC_CTO for each soil group, in the topsoil (0–30 cm), subsoil (30–60 cm), and over the entire soil profile (0–60 cm), is shown in Fig. 4. Contrary to our hypothesis, slope position was not a significant factor explaining soil PyC stocks in this watershed. On average PyC_CTO amounted to 2.5 ± 0.22 t ha⁻¹. PyC_CTO stocks were only slightly lower in a Hyperdystric Cambisol under woodland at the summit (1.9 ± 0.22 t ha⁻¹) and higher in a Solimovic Cambisol at the footslope (3.3 ± 0.26 t ha⁻¹). The proportion of PyC_CTO in total SOC stock varied little, from 1.5 ± 0.35 % to 3.1 ± 0.48 %. PyC in the subsoil represented between 37 % and 51 % of the PyC_CTO stock in the hyperdystric Cambisol under woodland at the summit and the Solimovic Cambisol at the footslope, respectively, and on average 44 % for the other soil groups. SOC stocks amounted to 100 ± 10 t ha⁻¹ on average (min = 68, max = 141 t ha⁻¹), with 12 % to 49 % of the stock below 30 cm (Fig. S4).

Figure 4 PyC_CTO stocks in the soil profile down to 60 cm (darker shade), topsoil (0–30 cm, intermediate shade), and subsoil (30–60 cm, lighter shade) by soil group according to the gradient of soil evolution with slope position. CM – Cambisol, LV – Luvisol, ST – Stagnosol, jd – hyperdystric, ha – haplic, lv – luvic, rw – relectistagnic, st – stagnic, ab – albic, fv – fluvic, sv – solimovic. The numbers above the bars are the proportion of the topsoil (resp. subsoil) PyC stock in the total stock, the numbers inside the bars are the proportion of PyC stock in the total SOC stocks. Error bars represent the propagated standard error.

3.4 Radiocarbon ages of PyC and SOC

There was no significant difference in F¹⁴C with slope position for either SOC or the PyC_HyPy fraction in the topsoil (Fig. 5, panel a) and the subsoil (not shown). PyC_HyPy in the topsoil was significantly older than SOC independent of slope position (Fig. 5a, F¹⁴C range 0.72–0.83 and 0.99–1.14, respectively, which corresponds to 2630–1530 ¹⁴C years BP and 50 BP – modern), and the difference increased with depth (Table S1 in the Supplement). In addition, both the PyC_HyPy fraction (Fig. 5b) and the total SOC (Table S1) were systematically older in the subsoil relative to the topsoil (F¹⁴C PyC_HyPy subsoil median = 0.60, IQR = 0.039, topsoil median = 0.79, IQR = 0.045) except for one profile at the footslope that showed a reverse “stratigraphy”. We found no trend in PyC_HyPy ¹⁴C values with altitude relative to the toeslope and/or soil type except for the Solimovic Cambisol which had a young PyC_HyPy over the whole profile (F¹⁴C 0.73–0.81). All ¹⁴C data can be found in Table S1.

Figure 5 (a) F¹⁴C of bulk SOC (white fill) and PyC_HyPy (grey fill) in the topsoil (0–10 cm) for upslope positions (orange) and downslope positions (purple). Two-sided Wilcoxon–Mann–Whitney test for the difference between C pools (${}^{\ast \ast \ast }$ – p<0.01) and between slope positions within each pool (ns – not significant, p>0.1). (b) F¹⁴C of PyC_HyPy in topsoil (0–10 cm, circles), intermediate (20 to 40 cm depth, triangles) and subsoil (50 to 70 cm depth, squares) horizons in upslope and dowslope positions. brown – Hyperdystric Cambisol at the summit, dark orange – Haplic Cambisols at intermediate flats and shoulders, aquamarine – Haplic Cambisol and relictistagnic Luvisols and the footslope, green – Stagnic Albic Luvisol and Albic Fluvic Stagnosol at the toeslope and pink – Solimovic Cambisol at the footslope. Note the different x-axis scale between the two panels.

3.5 Correlations between soil characteristics and PyC

Both PyC_CTO and PyC_HyPy were positively correlated with SOC (τ = 0.57, Fig. 6). This correlation was weak and insignificant in the topsoil (0–10 cm) and of intermediate strength at 50–60 cm depth (τ = 0.47 and 0.45 for PyC_CTO and PyC_HyPy, respectively). On the contrary, the proportion of PyC_CTO and PyC_HyPy in SOC was strongly negatively correlated with SOC (τ = −0.78 and −0.62, respectively) and the C:N ratio (τ = −0.66 and −0.65, respectively) and positively with δ¹³C (τ = 0.49 and 0.67, respectively). Overall, with the exception of CEC_Co and WI, correlations between PyC and variables reflecting soil chemistry (pH, CIA) and soil mineralogy and texture (Fe_oxalate and Fe_crystalline, clay, silt and sand content) were weaker than those with variables representing SOC quantity and quality, and when significant did not remain so in the 0–10 and 50–60 cm depth intervals separately, indicating that they may result from co-variation with depth rather than actual mechanistic relationship. They were however stronger for PyC_HyPy concentration, in particular the correlation with Fe_oxalate (τ = 0.47). F¹⁴C of PyC_HyPy was positively correlated with CEC and Fe_oxalate (τ = 0.66 and 0.37, respectively) and negatively correlated with clay content (τ = −0.40).

Figure 6 Kendall's rank correlation between PyC content (CTO and HyPy), PyC proportion in SOC (CTO and HyPy), F¹⁴C (PyC_HyPy and bulk SOC) and selected soil characteristics, including SOC and indicators of SOC quality (δ¹³C, C:N), soil chemistry (pH_H2O, CEC_Co, Fe_oxalate and Fe_crystalline), soil texture (Clay to Sand) and soil evolution (Silt:Clay, weathering index WI and chemical index of alteration CIA), across depths. Significant correlations at the 95 % confidence level are in bold. Correlations that remain significant and of the same sign in at least one depth increment (0–10 or 50–60 cm) are outlined in black.

4 Discussion

4.1 Cropland soils store old PyC after centuries of cultivation

Our results attest to the durable presence of PyC in cropland soils even after centuries of cultivation. The transition to continuous cropping may not have completely suppressed fire from the landscape, as stubble burning was practised in Western Europe, but we assumed that ensuing PyC inputs were either small (between 25 % and 50 % of cereal crop biomass is harvested and crop residues have many other uses besides burning, Smil, 1999), or/and of more labile nature due to low (< 350 °C) fire temperatures (Rasmussen et al., 1986). PyC_HyPy in the first 10 cm of the soil had substantially lower (i.e. older) F¹⁴C than SOC (Fig. 5a), highlighting their longer residence time of PyC in soils, up to millennia. These low F¹⁴C values (0.72–0.83, 2630–1530 years BP) further indicate the presence of a non-negligible proportion of old PyC remaining in the topsoil, consistent with the assumption of limited to no PyC inputs in the last two centuries. They are similar to F¹⁴C found at 0–15 cm depth in boreal soils with discontinuous permafrost in Canada (0.76–0.86, Schiedung et al., 2024) but significantly higher (i.e. younger) than values for continuous permafrost soils from the same study (0.44–0.56). They are lower (i.e. older) than in the 0–15 cm layer of tropical soils under dry shrubland and forest (0.85–0.95 and 0.89–1.05, respectively, Bellè, 2023). These differences could be explained by a combination of different PyC mineralisation rates (related to soil temperature and PyC quality) and different fire regimes (more frequent and recent PyC inputs in both climatic zones).

PyC_CTO concentrations in the topsoil under cropland (0.4 to 0.7 g C kg⁻¹ at 0–10 cm, Figs. 3 and S2) were very similar to that found in soils of the Région Centre in France (mostly croplands, Q1 = 0.4 to Q3 = 1 g C kg⁻¹ at 0–30 cm, Paroissien et al., 2012) and in Swiss cropland soils (0.4 to 0.8 g C kg⁻¹ at 0–20 cm, Agarwal and Bucheli, 2011). They were slightly lower than that measured at 0–20 cm in grassland (0.8 to 1.3 g C kg⁻¹) and forest (0.6 to 2.7 g C kg⁻¹) soils across Switzerland (Agarwal and Bucheli, 2011) and at 0–5 cm in remote and rural areas of the UK and Norway (range 0.45 to 3.13 g C kg⁻¹, average 1.27 g C kg⁻¹, Nam et al., 2008). These comparisons point to a potential loss of PyC_CTO from the topsoil in cultivated soils.

However, in our study area, PyC_CTO stocks over 60 cm were lower in a soil under woodland (1.9 ± 0.22 t ha⁻¹, Summit – CM-jd in Fig. 4) than under croplands in similar topographic positions (2.6 ± 0.22 t ha⁻¹, Flat and Shoulder – CM-ha in Fig. 4), after more than 150 years of cultivation. To the best of our knowledge, there exist no other study that measured PyC_CTO stocks in a similar pedoclimatic context. Silva et al. (2023) found an average PyC_CTO stock to 1 m depth of 1.37 ± 0.65 t ha⁻¹ in Amazon forest fragments in Brazil. This is similar to our stocks to just 30 cm depth (1.2 to 1.6 t ha⁻¹, Fig. 4). Santos et al. (2017) and Soucémarianadin et al. (2014) found PyC stocks to the lower limit of the B horizon between 1.75 and 4.9 t ha⁻¹ in Podzols of North American forests using a PyC isolation method based on chemical oxidation (KMD method). These higher values could be attributed to different fire regimes, soil properties or a larger coverage of the PyC continuum by the KMD method compared to the CTO-375 method which only detects the most condensed, soot-like PyC (Hammes et al., 2007; Gustafsson et al., 2001; Elmquist et al., 2004). Across pedoclimatic contexts, PyC stocks determined with the BPCA method (Lehndorff et al., 2016; Rodionov et al., 2010) are 2 to 20 times larger than PyC_CTO stocks (Schiedung et al., 2024; Bellè, 2023, and this study), which points to method differences as the most likely explanation and prevents meaningful comparisons between the data obtained from these two different methods.

The subsoil (30–60 cm) represented a important share of the total PyC_CTO stock (44 %, Fig. 4) in our study area. This proportion was slightly lower (not significant) in the Hyperdystric Cambisol under woodland (37 ± 6.2 %) than in the Cambisols under cropland (44 % to 46 %). Soucémarianadin et al. (2014) and Silva et al. (2023) found only 17 % and 30 % of the PyC stock below 30 cm in Podzols and Ferralsols under forest, respectively, whereas under cropland Lehndorff (2016) found more than 50 % of the PyC stock in the subsoil (below the A horizon) for an Andosol, an Alisol, a Cambisol and a Vertisol and less than 50 % for a Fluvisol and another Alisol, when not under paddy rice cultivation. Redistribution of PyC to the subsoil in croplands could be favoured by mixing and fragmentation of PyC particles due to tillage. All soil profiles in recently tilled fields showed a net homogenisation of PyC_CTO (Fig. S2), SOC and N content (not shown) in the first 20 to 30 cm of the soil compared to the woodland and grassland sites (summit and toeslope positions in transect 1, respectively). Rodionov et al. (2006) showed that the fine sand and silt fractions (particle size 2–250 µm) at 0–10 cm depth were enriched in PyC in arable fields compared to adjacent native grassland, which they attributed to the crushing of bigger PyC particles by cultivation practices. Using a soil physics model, Hobley (2019) demonstrated that smaller PyC particles were more likely to be transported to depth, especially when eluviation rate is high relative to decomposition rate, which is likely the case for PyC particles in well drained soils. PyC redistribution to the subsoil, where decomposition rates are usually lower and erosion inexistent, could contribute to its preservation in cropland soils.

Altogether, these results show that in temperate regions and on shallow slopes, PyC can persist in cropland soils for centuries, bearing traces of past fire regimes.

4.2 SOC becomes enriched in PyC with depth

We observed an increasing proportion of PyC in SOC with depth in all profiles and for both PyC determination methods (Fig. 2). This is common for PyC_CTO (Qi et al., 2017; Agarwal and Bucheli, 2011; Koele et al., 2017; Silva et al., 2023) and PyC measured using other quantification methods (Andreeva et al., 2011; Brodowski et al., 2005; Butnor et al., 2017; Dai et al., 2005; Hammes et al., 2008a; Matosziuk et al., 2020; Soucémarianadin et al., 2019; Güereña et al., 2015; Velasco-Molina et al., 2016), across soil type, vegetation and climate, although constant or decreasing proportion of PyC in SOC with depth have been observed under grasslands (Rodionov et al., 2006), in Andosols and in some cases under paddy rice cultivation (Lehndorff et al., 2014), in Podzols (Soucémarianadin et al., 2014), and where there are frequent PyC inputs at the surface from stubble burning (Rumpel et al., 2006a). At the same time, F¹⁴C of PyC_HyPy was always lower (older) than F¹⁴C of SOC (Table S1) and decreased sharply with depth in all profiles but two (Fig. 5b). We are the first to report F¹⁴C of PyC_HyPy in temperate cropland soils but Butnor et al. (2017) found a regular decrease in the F¹⁴C of the chemically resistant SOC fraction (which can be equated to PyC) from 0–20 cm (0.88–0.66) down to 100 cm (0.28–0.54) in forest soils in the southern US, in close agreement with our findings.

Several hypotheses have been brought forward by the above-cited authors to explain the increasing proportion of PyC in SOC with depth: (1) PyC is preferentially transported to depth relative to non-PyC, (2) PyC is preferentially preserved at depth (e.g. via retention and stabilisation on mineral surfaces) relative to non-PyC, or (3) PyC has an intrinsically higher biological and/or chemical stability relative to non-PyC, which allows a larger proportion of the inputs to reach and stay in the subsoil before being decomposed by microorganisms or abiotic reactions. Although our study was not designed to answer this question, we examined the three mechanisms of PyC enrichment with depth and argue that higher stability of PyC relative to non-PyC is sufficient to explain the profiles we observe (Fig. 7).

Figure 7 Conceptual, 1D model of PyC and non-PyC persistence in soils based on differences in biological and chemical stability. After a fire some PyC is incorporated in the topsoil (a) and starts to migrate downward into the soil profile by leaching, eluviation and turbation (b). Vegetation regrowth provides non-PyC inputs continuously throughout the post-fire succession (b, c), natural fire – roman type, slash and burn – italic type). Over time, transport and decomposition create a concentration gradient with depth (b, c). The transport rate of PyC may be slower than that of non-PyC, but because it is more stable, a larger proportion of it remains available to be transported and its concentration gradient is less pronounced. On the contrary, a large part of non-PyC is lost via microbial respiration before it can reach the altered bedrock. After several fire cycles and millennia of soil evolution, it results in the observed concentration profiles, under the current agricultural land-use (d).

(1) Estimates of downward vertical transport rate of particulate and dissolved PyC range from 2 to 30 mm yr⁻¹ (Rumpel et al., 2015). Such fast transport of PyC to depth is not compatible with our data (see also Fig. 5b in Hobley, 2019): at 2 mm yr⁻¹ it would take only 250 years for PyC to migrate from the 0–10 to 50–60 cm layer when the age difference we observe is between 1740 and 7410 ¹⁴C years. The sharp increase of PyC proportion in SOC at the E horizon of the Stagnic Albic Luvisol and Albic Fluvic Stagnosol (Fig. 2), where SOC eluviation is the strongest, is also contrary to the hypothesis of preferential vertical transport of PyC relative to SOC in the pedoclimatic context of our study area.

Fast vertical transport in the soil profile may only apply to a small portion of the PyC pool. Indeed, dissolution and leaching of PyC is usually limited relative to SOC (Abiven et al., 2011; Maestrini et al., 2014; Major et al., 2010; Hilscher and Knicker, 2011; Schiedung et al., 2020; Abney et al., 2024). Alternatively, the elevated rates measured in previous studies may apply only to soils with high porosity (Leifeld et al., 2007), coarse texture (Schiedung et al., 2020) and/or where annual rainfall is important (Alexis et al., 2012; Major et al., 2010; Nguyen et al., 2009), conditions that are not met in our study area.

Soil mixing processes (bioturbation, peloturbation, cryoturbation and that induced by windthrow and agricultural practices) can transfer particulate PyC to depth (Lehndorff et al., 2016; Schiedung et al., 2020; Carcaillet, 2001) but are not selective for PyC. As demonstrated by Hobley (2019, see figure 1d), these spatially heterogeneous processes are also inconsistent with the nearly systematic increase of PyC_HyPy age with depth that we observe here (Fig. 5) and in other studies (Carcaillet, 2001; Hajdas et al., 2007).

(2) Both PyC and SOC in the subsoil are effectively protected from erosion and usually experience lower microbial activity. The proportion of mineral-associated SOC, a fraction which typically has higher stability, is often found to increase with depth (e.g. Sanderman et al., 2021). PyC was found to form mineral-associations potentially mediated by iron (Solomon et al., 2012; Brodowski et al., 2005; Schiedung et al., 2023). Schiedung et al. (2020) have shown that PyC migrating in a saturated soil column was preferentially found in the mineral associated fraction in a finer textured soil with higher clay and Fe_oxalate contents, and interaction of PyC with iron (hydr)oxides was suggested as the reason for PyC retention in podzolic Bs horizons (Soucémarianadin et al., 2014; Santos et al., 2017), but there is no compelling evidence that PyC is preferentially preserved relative to non-PyC.

In our study, Fe_oxalate (iron in organo-metallic complexes + poorly crystalline oxides) showed weak, negative correlations with the proportion of PyC_CTO and PyC_HyPy in SOC and a weak positive correlation with F¹⁴C of PyC_HyPy (Fig. 6) indicating no preferential preservation of PyC with increasing presence of reactive iron phases. The proportion of PyC_HyPy in SOC presented an intermediate, positive correlation with iron in crystalline oxides overall (Fig. 6) but not in the 50–60 cm layer (τ = −0.09, p= 0.73), indicating that this form of iron was not a factor of PyC stabilisation in the subsoil. In previous studies, PyC content correlated positively with pyrophosphate and Meilich-3 solution extractable iron (iron in organo-metallic complexes, Eckmeier et al., 2010; Soucémarianadin et al., 2014; Butnor et al., 2017) but not with Fe_oxalate, Fe_DCB and total iron (Qi et al., 2017; Lehndorff et al., 2016; Velasco-Molina et al., 2016), warranting further investigation of PyC interaction with iron in soils.

We also found no significant correlation between the proportion of PyC_CTO and PyC_HyPy in SOC and clay content (Fig. 6). Their was a negative correlation of intermediate strength between PyC_HyPy F¹⁴C and clay content overall but it disappeared when looking at each depth interval separately, which points to a co-variation with depth rather than a mechanistic relationship. This is not inconsistent with the sparse evidence of correlations between PyC and clay content in the literature (Qi et al., 2017; Velasco-Molina et al., 2016; Butnor et al., 2017; Paroissien et al., 2012; Soucémarianadin et al., 2014; Schiedung et al., 2024; Sanderman et al., 2021). In addition, uncertainties remain about whether or not PyC in the clay fraction is actually mineral-associated (Lutfalla et al., 2017; Hilscher and Knicker, 2011; Brodowski et al., 2007; Skjemstad et al., 1996) and whether or not mineral-associated PyC is more stable than its particulate counterpart (Hilscher and Knicker, 2011; Chassé et al., 2021; Vasilyeva et al., 2011).

(3) Millennial persistence of at least some portion of PyC is attested by radiocarbon dating of charcoal fragments (Carcaillet, 2001; Hajdas et al., 2007; Gavin et al., 2023; Ohlson et al., 2009; Liang et al., 2008) and thermally or chemically isolated PyC (Butnor et al., 2017; Schiedung et al., 2024). Long-term field experiments have shown that PyC tends to be more stable than the rest of SOC (Lutfalla et al., 2017; Vasilyeva et al., 2011; Hammes et al., 2008b), which turns over on decadal to centennial time scales (Schmidt et al., 2011; Barré et al., 2016). Our observations support this difference, with modern to ≈ 1000 ¹⁴C years old total SOC and ≈ 1500 to ≈ 10 000 ¹⁴C years old PyC_HyPy. Higher stability of PyC in soil would give it time to be transported down the soil profile from the surface where it is produced, and explain its higher proportion in the subsoil relative to SOC, despite potentially large inputs of relatively stable SOC from roots directly to the subsoil (e.g. Rasse et al., 2006; Sokol et al., 2019). During ageing, PyC is fragmented and oxidised (Hockaday et al., 2006; Sorrenti et al., 2016; Pignatello et al., 2015), its solubility increases (Abiven et al., 2011) and colloids are formed which could move more easily through the soil porosity. These processes could favour PyC transport to depth over time (Carcaillet, 2001). In addition, small, oxidised PyC particles are more prone to form mineral association (Schiedung et al., 2020), which could favour their retention once in the subsoil and contribute to the overall persistence of PyC.

Thus we argue that intrinsically higher stability of PyC relative to SOC, not preferential transport or preservation, explains the large proportion of very old PyC in subsoil SOC stocks in most soils in the Kervidy-Naizin watershed (Fig. 7). Despite good drainage in the upslope domain and noticeable eluviation in the downslope domain (Walter and Curmi, 1998), vertical transport of PyC likely remains slow, consistent with the large F¹⁴C difference between topsoil and subsoil PyC and the lingering presence of old PyC in the topsoil.

4.3 Slope position does not control stocks and age of PyC

4.3.1 Toeslopes do not accumulate PyC in the long term

The hypothesis that in this watershed, toeslopes would represent large stocks of old PyC, is refuted. On the contrary, PyC_CTO stocks were relatively homogeneous along the gradient of soil evolution with slope position (Fig. 4), around 2.5 t ha⁻¹. In the topsoil, PyC content and proportion in SOC even decreased with elevation relative to the toeslope for PyC_HyPy, the fraction of PyC that is 7–14 aromatic rings or bigger, and showed a minimum at midslope for PyC_CTO, at the most condensed end of the PyC continuum (Fig. 3). Our hypothesis was based on the assumption that after a fire, PyC characteristics (low bulk density, hydrophobicity) and favourable conditions for erosion (lesser vegetation and litter cover) would lead to the redistribution of PyC to the toeslope, where the shallow slope and riparian vegetation would stop it from reaching the stream. Its incorporation into the soil profile would then protect it from further transport and decomposition. To explain our results, we suggest that PyC produced by past fires was either not significantly redistributed along the slope, washed out to the stream, or lost from the toeslope since then.

If evidence of (preferential) PyC erosion post fire is numerous (Rumpel et al., 2009; Cotrufo et al., 2016a; McGuire et al., 2021), most PyC erosional studies took place on steep slopes (Cotrufo et al., 2016a; Rumpel et al., 2006b; Güereña et al., 2015), or under high rainfall intensity on bare soil (Rumpel et al., 2009; Bellè et al., 2021). In case of high severity fires, PyC can be eroded independent of slope if sufficient rainfall occurs (McGuire et al., 2021; Boot et al., 2015). On shallow slopes and flat areas or where burn severity is lower or patchy, PyC may only be redistributed locally (Pyle et al., 2017; Galanter et al., 2018). Little is known about past fire regimes in the region, let alone fire severity (Appendix A1). However, under the current agricultural land-use which leaves the soil bare a large part of the year, erosion proceeds in successive steps, with overland flow infiltration several times along the slope and without formation of rills and gullies able to carry material far away (Cros-Cayot, 1996). This suggest that post-fire transport of PyC to the toeslope could have been limited.

Even where there was transport of PyC to the toeslope, several studies have failed to detect long term accumulation of PyC in low topographic positions (Rumpel et al., 2006a; Galanter et al., 2018; Güereña et al., 2015; Abney et al., 2017). Galanter et al. (2018) suggested that PyC directly flowed over the riparian area and into the stream, whereas Güereña et al. (2015) advocated that PyC was only temporarily stored in the lower part of a convex slope before being picked up by the river in later erosion/flood events. The same process could explain the disappearance of PyC at the bottom of a concave slope 10 years post-fire observed by Abney et al. (2017), although the absence of visible erosion features led them to favour the burial or decomposition hypotheses. Here, the presence of a dense albic eluvial layer under the A horizon in the Stagnic Luvisols and Stagnosols near the stream (Fig. 1) reduces infiltration (Curmi et al., 1998; Cros-Cayot, 1996), favouring surface flow that could take away light, freshly deposited PyC to the stream and prevent its incorporation deeper into the soil profile. Indeed, in another study based in the same watershed, Vongvixay et al. (2018) have shown that if riparian trees and grass buffer strips are efficient at intercepting suspended solids transported by runoff during most of the year, they seem to loose this ability during winter, maybe because the vegetation is submerged. Rumpel et al. (2006b) suggested that because it is preferentially found in the light fraction, PyC is easily eroded and transported out of the watershed. Studies considering different spatial scales showed that PyC was enriched in sediment eroded from the headwater catchment but depleted in the larger order catchment (Chaplot et al., 2005), indicating that PyC was transported away from its production site and redeposited further downstream, for instance at river banks (Cotrufo et al., 2016a).

Finally, PyC at the toeslope may have been lost overtime at a faster rate than in the upslope domain. Studies on SOC and DOC in porewater and in the stream in the Kervidy-Naizin catchment have shown that the fluctuation of the water table in the soil profile affects OC dynamics (Lambert et al., 2013, 2014; Jeanneau et al., 2014). In the area close to the stream where the water table is present to the surface a large part of the year, SOC mineralisation was reduced relative to midslope areas but dissolution was increased which resulted in export of OC from the subsoil. The dissolution of SOC was related to the destabilisation of OC-retaining iron oxides under reducing conditions (Lambert et al., 2013; Jeanneau et al., 2014). Indeed, the Stagnic Luvisols and Stagnosols at the toeslope were depleted in free and amorphous iron which partly re-precipitated as more crystalline iron oxides in mottles in the illuvial B horizon. PyC_HyPy content was significantly positively correlated with intermediate strength to Fe_oxalate (free and amorphous iron) and strongly to Fe_crystalline in the topsoil (τ = 0.47 and 0.67, respectively, p<0.01) but not in the subsoil, indicating that the same dissolution processes may have affected PyC_HyPy, preventing its accumulation in toeslope positions.

The difference between PyC_HyPy and PyC_CTO along the toposequence may result from differences in PyC quality. PyC_HyPy represents PyC from 7 aromatic rings in size while PyC_CTO only accounts for the most condensed, soot-like PyC. As such, PyC_HyPy is likely to be more prone to solubilization and more labile, and thus to depend more on mineral interactions (e.g. with iron oxides, see above) for its preservation, while PyC_CTO would behave more like an inert soil particle affected primarily by erosion (see Sect. 4.3.2) and eluviation (see Sect. 4.2).

4.3.2 Recent erosion caused redistribution of PyC with the soil matrix to specific locations

PyC_CTO and SOC stocks were higher in the Solimovic Cambisol at the toeslope of transect 3 relative to all other soils along the gradient of soil evolution with slope position (Fig. 4), whereas cross transects, PyC_CTO content in the topsoil decreased with increasing slope (Fig. S3). We interpret these observations as signs of recent (relative to PyC residence time) erosion processes that have redistributed topsoil PyC in and out of the watershed.

The soil profile that we classified as a Solimovic Cambisol presented a thick A horizon (55 to 60 cm). The composition of this horizon did not markedly differ from the A horizon of soils higher in the same transect, which supports the local origin of this material. PyC and SOC F¹⁴C did not decrease with depth like in the other undisturbed soil profiles but were instead slightly higher at 40–50 cm depth than at 0–10 cm (Table S1), a feature consistent with depositional processes (see also Fig. 2d in Hobley, 2019). Since PyC content in the soil decreased with depth (Fig. S2), a change in the relative soil surface due to erosion at steeper, midslope positions would explain the lower PyC_CTO concentrations relative to less eroding, shallower summit and toeslope positions.

Despite shallow slopes, erosion in this landscape is likely favoured by extended periods without vegetation cover, and cultivation in the direction of the slope. Based on the extent, thickness and density of deposits, Cros-Cayot (1996) estimated that within one year of conversion from grassland to maize, 15 t of topsoil were eroded and deposited at the concavity in a 3 ha field located in the same catchment as our transects. Using sediment traps along several toposequences in another catchment in the region (similar geology, climate, topography and soil types), she estimated that 152.5 g m⁻² of solids were eroded yearly, of which 17 % (25.9 g m⁻²) was exported to the river while the rest was redeposited along the slope. These values were considered as an upper limit for the Kervidy-Naizin catchment where structural stability of the soil was higher (Cros-Cayot, 1996). This is in line with the results of Vongvixay et al. (2018), who measured discharge and turbidity in the stream at the outlet of the catchment over 9 years, and calculated annual suspended sediment yields ranging from 3 to 22 g m⁻².

Another similar catchment in the region had higher suspended sediment yields despite a greater proportion of woodland and grassland (smaller sediment source area) and a higher density of hedgerows (more obstacles to sediment fluxes), which the authors attributed in part to the retention of sediments by riparian trees and grass buffer strips in the Kervidy-Naizin catchment (Vongvixay et al., 2018). The location of the Solimovic Cambisol we sampled, at the toeslope and about 15 m away from the riparian tree line, is consistent with this interpretation. However, we did not find noticeable amounts of solimovic material at the lowermost sampling site of T1 and T2, which suggest that such accumulation is localised. It may be related to a slightly raised bank providing a barrier between the slope and the stream and/or to the short toeslope at T3. Where the toeslope is longer, material eroded at midslope may be redistributed in a diffuse way, consistent with the stepwise erosion model proposed by Cros-Cayot (1996) for these watersheds. Alternatively, sediments deposited at the toeslope throughout the year may be washed out during winter floods.

In the Solimovic Cambisol at the toeslope, SOC had an F¹⁴C of 1.03 at 40–50 cm depth (Table S1), indicative of the presence of a non negligible proportion in the subsoil of organic carbon formed following the bomb peak (circa 1963). Throughout the A horizon, soil chemistry (pH, CEC, $\mathrm{C}/\mathrm{N}$ ratio, CIA, not shown), SOC and PyC content (Fig. S2) were almost constant, showing that the deposited material did not have time to evolve. On shallow slopes and structurally stable soils, erosion is likely negligible under permanent vegetation cover. Agricultural activity in the area may date back as far as 5000 BP but continuous cultivation is likely much more recent (Astill and Davies, 1997; Tonnerre, 1992). Slope continuity along transect 3 was only established in the 1970's when the hedgerows at the former plot boundaries were removed throughout the Kervidy-Naizin catchment (Appendix Fig. A1). In the following decades, this catchment saw an intensification of agricultural practices (see Appendix A2) with an increasing frequency of maize in the rotation (Bordenave and Merceron, 1999). When maize is planted after a winter crop (e.g. wheat), a common occurrence in the fields where our transects are located, the soil can remain bare for more than 6 months, which favours erosion (Bordenave and Merceron, 1999). This concurs to indicate that the establishment of the Solimovic Cambisol at the toeslope of T3 is recent.

There was no PyC enrichment or depletion during erosion and deposition at T3 (Fig. 3). PyC accumulation at this site was only related to accumulation of topsoil material, and would have gone undetected if we did not consider the entire soil profile. Aggregate stability measured in previous studies was high (Cros-Cayot, 1996) and soil was likely transported in aggregated form, as evidenced by the absence of textural differentiation in the solimovic material (not shown). Under these conditions, centuries to millennia old PyC was not preferentially eroded, as opposed to what was observed for fresh PyC (Bellè et al., 2021; Cotrufo et al., 2016b; Rumpel et al., 2009; Chaplot et al., 2005). It suggests that during aging, PyC becomes involved in aggregates and mineral interactions which stabilise it with regard to erosion.

5 Conclusions and perspectives

We showed that topographic position was not the main driver of PyC stocks in the landscape, in a shallow sloping watershed dominated by agricultural land-use under temperate oceanic climate, and after hundreds of years of cultivation. If it existed, PyC enrichment at the toeslope would have been only temporary and probably reversed in the decades following fire due to subsequent erosion or conditions favouring dissolution and leaching of PyC. Changes in erosion dynamics related to land-use redistributed already aged PyC to a localised area of soil accumulation without enrichment or depletion of PyC. This could make PyC an interesting indicator to trace erosion in temperate agricultural landscape where fire is rare. More studies are needed to assess whether different erosion modalities (sheet vs. rills and gullies) and aggregate stability (i.e. whether or not the soil disaggregates during erosion) affect the fate of aged PyC. Large stocks of old PyC in the subsoil supports the long persistence of PyC in soils and slow advection towards soil depth under the pedoclimatic conditions of our study area. More than 150 years of cultivation did not seem to deplete PyC stocks relative to an adjacent wooded area but may have favoured PyC redistribution below 30 cm depth, although replication is needed to confirm these findings. We call attention to the importance of subsoil sampling when studying persistent forms of organic carbon and estimating their stocks. These results are of interest in the context of biochar application to agricultural soils, as losses of biochar from the topsoil may be accounted for by lateral and/or vertical transport, and not only by mineralisation, which impacts the carbon budget of the system. Identifying the proportion of PyC produced which is quickly transported away from the watershed and that which remains in soils for millennia after a fire is an important knowledge gap that still needs to be investigated to close the terrestrial PyC budget, but it necessitates a better estimate of PyC production from fires and knowledge of fire history and past fire regimes. Pools of PyC of different quality measured by different PyC quantification methods showed different patterns of concentration in the topsoil along the gradient of soil evolution with slope position. The most condensed, soot-like form of PyC measured by CTO tended to be depleted at midslope, possibly as a results of recent erosion that exposed former subsoils less rich in PyC, whereas the intermediate to highly condensed PyC measured by HyPy was increasingly depleted with proximity to the river, which may be related to changes in iron-mediated mineral interactions favouring its dissolution and leaching. These interactions need to be studied in more details. Multi-method PyC characterisation is an opportunity to study the interplay between PyC quality and soil processes.

Appendix A: Land-use and fire history at the study site

A1 Land-use and fire history in Brittany over the Holocene

The loess that form the parent material of the soils in the Kervidy-Naizin watershed are derived from both local and remote sources and were likely deposited on the periglacial landscape during the last glacial period (∼ 115 000 to 11 700 BP, Pellerin and Van Vliet-Lanoe, 1998; Van Vliet-Lanoé et al., 1998). More specifically it is estimated that the loess deposition in the Northern European Loess Belt started at about 32 000 BP and that maximum deposition occurred around 21 800 BP Bosq et al. (2023), and loess deposits in Northern Brittany have been dated between 15 900 and 67 000 BP Van Vliet-Lanoe et al. (1997) but no datation exists at our study site. To the best of our knowledge, there are little indications as to the vegetation dynamics in the region at the Pleistocene to Holocene transition and into the early Holocene. Paleoenvironmental and archaeological studies have shown that the Atlantic oak forest that was likely predominant in Brittany before 7500 BP started to be opened by early human population during the neolithic period, from the coast inwards. The forest was first cleared for pasture, with only localised cropping. Abandoned land was re-colonised by heath and heliophytes tree species that could be maintained by grazing or fire (Gebhardt and Marguerie, 2006; Briard et al., 1989). It is estimated that most of the primary Atlantic forest had been cleared by the end of the Iron Age (Gaudin et al., 2014). Heath and cropland increased globally from the Bronze Age to the end of the Middle Ages with periods of acceleration or stability (Gaudin et al., 2014). Tonnerre (1992) points out that in the Early Middle Ages (1500–1000 BP), small pieces of land were cropped probably for short periods and rotated often, which could have led to frequent slash-and-burn practices. Blank spaces (small woods, some wetlands) were taken into cultivation in the late mediaeval – early modern period but heath and secondary growth woodland dominated the landscape until the end of the 19th century (Astill and Davies, 1997). However, heath was not synonymous of land abandonment as part of the heathland was only temporary (a form of fallow), part was grazed, and part was used as a source of plant material for animal bedding and fertilization of the fields after “composting” or burning the biomass to ashes (Antoine, 2000, see the definitions of Étrepage and Écobuage). The division in small plots, delimited by fences, ditches and/or hedgerows, originated in the 16th–17th century and lasted until the 1970's when plots were pooled together to allow the use of heavier machinery on larger surfaces (Astill and Davies, 1997).

There are no good nowadays equivalent of the Atlantic forest that would allow us to estimate pre-historical fire regimes in the area. Fires during the period of human presence after the first deforestation may have been mostly related to slash and burn practices to maintain the heather or clear the secondary growth forest (Gebhardt and Marguerie, 2006; Briard et al., 1989). Heath and moorland fires can be very intense, although very variable (Hobbs and Gimingham, 1987; Davies et al., 2022), and burn into the moss and litter layer (Davies et al., 2022). However, as it is a multi-millennial landscape shaping practice, users have developed a traditional knowledge on how to minimise soil impacts (in particular erosion) while achieving management goals (Davies et al., 2022). This could have limited the severity of heath fires. In croplands, the addition of ashes resulting from heathland vegetation burning, a practice described by Antoine (2000) for the modern period (late 15th to early 19th century), contributed an unknown amount of PyC to the soil.

A2 Land-use and fire in the study area in the last 150 years

The reconstruction of land-use and plot borders in the catchment area over the last 150 years was based on maps, aerial photographs and farmer's declarations. You can visualise most of the cited maps on Géoportail (https://www.geoportail.gouv.fr/carte?c=-2.8423422520194026,48.009738098366256&z=15&l0=GEOGRAPHICALGRIDSYSTEMS.ETATMAJOR40::GEOPORTAIL:OGC:WMTS(1)&l1=GEOGRAPHICALGRIDSYSTEMS.MAPS.SCAN50.1950::GEOPORTAIL:OGC:WMTS(1;h)&l2=ORTHOIMAGERY.ORTHOPHOTOS.1950-1965::GEOPORTAIL:OGC:WMTS(1;h)&l3=GEOGRAPHICALGRIDSYSTEMS.PLANIGNV2::GEOPORTAIL:OGC:WMTS(1;h)&l4=ORTHOIMAGERY.ORTHOPHOTOS::GEOPORTAIL:OGC:WMTS(0.25;h)&d5=3573961(1)&d6=3573952(1)&permalink=yes, last access: November 2025). In the following paragraphs we highlight the main changes and their potential impacts on the watershed hydrology.

The oldest document available for our study area is the land register established in 1833, commonly referred to as “adastre Napoléonien” (Ministère des Finances, Direction départementale des contributions directs, Bureau du cadastre, 2007). This register consists of two layers, one that displays the limits of plots of land based on ownership and one that represents land use. The later has unfortunately been lost for most municipalities in the “département” (french administrative division, similar to a county) where our study area is located. However, the military survey map established between 1820 and 1866 (IGN, 2021) indicated that most of the study area was open land (cropland, grassland or heath), apart from a wooded area at the south-west border of the catchment (most upslope sampling location in T1) and localised patches of trees, notably around the upslope area of T3. The stream was surrounded by small wetlands in most of the study area. Comparing the 1833 land register and its updated version of 1952 (Direction départementale des finances publiques, 2020), we observed that plot borders have remained the same over the period. Both registers are publicly available on the website https://patrimoines-archives.morbihan.fr/ (last access: November 2025) of the département du Morbihan.

The aerial photograph taken in 1952 (IGN and GEOPAL, 2011) revealed the presence of hedgerows at most plot borders in the area of T1 and T3 whereas only a sparse hedge separated the top-most sampling site in T2 from the rest of the transect. The wetlands around the stream appeared to have been taken into cultivation whereas the wood at the south-west border was still present (and still is to this day).

Figure A1 Elevation profiles derived from the 1 m resolution digital elevation model RGE ALTI^® 1M (IGN, 2013) along transects 1 to 3, with location of the sampling points (red crosses), hedges (circles), and roads or paths (diamonds). The evolution of obstacles (roads and hedges) through time is indicated by colours: green – 2021 (observed), pink – 1952 (derived from IGN and GEOPAL, 2011 and Direction départementale des finances publiques, 2020), orange – 1850's (derived from IGN, 2021, and Ministère des Finances, Direction départementale des contributions directs, Bureau du cadastre (2007)).

For the most part, hedgerows were removed in the 1970's as part of a government driven process to modernise agriculture via the use of heavier machinery on larger surfaces. Several path cutting through T1 and T3 were abandoned in the process but the roads separating the upslope sites in T2 and T3 from the summit remained in place (Fig. A1). These changes likely increased the possible path length for erosion fluxes along these transects. Drains have been installed in several of the fields cross-cut by T2 in the 1990's, significantly reducing the waterlogging in the midslope to toeslope section of the transect (C. Walter, personal communication, 2021). Agricultural practices in most of the study area since the 90's consist of a rotation of maize, cereals and sometimes vegetables. Aerial photographs from 1993, 1999, 2004, 2009, 2010, 2013, 2016 and 2019 show that cultivation was systematically carried in the direction of the slope at the plots where transect 2 and 3 are located, whereas the field in the upper part of transect 1 is cultivated perpendicular to the slope. The three sampling sites in the lower part of transect 1 are located in a field that had been a cultivated grassland since at least 1993 before being turned to maize in 2020. The wetland around the stream has been left fallow or used as grassland, with recent tree plantation in part of the buffer zone.

Hedgerow removal associated with in-situ residue and/or stump burning could have generated localised PyC inputs. Although some of our sampling points are close to former plot boundaries, we have found no charcoal accumulation layer that could corroborate such practice. Fire has been close to inexistent in croplands of north-western Europe (and in particular in Brittany) in the 21st century Forrest et al. (period 2002–2014,  2024, Fig. 4). Stubble burning was likely a common practice in Western Europe until the 1980's but we found no data that would allow us to characterize it and estimate the PyC inputs it represented. We assumed that, if it took place in our study area, stubble burning produced small amounts of relatively labile PyC due to low fire temperature (< 350 °C,  Rasmussen et al., 1986) and the small amounts of biomass involved (part of the biomass was exported before burning). This hypothesis is supported by the old radiocarbon age of the PyC^HyPy in our study area (see main text Sect. 3.4 and 4.1).