biochar/biochar.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img width=\"50\" src=\"https://carbonplan-assets.s3.amazonaws.com/monogram/dark-small.png\" style=\"margin-left:0px;margin-top:20px\"/>\n",
    "\n",
    "# Biochar lifetime analysis\n",
    "\n",
    "_by Jeremy Freeman (CarbonPlan), Created May 17, 2020, Last Updated May 24,\n",
    "2021_\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we present a simple toy model for evaluating the carbon removal and\n",
    "permanence of biochar projects. The data and analysis method is based directly\n",
    "on two publications\n",
    "\n",
    "- Spokas (2010) Review of the stability of biochar in soils: predictability of\n",
    "  O:C molar ratios, Carbon Management, doi: 10.4155/CMT.10.32\n",
    "\n",
    "- Campbell et al. (2018) Potential carbon storage in biochar made from logging\n",
    "  residue: Basic principles and Southern Oregon case studies, PLOS One, doi:\n",
    "  10.1371/journal.pone.0203475\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Notebook setup\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "from carbonplan_styles.colors import colors\n",
    "from carbonplan_styles.mpl import set_theme\n",
    "\n",
    "set_theme(style=\"carbonplan_light\")\n",
    "c = colors(\"carbonplan_light\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The basic model\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Campbell et al. (2018) present a simple model for biochar carbon dynamics by\n",
    "comparing the carbon content of biomass after biocharing to the carbon content\n",
    "that would have resided in the form of the source feedstock (e.g. logging\n",
    "residues).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The difference is\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$∆ = C_{biochar} - C_{biochar}$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the mass of carbon in both is modeled using a first-order differential\n",
    "equation\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$C_t = C_{t-1}e^{-k} + C_{input}$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll write a function that generates a complete carbon curve as a function of\n",
    "the input and the parameter k over 1000 years\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def model(t, initial, k):\n",
    "    return initial * np.exp(-k * t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can now plot carbon curves for both unmodified residue and biochar over a\n",
    "fixed duration, assuming an initial carbon content of 20 tC for the residue and\n",
    "12 tC for the biochar (which would be achieved through a pyrolysis process with\n",
    "60% efficiency).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.arange(1000)\n",
    "residue = model(t, 20, 0.03)\n",
    "biochar = model(t, 12, 0.003)\n",
    "plt.plot(t, residue)\n",
    "plt.plot(t, biochar)\n",
    "plt.xlim([0, 200])\n",
    "plt.ylim([0, 20])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This precisely matches Figure 1A from Campbell et al. (2019)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These curves makes clear that biochar is not removing carbon per se, but rather\n",
    "avoiding the emissions that would have been associated with the corresponding\n",
    "feedstock. For that reason, the appropriate quantity is the difference between\n",
    "the two curves.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plt.plot(t, biochar - residue)\n",
    "plt.xlim([0, 200])\n",
    "plt.ylim([-10, 10])\n",
    "plt.hlines(0, 0, 200, color=c[\"secondary\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As this curve makes clear, the cumulative effective carbon removal is initially\n",
    "negative, quickly reaches a compensation point, and then reaches a point termed\n",
    "by Campbell et al. (2018) as \"climate parity\" where the storage\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mapping O:C ratios to half life\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A key parameter in the above model is the decay rate (also referred to as the\n",
    "biochar's recalcitrance). Campbell et al. (2018) find that this parameter has\n",
    "little effect on the time at which climate parity is achieved, so long as it 10\n",
    "times greater than the decay rate feedstock. But it is also importantly related\n",
    "to the permanence, or time scale over which the carbon stored in the biochar\n",
    "remain.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use data digitized from a meta-analysis by Spokas (2010) that relate the\n",
    "oxygen to carbon (O:C) molar ratio to the predicted half-life of synthetic\n",
    "biochar in various laboratory conditions.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_csv(\"biochar.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(data.ratio, data.halflife, \".\", color=c[\"primary\"])\n",
    "plt.xlim([0, 0.8])\n",
    "plt.ylim([1, 10**8])\n",
    "plt.yscale(\"log\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We fit a simple linear model in log space so we can predict half-life as a\n",
    "function of ratio. In order to put bounds on our estimates, we use a simple\n",
    "bootstrap to fit the model for each of 1000 random samples (with replacement)\n",
    "from the data. We store the parameter estimates from each sample, and plot a\n",
    "regression line.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 10000\n",
    "plt.plot(data.ratio, data.halflife, \".\", color=c[\"primary\"])\n",
    "xhat = np.arange(0, 1, 0.1)\n",
    "indices = np.arange(34)\n",
    "alpha = np.zeros(k)\n",
    "beta = np.zeros(k)\n",
    "for i in range(k):\n",
    "    samples = np.random.choice(indices, 34)\n",
    "    mod = sm.OLS(\n",
    "        np.log(data.halflife[samples]),\n",
    "        sm.add_constant(data.ratio[samples], prepend=False),\n",
    "    )\n",
    "    res = mod.fit()\n",
    "    alpha[i] = res.params[1]\n",
    "    beta[i] = res.params[0]\n",
    "    yhat = res.predict(sm.add_constant(xhat, prepend=False))\n",
    "    if i % 10 == 0:\n",
    "        plt.plot(xhat, np.exp(yhat), \"-\", color=\"red\", alpha=0.005)\n",
    "plt.yscale(\"log\")\n",
    "plt.xlim([0, 0.8])\n",
    "plt.ylim([1, 10**8])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we write a simple function that, for a given ratio, returns a prediction\n",
    "from the bootstrapped distribution at a given percentile.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def predict(ratio, prctile):\n",
    "    dist = np.exp(alpha + beta * ratio)\n",
    "    return np.percentile(dist, [prctile])[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Project evaluation\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fixed fraction permanence\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now use the above to evaluate some aspects of a biochar projects. If we\n",
    "assume a project reports an O:C ratio of 0.08, we can use the simple linear\n",
    "model above to we can compute a half-life. We use the 2.5th percentile of the\n",
    "posterior predictive distribution as a crude, highly conservative estimate,\n",
    "given that permanence is only weakly correlated with composition, and likely\n",
    "depends as much or more so on the decay environment, which is often unknown.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ratio = 0.09\n",
    "halflife = predict(ratio, 2.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Still, given the decay kinetics assumed by our toy model, we can compute a decay\n",
    "constant from the half-life\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = np.log(2) / halflife"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now determine the duration after which a fixed percent of the biochar\n",
    "remains. For a target of 90% for example, we get the following number of years.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fraction = 0.9\n",
    "years = -np.log(fraction) / k"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can summarize our parameters\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"summary\")\n",
    "print(\"-------\")\n",
    "print(\"ratio: \" + str(ratio))\n",
    "print(\"half-life: \" + str(halflife) + \" years\")\n",
    "print(\"fraction: \" + str(fraction))\n",
    "print(\"k: \" + str(k))\n",
    "print(\"years: \" + str(years))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can plot this on the decay curve from above, assuming a initial volume of\n",
    "carbon storage in the biochar (tC).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "initial = 100\n",
    "t = np.arange(0, 20000)\n",
    "biochar = model(t, initial, k)\n",
    "plt.plot(t, biochar)\n",
    "plt.ylim([0, initial])\n",
    "plt.xlim([0, 2000])\n",
    "plt.vlines(years, 0, initial)\n",
    "plt.hlines(initial * fraction, 0, 20000, color=c[\"secondary\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In general, validating the volume and permanence for an actual biochar project\n",
    "requires knowing the composition (and thus recalcitrance), but perhaps more\n",
    "importantly, also requires knowing the conversion efficiency (the fraction of\n",
    "initial feedstock carbon retained in biochar after pyrolysis) and the decay rate\n",
    "of the feedstock. That said, simply by knowing the recalcitrance, and making\n",
    "some assumptions, we can approximate a permanence over which a fixed fraction of\n",
    "volume is likely to remain.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Counterfactual feedstock decay\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Per the Campell et al. (2018), biochar acheives carbon storage by decaying more\n",
    "slowly than its feedstock.\n",
    "\n",
    "Using the approximate permanence horizon calculated in the section above, we can\n",
    "ask how quickly the feedstock would have had to decay for the counterfactual\n",
    "carbon storage in the feedstock to be considered negligible.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assuming a pyrolysis efficiency (e.g. 60%), we can estimate the starting carbon\n",
    "storage of the feedstock relative to the biochar.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "efficiency = 0.6\n",
    "feedstock_start = initial / efficiency"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By setting a bar for \"negligible impact\" (e.g. feedstock carbon storage must be\n",
    "<0.5% of biochar carbon storage at the end of the permanence period), we can\n",
    "calculate an upper bound for feedstock carbon storage.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "negligible = 0.005\n",
    "feedstock_end = (initial * fraction) * negligible"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now determine a minimum decay constant for a feedstock's counterfactual\n",
    "carbon storage to be considered negligible. (As a reminder, lower decay constant\n",
    "means slower decay!)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k_feedstock = -np.log(feedstock_end / feedstock_start) / years"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot the feedstock and biochar decay curves over the permanence period\n",
    "calculated above.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.arange(0, years)\n",
    "biochar = model(t, initial, k)\n",
    "feedstock = model(t, feedstock_start, k_feedstock)\n",
    "plt.plot(t, biochar)\n",
    "plt.plot(t, feedstock)\n",
    "plt.ylim([0, initial / efficiency])\n",
    "plt.xlim([0, years])\n",
    "plt.vlines(years, 0, initial * fraction, color=c[\"secondary\"])\n",
    "plt.hlines(initial * 0.90, 0, 20000, color=c[\"secondary\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The net carbon storage at the end of the permanence period is the difference\n",
    "between the mass of biochar carbon storage year and the counterfactual feedstock\n",
    "carbon storage.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NCS = (biochar - feedstock)[years.astype(int)]\n",
    "print(str(np.round(NCS)) + \" tC\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can summarize our parameters and outputs:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"summary\")\n",
    "print(\"-------\")\n",
    "print(\"years: \" + str(np.round(years)) + \" years\")\n",
    "print(\"efficiency: \" + str(efficiency * 100) + \"%\")\n",
    "print(\"negligible impact threshold: \" + str(negligible * 100) + \"%\")\n",
    "print(\"min feedstock k: \" + str(np.round(k_feedstock, 3)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compare this feedstock decay rate bound against values found in\n",
    "literature to gain intuition about how important it is to take into account the\n",
    "feedstock counterfactual when crediting biochar.\n",
    "\n",
    "Publications we have queried for this information include:\n",
    "\n",
    "- Harmon et al. (2020) Release of coarse woody detritus-related carbon: a\n",
    "  synthesis across forest biomes, Carbon Balance Management, doi:\n",
    "  10.1186/s13021-019-0136-6\n",
    "\n",
    "- Ximenes et al. (2017) The decay of engineered wood products and paper\n",
    "  excavated from landfills in Australia, Waste Management, doi:\n",
    "  10.1016/j.wasman.2017.11.035\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}