Soil Carbon Sequestration

import matplotlib.pyplot as plt
import numpy as np
from scipy import stats

np.random.seed(999)

# Soil carbon data
depth = np.random.uniform(0, 100, 50)  # cm
carbon = 45 * np.exp(-0.025 * depth) + np.random.normal(0, 3, 50)

# Exponential fit on log scale
log_carbon = np.log(carbon)
slope, intercept, r_value, _, _ = stats.linregress(depth, log_carbon)

d_fit = np.linspace(0, 100, 100)
c_fit = np.exp(intercept) * np.exp(slope * d_fit)

fig, ax = plt.subplots(figsize=(10, 7), facecolor='#ffffff')
ax.set_facecolor('#ffffff')

ax.yaxis.grid(True, color='#f0f0f0', linewidth=1, zorder=1)
ax.xaxis.grid(True, color='#f0f0f0', linewidth=1, zorder=1)

# Soil horizon zones
ax.axvspan(0, 20, alpha=0.08, color='#9C2007', zorder=1)
ax.axvspan(20, 50, alpha=0.05, color='#F5B027', zorder=1)
ax.text(10, 5, 'O/A', fontsize=9, color='#9C2007', ha='center')
ax.text(35, 5, 'B', fontsize=9, color='#F5B027', ha='center')

ax.fill_between(d_fit, c_fit * 0.85, c_fit * 1.15, color='#27F5B0', alpha=0.15, linewidth=0, zorder=2)
ax.plot(d_fit, c_fit, color='#27F5B0', linewidth=2.5, zorder=3)
ax.scatter(depth, carbon, c='#9C2007', s=55, alpha=0.85, edgecolors='white', linewidths=0.6, zorder=4)

# Decay constant
ax.text(0.97, 0.95, f'C(z) = {np.exp(intercept):.1f}e^({slope:.3f}z)\nτ = {-1/slope:.1f} cm',
        transform=ax.transAxes, fontsize=10, color='#555555', ha='right', va='top', family='monospace')

for spine in ['top', 'right']:
    ax.spines[spine].set_visible(False)
for spine in ['bottom', 'left']:
    ax.spines[spine].set_color('#cccccc')

ax.set_xlabel('Soil Depth (cm)', fontsize=12, color='#333333', fontweight='500', labelpad=10)
ax.set_ylabel('Organic Carbon (g/kg)', fontsize=12, color='#333333', fontweight='500', labelpad=10)
ax.set_title('Soil Carbon Profile', fontsize=15, color='#1a1a1a', fontweight='bold', pad=20, loc='left')
ax.tick_params(colors='#666666', labelsize=10, length=0)

plt.tight_layout()
plt.show()