Global soil organic carbon changes and economic revenues with biochar application

Biochar has been proposed as a promising negative CO2 emission technology to mitigate future climate change with the additional benefit of increasing agricultural production. However, the spatial responses of soil organic carbon (SOC) to biochar addition in cropland are still uncertain, and the economic feasibility of large‐scale biochar implementation remains unclear. Here, we analyzed the response of SOC to biochar addition using 389 paired field measurements. The results show that biochar addition significantly increased SOC by 45.8% on average with large regional variations. Using a random forest model trained with soil, climate, biotic, biochar, and management factors, we found that the response of SOC to biochar addition was mainly dependent on biochar application rates, initial SOC, edaphic (e.g., pH), and climatic (e.g., mean annual precipitation) variables. Combined with the predicted SOC changes to biochar addition on the global cropland, we assessed the revenue of the biochar system based on the current and potential pyrolysis plants in the world using the life‐cycle analysis. Net revenue of the currently existing 144 pyrolysis plants increases with larger plant capacity and higher carbon price. Potential revenue of building new plants is high in regions like America and Europe but low in regions with infertile soil, low crop residues availability, and inconvenient transportation. The global CO2 removal of biochar application is 6.6 Tg CO2e (CO2 equivalent) year−1 with a net revenue of $ 177 million dollars at a carbon price of $ 50 t−1 CO2 for current pyrolysis plants with a biomass‐processing capacity of 20,000 t year−1. Our study provides a full economic assessment of idealized biochar addition scenarios and identifies the locations with maximal potential revenues with new pyrolysis plants.


| INTRODUCTION
To meet the Paris Agreement's goal of limiting global warming to below 2°C, negative emission technologies will be required (Minx et al., 2018;Smith, 2016;UNFCCC, 2015), such as increasing soil carbon storage with the help of biochar. The production of biochar starts from carbon removed from the atmosphere and incorporated in plant material and transforms it via pyrolysis onto a very longlived product that remains in the soil, being recalcitrant to microbial decomposition (Smith, 2016;Woolf et al., 2010). The CO 2 sequestration potential of biochar, when applied to soil, is estimated to be 0.5-2 GtCO 2e year −1 (Minx et al., 2018;Wu et al., 2019). In addition, biochar has potential co-benefits, such as increasing crop yield, reducing soil nutrient loss, and improving soil water-holding capacity (Jeffery et al., 2011;Jones et al., 2012;Liu et al., 2016).
Biochar amendments to soil affect the stability and turnover of SOC by changing pH, bulk density, and moisture conditions, and these effects vary with soil properties, biochar application rate, biochar type, and pyrolysis temperature (Liu et al., 2016). For example, biochar application may increase the priming of SOC by increasing soil microbial activity, accompanied with a decrease of soil bulk density and increased soil pH, moisture, and porosity (Luo, Wang, et al., 2017;Luo, Zang, et al., 2017). Biochar can also enhance the formation of soil aggregates through the associations of soil minerals and biochar by blocking contacts between microorganisms and SOC particles (Zheng et al., 2018). A greater microbial carbon use efficiency and carbon sequestration potential were found with biochar amendment compared to crop straw . Liu et al. (2016) analyzed the responses of soil CO 2 flux, SOC, and soil microbial biomass carbon (MBC) to biochar additions using 148 paired field SOC data from 50 publications and found that SOC was increased significantly by 40% (95% confidence interval [CI] = 32-51) within 4 years compared to the non-biochar addition treatment. Bai et al. (2019) compared the effects of biochar on SOC sequestration with other climate-smart agriculture management practices such as conservation tillage and cover crop and found that SOC content was increased by 28% (95% CI = 23-32) with field biochar application (96.7% sites <5 years). These meta-analyses summarized the direction and magnitude of SOC change with biochar addition in different regions, soil types (e.g., pH, texture), applied biochar types (e.g., feedstock, pyrolysis temperature), and management practices (e.g., fertilization, residue management). However, the spatial pattern of SOC change due to large-scale biochar implementation at the global scale remains unclear, which is valuable to identify the optimal regions for biochar application to gain the largest SOC accumulation and assess the carbon removal potential of biochar application as a natural climate solution.
In addition to the ecological benefits of SOC sequestration by biochar, the feasibility of biochar implementation depends on economic costs and revenues of the biochar supply chain. Cost-revenue analysis of the whole biochar system is thus required. The widely used method is life cycle analysis (LCA) that assesses the economic, energetic, and greenhouse gases (GHGs) emissions from feedstock collection to biochar application in biochar system. Previous studies assessed the revenue or climate change mitigation potential in the full life cycle of biochar implementation (Azzi et al., 2019;Hammond et al., 2011;Roberts et al., 2010) and showed that the economic feasibility of biochar-pyrolysis system largely depends on the cost of biochar feedstock production, pyrolysis system, and carbon offset value (Homagain et al., 2016;Roberts et al., 2010). However, these studies mainly focused on a small scale, such as a certain region or a single pyrolysis plant with specific conditions like energy structure, feedstock availability, and pyrolysis facility processing system on a small scale (Homagain et al., 2016;Roberts et al., 2010). In addition, biochar can not only increase SOC through the stable fraction of itself, but also change the original (native) SOC through changes, e.g., in soil properties and microbial activities. These processes of the native SOC changes with biochar addition were largely ignored in previous LCA analysis (Hammond et al., 2011). Therefore, it is important to analyze the native SOC changes with biochar addition at the global scale using the existing field observations, in order to fully assess the potential and feasibility of biochar application as a climate mitigation option. biochar application is 6.6 Tg CO 2 e (CO 2 equivalent) year −1 with a net revenue of $ 177 million dollars at a carbon price of $ 50 t −1 CO 2 for current pyrolysis plants with a biomass-processing capacity of 20,000 t year −1 . Our study provides a full economic assessment of idealized biochar addition scenarios and identifies the locations with maximal potential revenues with new pyrolysis plants.

K E Y W O R D S
biochar, climate change, economic revenue, life cycle analysis, random forest, SOC In this study, we first analyzed the responses of SOC and other soil variables to biochar application based on a meta-analysis of the 389 paired experiment data with and without biochar addition. We then trained a random forest (RF) model using the observation data of SOC responses to biochar addition and the corresponding edaphic, climatic, biotic, and management variables. The trained RF model was further used to predict SOC change due to biochar addition on the global cropland at 1 km × 1 km. Finally, we adopted the LCA method (Figures S1 and S2) to estimate the net revenue of the biochar system, in order to analyze the regional differences of cost-revenue contributions of current pyrolysis plants and to identify the potential locations with the maximum net revenue for building new plants.

| Data collection
We conducted a literature search in the Web of Science (http://apps.webof knowl edge.com/) and China Knowledge Resource Integrated (CNKI) databases (http:// www. cnki.net/) with the keywords of "soil organic carbon" or "SOC" and "biochar." Publications were selected using the following criteria: (1) Biochar addition experiment must be conducted in the field cropland with explicit locations; (2) at least one pair of SOC data (control and biochar-amended treatment) is available. A total of 389 paired field SOC data with control and biochar-amended treatments from 70 peer-reviewed publications were extracted, and their spatial locations are shown in Figure  S3. The experiment duration for most data is short-term (83.8% sites within 3 years).
Six attributes were taken from the selected publications of biochar application experiments when available (Table S1): (1) SOC concentration (SOC with control [SOC_ini] and biochar addition) and soil properties, including clay (CLAY), sand (SAND), silt (SILT), bulk density (BD), base saturation (BS), pH, cation exchange capacity (CEC), average annual soil moisture (SM), total N (TN), total P (TP), total K (TK), and C/N ratio ( Figure  S4); (2) climate, including the mean annual temperature (MAT), mean annual precipitation (MAP), mean of the minimum temperature of each month (T min ) and mean of the maximum temperature of each month (T max ); (3) biological variable (net primary productivity [NPP] of the present field crop); (4) biochar characteristics, including biochar feedstock types (i.e., biomass source type for biochar production, FS), pyrolysis temperature (Temp_BC), biochar application rate (Rate_BC), and aging time of biochar (Age_BC, i.e., interval between biochar application and soil sampling); (5) management practices, including irrigation and fertilizer rate (Fertilizer) ( Figure S5); (6) other auxiliary variables, including coordinates (latitude and longitude), experimental duration, crop types, rotation (whether or not) and the month of crop harvested. Mean and standard deviation (or error) of SOC concentrations and soil properties variables in each treatment were also extracted from the contexts or Supplementary Materials for the meta-analysis (Section 2.2). Base fertilizer was widely applied to improve crop yield in the field experiments, and thus some studies explored the effects of the combination of biochar and fertilizer on SOC. Such data were documented as fertilizer treatment and "fertilizer + biochar". For sites with more than one experiment available, each experiment was considered as independent paired data.

| Meta-analysis
Following the meta-analysis method described by Hedges et al. (1999), the natural log-transformed response ratio (RR) was used to calculate the effect size: where X t and X c are the means in the treatment and control groups, respectively.
The variance (v) of each individual RR is estimated as: where n t and n c refer to the number of replicates of the variable in the treatment and control groups; S t and S c are the standard deviations in the treatment and control groups. The variance (v) of each RR was calculated based on the study of Hedges et al. (1999). The effect size of each observation is weighted by the inverse of v, and then it is used to calculate the 95% CI. The mean effect size was quantified by the weighted RR (RR ++ ), RR ++ , and 95% CI were generated from the random-effects model using the metacont package in R (Adams et al., 1997). If the 95% CI do not overlap with zero, biochar addition is considered to pose a significant effect on SOC change. The effect size was transformed into a percent relative change (RC SOC ): A significance level of p < 0.05 was used in all the statistical analyses.

| Model training
We adopted RF to capture the complex nonlinear relationship between SOC changes and the biochar addition as well as the interactions with other explanatory variables. RF is a machine learning method with robust performance and low bias for prediction by building numerous trees based on the original observation dataset (Breiman, 2001) and has been widely applied for predicting spatial-temporal variation of soil carbon stock , soil respiration (Huang et al., 2020), GHGs emissions , and crop yields (Everingham et al., 2016;. The observed paired data were randomly split into two parts in this study: 291 sites (75% of the total) as the training subset and 98 sites (25%) as the test subset. We used the scikitlearn module in Python (Pedregosa et al., 2011) to perform the RF regression with bootstrapped 1000 trees and a maximum tree depth of 15.
In RF training, the natural logarithm-transformed RR (Equation 1) is used as the response variable, and the explanatory variables include edaphic variables (CLAY, SAND, SILT, BD, BS, pH, CEC, TN, TP and TK, SM, SOC_ini), climatic variables (MAT, MAP, T min , T max ), biological variable (NPP), biochar-related variables (Rate_ BC, FS, Temp_BC, Age_BC), and management variable (Fertilizer), of which all variables are continuous variables except the biochar feedstock type. The detailed information of all variables is listed in Table S1. NPP, as the input of SOC, is the main biological factor that controls SOC, and other biological-related factors (e.g., microbial activity) were not included due to a lack of relevant data. Management practice-related variables only include the fertilizer rate because very limited information on irrigation was available ( Figure S5).
The coefficient of determination (R 2 ) between predictions and observations in the training data and test data were derived to evaluate the RF model performance. The corresponding R 2 is 0.95 and 0.81 for the training and test data, respectively ( Figure S6). As a comparison, we also tried linear regression models, including multiple linear regression (MLR) and the least absolute shrinkage and selection operator (LASSO) regression model, and the R 2 (0.53 and 0.46, respectively, Figures S7b and S8b) is lower than the R 2 from the RF model.
The importance of variables can be directly obtained from the RF model training based on their contributions to the tree splits using the sklearn package in Python (Louppe et al., 2013). The partial dependence between the explanatory variables and a set of dependent features were conducted to analyze the response of SOC change with biochar addition to single important variables.

| Model prediction
The trained RF model was then used to project SOC change with idealized biochar addition on the global cropland at 1 km × 1 km resolution (same biochar application rate on each cropland grid cell). The RC SOC with biochar addition is assumed consistent in topsoil (0-30 cm). The base SOC stock (soil depth = 28.9 cm) derived from the gridded global soil dataset for the use in earth system models (GSDE) (Shangguan et al., 2014) with biochar application for the LCA. Cropland distribution at 1 km × 1 km was derived from global land cover map database (GLCM) (Ramankutty et al., 2008). The global explanatory variables include 12 edaphic variables, 4 climate variables, 1 biotic variable, 1 management variable, and 4 user-specified biochar-related variables, and the detailed information (e.g., source and spatial resolution) of these variables are listed in Table S2. For the soil moisture and nitrogen fertilizer with coarser resolutions, nearest-neighbor resampling was used to obtain the data at 1 km × 1 km resolution. We made a total of 21 predictions (Table S3) by combining three biochar feedstock types (straw, wood, and manure) and seven biochar application rates (5, 10, 20, 40, 60, 80, and 100 t ha −1 ). The seven application rates are within the range of application rates in the training dataset (1.5-120 t ha −1 , Figure S9a). We also made four sensitivity tests using four pyrolysis temperatures with straw as the feedstock type and a biochar application rate of 20 t ha −1 (Table S3).

| Life cycle analysis
Life cycle analysis is an assessment method based on processes from "cradle" to "grave," which has been widely used for the assessment of environmental consequences in a system in the field of energetic, economic, and carbon abasement (Guenet et al., 2021;Sykes et al., 2020). In this study, LCA was modified from Roberts et al. (2010). The LCA ( Figure S2) includes the processes of: (a) feedstock collection, (b) feedstock transportation from cropland to biochar pyrolysis plant, (c) biomass pyrolysis, (d) the return transportation from plant to cropland, (e) biochar applied, and its derivational influence; (f) energy generation from (c); (g) changes in SOC, N 2 O emissions and NPK fertilizer use efficiency from (e). Cost and revenue as well as carbon emission and reduced emission in each process were calculated for the economic and carbon balance analysis. The SOC changes due to biochar addition (20 t ha −1 ) were taken from the RF projection outputs (i.e., AC SOC ). Three carbon prices ($ 20, $ 50, $ 100 t −1 CO 2 ) were used in the LCA for further analysis.
The transport distance between cropland and pyrolysis plant is related to the plant capacity, the amount of available agricultural residue around the plant, and the transportation accessibility (Han et al., 2013). In this study, the plant capacity and the available agricultural residues were used to calculate the maximum base service distance (D max ) that satisfies the plant capacity based on the road network (https://www.natur alear thdata. com/) with ArcGIS 12.0 ( Figure S10). The distances thus varied with available agricultural feedstock and the plant capacity. Note different from the resolution for the RF (1 km × 1 km), the LCA analysis was conducted at a resolution of 10 km × 10 km due to the complexity in the road network and the corresponding computing cost of distance at a high resolution. The available residue on cropland within D max was first transported to the pyrolysis plant for biochar production, and then the produced biochar was sent back to the corresponding cropland grid cell for application. Because the biochar application rate is 20 t ha −1 , much higher than the biomass density in the cropland, it means only a part of cropland in each 10 km × 10 km grid cell is added with biochar. The net revenue of each biochar pyrolysis plant was summed up within D max (see Figure S10 for details). It should be noted that we did not explicitly consider the biomass competition between nearby plants.
We surveyed data from 144 companies in the world for biochar production in the State of the Biochar Industry 2015 from International Biochar Initiative (IBI, https:// bioch ar-inter natio nal.org/), and the coordinates of these plants were acquired using Google earth by relative information ( Figure S11). Due to the lack of information on the plant capacity, we tested three assumed biomass processing capacities: 2000, 20,000, and 100,000 t dry biomass (DM) year −1 , representing small-, medium-, and largesized plants, respectively (Hammond et al., 2011). The agricultural residues were estimated from satellite-based NPP after removing the fraction of root, grain, and the other usage for poultry feed (Yang et al., 2010) (24%) and industrial materials (Yang et al., 2010) (3%) using the root fraction (Bolinder et al., 2007) (20%) and the ratio of residue to grain (NDRC, 2015) (0.89). In addition, collection efficiency (Yang et al., 2010) (0.83) was used to account for the loss of uncollected parts of crop residues (default).
Finally, in order to identify possible locations for new pyrolysis plants, we estimated the revenues of plants with a biomass processing capacity of 2000, 20,000, or 100,000 t year −1 in each cropland grid cell at a carbon price of $ 20, $ 50, or $ 100 t −1 CO 2 . We thus provided global revenue maps to show the revenue if a new plant will be built in each grid cell.

F I G U R E 1
The effect of biochar addition on soil organic carbon (SOC) and other soil properties in croplands. The circle and error bar indicate the mean value and 95% confidence interval (CI). The acronyms SOC, microbial biomass carbon (MBC), dissolved organic carbon (DOC), bulk density (BD), total N (TN), C/N (ratio of soil carbon to nitrogen). If the error bar does not overlap with zero, the response is considered significant ('*'). The numbers of observations are also shown

| Importance of explanatory variables
Among all the explanatory variables, Rate_BC is the most important variable that contributes more than half (51.3%) of the tree splits in the RF model (Figure 2a). SOC response to Rate_BC is positive as indicated by the partial dependence of the SOC RR on Rate_BC (Figure 2b). This is expected since most applied biochar stays stable in soils and contributes to the total soil carbon content. The other biochar-related variables (Age_BC, Temp_BC, and FS) make a contribution of 9.3% in total (Figure 2a). Especially for Age_BC (i.e., the interval between biochar application and soil sampling), SOC RR with Age_BC ≥5 year is lower than Age_BC <5 year ( Figure S13). The acronym Temp_ BC (biochar pyrolysis temperature) is the third important variable and has a negative partial relation with SOC RR (Figure 2b). The total importance of biochar-related variables reaches 60.6%, which is the largest one among the five categories ( Figure 2a, the inset pie plot).
The second most important variable is SOC_ini (initial SOC content) with a contribution of 18.7% (Figure 2a), and the SOC RR responds negatively to SOC_ini (Figure 2b). When SOC_ini <10 g kg −1 , the decreasing trend of SOC RR is very sharp (Figure 2b). SOC RR is positively related to the other edaphic variable, pH (Figure 2b), with a contribution of 2.4% to the RF model (Figure 2a). The total contribution of edaphic variables (SOC_ini, CLAY, SAND, SILT, BD, BS, pH, CEC, SM, TN, TP, TK) is 31.4%, which is the second most important category next to the biochar-related variables (Figure 2a, the inset pie plot). The fourth important variable is MAP with a contribution of 2.8% (Figure 2a). However, it should be noted that the correlation between MAP and SM is high (R 2 =0.97; Figure S14), and thus the contributions between these two variables may not be well separated. We found that MAP has a negative impact on the SOC RR, indicated as a lower MAP associated with a higher SOC RR (Figure 2b). The total climate-related variables (MAP, MAT, T min , and T max ) contribute 5.0% of all explanatory variables (Figure 2a). Management (only fertilization included) has a minor contribution of 1.8% (Figure 2a), probably due to the lack of other management information (e.g., irrigation, rotation) and difficulties in quantifying management practices. The biological variable (NPP of the present live crop on the experimental field, not NPP for the crop whose residue is used for biochar production) has minimum importance F I G U R E 2 Variables with an importance ≥1% in the trained random forest (RF) model (a) and partial dependence of soil organic carbon (SOC) changes with biochar addition (i.e., response ratio [RR] in Equation 1) to the six most important variables (b). Variables in (a) are further grouped into five broad categories in the inset pie plot (Biochar-, Edaphic-, Climate-, Management-, and Biological-related variable). Variables with an importance <1% in the trained RF model are shown in Figure S15 of 1.1% of selected explanatory variables (Figure 2a). The response of RR to Age_BC is rather uncertain, due to the lack of long-term experiment data (Section 4.2).

| Predicted global SOC change due to biochar addition
Using the RF model trained by field observations, we estimated SOC RR to biochar addition on global croplands at 1 km × 1 km and converted to the relative changes of SOC (RC) (Figure S16) using Equation (3). We further derived the absolute changes of SOC (AC SOC ) based on the present cropland SOC stock ( Figure S17b). The spatial AC SOC patterns on global croplands with three biochar application rates (20, 60, and 100 t ha −1 ), Temp_BC = 450°C (median value of all field data) and Age_BC = 1 year (i.e., biochar addition to soil for 1 year, the median value of all field data = 1.2 year, Table S3) are shown in Figure 3, and maps for other biochar addition rates and other biochar are shown in Figures S18-S20.
Globally, AC SOC with biochar addition ranges from 12.6 to 63.8 t ha −1 , and the response is generally linear to the biochar addition rate (Figure 3g). Larger AC SOC is predicted in the temperate regions, such as North America, Europe, and Northern China with larger F I G U R E 3 Predicted absolute soil organic carbon (SOC) changes (AC SOC ) and absolute native SOC change (AC SOC-native ) with different biochar application rates on the global croplands. Predicted absolute SOC changes (AC SOC ) on the global croplands with (a) 20 t ha −1 , (c) 60 t ha −1 , (e) 100 t ha −1 straw biochar addition; and absolute native SOC change (AC SOC-native ) with (b) 20 t ha −1 , (d) 60 t ha −1 , (f) 100 t ha −1 straw biochar addition. (g) Predicted AC SOC (left y-axis) and AC SOC-native (right y-axis) across all cropland grid cells with seven different biochar application rates from different biochar types (straw-, wood-, and manure-derived; only straw biochar is shown in purple line for AC SOC-native ). Error bars represent 95% confidence intervals initial SOC content (Figure 3a,c,e; Figure S17b). Smaller AC SOC response occurs in the tropical and subtropical regions (e.g., Latin America) with warm and moist climate, smaller initial SOC content, and lower pH ( Figure  3a,c,e; Figure S17b,e,g,h). The AC SOC with straw-, wood-, and manure-derived biochar addition have similar spatial patterns ( Figures S18-S20), but the AC SOC of manure biochar is greater than the other two types (Figure 3g).
The SOC changes due to biochar addition comprise two parts. The first part is the stable fraction of biochar with a long turnover time, and thus it increases SOC. The second part is changes in the native SOC due to the impacts of biochar addition, including direct impacts on soil properties and microbial activities and indirect impacts on plant growth which further influence the return of plant residues into SOC. Therefore, we also represent the changes in the native SOC (AC SOC-native ) after 1 year since the biochar addition by excluding the stable biochar (assuming a fixed fraction (91%) (Wang et al., 2016) of the application rate with a biochar carbon content of 0.68 (Roberts et al., 2010) from the total SOC changes (Figure 3b,d,f). The global median AC SOC-native varies with the biochar application rates (an increase of 6.0 t ha −1 for 20 t ha −1 biochar application rate, and a reduction of 6.1 and 15.5 t ha −1 for 60 and 100 t ha −1 biochar application rates, respectively, Figure 3g). Spatially, biochar addition generally increases the native SOC in most regions at low addition levels but decreases native SOC at higher addition levels in some areas with lower initial SOC content, such as South Asia, southern Africa, and Australia ( Figure  3b,d,f; Figure S17b).

| Revenue of current pyrolysis plants
Globally, the net revenue of a pyrolysis plant increases with the plant capacity and the carbon price (Figure 4a). After grouping current pyrolysis plants into different regions ( Figure S21), the regional average net revenue patterns ( Figure S22) are similar to the global one ( Figure  4a). We further calculated the threshold of carbon price at the global scale and in each region, above which the average net revenue changes from negative to positive. The thresholds decrease with the increase in plant capacity (Figure 4b). High carbon price thresholds were found in Africa, East Asia, and Oceania, which are above the global level (Figure 4b). By contrast, Europe and Latin America show low carbon price thresholds (Figure 4b).
Take three-carbon prices ($ 20, $ 50, or $ 100 t −1 CO 2 ) and three biomass-processing capacities of pyrolysis plants (2000, 20,000, or 100,000 t year −1 ) as examples. Europe and Southeast Asia are profitable at a low carbon price of $ 20 t −1 CO 2 in the case of medium and large biomassprocessing capacity ( Figure S23a,d,g; red point, i.e., revenue minus cost). At the carbon price of $ 50 t −1 CO 2 , most regions are profitable except East Asia, Oceania, and F I G U R E 4 Net revenues averaged over current pyrolysis plants with a range of plant sizes and carbon prices (a) and threshold of carbon prices for each region and at the global scale, beyond which the net revenue turns from negative to positive (b). Note that the contour line with zero net revenue in (a) is the black line in (b) for the global carbon price threshold. The calculation of net revenue is shown in Figure S2 Africa in the case of 2000 t year −1 plant capacity ( Figure  S23b,e,h). At a carbon price of $ 100 t −1 CO 2 , all regions are profitable in three cases of biomass-processing capacity ( Figure S23c,f,i).
The regions with the largest average net revenue are Southeast Asia ( Figure S23b) and Europe ( Figure S23c-i) in nine cases. The revenue mainly comes from the carbon value (offsetting GHGs emissions by sequestrating carbon with biochar and renewable energy generated to replace fossil fuels, Figure S2) and energy generation ( Figure S23). Africa, in most cases, has the lowest net revenue ( Figure  S23b-g), and Latin America ( Figure S23a) and East Asia ( Figure S23h,i) also show a lower revenue in some cases, as a result of high cost in operation, capital, and feedstock collection. The transport cost is relatively low (<16.1%, the ratio of transport cost to the sum of all four cost components). The average D max is 33.9, 45.4, and 66.6 km for small-, medium-, and large-size plants, respectively ( Figure  S24), and it generally increases with the plant biomassprocessing capacities due to more feedstock fetched.
We also calculated the total net revenue of all pyrolysis plants in each region ( Figure S25), and the pattern is similar to the average net average shown above. North America and Europe have the highest total net revenues in most cases, not only due to the high carbon value and the efficient energy generation ( Figure S26), but also due to the great (22.2% and 45.8% of global total plants) number of pyrolysis plants in these two regions ( Figure S21). In contrast, Latin America, South Asia, and Africa with few plants (<2.8%) have relatively low total net revenues ( Figures S26 and S27).
The carbon-dioxide removal (CDR) of current pyrolysis plants is greater than the CO 2 emissions in the life cycle of the biochar system ( Figures S28 and S29), indicating that the great CDR potential of biochar is used as a negative carbon emission technology. The CDR is mainly contributed by the sequestrated stable carbon in the biochar and the byproducts of biomass pyrolysis (e.g., bio-oil) as substitutes of fossil fuels for combustion ( Figure S28). The net CDR averaged over all plants in each region is high in Latin America, Europe, and Southeast Asia ( Figure S28). Summing up the CDR of all plants in each region, however, the total CDR is high in North America, Europe, and Oceania ( Figure S29a-c). As expected, the global total net CDR increases with plant biomass-processing capacities ( Figure S29d). In a medium case with a plant capacity of 20,000 t year −1 and a carbon price of $ 50 t −1 CO 2 , the global net CDR potential is 6.6 Tg CO 2 e year −1 with a total net revenue of $ 177 million dollars (Figures S26j and  S29d). The maximum global net CDR and revenue potential reach 25.5 Tg CO 2 e year −1 and 2.1 × 10 3 million dollars with a high carbon price of $ 100 t −1 CO 2 and a large-size plant capacity ( Figure S29d).

| Potential revenue of building new plants
We further calculated the net revenue of a hypothetical pyrolysis plant in each cropland grid cell (10 km × 10 km) to identify the possible locations for new plants with maximum profits. The global average D max is increased with the increased plant capacities (20.0, 47.8, and 83.2 km for 2000, 20,000, and 100,000 t year −1 plant capacities, respectively) ( Figure S30). The nine potential revenue maps (3 biomass-processing capacities × 3 carbon prices) are shown in Figure S31. Obviously, the net revenue increases with carbon price since most revenues are from the carbon value (i.e., net mitigated CO 2 e amount × carbon price; Figure S2). We choose a median scenario (i.e., carbon price = $ 50 t −1 CO 2 and biomass-processing capacity = 20,000 t year −1 ) as an example for further analysis ( Figure 5, see all cases in Figure S31). Most land cells with larger profitable new plants are located in North America, Europe, Northeast China, and Southeast Asia, where a large area of cropland with the significant increase in SOC in response to biochar addition exists (Figures 3a and 5; Figure S17a). On the contrary, India has a large cropland area ( Figure S17a), but the potential revenue of new pyrolysis plants is low (Figure 5), because of the relatively small increase in AC SOC in response to biochar addition (Figure 3a). In most regions in Africa, South America, and Australia, the potential revenue is low due to the low availability of feedstock, a small increase in AC SOC , and the high cost of transportation (Figures 3a and 5; Figure  S32).

| Effects of biochar on SOC
Biochar can sequestrate carbon due to its recalcitrant carbon content, and it can also alter the SOC fate by influencing native SOC formation and mineralization when applied to the soil. The effects of biochar on SOC were widely studied under different conditions in the field experiments (Blanco-Canqui et al., 2020;Liu, Zhu, et al., 2019;Qi et al., 2020), which provided valuable field data for our meta-analysis here. 11.2% of our collected field data were also used in a previous meta-analysis by Liu et al. (2016), the RC SOC (45.8%) with biochar addition in this study is slightly larger than the RC SOC (40%, 95% CI = 32%-51%) from Liu et al. (2016). Our derived RC SOC , however, is larger than the RC SOC (28%, 95% CI = 23-32) from another meta-analysis by Bai et al. (2019). The response of MBC to biochar application (an increase of 17.8%, Figure 1) is consistent with an 18% (95% CI = 12-23) increase in MBC in Liu et al. (2016), but lower than an average increase of 29% (95% CI = 21%-37) in the field experiments (Zhou et al., 2017).
The impact of biochar addition on SOC is affected by multiple factors (e.g., soil properties, climate conditions, and biochar characteristics) through complex mechanisms. Stable carbon contained in biochar was a large contributor to the short-term SOC increase. In addition, biochar addition can trigger positive or negative priming effects on the native SOC mineralization through its impacts on soil properties (e.g., pH, Clay) indirectly (Sheng et al., 2016;Wang et al., 2016). Biochar can promote larger CO 2 emission because of greater biochar degradation, through increasing bioavailability of SOC and copiotrophic bacteria abundance in acid soil, which is contrary to that in alkaline soil with negative priming effects (Sheng & Zhu, 2018). Native SOC mineralization may be also inhibited by biochar due to its porous structure (Purakayastha et al., 2019) and enhancing microaggregate formation (Purakayastha et al., 2019;Zheng et al., 2018), which increases the physical protection of SOC and reduces the contact between microorganisms and SOC. These mechanisms may partly explain the negative AC SOC-native in the tropical regions like Africa, South Asia, and the positive AC SOC-native in North America, Europe (Figure 3b,d,f).

| Sensitivity tests and uncertainties in the RF model
The distribution of global cropland input data is within the range of raining data of our RF model building, indicating that our extrapolation to global scale is robust ( Figure S33). We used the absolute SOC change (experimental -control) and the relative SOC change ([experimental -control)/control] as the independent variable directly to train the RF model for sensitivity test, and the validation R 2 are 0.85 and 0.74, respectively ( Figure  S34a,b). Since our datasets used for RF model training is experimental dependence (e.g., one or more experiments were conducted on one site), we considered the effects of sites index on SOC RR with biochar addition using the site name as a categorical variable for model training and predictions and found no improvement (same as Figure  S6). We spilt data as training and test data with rates of 8:2 and 9:1 for uncertainties analysis and the RF models developed could explain 82% and 78% of the variance in the biochar impacts on SOC change.
We also tested the cultivated crop type on the cropland as an explanatory variable in the RF model, but its contribution is low (<4%, Figure S35). However, different microbial activities in the soil applied with biochar were reported among different crop types including soybean, radish, and wheat crops (Van Zwieten et al., 2010). We thus trained one RF model for each individual main crop type (maize, rice, and wheat), but the validation R 2 is relatively low (0.29, 0.50, and 0.60 Figure S34c-e). We also made predictions using four pyrolysis temperatures (350°C, 450°C, 500°C, 600°C) for biochar production because of the relatively high importance of biochar pyrolysis temperature (5.7%, Figure 2a) and the patterns are similar with mean change of 22.9-23.5 t ha −1 ( Figure S36). The pyrolysis temperature can impact the biochar stability. Higher pyrolysis temperature leads to a smaller O/C ratio, indicating higher stability (Ippolito et al., 2020). The F I G U R E 5 Revenue map for a medium scenario with a plant biomass-processing capacity of 20,000 t year −1 and a carbon price of $ 50 t −1 CO 2 by assuming a new plant existing in each cropland grid cell biochar yields from biomass, however, are lower under a higher pyrolysis temperature (Yang et al., 2021).
The interval between biochar application and soil sampling (Age_BC) explained a 2.2% variation in the RF model ( Figure 2). Our meta-analysis results show a lower SOC increase with Age_BC ≥5 year than Age_BC <5 year on average ( Figure S13). This may be caused by the continuous decomposition of biochar in the fields, and Dong et al. (2017) found that 40% of the biochar carbon was lost after 5 years since application. The mean decomposition rate of 0.025% day −1 for crop-derived biochar was reported (Wang et al., 2016). The SOC changes to biochar application are, therefore, time-dependent. However, because most observations in this study were collected from short-term experiments with a median of 1.2 year ( Figure  S9c), our RF model may be less effective for predicting the long-term impacts of biochar on SOC changes and capturing the SOC temporal dynamics. The long-term trials of biochar applications are thus needed to investigate the mechanism of long-term effects of biochar addition on SOC dynamics.
In addition, we also tried the recursive feature elimination method with fivefold cross-validation ( Figure  S37) to distinguish the most important variables. The derived important feature number is eight, including Rate_BC, SOC_ini, Temp_BC, MAP, pH, BD, Fertilizer, Age_BC ( Figure S37). We then trained a new RF model using these eight selected variables, and the variable importance ( Figure S38) is very similar to that with all variables considered (Figure 2a). R 2 in the new RF model (0.81, Figure S39b) is the same to the original R 2 (0.81, Figure S6b).

| Uncertainties in LCA
We calculated the revenue and identified potential locations for new plants based on available crop residues for biochar production. Instead of using the default fraction (23.1%, calculated from Section 2.4) of crop residues for biochar production, we used the theoretical maximum fraction (37.7%, calculated from Section 2.4) by only removing the root and grain fraction from NPP. As expected, the estimated average transportation distance to collect crop residues for the pyrolysis plants (73.6 km) is smaller than the original value (83.2 km, Figure S30c) because of more biomass residue available around the plants. The spatial pattern of revenues for potential new plants is similar to the original result ( Figure S40). More revenues were found in North America, Europe, and Southeast Asia than using the default biomass residue fraction (Figure S40), where a larger AC SOC increase with biochar application occurred (Figure 3).
In addition to the biochar technology with low-cost and great carbon sequestration potential, other usages of agricultural residues (e.g., return of residues to soil) were also adopted widely as soil sequestration options. Crop residue retention is considered as an economical method both increasing soil fertility and increasing SOC (Zhang et al., 2015). Zhao et al. (2020) and Liu et al. (2014) found that SOC increased by 12.3% and 7.6% with agricultural residues retention to cropland, respectively. Liu et al. (2020) suggested that biochar can outperform straw retention to increase biological carbon sequestration potential due to greater microbial carbon use efficiency in biochar-amended soil. Therefore, if the crop residues are used for biochar production, the SOC sequestration from crop residue retention will not be achieved. Taking an increase of 10% SOC with crop residue retention (average of the values from Zhao et al., 2020 andLiu, He, et al., 2019) as a baseline, the global total net revenue of current pyrolysis plants decreases by 45.4% in the medium scenario ( Figure S41).
For the moment, the revenue from the LCA analysis is based on one year of biochar production and application. When the pyrolysis plants are in operation, however, biochar would be produced and applied every year. The native SOC changes may respond differently to multiple applications year after year, and the soil carbon sink may saturate after many years. The carbon contained in biochar may also decay with time. Indeed, our analyses also show the importance of initial SOC to the RR (Figure 2b) and the lower SOC increase with biochar application ≥5 years than <5 years ( Figure S13). The lack of long-term experiment data and experiments with multiple biochar applications, however, precludes us from taking these aspects into account.
In this study, we provided a global-scale LCA of biochar system with spatially-explicit SOC changes with biochar addition and transportation distances, but there remain uncertainties. Some parameters were assumed to be spatially invariant, such as fractions of crop residues used for biochar production, prices of fertilizer, labor and transportation, and energy generation, which likely vary in space and time. For example, we used a constant reduction ratio of 25% in soil N 2 O emissions with biochar addition, but one study found that an inverted U-shaped curve of soil N 2 O emissions with increasing biochar application rates . Literature data show that biochar addition can inhibit N 2 O emissions by 20%-80% (Wang et al., 2014). Although most direct processes of the biochar systems were included in the LCA, some indirectly related processes are missing, such as the heat and energy consumption of ancillary pyrolysis equipment, and other environmental benefits with biochar addition. In addition, the diverse usages of biochar other than applications in cropland were not considered in this study. Further work is needed to spatialize parameters and to extend processes in the LCA, in order to comprehensively assess the feasibility of global biochar application.

| CONCLUSION
Our study explicitly estimated the SOC changes with biochar addition on the global cropland at 1 km × 1 km resolution and conducted a preliminary LCA of the biochar system for current pyrolysis plants and the potential new plants. Long-term field experiments with biochar application will be valuable to improve the RF prediction of long-term effects and to understand the mechanisms of biochar on native SOC changes. It would also help implementing the biochar effects in process-based models [e.g., MIMICS (Wieder et al., 2014)] that can be further coupled into a land surface model to simulate the interactions of biochar-induced SOC changes with future climate change. For the LCA, more spatially explicit parameters are needed to accurately calculate the cost and revenue of the biochar system and evaluate the feasibility of biochar addition globally. The potential economic benefits of the biochar system are not high because of the high feedstock collection and pyrolysis cost, and the low carbon trading prices (Laird et al., 2009) (Figure S22). The cost of feedstock supply will further increase with the diversified utilization of agricultural waste resources under various climate mitigation options (e.g., using biomass residue for bioenergy capture and storage). More research is thus needed to analyze the trade-off among these climate mitigation options.