-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathMagOxyDriver.R
148 lines (107 loc) · 4.83 KB
/
MagOxyDriver.R
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
# Driver for time series data inversion using forward MagOxyPSM
#########################################################################################
# Load libraries
library(rjags)
library(R2jags)
############################################################################################
############################################################################################
# INPUT TO PASS TO JAGS
############################################################################################
# Input approximate % of calcite non-primary (recrystallized), with 1sd uncertainty,
# and estimated Dd18O of primary-inorganic (secondary) calcite
seccal = 20 # Percent recrystallized
seccal.u = 2.5 # 1sd % recrystallized
d18Oseccal = 0.85 # Calculated following Edgar et al. (2015)
############################################################################################
# Mg/Ca SST proxy vital effects and calibration parameters
# Nonlinearity of the relationship b/w shell and Mg/Casw (Haynes et al. 2023, T. sacculifer)
Hp.mean = 0.74
Hp.sd = 0.05
# Modern (pre-corrected) pre-exponential constant in Mg/Ca-SST calibration (Anand et al., 2003; Evans et al., 2016)
Bmod.mean = 0.38
Bmod.sd = 0.02
# Exponential constant in Mg/Ca-SST calibration (Evans et al., 2016) for low (P-E) Mg/Casw
A.mean = 0.061
A.sd = 0.005
############################################################################################
# INVERSION DRIVER
############################################################################################
# These parameters will be recorded in the output
parms <- c("sal", "tempC", "xca", "xmg", "d18Osw", "d18Osw.sc", "d18Of", "mgcaf")
# Read in proxy time series data
prox.in <- read.csv('1258PETM/data/1258PETM.csv')
prox.in <- prox.in[,c(1:5)]
names(prox.in) <- c("depth", "age", "d13C", "d18O", "MgCa")
# Setup age range and bins
ages.prox <- unique(round(prox.in$age))
ages.prox.max <- max(ages.prox)
dt <- abs(diff(ages.prox, lag=1))
ages.prox.ai <- seq(1,length(ages.prox), by=1)
# Age index proxy data
prox.in <- transform(prox.in,ai=as.numeric(factor(round(age*-1))))
clean.mgca <- prox.in[complete.cases(prox.in$MgCa), ]
mgcafu <- clean.mgca$MgCa*0.015
clean.d18O <- prox.in[complete.cases(prox.in$d18O), ]
# Vector of age indexes that contain Mg/Ca proxy data
ai.mgca <- c(clean.mgca$ai)
# Vector of age indexes that contain d18O proxy data
ai.d18O <- c(clean.d18O$ai)
ai.all <- c(ai.mgca, ai.d18O)
# Index vector which contains each environmental time step that has one or more proxy data
ai.prox <- unique(ai.all)
ai.prox <- sort(ai.prox, decreasing = FALSE)
n.steps <- length(ai.prox)
# Prior time bin vectors for which there are proxy data (includes duplicates)
ai.mgca <- match(ai.mgca, ai.prox)
ai.d18O <- match(ai.d18O, ai.prox)
############################################################################################
# Input prior mean and precision estimates for environmental parameters
sal.m = 35
sal.p = 1/0.5^2
tempC.m = 33
tempC.p = 1/2^2
xca.m = 21
xca.p = 1/1^2
xmg.m = 52
xmg.p = 1/1^2
d18Osw.m = -0.4
d18Osw.p = 1/0.1^2
############################################################################################
# Data to pass to jags
data <- list("d18Of.data" = clean.d18O$d18O,
"mgcaf.data" = clean.mgca$MgCa,
"mgcafu.data" = mgcafu,
"n.steps" = n.steps,
"dt" = dt,
"ages.prox" = ages.prox,
"ai.prox" = ai.prox,
"ai.d18O" = ai.d18O,
"ai.mgca" = ai.mgca,
"seccal" = seccal,
"seccal.u" = seccal.u,
"d18Oseccal" = d18Oseccal,
"Hp.mean" = Hp.mean,
"Hp.sd" = Hp.sd,
"Bmod.mean" = Bmod.mean,
"Bmod.sd" = Bmod.sd,
"A.mean" = A.mean,
"A.sd" = A.sd,
"sal.m" = sal.m,
"sal.p" = sal.p,
"tempC.m" = tempC.m,
"tempC.p" = tempC.p,
"xca.m" = xca.m,
"xca.p" = xca.p,
"xmg.m" = xmg.m,
"xmg.p" = xmg.p,
"d18Osw.m" = d18Osw.m,
"d18Osw.p" = d18Osw.p)
############################################################################################
# Run the inversion
jout = jags.parallel(model.file = "1258PETM/MagOxyPSM.R", parameters.to.save = parms,
data = data, inits = NULL, n.chains = 9, n.iter = 150000,
n.burnin = 50000, n.thin = 100)
############################################################################################
write.csv(jout$BUGSoutput$summary, "1258PETM/inversion.sum.csv")
PETM.DSST <- rowMeans(jout$BUGSoutput$sims.list$tempC[,10:14]) - rowMeans(jout$BUGSoutput$sims.list$tempC[,6:7])
print(c(mean(PETM.DSST), quantile(PETM.DSST, 0.025), quantile(PETM.DSST, 0.975)))