Numerical modelling was performed in R Core Team (2021) and RStudio Team (2020) by adapting the script developed by Yu and Flury (2021). In there, the three-dimensional soil matrix is imagined as two-dimensional space, biochar particles were aligned in clusters representing row or tree circumference application on a 50 × 50 m field (0.25 ha). Each particle had a specific location and expansion (surface area) on the field. Biochar particles were allocated to the clusters, e.g., a row, and located within these following a random uniform distribution (Figure 2). Soil cores (diameter = 80 mm) were allocated numerically to represent the three sampling approaches as detailed below (chapter 2.3).

To assess whether a biochar particle is numerically sampled, the coordinates of the soil core expansion and the biochar particle were compared, if the center of the biochar particle was within the boundaries of the soil core, the particle was regarded as within the sample. For each soil core, the number (n) of individual biochar particles (i = [1, n]) and the sum of their surface areas ( A = ∑ i = 1 n π r i 2 ) was calculated, considering the radius (r_i) of the individual biochar particle. The mean biochar concentration for a random sequence of a defined number of cores was calculated. The order of soil cores in the sequence was randomized. Each simulation was carried out 100 times. It was assessed whether the mean biochar concentration of the sequence of soil cores represented the biochar concentration of the field (± 5% error). Within the 2D field model, biochar concentrations are calculated as surface area (biochar) per surface area of the field. The true biochar concentration of the field is given by the sum of the areas of all biochar particles ( A = ∑ i = 1 n π r i 2 in m²) divided by the field surface area (2,500 m²). Within the scope of this paper, we consider representative sampling was achieved if the mean SOC concentration, defined as the sum of the biochar-C concentration and the native SOC concentration (see 2.4) of the soil samples (sequence of soil cores) was within a margin of ±5% of the true SOC concentration of the field (sum of biochar-C and native SOC concentration).

The biochar application rate was 5 t ha⁻¹ for all simulations. Biochar particles were assumed to be of spherical shape with diameters between 2 mm to 4 mm (following a uniform distribution), with a porosity (ϕ) of 0.4 m³ m⁻³ (Hagemann et al., 2020) and a true biochar density (ρ_T) of 2 t m⁻³ (Brewer et al., 2014). Total volume of biochar per granule without porosity ( V = 4 3 π r 3 × ( 1 − ϕ )) and its respective mass ( m = V × ρ T ) was calculated. Further the number of biochar particles and the surface area of the biochar circles representing the biochar application rate in a two-dimensional space was calculated. In total the application rate was represented by 73.7 million biochar particles with a mean surface area of 7.33 mm² per particle on the 50 × 50 m field. In total, in this 2D projection, biochar covered 166.67 m² of the 2,500 m² (0.25 ha) plot (6.67%).

Biochar application was either modelled as row application or tree circumference application (Figure 2). For row application, the row width was 5 cm with a row distance of 75 cm between the centers of the rows, unless specified otherwise. For tree circumference application, a planting density of 1,667 trees ha⁻¹ following a diamond (quincunx) planting pattern was assumed. For circumference application, biochar was applied in a 15 cm wide circle around the tree (radius = 1 m, i.e., biochar was applied within 0.925–1.075 m distance to the center of the tree). As controls, models were run both without biochar as well as assuming a random uniform (even) biochar application over the whole plot.

Three sampling strategies were modelled:

Biochar carbon content (C_BC) was assumed to be 0.7 kg kg⁻¹ and biochar carbon per soil core (C_BC,core) was calculated based on the number and surface area of the particles. Soil characteristics were assumed as detailed in Table 1. For non-biochar SOC, each core contained a SOC of SOC_soil1 = 6.1 g kg⁻¹ and SOC_soil2 = 22.2 g kg⁻¹ representing a low and high C-stocks, respectively, assuming a perfectly homogenous distribution of SOC. To account for natural heterogeneity, non-biochar SOC (SOC_stock) content was calculated as follows:

Here SOC_hom, is the respective (homogenous) non-biochar SOC content of each soil (Table 1), CV_SOC, the normally distributed coefficient of variation (mean = 0.088, sd = 0.048) of the SOC-content, FSD_soil, the respective fine soil density (Table 1) and CV_FSD the normally distributed coefficient of variation (mean = 0.065, sd = 0.015) of the FSD. Soil characteristics and variation coefficients were taken from an assessment of potential resampling errors of SOC in German croplands by Poeplau et al. (2022).

According to our model, the number of soil cores required for a representative sampling of native SOC-stocks, without biochar amendment (“zero control”), depends on the sampling strategy. Following the random soil sampling strategy, at least 15 soil cores per plot (0.25 ha) were required for a representative sampling in at least 90% of repeated sampling simulation cases, for both soil types. To achieve a representative sampling of at least 95 and 100% of simulations, sampling density of 17 and 61 randomly distributed soil cores per plot, respectively, were required (Figures 3A,B).

Following the AFRS sampling, one single composite sample (comprising 25 soil cores) per 0.25 ha led to a representative SOC-stock assessment in 99% of repeated sampling simulation cases for both soil types (Figures 4A,B). Noteworthy, we calculated the native SOC-stock based on Eq. (1) for each soil cores independently, spatial factors were not considered, thus, factors such as expected lower heterogeneity between cores spatially close to each other were neglected. These results were basically defined by the statistical data, particularly the coefficients of variation, we adopted from Poeplau et al. (2022). Similar coefficients of variation were also described by other authors (Saby et al., 2008; Conant et al., 2011) for agricultural soils elsewhere, underlining that the data utilized in the model represents the SOC heterogeneity of many agricultural sites adequately.

Uniformly applied biochar (following a uniform distribution) of 5 t ha⁻¹ (“positive control”) reduced the amount of soil cores required for a representative sampling in at least 90% of the repeated sampling simulations to 13 for soil 1 (low native SOC), while it remained at 15 soil cores for soil 2 (the high native SOC soil type) for the random sampling. The reduction of the threshold sample size for a representative SOC assessment by uniformly applied biochar in the soil with low non-biochar SOC-stocks (soil 1) was the result of the added biochar-C (3.5 t ha⁻¹) accounting for 11.7% of the total SOC (native SOC+ biochar-C). As this (substantial) share was distributed uniformly, it de facto facilitated representative sampling in the model. Regarding soil 2, the evenly applied biochar did not affect the number of randomly distributed soil cores, as the biochar-C added (3.5 t ha⁻¹) only represented 3.8% of total SOC.

Taking the soil samples numerically based on the AFRS sampling method, one single composite sample was sufficient to represent the SOC-stock after uniform application of biochar in 99 and 100% of repeated sampling simulation cases for low and high baseline SOC-stocks, respectively.

Soil sampling guidelines commonly propose to first calculate a mean and standard deviation based on few samples randomly drawn across the plot. Based on these and the difference in means that the sampling aims to detect, the minimum number of required samples could be calculated (Aynekulu et al., 2011; FAO, 2020). Following this approach, and taking ten random samples before the biochar application, we found a standard deviation of 4.7 t ha⁻¹ for soil 1 with a native SOC-stock 26 t ha⁻¹ and a standard deviation of 15.8 t ha⁻¹ for soil 2, with a native SOC-stock of 88 t ha⁻¹ (Table 1), corresponding to a relative standard deviation of 18% in both cases. Aiming to detect the 3.5 t ha⁻¹ of biochar-C applied this would lead to recommendations of 30 and 320 samples, respectively (considering a power of 0.8, cf. SI, section 2). However, native SOC-content heterogeneity can be considerably higher than assumed here, for example Heim et al. (2009) showed that sampling densities of 27 composite samples (of 25 cores each) were required to detect a variation in SOC of 2.6 t per hectare (equaling to about 3.9% of SOC-stocks).

Considering a baseline SOC-stock of 26 t ha⁻¹ (soil 1) and 88 t ha⁻¹ (soil 2), random sampling with a sampling density of 100 cores per plot after biochar row application led to a representative sampling of SOC-stocks (± 5% error) in 85 and 100%, respectively (Figures 3A,B). For the soil 2, the assumed biochar-C application (3.5 t ha⁻¹) accounted for 3.8% of SOC and it was smaller than the tolerated error margin of 5% (4.4 t ha⁻¹), thus, the biochar induced heterogeneity was consistently low in this soil (Table 2). In contrary, for the scenario with baseline soil 1 the biochar-C amendment (3.5 t ha⁻¹) accounted for 11.7% of SOC, hence the spatially discrete distribution of this share hampered representative sampling severely. Baseline native SOC-stock and relative biochar-C application rate are hence pivotal factors regarding the heterogenous biochar application’s impact on representative sampling. Here, we only quantify the implications for two contrasting SOC-stocks of agricultural sites in Germany. However, the SOC-contents given for soil 1 (6.1 g kg⁻¹) and soil 2 (22.2 g kg⁻¹) may represent other soils suitable for biochar application elsewhere. Martinsen et al. (2014) report successful biochar root-zone application for soils with SOC-contents (5.4 to 6.2 g kg⁻¹) similar to “soil 1” in Zambia. Yield increases for concentrated biochar root-zone application on Nepalese soils reported by Schmidt et al. (2015) and Schmidt et al. (2017) took place on soils with a wide range of SOC- from 15 to 54 g kg⁻¹), with many sites showing SOC-contents similar to “soil 2.”

In general, biochar induced heterogeneity (due to heterogenous biochar application) is larger if biochar accounts for a high share of total SOC, leading to high sampling densities required for representative assessments.

High sampling densities of 100 soil cores per 2,500 m² may be used in research, e.g., in long-term field experiments (Schrumpf et al., 2011). For lower, more economical sampling densities, e.g., 40 randomly distributed soil cores per hectare (10 soil cores per plot), 47 and 68% of repeated sampling simulation cases were representative for the low and the high baseline SOC, respectively (Table 2). For native SOC alone the same sampling density led to representative assessments of SOC in 75% of simulations (Table 2).

To achieve representative SOC sampling in at least 90% of repeated sampling simulation cases, a density of 100 randomly distributed soil cores was not sufficient for soil 1 (Table 2); at least 140 soil cores per plot were required. Thus, the number of required randomly taken samples increased by a factor of 9.3, compared to the model for the native SOC (from 15 to 140; Figure 3A). For soil 2 (Figure 3B), a sampling density of more than 23 randomly distributed soil cores per plot was sufficient for this purpose. The variability induced by biochar row application was bigger for low native SOC, as the biochar-C amendment represented a higher share of total SOC. For soil 2, where the biochar-C amendment represents a small fraction of total SOC (3.8%), soil sampling density had to be increased by 53% (from 15 to 23 soil cores per plot) to achieve representative results in 90% of repeated sampling simulations (Figures 3A,B).

To achieve representative sampling in 100% of repeated sampling simulation cases, the minimum distance of 3 m between soil sampling points had to be abandoned, as this would have limited the number of soil cores to a theoretical maximum of 256 soil cores per plot. For soil 1 and soil 2, at least 405 and 152 cores per plot (1,620 and 608 ha⁻¹), respectively, had to be numerically assessed to quantify the SOC-stock in 100% of simulations.

Ten composite samples distributed according to the AFRS sampling method, i.e., 250 soil cores, reflected the SOC stock after biochar row application correctly (5% error margin) in 97% of all cases for soil 1. The same SOC-stock, achieved by uniform biochar application could be representatively quantified with a single composite sample in 99% of simulations (positive control). Hence, for a soil with a low baseline SOC-content, the measuring effort might need to be increased tenfold, to gain equally representative results after concentrated, i.e., spatially discrete biochar applications (Figure 4A). For soils with high baseline SOC-contents, the number of composites samples required to be taken, following the AFRS sampling scheme with post biochar row application was tripled to achieve representative SOC quantification in 99% of simulations (Figure 4B). To assess the SOC content representatively in 90% of repeated sampling simulations, 7 composite samples (175 soil cores) and one composite sample (25 soil cores) were sufficient for soil 1 and soil 2, respectively.

Especially for SOC depleted soils, or soils with naturally low SOC contents, biochar C amendments can contribute significantly to the overall carbon stock and the biochar’s heterogenous distribution hampers representative sampling of overall SOC-stocks. It must be noted that such soils are very relevant in this context, as the use of biochar in marginal soils is particularly interesting from an agronomic point of view (Melo et al., 2022).

Comparing biochar row application and biochar tree circumference application, we found that the type of spatially heterogenous biochar application does not impact the biochar induced variability significantly (Table 2, Figures 3, 4). Similar to row application, the induced variability was higher for low native SOC-stock and minimal to absent in the soil with high native SOC-stock (Table 2).

We tested the effect of different distances between rows of biochar application (0.5 m, 0.75, 1 m, 1.5 m), which did not have a significant effect (based on Kruskal-Wallis test) on the assessed SOC concentrations (cf. SI, Supplementary Figure S4). A tendency towards more accurate SOC estimations with smaller row distances could be observed. The main factor for this effect is most probably that larger row distances lead to higher concentration of biochar particles in each row (more heterogenous distribution) as long as the per hectare application rate is kept constant. Thus, when biochar particles are present at a higher concentration on a smaller area within the plot, the increased heterogeneity further impedes representative sampling.

Furthermore, we tested the effect of biochar application being numerically carried out as multiple smaller biochar applications with conventional agricultural equipment (non-precision farming). This led effectively to a larger area covered with lower concentrations of biochar as we assumed biochar in the model to be applied with a random scatter (0 to 10 cm). Therefore, the overall area covered with biochar is larger and hence multiple application facilitated representative sampling (cf. SI, Supplementary Table S1; Supplementary Figure S3). Still, after 5 years of applying 1 t biochar ha⁻¹ yr.⁻¹ in rows, the number of cores in random sampling had to be increased by 260% (from 15 to 54 cores per plot) and by 27% (from 15 to 19 cores per plot) for soil 1 and 2, respectively, to achieve representative SOC quantification in at least 90% of simulations in our model. We understand precision farming, taking advantage of real-time kinematic positioning, could achieve repeated biochar root-zone application with accuracy at the cm-scale, but this technology is still not widespread. Also, non-accuracy in repeated root-zone application will create more soil volume that is in touch with biochar and thus might be beneficial for overall soil improvement in the long-term, however, we are not aware of any field study looking into this question.

Within the models we considered biochar particles to remain at the location of original application. Potential vertical mobility of biochar particles (Schiedung et al., 2020) through leaching or biological activity is captured indirectly as vertical mobility does not alter the results of the 2D model. Potential lateral mobility of biochar particles through biological activity, erosion and soil tilling were neglected here, as the extent of such dispersion could not be quantified reliably. Also, we assume the tillage is commonly reduced to, e.g., the upper 10 cm, while biochar for root-zone application might be placed below this zone.

To improve the sampling approach, i.e., to reduce the number of individual samples needed to achieve representative sampling, we tested if knowledge on the mode of biochar application could help to optimize the sampling strategy.

When the orientation and spacing of biochar row application is known, we suggest the following approach: In the scenario with 75 cm interrow space, aligning 10 soil cores (diameter = 80 mm) as described in 2.3. With 10 soil cores per plot distributed in such manner, the SOC-stock can be assessed representatively in 82 and 78% of repeated sampling simulation cases after biochar row application for soils with low and high native SOC-stocks, respectively (Figure 5). In our model the biochar distribution within the rows is uniform along the row (Y-Axis). In the real-world biochar distribution within the row might be less uniform. Furthermore, the native SOC’s heterogeneity was accounted for numerically without accounting for spatial effects (SOC was calculated for each core, irrespective of its location), such as potentially lower heterogeneity with closely aligned cores and higher potentially heterogeneity with larger distances. Selecting multiple rows, instead of a single row might help to account for these effects, that were neglected in the model. While a higher number of soil cores should be considered for unknown biochar patterns, if the biochar application pattern is well-known, a cost-effective highly representative SOC assessment with much lower biochar induced variability, especially in soils with low native SOC-stocks, can be achieved with the adapted soil sampling pattern.

To limit costs and effort of soil sampling without sacrificing accuracy in SOC-stock assessments after heterogenous biochar applications, such as biochar row applications, it is mandatory to precisely record where biochar (row) application was carried out, and which pattern was followed.

Biochar application may lead to positive priming (increased decomposition of non-pyrogenic SOC), though over time it overall may result in negative priming, meaning reduced decomposition and net increase of non-pyrogenic SOC (Ding et al., 2018; Schmidt et al., 2021). Positive priming is the result of increased mineralization due to enhanced microbial activity, e.g., due to biochar supporting microbial metabolisms as electron shuttle and geo-conductor (Kappler et al., 2014; Sun et al., 2017; Pokharel et al., 2020). In turn, this increased activity and accelerated nutrient cycling may increase plant productivity and input of plant derived organic carbon to the soil, e.g., by increased root growth and root exudates (Weng et al., 2017; Xiang et al., 2017). This leads to a higher SOC equilibrium over time, though the scale of biochar induced negative priming is not yet consistently determined (Schmidt et al., 2021), and the overall extent of the microbial response to biochar application and thus the relevance of biochar induced positive and negative priming strongly depend on the environmental factors and biochar properties (Pokharel et al., 2020). For instance Blanco-Canqui et al. (2020) found that 6 years after the application of 7 t biochar per hectare, the total SOC increased (i.e., negative priming) by 14 t ha⁻¹. The biochar carbon led thus to an additional accumulation of 7 t SOC ha⁻¹ from plant derived organic carbon. The effect found by Weng et al. (2017) in an Australian field trial was smaller but still significant. Priming might occur in close proximity to biochar amendments, as, e.g., the microbes need direct interaction with biochar surfaces to achieve electron exchange (Kappler et al., 2014). In contrast, biochar stimulates root growth (Xiang et al., 2017) which is not expected to occur solely within the biochar amended area, but throughout the soil in unknown, but most likely heterogenous patterns. Although any priming effect induced by biochar are beyond the scope of this work, we expect these effects to further complicate an effective representative SOC sampling of biochar amended agricultural soils.