有機質資材からの窒素放出量を予測したい

有機質資材からの窒素放出

ほ場に有機質資材、例えば牛ふん堆肥であったり稲わらであったりというものを施用した場合、これらに含まれている窒素というのは必ずしもすぐには使いやすい形にならない。

有機質資材中の窒素は基本的には微生物の分解を受けることで植物が利用しやすい形態となる。有機質資材はそれぞれ特有の割合で炭素と窒素を含んでいて、これをC/N比と言う。そして、有機物を分解する微生物も炭素と窒素を体に含んでいて、微生物のC/N比はだいたい20くらいである。従って、C/N比が20より小さいときは、微生物が有機物を体内に取り込むと、窒素のほうが多いので窒素が土壌中に放出される。アバウトな説明だがこれが有機物の肥料効果が現れる流れである。従って、C/N比が低い有機物ほど肥料効果が高く、早い。

では20より大きいと炭素が放出されるかというとそうではなくて、微生物は土壌中にもともとあった窒素を利用して増殖をはじめる。従って土壌中の窒素分が吸われてしまい、負の肥料効果とでも言うべきものが現れる。C/N比が大きな有機物を投入することで土壌中の窒素がなくなってしまう現象を窒素飢餓と呼ぶ。

しかし、微生物に取り込まれた窒素も、微生物が死ねば土壌中に放出される。炭素は呼吸に伴って炭酸ガスとして減少していくので、基本的には有機物中の窒素はいつかは全て放出される(ガス化する場合もあるが)。

したがって、毎年有機物を入れ続ければ、土壌中に放出される窒素の量は年々増えていき、放出量は投入量に漸近していく。しかし、平衡状態に達するには通常長い年月がかかる。何年か連用したときの窒素放出量を知るすべはないだろうか。

土壌中における有機物の分解モデル

土壌中の有機物の分解モデルはいくつか考えられているものがあって、有機物が半減期の異なるいくつかの種類のコンパートメントの複合したものであるという前提でモデル化したもので比較的良好な結果がえられているものが多い。

前置きはこの程度にしておいて

今回はhttp://agriknowledge.affrc.go.jp/RN/2010322673のモデルを使う*1。

Rでやっていこう。

まず、paramというデータテーブルが次のように用意してあるものとする。

資材名	a	c	f
乾燥牛ふん	0.10	0.58	0.32
牛ふん堆肥	0.04	0.15	0.81
牛ふんオガクズ堆肥	0.39	0.11	0.85
豚ぷん堆肥	0.02	0.39	0.59
豚ぷんオガクズ堆肥	0.02	0.34	0.64
バーク堆肥	0.14	-0.08	0.96
未熟稲ワラ堆肥	-0.06	0.26	0.80
中熟稲ワラ堆肥	0.03	0.18	0.79
完熟稲ワラ堆肥	0.15	0.11	0.74
稲ワラ	-0.40	0.95	0.45
小麦ワラ	-1.83	1.63	1.20
モミガラ	0.19	-0.11	0.92
水稲根	-0.35	0.63	0.72
オガクズ	-0.03	-2.77	3.80
製紙カス	-1.15	0.49	1.66
発酵製紙カス	0.43	-0.37	0.94
余剰汚泥	0.65	0.28	0.07

パッケージ

こんな感じで。data.tableは気分。

library(dplyr)
library(data.table)
library(ggplot2)
library(ggrepel)

残存率、蓄積率、放出率を求める

毎年の有機物投入量を1として、t年後の残存率、累積率、放出率は、

## functions ---------------
residual_rate <- function(a, c, f, t){ # 残存率
  a * 0.01^t + c * 0.63^t + f * 0.955^t
}

accum_rate <- function(a, c, f, t){ # 累積率
  a * 0.01 * (1-0.01^t)/(1-0.01) + 
    c * 0.63 * (1-0.63^t)/(1-0.63) +
    f * 0.955 * (1-0.955^t)/(1-0.955)
}

release_rate <- function(a, c, f, t){ #放出率
  1 - (a * 0.01^t + c * 0.63^t + f * 0.955^t)
}

そしてtをベクトルで与えられるように関数を整備して、プロットに備える。

sim_fun_gen <- function(t, name, fun){
  tmp_p <- param[資材名==name]
  a <- tmp_p$a
  c <- tmp_p$c
  f <- tmp_p$f
  sapply(t, function(t) fun(a, c, f, t))
}

resid_rate_sim <- function(t, name){
  sim_fun_gen(t, name, residual_rate)
}
accum_rate_sim <- function(t, name){
  sim_fun_gen(t, name, accum_rate)
}
rel_rate_sim <- function(t, name){
  sim_fun_gen(t, name, release_rate)
}

窒素放出量を予測する

実際に一番知りたいのは窒素の放出量だ。連用何年目に何kg窒素が出てくるのかということだ。

放出量は、放出率に実際に資材が持っている窒素の量を乗ずれば出てくる。

t = 0:50
shizai <- c("豚ぷん堆肥", "豚ぷんオガクズ堆肥", "牛ふん堆肥", "牛ふんオガクズ堆肥",
            "稲ワラ", "中熟稲ワラ堆肥", "小麦ワラ", "オガクズ")
n_rate <- c(2.7, 1.4, 1.1, 0.8, 0.6, 0.5, 0.4, 0.1)
names(n_rate) <- shizai

result <- data.table(資材名=character(0), 連用年数=numeric(0), 窒素放出量=numeric(0))

for(i in 1:length(n_rate)){
  result <-
    rbind(result,
          data.table(資材名 = names(n_rate)[i],
                     連用年数 = t,
                     窒素放出量 = rel_rate_sim(t, names(n_rate)[i])* n_rate[i] * 10)
          )
}

プロットする

result %>%
  ggplot(aes(x = 連用年数, y = 窒素放出量, color = 資材名)) +
  geom_line() + 
  theme_bw(base_family = "Osaka") +
  geom_text_repel(
    data = subset(result, 連用年数==max(連用年数)),
    aes(label = 資材名),
    family = "Osaka",
    nudge_x = 10,
    size = 3,
    segment.size = 0.25,
    segment.alpha = 0.5
  ) +
  lims(x = c(min(result$連用年数), max(result$連用年数)*1.3)) + 
  theme(legend.position = "none") + 
  geom_hline(aes(yintercept = 0), lty = 2, cex = 0.1) +
  labs(title = "有機質資材を1t/10aずつ連用した水田での無機体窒素放出経過",
       x = "連用年数(年)",
       y = "無機体窒素放出量(kg/10a)")