Where should we apply biochar?

We built our database on top of the raw data previously collected from a meta-analysis [2]. The original dataset was built upon 40 studies (published up to 2013) with 17 variables and 685 observations [2]. New peer-reviewed studies from both pot and field studies were found using academic search engines (Google Scholar, Web of Science, Scopus) and the same variables were extracted and added to the database. These search engines were searched using 'biochar AND crop productivity' OR 'biochar AND crop yield' keywords and the latest date for inclusion of the studies was 31 April 2016.

Among all studies, publications which examined the effect of biochar application on crop production (grain/biomass) were selected. We did not include studies or observations in the database if either biochar application rate or grain/biomass yield data were missing. In total, from the 63 new studies which met our criteria, 575 more observations were added to the database.

Soil organic carbon (SOC), sand, silt, clay content, CEC and soil pH were extracted from all studies to describe the chemical and physical properties of soils. Biochar carbon, nitrogen, ash content, pH, carbon-to-nitrogen (C:N) ratio, highest pyrolysis temperature (HPT), feedstock, and thermochemical process were variables extracted to account for differences in biochar type. Biochar feedstock was classified into woody, non-woody, and manure, while pyrolysis type was characterized as fast and slow.

In addition to soil and biochar properties, our model also included latitude and both N fertilizer and biochar application rates. Given that absolute yield is not readily comparable among studies, response ratio (RR) was used as the target variable [28]. The RR is defined as:

$$\begin{eqnarray}&&{\rm{RR}}=\mathrm{ln}\left(\displaystyle \frac{{{Yield}}_{{Biochar}}}{{{Yield}}_{{control}}}\right).\end{eqnarray} \tag{ 1 }$$

A positive RR indicates a positive yield response to biochar application, whereas a RR of 0 shows no change from control treatment. This variable also can be easily transformed back to percentage relative increase (RI) using ${\rm{RI}}=({e}^{{\rm{RR}}}-1)\times 100$. Many key drivers of crop yields (e.g. soil N) particularly soil nutrient levels are not included in the database because nutrient availability does not perfectly correlate to total nutrient content in soil given spatially heterogeneous soil chemistry [2]. Therefore, it is expected to see smaller variability in the models estimation of RR compared to the observed RR.

A Bayesian network (BN) was used for modeling the yield response to biochar applications. BN models usually are made of qualitative and quantitative components. The qualitative component is a graphical model which represents how the variables are statistically dependent on each other; nodes indicate variables and arcs show dependencies. The quantitative component is the conditional probability distribution of a node x_i (specified in the graphical model) on its parents ${pa}({x}_{i})$. Taking into account the conditional independence assumption (Markov condition), the joint distribution over all the variables x_i for $i=1,\,...,\,n$ is equal to [25]:

$$\begin{eqnarray}&&p({x}_{1},\,...,{x}_{n})=\displaystyle \prod _{i=1}^{n}p({x}_{i}| {pa}({x}_{i}))\end{eqnarray} \tag{ 2 }$$

which is the product of conditional distributions defined for each variable. When new evidence is introduced, it propagates through the BN and the posterior probabilities are computed. This is called an inference and it allows for detecting the change in the probabilities of some variables given a value for other variables [27].

In the current study, a hybrid BN model (a model that includes both discrete and continuous variables) was developed within the R environment [29] and by using the Bayes Server software [30]. The final design of the model's structure resulted in 84 parameters. The heterogeneity among crop species and other key drivers of yield that are not included in the model, convinced us to use the probability of yield increase P(RR > 0) instead of directly using estimated mean of RR; this inherently accounts for the variability around the average estimate of RR. Thereby, the estimated mean and variance were used to estimate P(RR > 0) and identify places with high probability of yield increase. We setup a repeated cross validation procedure with 250 iterations for training and testing our BN. In each iteration, the model was trained with 80% of the observations in the dataset and tested against the remaining 20% of observations. Predictions made by the BN were then compared with the previously developed statistical model [2]. Model efficiency (EF) and mean absolute difference (MAD) were used to compare the performance of these models as follows:

$$\begin{eqnarray}&&{\rm{EF}}=1-\displaystyle \frac{\sum {\left({Y}_{m}-{Y}_{O}\right)}^{2}}{\sum {\left({Y}_{O}-{\bar{Y}}_{O}\right)}^{2}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{\rm{MAD}}=\displaystyle \frac{| {Y}_{m}-{Y}_{O}| }{n},\end{eqnarray} \tag{ 4 }$$

where Y_m is modeled yield, Y_O is observed yield, ${\bar{Y}}_{O}$ is the mean of the observed yield and n is the number of the observations. EF varies between $-\infty$ to 1; EF = 1 represent the perfect match between modeled and observed data, whereas ${\rm{EF}}\lt 0$ indicates an unsatisfactory model performance. Biochar properties were produced from common biomass feedstock materials such as corn (C), soybean (S), switchgrass (G) and hardwood (W) using both fast (F) and slow (S) pyrolysis engineering (hereafter designated as CS, CF, GS, SF, and WS). Details on biochars properties can be found in table 1.

The final BN model was then projected onto all cultivated lands in the US based on the 2016 Cropland Data Layer (CDL) with 30 m × 30 m resolution [31, 32] using the Gridded Soil Survey Geographic (gSSURGO) [33] database for the five biochar types (table 1) and two application rates, i.e. 5 and 15 Mg ha⁻¹. The gSSURGO database provides a wide range of soil properties in 10 m × 10 m resolution including all the variables required for our model [34]. A pixel in the map is considered cultivated if it is identified as cultivated in at least two out of the five years of CDL data. The gSSURGO database includes the heterogeneity of soil properties (as different compartments) for each pixel. In the current study, the most representative soil type was used to extract the basic soil properties.

The estimated average RR to biochar showed a 12% increase for all studies in our database. A large variability in RR was also observed, ranging from −24.4% to 98%, with the interquartile range ranging from 0 to 21% (figure 1). Among all soil properties, clay content, SOC, pH and CEC showed a significant negative correlation with RR, while sand and silt content were positively correlated with RR. Yield response was invariant with nitrogen application rate, biochar ash content, and biochar pH. The HPT, biochar N, and C:N ratio showed a significant negative correlation with RR. Higher biochar C content was significantly correlated with a higher RR (figure 1). A linear model analysis revealed a minor (not significant) association between feedstock and crop type with RR while no direct association was found between thermochemical technology and RR. Note that we do not make any assumptions about the functional form of yield response to biochar application. The model predicts that the biochar application rate has a diminishing marginal effect on the yield response, i.e. the yield effect is higher with low biochar application rates and flattens out at higher application rates.

Zoom In Zoom Out Reset image size

We compared our BN model and the only available statistical model for explaining the heterogeneity in yield response to soil and biochar properties [2]. The authors of that study used a generalized additive model (GAM) with 162 parameters to develop smooth functions mapping of independent variables to RR. In more than 250 iterations, the BN model consistently outperformed the GAM; as the BN model average EF and MAD values were 0.23 and 0.10, respectively, compared to −1.96 and 0.18 for the GAM. Negative EF implies that the observed mean is a better predictor than the model; and the positive average EF proves the merit of the BN and indicates that it can be used to explore different scenarios of soil and biochar properties within the scope of the training dataset. Given the proficiency of the BN model, spatially explicit analysis of response to biochar was explored for cropland in the US under different biochar scenarios. Figure 2 shows the estimated probability, magnitude, and expected yield increase for hardwood biochar and 15 Mg ha⁻¹ application rate scenarios. We focus on locations that have a positive expected yield increase because farmer's will not apply biochar if the expected yield change is negative. This assumes that the only incentive for farmers to apply biochar is the yield increase.

Zoom In Zoom Out Reset image size

Regions known to have high soil quality (e.g. Des Moines lobe in north central Iowa) showed a low probability of having a yield increase under all biochar scenarios. Our model indicates a high probability of a positive RR in areas with highly weathered soils (e.g. Eastern half of San Joaquin valley in California). Yield response was predicted to be the lowest in areas with very high SOC, CEC, or soil pH such as those found in north Texas and Minnesota (figure 2).

Our results show the highest increase in expected yields of 6.43% for biochar derived from hardwood with an application rate of 15 Mg ha⁻¹. The lowest increase increase of 4.69% was observed for biochar derived from soybeans at an application rate of 5 Mg ha⁻¹. Assuming farmers with the highest expected yield increase are the first to apply biochar to the field, then minimum yield increases ranging from 8.8% (soybeans biochar at 5 Mg ha⁻¹) to 11.4% (hardwood biochar 15 Mg ha⁻¹) are necessary to cover 10% of crop area. A minimum yield increases ranging from 6.1% (switchgrass biochar at 5 Mg ha⁻¹) to 7.9% (hardwood biochar 15 Mg ha⁻¹) are necessary to cover 25% of crop area.

We estimated the total area where a 75% (or higher) probability of yield increase is expected for each scenario as well as areas where this is expected to be lower than 25%. This helped us identify the most and least responsive regions to biochar applications across the United States. Total high probability areas for biochar application range from 8.4% to 30% of total cropland in the United States. The slow pyrolysis hardwood biochar with an application rate of 15 Mg ha⁻¹ was the scenario with the largest high probability area with 39.7 Mha, whereas fast pyrolysis soybean biochar with an application rate of 5 Mg ha⁻¹ resulted in the lowest high probability area with 11.2 million ha. In general, for similar application rates, biochar produced from hardwood and corn stover resulted in larger responsive areas than switchgrass, while fast pyrolysis soybean biochar showed the lowest response. In low-quality soils, the higher application rate (15 Mg ha⁻¹) resulted in a higher probability of a yield increase, whereas high quality soils showed no response to application rate.

We selected the Central Valley of California to assess the model's response to different soil properties under the slow pyrolysis from corn stover scenario with an application rate of 15 Mg ha⁻¹. The Central Valley is an agricultural region drained by the Sacramento and San Joaquin rivers. It is about 82 kilometers wide and extends 600 km northwest from the Tehachapi Mountains to Redding. The southern part of the valley, also known as San Joaquin valley, is well-known for having highly variable alluvial soils ranging from very acidic-low organic matter soils to very alkaline-high organic matter soils. The BN model successfully captured the essence of our general understanding of crop yield response to biochar on high and low-quality soils. Yield response was the weakest in west side of the San Joaquin valley, which is dominated by soils with high clay, pH, CEC, and organic matter (figure 3). On the other hand, the old, highly weathered soils on the east side of the valley, which have low organic matter, CEC, and pH, showed the greatest response to biochar applications. The model's response was not attributed to just one variable, as all soil properties contributed to the estimation of yield response (figure 3). The total area with high probability of yield response to biochar application varied from 0.01 Mha for fast pyrolysis soybean biochar (5 Mg ha⁻¹) to 0.6 Mha for slow pyrolysis hardwood biochar (15 Mg ha⁻¹) which is approximately 1% and 15% of total cropland in the Central Valley, respectively.

Zoom In Zoom Out Reset image size

There are various economic aspects of biochar application such as expected revenue from applying biochar, cost of biochar application, commodity price effects of biochar-induced yield increases in the long-run, or potential revenue from the provision of environmental services (e.g. increased carbon sequestration and/or reduction in nitrogen leaching). For this case study, we focus only on the increase in revenue at the county level and the long-term effects on commodity prices. The differences in the expected yield increase between the 5 Mg ha⁻¹ and the 15 Mg ha⁻¹ biochar application rate is not sufficient to justify the higher application rate and thus, we focus on the lower rate. We assume that the application of biochar occurs in the first year and lasts for the remainder of the projection period, i.e., 10 and 20 years.

To determine the farmers' willingness to pay per ton of biochar, we calculate the net present value (NPV) of the additional revenue triggered by the expected increase in yield over a period of 10 and 20 years (figure 4). Table 2 shows the crop area covered under three different biochar prices over 10 and 20 years. If the yield increase due to biochar is limited to 10 years, the crop area covered is much smaller because the timeframe to recuperate the initial cost is shorter. Because we calculated the NPV per ton of biochar applied, the values in figure 4 represent farmers' willingness to pay per ton of biochar. The highest NPV is observed for corn followed by soybeans and wheat. The highest revenue increases are observed in the southeast for corn and the east for soybeans. Note that those areas do not coincide with large corn and soybean acreage. For wheat and soybeans, the revenue gains are more moderate compared to corn suggesting that the best use of biochar is its application to corn area.

Zoom In Zoom Out Reset image size

Large-scale application of biochar on US agricultural land will have long-term commodity price effects because the increase in supply will decrease prices. The decrease in commodity prices from the expected yield increase is about 7.5% and 7.6% for soybeans and corn, respectively and 6.3% for wheat based on an economic simulation model for US land-use change [36, 37]. An increase in commodity prices is observed when biochar-induced yield increase is compensated by land being taken out of food and feed crop production. At a biomass price of approximately $70 per dry ton ($4 GJ⁻¹), we see an average increase in corn, soybean, and wheat prices of 6.3%, 10%, and 15.8%, respectively. At a biomass price of $4 GJ⁻¹, our model estimates a biomass production of 117–119 million tons per year from both, corn stover and switchgrass under a low biomass production cost scenario and 87–89 million tons at a biomass price of $5 GJ⁻¹ under high biomass production cost. The effect on commodity prices of biochar-induced yield increases is speculative but highlights one additional aspect of biochar application that can become important under large-scale deployment.

There is a wide range of biochar prices which is the reason why we report willingness-to-pay to avoid any cost assumptions. The biochar price is assumed to range from $87 to $351 per metric ton for farmers [38] but small scale manufacturing costs are estimated to be closer to $500 per metric ton [39].