Soil Carbon Analysis

1. Setup

1.1 Packages

library(tidyverse)   # dplyr, tidyr, ggplot2, readr
library(broom)       # tidy() model summaries
library(knitr)       # kable() tables
library(kableExtra)  # styled tables

1.2 Simulated soil survey data

The dataset is modelled on soil surveys from paired land-use plots in Southern California: native chaparral, annual grassland (invaded by non-native species), cultivated cropland, and restored native (10-year revegetation). Each land use has 5 replicate sites, each sampled at standard horizon mid-depths (5, 15, 30, 50, 75, and 100 cm).

set.seed(2024)

land_uses <- c("native_chaparral", "annual_grassland",
               "cultivated",       "restored_native")
n_sites   <- 5    # replicates per land use
depths_cm <- c(5, 15, 30, 50, 75, 100)

# Exponential decay parameters (A = surface SOC g/kg, k = decay rate /cm)
# Bulk density (g/cm³) increases with depth for all land uses
params <- tribble(
  ~land_use,          ~A,    ~k,      ~bd_surface, ~bd_deep,
  "native_chaparral",  62,   0.028,   0.82,        1.25,
  "annual_grassland",  34,   0.032,   0.95,        1.32,
  "cultivated",        22,   0.038,   1.10,        1.45,
  "restored_native",   48,   0.030,   0.88,        1.28
)

soil_raw <- params %>%
  rowwise() %>%
  reframe(
    land_use    = land_use,
    site_id     = paste0(substr(land_use, 1, 3), "_", str_pad(1:n_sites, 2, "left", "0")),
    depth_cm    = list(depths_cm)
  ) %>%
  unnest(depth_cm) %>%
  left_join(params, by = "land_use") %>%
  mutate(
    # SOC concentration: exponential decay + site-level random effect + noise
    site_effect = rnorm(n(), 0, 4),
    soc_g_kg    = pmax(0.5,
                    A * exp(-k * depth_cm) +
                    site_effect * exp(-0.02 * depth_cm) +
                    rnorm(n(), 0, 1.5)),
    # Bulk density: linear increase from surface to 100 cm
    bulk_density_g_cm3 = bd_surface +
                         (bd_deep - bd_surface) * depth_cm / 100 +
                         rnorm(n(), 0, 0.03),
    bulk_density_g_cm3 = pmax(0.6, bulk_density_g_cm3),
    # Horizon thickness (cm) for stock calculation
    horizon_thick_cm   = case_when(
      depth_cm ==   5 ~ 10,
      depth_cm ==  15 ~ 10,
      depth_cm ==  30 ~ 20,
      depth_cm ==  50 ~ 20,
      depth_cm ==  75 ~ 25,
      depth_cm == 100 ~ 25
    )
  ) %>%
  select(land_use, site_id, depth_cm, horizon_thick_cm,
         soc_g_kg, bulk_density_g_cm3)

glimpse(soil_raw)

## Rows: 120
## Columns: 6
## $ land_use           <chr> "native_chaparral", "native_chaparral", "native_cha…
## $ site_id            <chr> "nat_01", "nat_01", "nat_01", "nat_01", "nat_01", "…
## $ depth_cm           <dbl> 5, 15, 30, 50, 75, 100, 5, 15, 30, 50, 75, 100, 5, …
## $ horizon_thick_cm   <dbl> 10, 10, 20, 20, 25, 25, 10, 10, 20, 20, 25, 25, 10,…
## $ soc_g_kg           <dbl> 57.056398, 41.632982, 24.573514, 14.679824, 7.55876…
## $ bulk_density_g_cm3 <dbl> 0.8480238, 0.9172507, 0.9899496, 1.0690255, 1.11519…

📄 Using your own data. If you are bringing real soil survey data from the EA030 field campaign, swap in your CSV at the import step below. Your file needs at minimum: a site identifier, a depth column (cm), an SOC column (g/kg or %), and a bulk density column (g/cm³). If bulk density was not measured, use a published pedotransfer function — the Saxton-Rawls equations are a common choice for California soils.

2. Import and Clean

2.1 Write and re-read as a CSV (simulating a real import)

In a real project you would load your data with read_csv(). Here we write the simulated data to disk and read it back so the workflow is identical.

# Write to a temporary file — replace this path with your real data file
tmp <- tempfile(fileext = ".csv")
write_csv(soil_raw, tmp)

# Import exactly as you would in a real project
soil <- read_csv(tmp, show_col_types = FALSE)
glimpse(soil)

## Rows: 120
## Columns: 6
## $ land_use           <chr> "native_chaparral", "native_chaparral", "native_cha…
## $ site_id            <chr> "nat_01", "nat_01", "nat_01", "nat_01", "nat_01", "…
## $ depth_cm           <dbl> 5, 15, 30, 50, 75, 100, 5, 15, 30, 50, 75, 100, 5, …
## $ horizon_thick_cm   <dbl> 10, 10, 20, 20, 25, 25, 10, 10, 20, 20, 25, 25, 10,…
## $ soc_g_kg           <dbl> 57.056398, 41.632982, 24.573514, 14.679824, 7.55876…
## $ bulk_density_g_cm3 <dbl> 0.8480238, 0.9172507, 0.9899496, 1.0690255, 1.11519…

2.2 Quality checks

Before any analysis, verify that the data look physically reasonable.

cat("=== Basic dimensions ===\n")

## === Basic dimensions ===

cat("Rows:", nrow(soil), "| Sites:", n_distinct(soil$site_id),
    "| Land uses:", n_distinct(soil$land_use), "\n\n")

## Rows: 120 | Sites: 20 | Land uses: 4

cat("=== Missing values per column ===\n")

## === Missing values per column ===

soil %>% summarise(across(everything(), ~ sum(is.na(.)))) %>% print()

## # A tibble: 1 × 6
##   land_use site_id depth_cm horizon_thick_cm soc_g_kg bulk_density_g_cm3
##      <int>   <int>    <int>            <int>    <int>              <int>
## 1        0       0        0                0        0                  0

cat("\n=== SOC range (g/kg) ===\n")

## 
## === SOC range (g/kg) ===

summary(soil$soc_g_kg)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.500   4.612  10.982  15.605  23.845  57.056

cat("\n=== Bulk density range (g/cm³) ===\n")

## 
## === Bulk density range (g/cm³) ===

summary(soil$bulk_density_g_cm3)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.7653  0.9940  1.1243  1.1172  1.2448  1.4635

# Flag physically implausible values
soil <- soil %>%
  mutate(
    flag_soc = soc_g_kg < 0 | soc_g_kg > 200,
    flag_bd  = bulk_density_g_cm3 < 0.4 | bulk_density_g_cm3 > 2.0
  )

cat("\n=== Flagged rows ===\n")

## 
## === Flagged rows ===

soil %>% filter(flag_soc | flag_bd) %>% print()

## # A tibble: 0 × 8
## # ℹ 8 variables: land_use <chr>, site_id <chr>, depth_cm <dbl>,
## #   horizon_thick_cm <dbl>, soc_g_kg <dbl>, bulk_density_g_cm3 <dbl>,
## #   flag_soc <lgl>, flag_bd <lgl>

⚠️ Common real-data problems to watch for. SOC values above 200 g/kg suggest an organic horizon (O horizon) that may need separate treatment. Bulk density below 0.4 g/cm³ is physically unrealistic for mineral soil. Negative values in either column indicate data-entry errors. Check them before proceeding.

2.3 Standardise land-use labels and factor ordering

soil <- soil %>%
  filter(!flag_soc, !flag_bd) %>%         # remove any flagged rows
  select(-flag_soc, -flag_bd) %>%
  mutate(
    land_use = factor(land_use,
      levels = c("native_chaparral", "restored_native",
                 "annual_grassland",  "cultivated"),
      labels = c("Native chaparral", "Restored native",
                 "Annual grassland",  "Cultivated")
    )
  )

# Confirm sample size per land use × depth combination
soil %>%
  count(land_use, depth_cm) %>%
  pivot_wider(names_from = depth_cm, values_from = n) %>%
  kbl(caption = "Sample counts per land use and depth (cm)") %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE)

Sample counts per land use and depth (cm)
land_use	5	15	30	50	75	100
Native chaparral	5	5	5	5	5	5
Restored native	5	5	5	5	5	5
Annual grassland	5	5	5	5	5	5
Cultivated	5	5	5	5	5	5

3. Depth Profile Visualisation

3.1 Profile plot with depth on the vertical axis

In soil science, depth always increases downward. We achieve this in ggplot2 with scale_y_reverse() on a flipped-coordinate plot.

lu_colors <- c(
  "Native chaparral" = "#1e3a2f",
  "Restored native"  = "#4a7c5e",
  "Annual grassland" = "#c8832a",
  "Cultivated"       = "#8b4513"
)

ggplot(soil, aes(x = soc_g_kg, y = depth_cm,
                 color = land_use, group = interaction(land_use, site_id))) +
  geom_line(alpha = 0.35, linewidth = 0.6) +
  geom_point(alpha = 0.55, size = 1.8) +
  # Add land-use mean profile
  stat_summary(
    aes(group = land_use),
    fun   = mean,
    geom  = "line",
    linewidth = 1.4,
    alpha = 0.95
  ) +
  scale_y_reverse(breaks = depths_cm) +
  scale_color_manual(values = lu_colors, name = NULL) +
  labs(
    title    = "Soil organic carbon depth profiles by land use",
    subtitle = "Thin lines = individual sites; thick lines = land-use means",
    x        = "SOC concentration (g kg\u207b\u00b9)",
    y        = "Depth (cm)"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title     = element_text(family = "serif", face = "bold",
                                  color = "#1e3a2f", size = 13),
    plot.subtitle  = element_text(color = "#7a7a72", size = 10),
    legend.position = "bottom",
    panel.grid.minor = element_blank()
  )

📌 Reading a depth profile. SOC is always highest at the surface where plant litter inputs are greatest and declines with depth as microbial decomposition and root inputs diminish. The steepness of the decline captures how deeply carbon is stored — native vegetation with deep roots typically maintains higher SOC at depth than disturbed or cultivated soils.

3.2 Faceted profiles — one panel per land use

# Pre-compute mean ± 1 SD per land use × depth for the ribbon
soil_ribbon <- soil %>%
  group_by(land_use, depth_cm) %>%
  summarise(
    soc_mean = mean(soc_g_kg),
    soc_lo   = mean(soc_g_kg) - sd(soc_g_kg),
    soc_hi   = mean(soc_g_kg) + sd(soc_g_kg),
    .groups  = "drop"
  )

ggplot(soil, aes(x = soc_g_kg, y = depth_cm, color = land_use)) +
  geom_line(aes(group = site_id), alpha = 0.4, linewidth = 0.7) +
  geom_point(size = 2, alpha = 0.6) +
  # Mean ± 1 SD ribbon from pre-computed data
  geom_ribbon(data = soil_ribbon,
              aes(x = soc_mean, y = depth_cm,
                  xmin = soc_lo, xmax = soc_hi,
                  fill = land_use),
              alpha = 0.18, color = NA) +
  # Mean profile line from pre-computed data
  geom_line(data = soil_ribbon,
            aes(x = soc_mean, y = depth_cm, color = land_use),
            linewidth = 1.6, alpha = 0.95) +
  scale_y_reverse(breaks = depths_cm) +
  scale_color_manual(values = lu_colors, name = NULL) +
  scale_fill_manual(values  = lu_colors, guide = "none") +
  facet_wrap(~ land_use, nrow = 1) +
  labs(
    title    = "SOC depth profiles — mean ± 1 SD by land use",
    subtitle = "Shaded band = mean ± 1 SD across replicate sites",
    x        = "SOC (g kg\u207b\u00b9)",
    y        = "Depth (cm)"
  ) +
  theme_minimal(base_size = 11) +
  theme(
    plot.title      = element_text(family = "serif", face = "bold",
                                   color = "#1e3a2f", size = 13),
    plot.subtitle   = element_text(color = "#7a7a72", size = 9.5),
    legend.position = "none",
    strip.text      = element_text(face = "bold", size = 9.5),
    panel.grid.minor = element_blank()
  )

4. Fitting Depth-Profile Models

Two model families are commonly used for SOC depth profiles:

Model	Formula	When to use
Exponential decay	SOC = A · e^(−k · depth)	Physically interpretable; A = surface SOC, k = decay rate
Polynomial (log-linear)	log(SOC) = β₀ + β₁·depth	Flexible; easy to fit with `lm()`; less interpretable

We fit both, compare AIC, and visualise the fits with confidence intervals.

4.1 Exponential decay with `nls()`

nls() (non-linear least squares) fits models where parameters enter non-linearly. We must supply starting values — use the known true values as starting guesses.

# Fit one exponential model per land use
exp_fits <- soil %>%
  group_by(land_use) %>%
  group_modify(~ {
    tryCatch(
      broom::tidy(
        nls(soc_g_kg ~ A * exp(-k * depth_cm),
            data  = .x,
            start = list(A = 50, k = 0.03),
            control = nls.control(maxiter = 200))
      ) %>% mutate(.rows = nrow(.x)),
      error = function(e) tibble(term = "error", estimate = NA_real_,
                                 std.error = NA_real_, statistic = NA_real_,
                                 p.value = NA_real_, .rows = nrow(.x))
    )
  }) %>%
  ungroup()

exp_fits %>%
  select(land_use, term, estimate, std.error, p.value) %>%
  mutate(across(where(is.numeric), ~ round(.x, 4))) %>%
  kbl(caption = "Exponential decay parameters by land use  (SOC = A · exp(−k · depth))") %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE)

Exponential decay parameters by land use (SOC = A · exp(−k · depth))
land_use	term	estimate	std.error
Native chaparral	A	62.1523	1.3858
Native chaparral	k	0.0294	0.0012
Restored native	A	44.8325	1.2097
Restored native	k	0.0284	0.0014
Annual grassland	A	34.4527	1.8000
Annual grassland	k	0.0307	0.0029
Cultivated	A	23.1039	1.4244
Cultivated	k	0.0372	0.0040

⚙️ Interpreting A and k. A (g/kg) is the modelled SOC at depth = 0 (the soil surface). k (/cm) is the exponential decay rate — larger k means carbon declines more steeply with depth. Native chaparral should have the highest A and one of the lower k values (deep carbon storage), while cultivated soil should have low A and high k (carbon concentrated near the ploughed surface).

4.2 Log-linear model with `lm()`

lm_fits <- soil %>%
  group_by(land_use) %>%
  group_modify(~ broom::glance(
    lm(log(soc_g_kg) ~ depth_cm, data = .x)
  )) %>%
  ungroup()

lm_fits %>%
  select(land_use, r.squared, AIC, BIC) %>%
  mutate(across(where(is.numeric), ~ round(.x, 3))) %>%
  kbl(caption = "Log-linear model fit statistics by land use") %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE)

Log-linear model fit statistics by land use
land_use	r.squared	AIC	BIC
Native chaparral	0.959	-11.188	-6.984
Restored native	0.947	2.535	6.739
Annual grassland	0.926	27.574	31.777
Cultivated	0.823	46.968	51.172

4.3 Profile plot with model fits and confidence intervals

# Generate predicted values over a fine depth grid for each land use
depth_grid <- tibble(depth_cm = seq(2, 105, by = 1))

pred_curves <- soil %>%
  group_by(land_use) %>%
  group_modify(~ {
    fit <- tryCatch(
      nls(soc_g_kg ~ A * exp(-k * depth_cm),
          data  = .x,
          start = list(A = 50, k = 0.03),
          control = nls.control(maxiter = 200)),
      error = function(e) NULL
    )
    if (is.null(fit)) return(tibble())
    # Confidence interval via parametric bootstrap
    set.seed(42)
    pred_boot <- replicate(500, {
      coefs <- MASS::mvrnorm(1,
                             mu    = coef(fit),
                             Sigma = vcov(fit))
      coefs["A"] * exp(-coefs["k"] * depth_grid$depth_cm)
    })
    tibble(
      depth_cm  = depth_grid$depth_cm,
      fitted    = predict(fit, newdata = depth_grid),
      ci_lo     = apply(pred_boot, 1, quantile, 0.025),
      ci_hi     = apply(pred_boot, 1, quantile, 0.975)
    )
  }) %>%
  ungroup()

ggplot() +
  # Raw data points
  geom_point(data = soil,
             aes(x = soc_g_kg, y = depth_cm, color = land_use),
             alpha = 0.4, size = 1.8) +
  # CI ribbon
  geom_ribbon(data = pred_curves,
              aes(x = fitted, y = depth_cm,
                  xmin = ci_lo, xmax = ci_hi,
                  fill = land_use),
              alpha = 0.18) +
  # Fitted curve
  geom_line(data = pred_curves,
            aes(x = fitted, y = depth_cm, color = land_use),
            linewidth = 1.3) +
  scale_y_reverse(breaks = depths_cm) +
  scale_color_manual(values = lu_colors, name = NULL) +
  scale_fill_manual(values  = lu_colors, guide = "none") +
  labs(
    title    = "SOC depth profiles with exponential decay fits",
    subtitle = "Shaded bands = 95% parametric bootstrap CI",
    x        = "SOC (g kg\u207b\u00b9)",
    y        = "Depth (cm)"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title      = element_text(family = "serif", face = "bold",
                                   color = "#1e3a2f", size = 13),
    plot.subtitle   = element_text(color = "#7a7a72", size = 10),
    legend.position = "bottom",
    panel.grid.minor = element_blank()
  )

5. Calculating Soil Carbon Stocks

5.1 The stock calculation

Concentrations (g/kg) need to be converted to stocks (Mg C ha⁻¹) so we can compare how much carbon is actually stored per unit land area. The formula for each horizon is:

Stock (Mg C ha⁻¹) = SOC (g/kg) × BD (g/cm³) × thickness (cm) × 0.1

The factor 0.1 converts g/kg × g/cm³ × cm to Mg/ha.

soil_stocks <- soil %>%
  mutate(
    # Stock per horizon (Mg C ha⁻¹)
    stock_Mg_ha = soc_g_kg * bulk_density_g_cm3 * horizon_thick_cm * 0.1
  ) %>%
  group_by(land_use, site_id) %>%
  summarise(
    stock_0_30  = sum(stock_Mg_ha[depth_cm <= 30]),   # 0–30 cm
    stock_0_100 = sum(stock_Mg_ha),                    # 0–100 cm
    .groups = "drop"
  )

# Summary statistics by land use
soil_stocks %>%
  group_by(land_use) %>%
  summarise(
    n          = n(),
    mean_30    = round(mean(stock_0_30),  1),
    sd_30      = round(sd(stock_0_30),    1),
    mean_100   = round(mean(stock_0_100), 1),
    sd_100     = round(sd(stock_0_100),   1)
  ) %>%
  kbl(
    caption   = "SOC stocks (Mg C ha\u207b\u00b9) to 30 cm and 100 cm by land use",
    col.names = c("Land use", "n",
                  "Mean 0–30", "SD 0–30",
                  "Mean 0–100", "SD 0–100")
  ) %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE)

SOC stocks (Mg C ha⁻¹) to 30 cm and 100 cm by land use
Land use	n	Mean 0–30	SD 0–30	Mean 0–100	SD 0–100
Native chaparral	5	128.7	5.0	193.1	11.1
Restored native	5	100.0	4.2	149.6	9.2
Annual grassland	5	79.3	7.8	109.8	5.2
Cultivated	5	54.5	10.6	74.4	15.9

5.2 Visualise stocks

soil_stocks_long <- soil_stocks %>%
  pivot_longer(c(stock_0_30, stock_0_100),
               names_to  = "depth_range",
               values_to = "stock") %>%
  mutate(depth_range = recode(depth_range,
    "stock_0_30"  = "0–30 cm",
    "stock_0_100" = "0–100 cm"
  ))

ggplot(soil_stocks_long,
       aes(x = land_use, y = stock, fill = land_use)) +
  geom_boxplot(alpha = 0.75, outlier.shape = 21,
               outlier.size = 2.5, outlier.fill = "white") +
  geom_jitter(width = 0.12, size = 2.2, alpha = 0.65,
              aes(color = land_use)) +
  scale_fill_manual(values  = lu_colors, guide = "none") +
  scale_color_manual(values = lu_colors, guide = "none") +
  scale_x_discrete(guide = guide_axis(angle = 30)) +
  facet_wrap(~ depth_range, scales = "free_y") +
  labs(
    title    = "SOC stocks by land use and reference depth",
    x        = NULL,
    y        = "SOC stock (Mg C ha\u207b\u00b9)"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title   = element_text(family = "serif", face = "bold",
                                color = "#1e3a2f", size = 13),
    strip.text   = element_text(face = "bold"),
    panel.grid.major.x = element_blank()
  )

📌 0–30 cm vs 0–100 cm stocks. The 0–30 cm stock captures the biologically active topsoil where most carbon inputs and turnover occur. The 0–100 cm stock is the standard reference depth for national carbon inventories (IPCC 2006). The ratio of the two tells you how deeply carbon is stored: native vegetation has a lower ratio (more carbon at depth) while cultivated soils have a higher ratio (carbon concentrated near the surface by tillage).

6. Inferential Statistics: Do Land Uses Differ?

6.1 Check ANOVA assumptions

Before running ANOVA, verify that residuals are approximately normal and that group variances are similar.

# Focus on 0–100 cm stocks
lm_full <- lm(stock_0_100 ~ land_use, data = soil_stocks)

par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
qqnorm(residuals(lm_full), main = "Normal Q-Q — SOC stock residuals",
       col = "#4a7c5e", pch = 19)
qqline(residuals(lm_full), col = "#c8832a", lwd = 2)
plot(fitted(lm_full), residuals(lm_full),
     xlab = "Fitted values", ylab = "Residuals",
     main = "Residuals vs fitted",
     col = "#4a7c5e", pch = 19)
abline(h = 0, col = "#c8832a", lwd = 2, lty = 2)

par(mfrow = c(1, 1))

# Levene's test for homogeneity of variance
if (!requireNamespace("car", quietly = TRUE)) install.packages("car")
car::leveneTest(stock_0_100 ~ land_use, data = soil_stocks)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3   1.876 0.1743
##       16

💡 What to do if assumptions are violated. If residuals are strongly non-normal or variances are heterogeneous, two options: (1) log-transform SOC stocks before ANOVA — log stocks are often more normally distributed; (2) use the non-parametric Kruskal-Wallis test as a robust alternative. Both are shown below.

6.2 One-way ANOVA

aov_result <- aov(stock_0_100 ~ land_use, data = soil_stocks)
summary(aov_result)

##             Df Sum Sq Mean Sq F value  Pr(>F)    
## land_use     3  39261   13087   107.2 8.3e-11 ***
## Residuals   16   1953     122                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

6.3 Tukey HSD post-hoc test

A significant ANOVA F-test tells you that at least one pair of land uses differs. Tukey’s HSD identifies which pairs.

tukey_result <- TukeyHSD(aov_result)
print(tukey_result)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = stock_0_100 ~ land_use, data = soil_stocks)
## 
## $land_use
##                                         diff        lwr       upr     p adj
## Restored native-Native chaparral   -43.55669  -63.54673 -23.56664 0.0000642
## Annual grassland-Native chaparral  -83.35077 -103.34081 -63.36072 0.0000000
## Cultivated-Native chaparral       -118.68967 -138.67972 -98.69963 0.0000000
## Annual grassland-Restored native   -39.79408  -59.78412 -19.80404 0.0001761
## Cultivated-Restored native         -75.13299  -95.12303 -55.14295 0.0000001
## Cultivated-Annual grassland        -35.33891  -55.32895 -15.34887 0.0006080

# Tidy and select the columns we want for the table
tukey_tbl <- broom::tidy(tukey_result) %>%
  select(contrast, estimate, conf.low, conf.high, adj.p.value) %>%
  mutate(across(where(is.numeric), ~ round(.x, 3))) %>%
  arrange(adj.p.value)

tukey_tbl %>%
  kbl(
    caption   = "Tukey HSD pairwise comparisons — SOC stock 0–100 cm",
    col.names = c("Comparison", "Difference (Mg/ha)",
                  "Lower 95% CI", "Upper 95% CI", "Adj. p-value")
  ) %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) %>%
  row_spec(which(tukey_tbl$adj.p.value < 0.05),
           bold = TRUE, background = "#f1f8f1")

Tukey HSD pairwise comparisons — SOC stock 0–100 cm
Comparison	Difference (Mg/ha)	Lower 95% CI	Upper 95% CI	Adj. p-value
Restored native-Native chaparral	-43.557	-63.547	-23.567	0.000
Annual grassland-Native chaparral	-83.351	-103.341	-63.361	0.000
Cultivated-Native chaparral	-118.690	-138.680	-98.700	0.000
Annual grassland-Restored native	-39.794	-59.784	-19.804	0.000
Cultivated-Restored native	-75.133	-95.123	-55.143	0.000
Cultivated-Annual grassland	-35.339	-55.329	-15.349	0.001

6.4 Non-parametric alternative: Kruskal-Wallis

kruskal.test(stock_0_100 ~ land_use, data = soil_stocks)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  stock_0_100 by land_use
## Kruskal-Wallis chi-squared = 17.857, df = 3, p-value = 0.0004707

6.5 Effect size: η² (eta-squared)

A significant p-value does not tell you how large the difference is. η² is the proportion of total variance explained by the treatment.

ss_treatment <- sum(aov_result$effects[2:length(levels(soil_stocks$land_use))]^2)
ss_total     <- sum((soil_stocks$stock_0_100 - mean(soil_stocks$stock_0_100))^2)
eta_sq       <- ss_treatment / ss_total

cat(sprintf("eta-squared (\u03b7\u00b2) = %.3f\n", eta_sq))

## eta-squared (η²) = 0.953

cat(sprintf(
  "Interpretation: %.0f%% of variance in SOC stock is explained by land use.\n",
  eta_sq * 100
))

## Interpretation: 95% of variance in SOC stock is explained by land use.

7. Summary

💡 Key takeaways

Depth profiles always decrease with depth — flip the y-axis with scale_y_reverse() so depth increases downward, matching the physical reality.
Exponential decay (fit with nls()) is the standard SOC profile model. Parameter A gives surface SOC; k gives the decay rate. Larger k = steeper decline = less deep storage.
Carbon stocks require bulk density. Concentration × bulk density × thickness × 0.1 converts g/kg to Mg C ha⁻¹. Always report the reference depth (e.g., 0–30 cm or 0–100 cm).
Native and restored vegetation consistently stores more carbon than cultivated soils, particularly at depth where roots contribute to stable carbon pools.
Check ANOVA assumptions before testing. Q-Q plots, residual plots, and Levene’s test are a minimal checklist. Log-transform stocks or use Kruskal-Wallis when assumptions are violated.
Report effect size alongside p-values. η² tells you whether a statistically significant difference is also ecologically meaningful.

💬 Discussion questions

The exponential decay model fits well for mineral soils but often underestimates SOC in the top 5 cm where organic matter accumulates. What modification to the model structure would you make, and how would you test whether it improves fit?
Bulk density was not measured in many historical soil surveys, so it is often estimated from pedotransfer functions using texture and organic matter content. How would uncertainty in bulk density propagate into the stock estimates?
Restored native sites were planted 10 years ago. The model predicts they will eventually approach native chaparral SOC levels. Using the fitted k values from both land uses, estimate how many years it would take for the restored site to reach 90% of the native chaparral surface SOC (A parameter).
The Tukey test adjusts p-values for multiple comparisons. If you were testing all six pairwise combinations of four land uses, what is the familywise error rate without correction at α = 0.05, and why does this motivate the adjustment?
Carbon sequestration projects often report “additional carbon” relative to a baseline land use. Using the mean stocks from this analysis, calculate the additional carbon stored by native chaparral relative to cultivated land at 0–100 cm depth. Express your answer in Mg C ha⁻¹ and in tonnes CO₂-equivalent per hectare.

Soil Carbon Analysis

Marc Los Huertos

EA030 · Compiled 2026-05-26

Soil Carbon Analysis: Depth Profiles and Land-Use Comparisons

🎯 Learning Outcomes

The Analytical Pipeline

1. Setup

1.1 Packages

1.2 Simulated soil survey data

2. Import and Clean

2.1 Write and re-read as a CSV (simulating a real import)

2.2 Quality checks

2.3 Standardise land-use labels and factor ordering

3. Depth Profile Visualisation

3.1 Profile plot with depth on the vertical axis

3.2 Faceted profiles — one panel per land use

4. Fitting Depth-Profile Models

4.1 Exponential decay with `nls()`

4.2 Log-linear model with `lm()`

4.3 Profile plot with model fits and confidence intervals

5. Calculating Soil Carbon Stocks

5.1 The stock calculation

5.2 Visualise stocks

6. Inferential Statistics: Do Land Uses Differ?

6.1 Check ANOVA assumptions

6.2 One-way ANOVA

6.3 Tukey HSD post-hoc test

6.4 Non-parametric alternative: Kruskal-Wallis

6.5 Effect size: η² (eta-squared)

7. Summary

💡 Key takeaways

💬 Discussion questions

References

Soil Carbon Analysis

Marc Los Huertos

EA030 · Compiled 2026-05-26

Soil Carbon Analysis: Depth Profiles and Land-Use Comparisons

🎯 Learning Outcomes

The Analytical Pipeline

1. Setup

1.1 Packages

1.2 Simulated soil survey data

2. Import and Clean

2.1 Write and re-read as a CSV (simulating a real import)

2.2 Quality checks

2.3 Standardise land-use labels and factor ordering

3. Depth Profile Visualisation

3.1 Profile plot with depth on the vertical axis

3.2 Faceted profiles — one panel per land use

4. Fitting Depth-Profile Models

4.1 Exponential decay with nls()

4.2 Log-linear model with lm()

4.3 Profile plot with model fits and confidence intervals

5. Calculating Soil Carbon Stocks

5.1 The stock calculation

5.2 Visualise stocks

6. Inferential Statistics: Do Land Uses Differ?

6.1 Check ANOVA assumptions

6.2 One-way ANOVA

6.3 Tukey HSD post-hoc test

6.4 Non-parametric alternative: Kruskal-Wallis

6.5 Effect size: η² (eta-squared)

7. Summary

💡 Key takeaways

💬 Discussion questions

References

4.1 Exponential decay with `nls()`

4.2 Log-linear model with `lm()`