Quantifying inertinite carbon in biochar

Quantifying inertinite carbon in biochar Hamed Sanei a,* , Małgorzata Wojtaszek-Kalaitzidi b , Niels Hemmingsen Schovsbo c , Rasmus Stenshøj a,c , Zhiheng Zhou a , Hans-Peter Schmidt d , Nikolas Hagemann d , David Chiaramonti e , Tryfonas Kiaitsis a , Arka Rudra a , Anna J. Lehner f , Robert W. Brown g,h , Sophie Gill g , Erica Dorr i , Stavros Kalaitzidis j , Fariborz Goodarzi k , Henrik Ingermann Petersen c a Lithospheric Organic Carbon (LOC), Department of Geoscience, Aarhus University, Denmark b Instytut Technologii Paliw i Energii (ITPE), Poland c Geological Survey of Denmark and Greenland (GEUS), Denmark d Ithaka Institute for Carbon Intelligence, Switzerland e Politecnico di Torino and RE-CORD, Italy f Carbonfuture GmbH, Germany g Isometric, London, United Kingdom h School of Environmental and Natural Science, Bangor University, Wales, United Kingdom i Rainbow, Paris, France j Department of Geology, University of Patras, Greece k FG & Partner Ltd, Research Group, Calgary, Alberta, Canada ARTICLE INFO Keywords: Biochar permanence Inertinite benchmarking Random reflectance (R o ) Carbon dioxide removal (CDR) Sampling frequency ABSTRACT The carbon dioxide removal (CDR) potential of biochar is determined by the long-term stability of its biogenic carbon, derived from atmospheric CO₂ fixed by photosynthesis and stabilized in solid form. This stability (carbon permanence) is commonly assessed using decay models to evaluate resistance to re-emission as greenhouse gases. However, these models are limited, as they focus primarily on short-term degradation of labile carbon fractions and are not suited to project the behavior of the highly recalcitrant component of biochar over extended timescales. Inertinite represents highly aromatized and condensed carbon structures that are geochemically stable over millennia. This paper builds upon the Inertinite Benchmarking (IBR o 2) methodology, directly quantifying the stable carbon fraction in biochar rather than relying on modeling. The method combines thermochemical analysis and incident-light microscopy to measure the reactive (labile) component and solid carbonized macerals, respectively. Random reflectance analysis (R o ) provides a representative distribution of carbonization states, with R o values >2.0 % defining the inertinite fraction after discounting reactive organic carbon. The R o distribution is processed using kernel density estimation (KDE) and numerical integration to classify inertinite carbon with precision and statistical robustness. As CDR crediting can be linked to measured inertinite content, statistical validity is essential. A Monte Carlo simulation model evaluates uncertainties from sampling frequency and production variability. Results show that increased sampling reduces uncertainty and lowers the conservative safety margin needed for potential errors. This framework supports a justified safety margin applied to reported inertinite carbon and corresponding CDR values, enabling conservative and robust crediting. By combining direct quantification of inertinite carbon with probabilistic modeling of uncertainty, the IBR o 2 method offers a transparent and rigorous framework for assessing biochar permanence, aligned with emerging international certification and national inventory methodologies. * Corresponding author. E-mail address: [email protected] (H. Sanei). Contents lists available at ScienceDirect International Journal of Coal Geology journal homepage: www.elsevier.com/locate/coal https://doi.org/10.1016/j.coal.2025.104886 Received 7 August 2025; Received in revised form 24 September 2025; Accepted 2 October 2025 International Journal of Coal Geology 310 (2025) 104886 Available online 3 October 2025 0166-5162/© 2025 Published by Elsevier B.V. 1. Introduction The effectiveness of biochar as a means of carbon dioxide removal (CDR) depends on the permanence of its stored carbon. For soil-applied biochar, common assessment methods rely on closed-system laboratory incubations over short timeframes (often one to three years) (Budai et al., 2016; Dharmakeerthi et al., 2015; Fang et al., 2014; Herath et al., 2015; Kuzyakov et al., 2014; Major et al., 2010; Singh et al., 2012; Wu et al., 2016; Zimmerman, 2010; Zimmerman and Gao, 2013). While these experiments primarily capture the short-term decay of labile organic matter, the highly carbonized, stable fraction remains largely unreactive (Gross et al., 2024; Sanei et al., 2025). As a result, derived decay models reflect only initial degradation processes and offer no direct empirical insight into the long-term fate (>100 years) of the stable carbon fraction, making extrapolations highly uncertain (Azzi et al., 2024; Sanei et al., 2025). A complementary approach addresses this limitation by directly quantifying the fraction of stable carbon expected to persist over environmentally relevant timescales (>1000 years). This method draws on geological principles, recognizing that highly carbonized organic matter transforms into inertinite macerals, condensed and aromatized structures known for their resistance to degradation and long-term preservation in sedimentary rocks (Ascough et al., 2011; Ascough et al., 2010; Azzi et al., 2024; Hudspith et al., 2015; Petersen et al., 2023; Sanei et al., 2024). Most inertinite macerals represent terminal transformation state of organic matter within the sedimentary carbon cycle, analogous with the endpoint of inorganic carbon becoming carbonate rock. From a geological perspective, inertinite macerals and carbonate minerals represent two complementary pathways of the Earth’s natural long-term carbon sequestration system. Both serve as mechanisms by which atmospheric or biospheric CO 2 is transformed into rock-forming components. These transformations have over geological timescales contributed to the regulation of the Earth’s climate, acting as natural thermostats during periods of elevated atmospheric CO 2 . The carbonization of plant-derived organic matter into inertinite macerals, in particular, is a well-documented and geochemically validated mechanism by which the terrestrial biosphere has contributed to geologicalscale carbon storage throughout Earth’s history (Kroeger et al., 2011; Morga, 2011; Sanei et al., 2024; Tissot and Welte, 2013). The Inertinite Benchmarking methodology (IBR o 2) builds on this principle by proposing direct measurement of inertinite carbon content (C Inert ) in biochar, which represents the portion of organic carbon that can be assumed to remain stable over millennial timescales. The recognition and quantification of the inertinite fraction within organic matter is well established in the field of organic petrology (Diessel, 1983; International Committee for Coal and Organic Petrology (ICCP), 2001; Scott and Glasspool, 2007; Morga, 2011; Petersen et al., 2023; Sanei et al., 2024; Mastalerz et al., 2023, 2025). This methodology, as applied to biochar, was described by Sanei et al. (2024) and further developed in subsequent studies (Mastalerz et al., 2025; Petersen et al., 2025; Petersen and Sanei, 2025; Rudra et al., 2024). The method was also tested in the longest-running biochar field trial in the European Union, located at La Braccesca, Italy, where Chiaramonti et al. (2024) evaluated the inertinite fraction in topsoil after 15 years. These have contributed to the refinement and streamlining of the methodology, enabling its broader application in carbon crediting frameworks and its integration into standards used by carbon registries and monitoring, reporting, and verification (MRV) systems. While further fundamental work is needed to reconcile the maceral composition with the findings from two decades of research on carbon speciation in biochars into a coherent picture to further advance the topic of biochar persistence. However, the aim of this paper is to consolidate these developments and formally establish a standardized methodology for inertinite quantification in biochar. Building on prior research and practical application, the proposed method is presented as a best-practice approach for assessing the long-term carbon permanence of biochar in the context of CDR certification. The IBR o 2 method enables certification schemes and their certification bodies to move to an advanced analysis method that allows a more nuanced and fit-forpurpose characterization of individual biochars regarding their persistence in soil. 2. Inertinite Benchmarking (IBR o 2) The IBR o 2 method integrates two complementary analytical approaches, (i) thermochemical analysis and (ii) incident light microscopy, to quantify the reactive organic carbon (C React ) content of biochar (Fig. 1). This is because organic matter in biochar is present in two primary forms: (i) Reactive organic matter includes more labile compounds, consisting of secondary-generated or adsorbed material during pyrolysis such as condensates, semi-liquid tars or bituminous substances, and small carbonaceous compounds at the fringes of fused aromatic macromolecules as well as residual noncarbonized material. These are more prone to degradation but are not observable under incident light microscopy and cannot be quantified by light microscopy as introduced in the following paragraph. However, they are thermally reactive and can be volatilized through controlled re-pyrolysis. Reactive organic carbon content (C React ) can be quantified using thermochemical methods such as Rock-Eval 6 pyrolysis, thermal gravimetric methods, or other simmilar analysis (see Section 3). (ii) Macerals are remaining solid and semi-solid, particulate organic matter amenable to identification and quantification using incident light microscopy (visible light), where random reflectance (R o ) can be measured. Because R o measurements require mechanical polishing of the sample, only macerals with their solid and semi-solid nature are observable and quantifiable under the microscope (see Section 4). To accurately determine the C Inert , both carbon fractions must be quantified using separate but integrated methods ((Petersen et al., 2025; Petersen and Sanei, 2025; Sanei et al., 2024); see Section 2.3). The biochar sample is therefore divided into two representative subsamples; subsample A is subjected to thermochemical analysis to quantify the C React (dry wt%), and subsample B undergoes incident-light microscopy to determine the R o distribution of the carbonized macerals (Fig. 1). It is critical that the sample accurately represents the entire production batch to ensure analytical accuracy and precision. Detailed guidance on appropriate sampling procedures is provided in the guidelines of the European Biochar Certificate (Schmidt et al., 2024). The use of grab samples alone can introduce high variability (Bucheli et al., 2014), which increases the risk of elevated safety margin deductions during credit issuance (see Section 5). The procedure follows these steps (Fig. 1): 2.1. Quantify reactive organic carbon fraction Using subsample A, the C React can be measured thermochemically (see Section 3). This value, expressed as a fraction of organic carbon (C Org , obtained from elemental analysis, cf. EBC 2025 (Schmidt et al., 2024) present in biochar, is defined as: FReact =CReact COrg (1) Where C React and C Org are both expressed as weight percentages on a dry basis (dry wt%) and FReact is expressed as a fraction of reactive organic carbon (i.e., values between 0 and 1). H. Sanei et al. International Journal of Coal Geology 310 (2025) 104886 2 2.2. Calculate inertinite fraction In subsample B, reflectance can be measured at 500 points across the polished surface (see Section 3). The relative proportion of measurements with R o >2.0 % is denoted as FRo>2, representing the fraction of data points exceeding this threshold (Mastalerz et al., 2023, 2025; Petersen et al., 2023; Sanei et al., 2024). FInert =FRo>2×(1−FReact)(2) Where FInert and FRo>2 are both expressed as fractions (i.e., nondimensional values between 0 and 1). CInert on a dry weight basis is then calculated by, based on assumptions detailed in Section 2.3: CInert (dry wt%) = FInert ×COrg(dry wt%)(3) To quantify the CDR, the CInert is converted to its CO 2 -equivalents using the molar mass ratio of CO 2 /C (44.01/12.01): CDR (wt%) = CInert (dry wt%)×44.01 12.01 (4) This integrated methodology provides a direct, empirical measure of the mass fraction of C Org in biochar that has transformed into an inertinite-like structure. By isolating and quantifying this fraction, the IBR o 2 method offers a scientifically robust framework for assessing longterm biochar permanence in CDR applications. 2.3. Combining thermochemical and R o methods to estimate inertinite carbon To quantify CInert (dry wt%) using FReact and FRo>2requires harmonizing differences in the units. The FReact is estimated as the mass ratio of carbon released during thermal decomposition, normalized to the C org . In contrast, the measured FRo>2, derived from the volumetric distribution of R o values, is based on the well-established point-counting method in organic petrology and represents an equally probable sampling of macerals present throughout the sample (International Organization for Standardization (ISO), 2009b; Gordon et al., 2021; Sanei et al., 2024; Zhou and Sanei, 2025). Therefore, the resulting frequency histogram of R o values provides a statistically representative approximation of the volume fraction of biochar at different R o classes: (i) poorly carbonized: FRo≤1.2; (ii) semi-inertinite: F1.2<Ro≤2, and (iii) inertinite: FRo>2 (Sanei et al. (2024); see Section 4). Combining the datasets reported in different units requires conversion of R o -based volume fractions into weight percentages. Because direct measurements of maceral grain densities are not practical, the analysis assumes that all three carbonization classes (poorly carbonized, semi-inertinite, and inertinite) have identical densities. With this assumption, volume fractions are treated as equivalent to weight fractions, allowing for direct comparison without the need for density correction. In both formats, the total remains normalized to 100 %. In reality, inertinite macerals have a higher density than other carbonaceous components due to their elevated carbon content, which results from increased aromaticity and lower volatile matter (Wang et al., 2024; Wang et al., 2023). Consequently, the common assumption of equal density across all maceral types leads to an underestimation of the inertinite contribution. For example, if the average grain densities for poorly carbonized, semi-inertinite, and inertinite components are approximately 1.30, 1.40, and 1.50 g cm −3 (Wang et al., 2023, 2024), respectively, and each class occupies one third of the total volume, the resulting weight distribution would be 30.9 %, 33.4 %, and 35.7 %. Thus, the actual inertinite weight fraction exceeds the estimate based on equal-density assumptions by about 2.4 percentage points. This Fig. 1. Schematic protocol for quantifying the inertinite fraction of organic carbon in biochar using the inertinite benchmarking method (IBR o 2). The dry biochar sample is divided into two subsamples: (A) analyzed thermochemically to quantify the reactive organic carbon fraction (F React ) fraction, and (B) analyzed petrographically via reflected light microscopy to obtain random reflectance (R o ) values. Based on R o thresholds, carbonized particles are classified as poorly carbonized (R o ≤1.2 %), semi-inertinite (1.2 % <R o ≤2 %), and inertinite (R o >2 %). The inertinite fraction of organic carbon in biochar (F Inert ) is calculated as the proportion of macerals with R o >2 % (F Ro>2 ) multiplied by the non-reactive organic carbon fraction (1 – F React ). H. Sanei et al. International Journal of Coal Geology 310 (2025) 104886 3 simplification introduces a conservative bias into the carbon accounting framework, reducing the risk of overestimating the stable carbon fraction. 3. Thermochemical measurement of reactive organic carbon The quantification of C React content in biochar can be achieved by using thermochemical analysis such as Rock-Eval 6, thermal gravimetric methods, or other similar analysis (Bordenave et al., 1993; Buss and Maˇ sek, 2014; Lafargue et al., 1998; Petersen et al., 2023). These methods aim to re-pyrolyze the biochar sample under controlled conditions to volatilize labile carbon bonds. The mass or chemical signature of these volatile components is then used to estimate the fraction of reactive organic matter susceptible to thermal decomposition below a functionally defined temperature. In thermal gravimetric methods, the reactive organic matter fraction (or volatime matter) is estimated gravimetrically by monitoring the sample’s weight loss during pyrolysis. In contrast, Rock-Eval analysis directly measures the evolved hydrocarbons and oxygen-containing compounds using flame ionization detection (FID) and infrared (IR) spectroscopy, respectively. Since approximately 95 % of the organic molecules in biochar are composed of carbon, hydrogen, and oxygen, the sum of measured hydrocarbons, CO, and CO₂ released during pyrolysis accounts for the majority of the reactive organic matter in the biochar sample. From these evolved gases, the C React content can be stoichiometrically estimated. These methods differ in their potential biases. Rock-Eval method may underestimate C React , as it does not account for minor contributions from heteroatoms such as sulfur, nitrogen, phosphorus, and trace elements. In contrast, thermal gravimetric methods may overestimate the reactive fraction, since the gravimetric signal can include weight loss from thermally unstable mineral phases, such as siderite, that decompose within the pyrolysis temperature range. Furthermore, factors such as the heating rate, maximum temperature, and residence time at the maximum temperature can also influence the measured results. Rock-Eval 6 methodology is calibrated against Institut Français du P´ etrole (IFP) standards and validated to recover the total organic carbon (C Org ) content. The pyrolysis stage involves heating the sample isothermally at 300 ◦C for 3 min, followed by a linear temperature ramp at 25 ◦C per minute up to 650 ◦C. This step facilitates the thermal degradation of labile carbon structures, releasing volatile products derived from C – H and C – O bonds. The evolved hydrocarbons, CO, and CO 2 are then quantified, and their stoichiometric carbon equivalents are used to calculate the C React content, reported as a weight percentage on a dry basis (dry wt%) (Lafargue et al., 1998). Following pyrolysis, the remaining sample is subjected to an oxidation phase. It is transferred to a combustion furnace, purged with air, and heated from 150 ◦C to 850 ◦C at a rate of 25 ◦C per minute. During this phase, the remaining organic matter is oxidized, and the resulting CO and CO 2 are measured in real time. Their stoichiometric carbon contributions define the residual organic carbon, a fraction analogous to fixed carbon (Lafargue et al., 1998; Sanei et al., 2024). The sum of reactive and residual organic carbon yields the total C Org content of the biochar sample (dry wt%). However, the C Org should still be determined by elemental analysis according to International Organization for Standardization (ISO), 2010, followed by the arithmetic deduction of inorganic carbon (Bachmann et al., 2016; Schmidt et al., 2024). 4. Random reflectance analysis of biochar 4.1. Sampling and sample preparation The biochar sample is first dried at 40 ◦C, then crushed to a particle size below 0.2 mm, allowing internal cross sections of the particles to be exposed for optical measurement. After preparation, the resulting fragments are embedded in a 2.54-cm diameter cold-setting epoxy resin pellet (or similar cold-curing resin). The embedding process is carried out in two stages. In the first stage, a minimal amount of resin is mixed thoroughly with the biochar to concentrate the fragments at the base of the pellet. This is essential because biochar tends to float within the resin due to its relatively low density, and it is the base of the pellet that will later be ground and polished for microscopic analysis. Once the initial layer cures, a second layer of resin is poured to increase the pellet height for easier handling. After complete curing, the base of the pellet is ground and polished following standard procedures outlined in International Organization for Standardization (ISO) (2009a) to obtain a scratch-free, relief-free, highly polished surface. This ensures optimal exposure of biochar fragments in random orientation on the polished cross-section for reflectance measurements. 4.2. Reflectance measurement procedure Reflectance measurements are performed using incident-, whitelight microscopy. An enhanced-contrast, oil immersion 50×objective lens is recommended, and the use of a camera-based photometric system is advised due to its long-term calibration stability and precision. The microscope should be calibrated with reflectance standards near the sample’s expected R o range, preferably using standards with higher values. A preliminary estimate of the sample’s expected R o can be obtained from its reported pyrolysis temperature using published empirical relationships between carbonization temperature and R o (see Fig. 15 in Sanei et al., 2024). Reflectance measurements should be carried out according to International Organization for Standardization (ISO), 2009b, using the smallest available light-probe diameter to minimize the effect of surface imperfections such as scratches and micro-relief, but also to avoid excluding particles with a naturally more delicate, finer structure. The size of the light probe must not exceed the area intended for measurement, in order to avoid unintentional bias from surrounding material such as scratches, debris, or adjacent phases, and to ensure that only the target surface is analyzed. Only carbonized organic matter is targeted for measurement; reflectance values are therefore inherently on an ash-free basis. Proper polishing is essential to ensure the accuracy of R o readings and eliminate surface artefacts. 4.3. Measurement strategy and statistical requirements Each polished pellet should be systematically examined to obtain R o measurements from a statistically representative population of carbonized organic fragments of the biochar sample. It is recommended to collect up to three point measurements on three distinct fragments within each microscope field of view during scanning. Using the central crosshair as a guide, selections should be made from different quadrants of the field to maximize spatial randomness and avoid bias. Following the methodology outlined by Sanei et al. (2024), a total of 500 individual point measurements per sample is advised to ensure both statistical reliability and reproducibility (Fig. 2). Given a maximum of three points per frame, this corresponds to approximately 170 microscopic fields per sample. These fields should be evenly distributed across the entire polished surface to ensure representative coverage of the sample’s internal heterogeneity (Fig. 2). The mean R o can be typically estimated with as few as 100 measurements according to International Organization for Standardization (ISO), 2009b. This estimate is based on the typical standard deviation ( σ ) observed in R o datasets and a 5 % constraint on the uncertainty of the mean value. Applying the standard formula for the Standard Error of the Mean (SE), where SE represents the uncertainty in the mean estimate, the following relationship is used: SE =1.96 σ / (N) √=0.05. Solving this equation provides a rough estimate of the minimum number of R o measurements (N) required to achieve a 5 % uncertainty on the mean, assuming a normal distribution and a confidence level of approximately H. Sanei et al. International Journal of Coal Geology 310 (2025) 104886 4 95 % (Altman and Bland, 2005; Ecampus Ontario, 2022). However, reliable quantification of FRo>2 requires a representative frequency distribution of R o values. This distribution must reflect the volumetric abundance of biochar particles across different carbonization levels and is essential for determining the fraction of highly carbonized material (FRo>2) (Sanei et al., 2024). Recent work by Mastalerz et al. (2025) suggests that for low-ash, laboratory-produced lignocellulosic biochars, the number of R o Fig. 2. Schematic illustration of the recommended protocol for measuring random reflectance (R o ) in a biochar sample. The figure shows systematic scanning of the resin-embedded pellet using a 50×oil immersion objective (total magnification 500×). Up to three R o measurements are recommended to be taken per microscopic field of view (frames), each on a distinct maceral located in a different quadrant of the frame. Scanning continues until at least 500 individual measurements are obtained, typically requiring approximately 170 frames. Fig. 3. Example of a bimodal distribution of random reflectance (R o ) in a biochar sample, showing the classification of carbonization stages based on R o thresholds. The gray bars show the histogram of R o measurements, while the red curve represents the kernel density estimate (KDE) of the R o distribution. Color-shaded regions indicate the three carbonization classes: poorly carbonized material (R o ≤1.2 %) in light blue, semi-inertinite (1.2 % <R o ≤2.0 %) in light green, and inertinite (R o >2.0 %) in light red. The green vertical line indicates the inertinite benchmark at 2.0 % R o (as defined in the IBR o 2 method). The dashed blue line shows the mean R o value (in this case mean R o =3.28 %). Proportions of each class are labeled within the shaded zones (F Ro≤1.2 =6.4 %, F 1.2<Ro≤2 =12.6 %, and F Ro>2 =81.0 %), based on the integrated area under the KDE curve. The bimodal shape reflects heterogeneous thermal conditions during pyrolysis, with a dominant mode in the inertinite range and a secondary mode in the lower reflectance region. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) H. Sanei et al. International Journal of Coal Geology 310 (2025) 104886 5 measurements required to determine FRo>2 may in most cases be reduced to 200 without significant loss of precision. However, industrial biochars typically exhibit higher ash content and more complex, often multimodal, R o distributions. For such materials, no consensus has yet been reached on a reduced measurement count. A global round-robin study currently underway under the auspices of the International Committee for Coal and Organic Petrology (ICCP) is expected to provide further guidance. Until more definitive evidence becomes available, the 500-point measurement protocol remains a conservative and statistically robust standard for reproducible quantification of FRo>2, particularly in heterogeneous industrial biochars, providing the resolution necessary for batch-level crediting and the consistency required for certification purposes (Sanei et al., 2024). Reflectance values are plotted as frequency distributions to illustrate the spatial variability of carbonization within a biochar sample (Fig. 3). The arithmetic mean R o provides a single value proxy for the overall degree of carbonization; however, its interpretive value depends strongly on the characteristics of the underlying distribution. In cases where the distribution is bimodal or polymodal, reflecting heterogeneous carbonization and/or heterogeneity or blending of feedstock materials, the mean R o does not adequately represent the variability within the sample (Petersen and Sanei, 2025) (see the example in Fig. 3). 4.4. Random reflectance (R o ) profiling in biochar The spatial distribution of carbonization levels within a biochar sample is analyzed using the frequency distribution of R o values. The R o data are compiled into a frequency histogram with bin widths dynamically determined to ensure statistically meaningful resolution (Fig. 3). This distribution provides a quantitative basis for assessing the extent and uniformity of carbonization, which is critical for evaluating biochar carbon stability and measuring the inertinite fraction (Sanei et al., 2024). For instance, a critical use of the R o distribution is the determination of the FRo>2 value. This can be done directly using frequency counting. However, this approach can introduce reproducibility issues in samples with complex (bior multimodal) R o distributions due to inter-operator measurement variability. Operators may unintentionally select different populations of macerals and may be drawn to particles with differing reflectivity, introducing slight biases (e.g., favoring loweror higherreflecting particles). This leads to noticeable differences in the resulting frequency distributions and, consequently, significant differences in the calculated proportions of FRo>2. Such variability produces reproducibility issues for non-automated measurements that involve human operators, each with inherent differences in working style and judgment. To minimize the reproducibility impact related to R o distribution analysis, applying kernel density estimation (KDE) is critical. KDE smooths slight inter-operator variabilities, ensuring reproducibility and reliability of fraction estimates. 4.4.1. Kernel density estimation (KDE) of R o distribution To obtain a continuous approximation of the distribution of R o values, a KDE method is applied using a univariate Gaussian kernel. A Gaussian kernel is considered suitable for modeling the R o distribution, as each component within a biochar sample is expected to follow a normal distribution when measured individually. In homogeneous samples, this yields a unimodal distribution, while heterogeneous samples produce a polymodal distribution approximated as a sum of NGaussians, where N represents distinct components. Rather than estimating N beforehand, the KDE method captures the composite distribution directly. KDE is a nonparametric technique that estimates the underlying probability density function f(x)of a random variable based on a finite set of observations, without assuming any predefined distribution shape (Chen, 2017; Parzen, 1962). Given a set of measured R o values {x 1 ,x 2 , …, x n }, the KDE is computed as: f(x) = 1 nh∑n i=1K(x−xi) h(5) Where: •f(x)is the estimated probability density function at point x, •n is the number of R o measurements (i.e., the sample size), •x i are the individual measured R o values, •h is the bandwidth, a smoothing parameter that determines the width of the kernel and controls the balance between bias and variance in the estimate, •K(u) is the Gaussian kernel function, where the variable u is defined as the standardized distance u=x−xi h, and K(u) is given by the standard normal distribution: K(u) = 1  2 π √exp(−u2 2). The bandwidth h plays a crucial role in the shape of the resulting density estimate. A small h may lead to overfitting (a noisy estimate), while a large h may oversmooth the distribution and obscure important features. In this study, the bandwidth is selected using Silverman’s rule of thumb, which provides a data-driven, closed-form solution for optimal smoothing under the assumption of normally distributed data. Silverman’s rule is expressed as Silverman (2018): h=0.9min( σ ,IQR 1.34)n−1/5(6) Where: • σ is the standard deviation of the R o values, •IQR is the interquartile range of the dataset, •n is the number of observations. This formulation adaptively scales the bandwidth to reflect the dispersion of the data while minimizing estimation error in most practical cases. To ensure sufficient resolution across the R o range, the domain x ε [x min ,x max ] is discretized into 500 equally spaced intervals. The resulting KDE curve is then overlaid on the empirical frequency histogram, providing a smooth and continuous representation of the R o distribution (Fig. 3). This approach enhances visual interpretation of the carbonization profile and allows better detection of multimodality or subtle shifts in the distribution that may not be evident in the binned histogram alone. 4.4.2. Numerical integration of the KDE curve To quantify the distribution of R o values across different carbonization categories, the output of the KDE is numerically integrated using Simpson’s rule (Talvila and Wiersma, 2012). This higher-order integration method approximates the area under the KDE curve by fitting a second-order polynomial to each pair of adjacent intervals, providing a more accurate estimate than simpler methods such as the trapezoidal rule. Atotal =∫∞ 0f(x)dx ≈Δx 3[f(x0)+4∑ odd jf(xj)+2+∑ even j j∕=0,nf(xj)+f(xn)](7) In this expression: •A total is the total area under the KDE curve, representing the integral of the estimated probability density function over the entire reflectance domain, H. Sanei et al. International Journal of Coal Geology 310 (2025) 104886 6 •f(xj)are the KDE-estimated density values at the discretized points x j , •Δx is the spacing between consecutive x-values (uniform across the interval), •The summations are taken over odd and even indices j (excluding the endpoints for even j), •x 0 and x n are the lower and upper bounds of the discretized domain, respectively. The total area is then subdivided into three carbonization categories based on R o thresholds (Fig. 3): poorly carbonized (R o ≤1.2 %), semiinertinite (1.2 % < R o ≤2.0 %), and inertinite (R o >2.0 %). The classification thresholds follow Sanei et al. (2024), Mastalerz et al. (2023), Mastalerz et al. (2025), Petersen and Sanei (2025), Petersen et al. (2023) and Petersen et al. (2025), with inertinite representing the most aromatized and condensed fraction, essential for long-term carbon permanence in CDR applications. Each category’s partial area is calculated as: ARo≤1.2=∫1.2 0f(x)dx,A1.2<Ro<≤2.0=∫2.0 1.2f(x)dx,ARo>2=∫∞ 2.0f(x)dx (8) Where ARo≤1.2, A1.2<Ro≤2, and ARo>2 are the integrated areas corresponding to each R o range fraction. The relative proportions of each carbonization category are then computed by normalizing these partial areas by the total area under the KDE curve (Fig. 3): FRo≤1.2=ARo≤1.2 Atotal ;F1.2<Ro≤2=A1.2<Ro≤2 Atotal ;FRo>2=ARo>2 Atotal (9) Where FRo≤1.2, F1.2<Ro≤2, and FRo>2 represent the fraction of each carbonization class. Normalization ensures the sum of these proportions sums to 1 (100 %), accounting for rounding errors and numerical approximations inherent in the discretization and integration process. This methodology provides a robust and reproducible framework for quantifying carbonization distributions in biochar. KDE smoothing and Simpson-rule integration enable accurate estimation of inertinite fraction within each sample. This computational workflow supports standardization and certification by enabling transparent reporting of biochar maceral composition in accordance with the IBR o 2 protocol. 4.4.3. Example of biochar’s R o profiling A case study presented here demonstrates the applicability of the KDE method through comparing two approaches for quantifying carbonization class fractions in a biochar sample with a bimodal R o distribution (Table 1). Two operators independently measured 500 R o points on the same polished specimen using incident light microscopy. Carbonization classes were defined based on established reflectance thresholds: Ro≤1.2 % (poorly carbonized), 1.2% <Ro≤2.0 % (semiinertinite), and Ro>2.0 % (inertinite). The purpose was to assess how consistent the results are when derived using two different data processing methods. The fractions were calculated using two approaches: (1) direct frequency counting, where the number of R o points falling within each defined threshold range was counted and then normalized to the total number of measurements (500 points), and (2) kernel density estimation (KDE), in which a continuous probability density function was generated from the same dataset, and the proportion of each class was calculated by integrating the area under the KDE curve within the corresponding thresholds (Table 1). The results show that carbonization fractions derived directly from the frequency distribution exhibit considerable variability between the two operators. This is particularly evident in the intermediate class (semi-inertinite), where relative differences between operators reached up to 30 % (Table 1). These discrepancies likely arise from the discrete nature of the dataset and the sensitivity of direct counting to small shifts in the distribution near classification boundaries. In contrast, the KDEderived values show much closer agreement between the two operators, with relative differences consistently within <10 % for all carbonization fractions and <1.5 % for the inertinite fraction (Table 1). This highlights the advantage of KDE in smoothing out minor local variations in the dataset, thereby reducing operator-dependent variability and enhancing reproducibility. Overall, the KDE approach provides a more reliable and consistent method for determining carbonization fractions in biochar using reflectance data. 4.5. Identification and significance of inertinite Under incident white light, carbonized macerals in biochar are easily recognized by their gray tones (varying with reflectivity), their topographic relief compared to surrounding resin or mineral matter, and a distinctive black rim surrounding the polished surface (Petersen et al., 2023, 2025; Sanei et al., 2024). This rim is the portion of the carbonized material that remains unpolished and lies beneath the exposed cross section, still embedded within the resin matrix. Inertinite represents a chemically inert maceral group composed of highly fused, polyaromatic carbon structures. It is characterized by high reflectivity, non-fluorescent behavior, and distinct morphological features such as vacuoles, which result from devolatilization during thermochemical conversions and from anatomical structures inherited from the original plant tissue (Diessel, 1983; International Committee for Coal and Organic Petrology (ICCP), 2001; Morga, 2011). The predominant inertinite maceral is fusinite, typically derived from lignocellulosic biomass, although other forms may originate from liptinitic matter (Hower et al., 2009). In the context of biochar, inertinite constitutes the most chemically stable fraction of the organic matter. Reflectance measurements (R o ) exceeding 2.0 % mark the beginning of the inertinite domain in the reflectance distribution histogram (Mastalerz et al., 2025; Sanei et al., 2024). These elevated R o values indicate extensive carbon condensation and are associated with high resistance to microbial degradation and abiotic oxidation (Petersen et al., 2025; Sanei et al., 2024). Quantifying the proportion of biochar fragments with R o >2.0 % is therefore critical for evaluating long-term carbon stability and permanence. 4.6. Inertinite R o threshold As discussed in the preceding sections, R o values greater than 2.0 % have been proposed as the inertinite benchmark for biochar (Sanei et al., 2024). A common misconception is to interpret this threshold as referring to the mean R o value. However, the R o >2.0 % benchmark denotes the lower boundary of the inertinite R o range, not its mean (see Fig. 3). In other words, the inertinite benchmark is defined by the onset of the Table 1 Comparison of carbonization fractions in a biochar sample, independently measured by two operators using two methods: (1) direct frequency counting based on fixed R o thresholds (Histo) and (2) kernel density estimation (KDE). Each operator analyzed 500 R o measurements on the same polished specimen. Carbonization classes were defined as R o ≤1.2 % (poorly carbonized), 1.2 % < R o ≤2.0 % (semi-inertinite), and R o >2.0 % (inertinite). Results from direct frequency counting show higher inter-operator variability, while KDE-derived values exhibit improved agreement, demonstrating the enhanced reproducibility of KDE. Fractions Method Operator 1 Operator 2 Relative Difference F Ro≤1.2 Histo: 5.6 % 6.6 % 17 % KDE: 5.7 % 5.6 % 1.8 % F 1.2<Ro≤2 Histo: 15 % 11 % 30 % KDE: 13 % 12 % 10 % F Ro>2 Histo: 79 % 82 % 3.8 % KDE: 81 % 82 % 1.5 % H. Sanei et al. International Journal of Coal Geology 310 (2025) 104886 7 reflectance distribution, not its central tendency. This means that a sample with a mean R o just above 2.0 % (e.g., 2.1 %) would still contain a substantial portion of reflectance values below the 2.0 % threshold. Assuming a perfectly Gaussian (normal) distribution of R o values, nearly half of the population could fall below the inertinite threshold. Consequently, the mean R o alone is insufficient for classifying a sample as inertinite-rich biochar. Achieving a reflectance distribution entirely above the 2.0 % threshold under normal distribution assumptions requires a mean R o greater than approximately 2.5 %. This distributional approach is consistent with the reflectance characteristics of inertinite-rich macerals observed in coal. Petersen et al. (2025) demonstrated that geological fusinite consistently exhibits mean R o values exceeding 2.5 %, with structural and chemical characteristics comparable to those of highly carbonized biochars. 4.7. Advantages of reflectance over bulk chemical proxies Unlike bulk chemical indices such as the molar H/C and O/C ratios, which can be influenced by sample heterogeneity and mineral contamination, the R o method specifically targets carbonized organic matter in a spatially resolved manner. Because R o values are measured directly on individual organic fragments, the results are not affected by the presence of ash or extraneous inorganic material, making the method particularly robust for high-ash biochars (Enders et al., 2012; Petersen and Sanei, 2025; Sanei et al., 2024). In addition, reflectance analysis captures the inherent variability in carbonization across the biochar sample by providing a distribution of results rather than a bulk average, providing a spatially resolved representation of the degree of carbonization within a sample. The resulting reflectance histogram reflects the full distribution of low-, intermediate-, and highly-carbonized components, enabling a more accurate and nuanced assessment of the sample’s carbonization. This is especially important for heterogeneous biochars, where carbonization may vary widely. 5. Uncertainty in inertinite quantification and its impact on CDR This section presents a Monte Carlo simulation model developed to quantify how sampling frequency and variance in measured CInert affect the uncertainty in estimating CDR from biochar. Based on the measured variance, the model subsequently derives a statistically justified safety margin with regard to CInert quantification to be applied to reported values, ensuring conservative and reliable crediting. In the IBR o 2 method, the amount of CDR that can be credited is directly determined by the measured content of CInert in the biochar. However, the reported CInert value can vary due to two main sources: measurement uncertainty (from analytical variability) and sampling variability (differences between individual samples). Sampling variability is particularly important, as biochar often shows considerable heterogeneity within a single production batch (Bucheli et al., 2014), such as the annual output of a continuously operated pyrolysis facility, or between batches produced under the same conditions in discontinuously operated systems. Preparing composite sampling, where multiple subsamples are combined into one mixed sample, can reduce the effect of small-scale variations and give a more representative sample. However, it also has a drawback in smoothing out or hiding important differences between production runs. These differences may reflect real changes in carbon stability and, if not properly accounted for, could lead to inaccurate CDR estimates. To quantify this uncertainty, it is crucial to establish the coefficient of variation (CV) from independently collected samples of the same production batch. It is recommended to collect at least three individual samples representing either an entire batch or separate runs of a discontinuous pyrolysis unit (when assessing variability within a production process). The mean and standard deviation ( σ ) of CInert (dry wt%) from these samples are then used to calculate the CV, a key input parameter for the Monte Carlo simulation. This allows statistical uncertainty to be embedded into credit issuance rules, consistent with quality criteria set forth by frameworks such as International Organization for Standardization (ISO), 2019, the EU CRCF regulation (European Union, 2024), UNFCCC guidance (United Nations Framework Convention on Climate Change (UNFCCC), 2024), or the ICVCM assessment framework (Integrity Council for the Voluntary Carbon Market (ICVCM), 2024). This framework makes a clear distinction between two sources of variability: (1) natural variability within a single production system or batch, including analytical and sampling uncertainty; and (2) deliberate changes made by the producer to optimize or modify the production process. The primary purpose is to assess and certify the level of variability within a defined production batch. If variability increases significantly, the framework helps determine whether this is due to random fluctuations or a modification of the production parameters. In the latter case, producers are expected to register a new batch and submit new baseline samples for analysis. Certification bodies can use the submitted data to monitor such changes. The framework allows producers to manage the homogeneity of their production and to clearly define production batches. If variability rises unexpectedly, a new batch should be started, or the system will automatically increase the safety margin applied to CDR credits. This safety margin can only be adjusted downward once additional data confirm a more consistent output. It is also important to clarify what this method does and does not provide. The presented statistical approach does not inherently ensure the representativeness of individual samples across an extended production period or large production volumes. The model is valid when supported by adherence to representative sampling at the production premises (see Schmidt et al. (2024)). Sampling protocols must be specifically designed for the type of production, and implemented consistently using a validated procedure, whether samples are collected cross flow or from storage piles. It offers, however, a transparent, reproducible, and operator-independent means to quantify sampling-induced uncertainty at any given point in time or production interval. Producers and certification bodies should therefore interpret the resulting safety margins and CDR risk estimates as tools guiding operational decision-making and biochar carbon certification. 5.1. Variance estimation of inertinite carbon in biochar samples Let ψ iwith i=1,2,3... denote the reported CInert (dry wt%) for isubmitted biochar samples from a production site, then the CV can be calculated using arithmetic mean ( ψ )and standard deviation ( σ )from the sample: CV = σ ψ =  (1 N−1∑N i=1( ψ i− ψ )2) √ √ √ √1 N∑N i=1 ψ i (10) The CV provides a normalized measure of dispersion that is independent of the magnitude of the mean and is used as a key input for risk modeling. For a biochar production site, the total amount of inertinite carbon produced annually is estimated by multiplying the annual production (tonnes per year), by the mean CInert (dry wt%). The mass of inertinite carbon is then converted to its equivalent mass of removed CO 2 using the molecular mass ratio of carbon dioxide (44.01 g/mol) to elemental carbon (12.01 g/mol). This yields the baseline estimate of CDR in tonnes per year: CDRi=Pi×(Cmean Inert,i 100 )×44.01 12.01 (11) H. Sanei et al. International Journal of Coal Geology 310 (2025) 104886 8 Where P i is the production volume (tonnes per year; or tonnes for a biochar batch) and Cmean Inert,i is the measured CInert (dry wt%). However, in practical applications, the reported mean CInert is subject to analytical and spatial variability, particularly when the sample frequency is low. To model this uncertainty, the simulation assumes that the measured CInert follows a normal distribution with a standard deviation ( σ ) derived from the coefficient of variation: σ i=(CVi 100)×Cmean inert,i(12) 5.2. Predicting CDR variability through Monte Carlo simulation A Monte Carlo simulation is then used to generate synthetic sample means for different sampling frequencies. For each sampling frequency n (e.g., weekly, monthly), the model generates 1000 random sample sets. Each set consists of n individual samples drawn from the specified normal distribution. The mean of each sample set is then converted into an annual CDR estimate, yielding a distribution of 1000 possible outcomes for that sampling frequency: CDR(n) j=Xj n×44.01 12.01 ×Pi 100 (13) From this distribution, the model calculates the expected (mean) annual CDR and the 5th percentiles of the simulated values. 5.2.1. Quantifying CDR risk and estimating safety margins The lower 5th percentile is interpreted as the conservative bound on creditable CDR, representing the minimum quantity of CO 2 removal that can be claimed with 95 % confidence given the sampling variability. The difference between the expected CDR and this lower bound defines the CDR risk: Risk(n) CDR = μ (n) CDR −CI(n) 5% (14) Where μ (n) CDR is the mean CDR value based on n independent samples, and CI(n) 5% is the 5th percentile of the corresponding confidence interval for the same sample size. The relative risk (%), which is also regarded as the percentage safety margin, is computed as: Relative Risk(n) CDR =(Risk(n) CDR μ (n) CDR )×100% (15) Where Relative Risk(n) CDR expresses the magnitude of the risk as a percentage of the expected CDR value. This formulation allows for quantifying the safety margin by capturing the extent to which low frequency sampling increases uncertainty in reported CDR and determining the percentage that must be deducted from the reported mean to ensure conservative and statistically robust crediting. The matrix chart shown in Fig. 4 is derived from the Monte Carlo simulation model and displays expected safety margin percentages across a range of CV values and sampling frequencies. It enables producers and certifiers to determine deduction levels appropriate to the degree of variability in measured CInert within a given production system. To apply the matrix, the CV must first be calculated from CInert values measured in at least three independent samples submitted for certification (see Sections 5 and 5.1). Once the CV is established, the matrix provides the corresponding safety margin for the intended sampling frequency (e.g., weekly, monthly, or annually). Alternatively, when regulatory requirements do not demand exact decimal precision but only a reliable quick estimate, the reporting safety margin can be obtained from a simple closed-form formula rather than running full Monte Carlo simulations. In this formulation, the safety margin is closely reproduced by a linear relationship between CV and the inverse square root of the sampling frequency: Saftey margin (%) ≈ kCV(%) n √(16) where n is the number of samples per year and k is a scaling constant. Across a wide range of simulations, the best fit was consistently obtained for k≈1.65, which corresponds to the one-sided 95 % confidence limit. Fig. 5a-c compare Monte Carlo simulations with the above closed-form equation using k=1.65, shown as dashed linear fits across three representative mean CInert values of 10, 40, and 80 wt%. With this choice of k, Fig. 4. Matrix illustrating the reporting safety margin (%) required to account for sampling uncertainty in measured inertinite carbon content (C Iner ) across a range of coefficient of variation (CV) values and sampling frequencies. Sampling frequencies are expressed as number of samples per year: 96 (twice per week), 48 (weekly), 24 (biweekly), 12 (monthly), 4 (quarterly), 2 (semi-annually), and 1 (annually). Higher CV values lead to increased levels of safety margin (%) deductions from the carbon dioxide removal (CDR) credit issuance due to greater sampling variability. Conversely, increasing the number of samples per year reduces uncertainty, thereby lowering the required safety margin and improving the accuracy, credibility, and marketability of issued CDR credits for biochar facilities. H. Sanei et al. International Journal of Coal Geology 310 (2025) 104886 9