-
Notifications
You must be signed in to change notification settings - Fork 0
/
N2O_SoilEmissions.qmd
167 lines (129 loc) · 6.29 KB
/
N2O_SoilEmissions.qmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
---
title: "Soil N2O Emissions-Regression Model"
---
This project contains two parts. The firs part estimates carbon stock change factors (i.e., emission factors) for tillage management and a developed Tier 3 regression model for soil N~2~O emissions through an
analysis of measurement data. I derived Tier 2 emission factors using simple averages and linear regression models.
The second part of this project fits the Tier 3 regression model from part 1 based on measured N2O emissions. An R script function was created to estimate emissions from soil carbon stock changes.
###
Part 1: Emission Factor Development R Script
```{r eval=FALSE}
library(nlme)
###### Management Factor
# Read data from csv file
management.data<-read.csv("SoilCManagement.csv", header=TRUE)
# test for correlation in predictor variables
cor(management.data[,c("years", "dep1", "dep2")])
# fit full model with all variables as main effects for input
test.fit<-lme(ch.cstock~ch.till+years+years2+dep1+dep2+moisture+temp, random=~1 | ran.exp/ran.yrexp,
data=management.data, method="ML", na.action=na.omit)
summary(test.fit)
# Diagnostic Plots, Residual Plot
resid<-residuals(test.fit)
plot(fitted(test.fit), resid)
abline(0,0)
# QQ normal plot
qqnorm(resid)
qqline(resid)
# Full model
test.fit<-lme(ch.cstock~ch.till+years+years2+dep1+dep2+moisture+temp, random=~1 | ran.exp/ran.yrexp,
data=management.data, method="ML", na.action=na.omit)
# Backward stepwise fit method
test.fit<-lme(ch.cstock~ch.till+years+dep1+dep2+moisture+temp, random=~1 | ran.exp/ran.yrexp,
data=management.data, method="ML", na.action=na.omit)
test.fit<-lme(ch.cstock~ch.till+years+dep1+dep2+moisture, random=~1 | ran.exp/ran.yrexp,
data=management.data, method="ML", na.action=na.omit)
summary(test.fit)
# Test interactions
test.fit<-lme(ch.cstock~ch.till+years+dep1+dep2+moisture+ ch.till*years+ch.till*dep1+
ch.till*dep2+ch.till*moisture+years*dep1+years*dep2+years*moisture+
dep1*moisture+dep2*moisture, random=~1 | ran.exp/ran.yrexp,
data=management.data, method="ML", na.action=na.omit)
test.fit<-lme(ch.cstock~ch.till+years+dep1+dep2+moisture+ ch.till*years+ch.till*dep1+
ch.till*dep2+ch.till*moisture+years*dep1+years*dep2+
dep1*moisture+dep2*moisture, random=~1 | ran.exp/ran.yrexp,
data=management.data, method="ML", na.action=na.omit)
test.fit.management<-lme(ch.cstock~ch.till+years+dep1+dep2+moisture+ ch.till*years+ch.till*dep1+
ch.till*dep2+years*dep1+years*dep2+
dep1*moisture+dep2*moisture, random=~1 | ran.exp/ran.yrexp,
data=management.data, method="REML", na.action=na.omit)
summary(test.fit.management)
# Derive PDF
fixed.management<-fixed.effects(test.fit.management)
management.cov<-test.fit.management$varFix
# Variables
x.rt.wet<-c(1,1,20,15,300,1,20,15,300,300,6000,15,300)
x.rt.dry<-c(1,1,20,15,300,0,20,15,300,300,6000,0,0)
x.nt.wet<-c(1,0,20,15,300,1,0,0,0,300,6000,15,300)
x.nt.dry<-c(1,0,20,15,300,0,0,0,0,300,6000,0,0)
# Estimates
t(x.rt.wet)%*%fixed.management
t(x.rt.dry)%*%fixed.management
t(x.nt.wet)%*%fixed.management
t(x.nt.dry)%*%fixed.management
# Variance
v.rt.wet<-(t(x.rt.wet)%*%management.cov%*%x.rt.wet)
v.rt.dry<-(t(x.rt.dry)%*%management.cov%*%x.rt.dry)
v.nt.wet<-(t(x.nt.wet)%*%management.cov%*%x.nt.wet)
v.nt.dry<-(t(x.nt.dry)%*%management.cov%*%x.nt.dry)
# Standard Deviation
sqrt(v.rt.wet)
sqrt(v.rt.dry)
sqrt(v.nt.wet)
sqrt(v.nt.dry)
```
### Part 2: Soil N2O Emissions Estimation R Script
```{r eval=FALSE}
"SynFert.N2O.emissions.Regression"<-
function(mineralN.amount=75, mineralN.amount.sd=5, beta= a, cov.beta=b, MAPPET=1, nreps=10000, iseed=230984, return.option=1)
# Script developed by Natalie Wiley
# Originally developed: 3/2/2022
# Last updated:
# Script estmates N2O emissions (kg CO2 eq. per ha per yr) from mineral N fertilization
### Arguments
# mineralN.amount The amount of mineral N fertilizer added to soil (Kg N per HA)
# mineralN.amount.sd Standard deviation of the N mineral fertilizer amount
# MAPPET Mean annual precipitation to potential evapotranspiration ratio
# beta R object with betas from LME model
# cov.beta R object with covriance matrix for betas from LME model
# nreps Number of monte carlo simulations
# iseed Initial seed for random draws
# return.option 1) list object with the emission mean and confidence intervals for N2O emissions, and
# 2) the full vector of all Monte Carlo simulations
##
##### Begin Script
{
##### Set seed
set.seed(iseed)
##### Check validity of input variables
# values equal to or greater than 0 are valid
check.mineralN.amount<-mineralN.amount>=0&mineralN.amount<=880
if(!check.mineralN.amount){stop("MineralN amount is not valid.")
} else {cat("NOTE: Mineral N amount is valid")}
check.mineralN.amount.sd<-mineralN.amount.sd>=0
if(!check.mineralN.amount.sd){stop("MineralN sd is not valid.")
} else {cat("NOTE: Mineral N sd is valid")}
check.MAPPET<-MAPPET>=0.7&MAPPET<=3.3
if(!check.MAPPET){stop("MAP:PET ratio is not valid.")
} else {cat("NOTE: MAP:PET is valid")}
# End Validity Checks
##
##### Estimate direct N2O emissions using linear mixed effect model
# Deterministic Calculation
# estimate emissions and backtransform (since we had to transform our data to be homogenous and use a log transformation, we now have to backtransform a log # into regular units)
direct.emission.deterministic.ln<-beta%*%t(cbind(1, mineralN.amount,MAPPET))
direct.emission.deterministic<- (exp(direct.emission.deterministic.ln)) * (44/28) * 298
# Probabilistic calculation
# Simulate nreps of fertilizer amounts
mineralN.amount.sim <- rnorm(nreps, mean = mineralN.amount, sd = mineralN.amount.sd)
# Simulate nreps of beta parameters based on LME model
# determine number of parameters
numpar<-length(beta) #number of parameters in the model
# computer choleski decomposition
M<-t(chol(cov.beta))
#
# generate random normals
z<-matrix(rnorm(nreps*numpar), numpar, nreps)
# produce simulated betas
sim.beta<-M%*%z+beta
}
```