-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathREADME.Rmd
297 lines (230 loc) · 12.4 KB
/
README.Rmd
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
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
---
output: github_document
---
<!-- README.md is generated from README.Rmd. Please edit that file -->
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
message = FALSE,
eval = TRUE,
fig.path = "man/figures/README-",
fig.height = 5,
fig.width = 5,
dpi = 150,
out.width = "60%"
)
knitr::opts_knit$set(global.par = TRUE)
r = getOption("repos")
r["CRAN"] = "http://cran.us.r-project.org"
options(repos = r)
```
```{r, include = FALSE}
#> Set global graphic parameters. This code has to be in its own chunk.
par(mfrow = c(1, 1), mar = c(4.5, 4.5, 0.3, 0.3), bg = "white")
```
# isogeochem: <img src="man/figures/logo.png" align="right" width="140"/> <br/> Tools for Stable Isotope Geochemistry
<!-- badges: start -->
[](https://github.com/davidbajnai/isogeochem/actions/workflows/R-CMD-check.yaml)
[](https://app.codecov.io/gh/davidbajnai/isogeochem)
[](https://www.codefactor.io/repository/github/davidbajnai/isogeochem)
[](https://www.repostatus.org/#active)
[](https://CRAN.R-project.org/package=isogeochem)
[](https://cranlogs.r-pkg.org/badges/grand-total/isogeochem)
[](https://zenodo.org/badge/latestdoi/401782303)
<!-- badges: end -->
`isogeochem` makes working with stable oxygen, carbon, and clumped isotope data straightforward and reproducible. It offers tools to quickly calculate:
* carbonate *δ*<sup>18</sup>O, *δ*<sup>17</sup>O, *∆*<sub>47</sub>, and *∆*<sub>48</sub> values at a given temperature
* carbonate growth temperatures from *δ*<sup>18</sup>O, *∆*<sub>47</sub>, and *∆*<sub>48</sub> values
* isotope fractionation factors, e.g., carbonate/water, quartz/water
* model DIC speciation as a function of temperature, pH, and salinity
* convert between the VSMOW and VPDB scales
The list of available proxy–temperature calibrations is growing with each new released version. </br>
<b>Please get in touch if you have suggestions to include!</b>
## Getting started
### Installation
Install the released version of `isogeochem` from CRAN.
``` {r, include = FALSE}
if (!require("isogeochem")) install.packages("isogeochem")
```
``` {r, eval = FALSE}
install.packages("isogeochem")
```
Install the development version of `isogeochem` from GitHub.
``` {r, eval = FALSE}
if (!require("devtools")) install.packages("devtools")
if (!require("rmarkdown")) install.packages("rmarkdown")
devtools::install_github("davidbajnai/isogeochem", build_vignettes = TRUE)
```
``` {r, include = FALSE}
library("isogeochem")
```
### Vignettes
Case studies demonstrating the use and scope of the functions in `isogeochem` are
available as vignettes.
``` {r, eval = FALSE}
browseVignettes("isogeochem")
```
## Dual clumped isotope thermometry
Use `D47()` and `D48()` to calculate equilibrium carbonate clumped isotope values (*∆*<sub>47</sub>, *∆*<sub>48</sub>) for a given temperature. `temp_D47()` calculates carbonate growth temperatures from *∆*<sub>47</sub> values, while `temp_D48()` calculates growth temperature corrected for kinetic effects considering both the *∆*<sub>47</sub> and the *∆*<sub>48</sub> value.
``` {r, label = "Figure1"}
if (!require("shades")) install.packages("shades")
# Model equilibrium carbonate ∆47 and ∆48 values
temp = seq(0, 100, 10) # temperature range: 0—100 °C
D47eq = D47(temp, eq = "Fiebig21")
D48eq = D48(temp, eq = "Fiebig21")
# Sample data
D47_coral = 0.617; D47_coral_err = 0.006
D48_coral = 0.139; D48_coral_err = 0.022
D47_speleo = 0.546; D47_speleo_err = 0.007
D48_speleo = 0.277; D48_speleo_err = 0.029
## Plot in ∆47 vs ∆48 space ##
plot(0, type = "l", axes = TRUE, ylim = c(0.4, 0.7), xlim = c(0.1, 0.3),
ylab = expression(Delta[47] * " (CDES90, ‰)"),
xlab = expression(Delta[48] * " (CDES90, ‰)"),
lty = 0, font = 1, cex.lab = 1, las = 1)
# Plot the equilibrium curve and points
lines (D48eq, D47eq, col = "purple", lwd = 2)
points(D48eq, D47eq, col = shades::gradient(c("blue", "red"), length(temp)),
pch = 19, cex = 1.2)
# Plot the sample data,
# ... the kinetic slopes,
# ... and calculate growth temperatures corrected for kinetic effects
# ... using a single function!
temp_D48(D47_coral, D48_coral, D47_coral_err, D48_coral_err, ks = -0.6,
add = TRUE, col = "seagreen", pch = 15)
temp_D48(D47_speleo, D48_speleo, D47_speleo_err, D48_speleo_err, ks = -1,
add = TRUE, col = "darkorange", pch = 17)
# Add labels to the plot
text(D48(temp, eq = "Fiebig21"), D47(temp, eq = "Fiebig21"), paste(temp, "°C"),
col = shades::gradient(c("blue", "red"), length(temp)), pos = 4, cex = 0.8)
```
## Triple oxygen isotopes
`d17O_c()` calculates equilibrium carbonate oxygen isotope values (*δ*<sup>18</sup>O, *δ*<sup>17</sup>O, *∆*<sup>17</sup>O) for a given temperature and ambient water composition. Use the `mix_d17O()` function to calculate mixing curves in triple oxygen isotope space, e.g., for modeling diagenesis.
``` {r, label = "Figure2"}
if (!require("shades")) install.packages("shades")
# Model equilibrium *calcite* precipitating from seawater
temp_sw = seq(0, 50, 10) # temperature range: 0—50 °C
d18O_sw = 0 # d18O value of seawater
D17O_sw = -0.004 # D17O value of seawater
d18Op = prime(d17O_c(temp_sw, d18O_sw, D17O_sw, eq18 = "Daeron19")[, 1])
D17O = d17O_c(temp_sw, d18O_sw, D17O_sw, eq17 = "Wostbrock20", eq18 = "Daeron19")[, 3]
# Model progressing meteoric diagenetic alteration
d18O_ds = -8 # d18O value of diagenetic fluid
D17O_ds = 0.020 # D17O value of diagenetic fluid
em_equi = d17O_c(temp = 10, d18O_H2O = d18O_sw, D17O_H2O = D17O_sw,
eq17 = "Wostbrock20", eq18 = "Daeron19") # equilibrium endmember
em_diag = d17O_c(temp = 80, d18O_H2O = d18O_ds, D17O_H2O = D17O_ds,
eq17 = "Wostbrock20", eq18 = "Daeron19") # diagenetic endmember
mix = mix_d17O(d18O_A = em_equi[1], d17O_A = em_equi[2],
d18O_B = em_diag[1], d17O_B = em_diag[2], step = 25)
## Plot in ∆17O vs d'18O space ##
plot(0, type = "l", ylim = c(-0.1, 0.05), xlim = c(-10, 40),
xlab = expression(delta * "'" ^ 18 * "O (‰, VSMOW)"),
ylab = expression(Delta ^ 17 * "O (‰, VSMOW)"),
lty = 0, font = 1, cex.lab = 1, las = 1)
# Plot meteoric waters from the build-in dataset
points(prime(meteoric_water$d18O), D17O(meteoric_water$d18O, meteoric_water$d17O),
col = "lightblue1", pch = 20)
text(-4, 0.05, "meteoric water", pos = 4, col = "lightblue1")
# Plot the composition of the fluids
points(prime(d18O_sw), D17O_sw, col = "darkmagenta", pch = 8) # seawater
text(prime(d18O_sw), D17O_sw, "seawater", pos = 4, col = "darkmagenta")
points(prime(d18O_ds), D17O_ds, col = "deeppink", pch = 8) # diagenetic fluid
text(prime(d18O_ds), D17O_ds, "diagenetic fluid", pos = 4, col = "deeppink")
# Plot the equilibrium curve and points
lines(d18Op, D17O, col = "darkmagenta", lwd = 2)
points(d18Op, D17O, pch = 19, cex = 1.2,
col = shades::gradient(c("blue", "red"), length(temp_sw)))
text(d18Op, D17O, paste(temp_sw, "°C"), pos = 4, cex = 0.8,
col = shades::gradient(c("blue", "red"), length(temp_sw)))
text(30, -0.05, paste("equilibrium calcite \nfrom seawater"),
pos = 3, col = "darkmagenta")
# Plot the mixing model between the equilibrium and diagenetic endmembers
lines(prime(mix[, 1]), mix[, 2], col = "deeppink", lty = 3, lwd = 2)
points(prime(mix[, 1]), mix[, 2], pch = 18, cex = 1.5,
col = shades::gradient(c("#3300CC", "deeppink"), length(mix[, 2])))
text(prime(mix[, 1]), mix[, 2], paste(mix[, 3], "%", sep = ""), pos = 2, cex = 0.8,
col = shades::gradient(c("#3300CC", "deeppink"), length(mix[, 3])))
text(22, -0.09, paste("progressing", "\ndiagenetic alteration", "\n(recrystallisation) at 80°C", sep =""),
pos = 2, col = "deeppink")
```
## Thermometry
Use `isogeochem` to calculate crystallization temperatures from carbonate *δ*<sup>18</sup>O and *∆*<sub>47</sub> values.
``` {r}
# Temperature from D47 with or without errors
temp_D47(D47_CDES90 = 0.601, eq = "Petersen19")
temp_D47(D47_CDES90 = 0.601,
D47_error = 0.008 ,
eq = "Anderson21")
# Temperature from d18O
temp_d18O(
d18O_c_VSMOW = 30,
d18O_H2O_VSMOW = 0,
min = "calcite",
eq = "Watkins13")
```
## Fractionation factors
Use `isogeochem` to calculate <sup>16</sup>O/<sup>18</sup>O fractionation factors at given temperatures.
``` {r, label = "Figure3"}
if (!require("viridisLite")) install.packages("viridisLite")
plot(0, type = "l", las = 1, yaxt = "n",
xlim = c(10, 30), ylim = c(-30, 50),
xlab = "Temperature (°C)",
ylab = expression("Equilibrium enrichment in "^18*"O relative to H"[2]*"O (‰)"))
axis(2, seq(-30, 50, 10), las = 1)
temps = seq(10, 30, 1)
d18O_H2O_VSMOW = 0
cols = viridisLite::viridis(7, option = "C")
text(10, 45, expression("CO"[2]*" (aq)"), col = cols[1], adj = c(0, 0))
lines(temps, A_from_a(a18_CO2aq_H2O(temps), d18O_H2O_VSMOW),
lwd = 2, lty = 2, col = cols[1])
text(10, 35, expression("HCO"[3]^"–"), col = cols[2], adj = c(0, 0))
lines(temps, A_from_a(a18_HCO3_H2O(temps), d18O_H2O_VSMOW),
lwd = 2, lty = 2, col = cols[2])
text(10, 30, "calcite", col = cols[3], adj = c(0, 0))
lines(temps, A_from_a(a18_c_H2O(temps, "calcite", "Daeron19"), d18O_H2O_VSMOW),
lwd = 2, lty = 1, col = cols[3])
text(10, 21, expression("CO"[3]^"2–"), col = cols[4], adj = c(0, 0))
lines(temps, A_from_a(a18_CO3_H2O(temps), d18O_H2O_VSMOW),
lwd = 2, lty = 2, col = cols[4])
text(10, 1, expression("H"[2]*"O"), col = cols[5], adj = c(0, 0))
lines(temps, rep(d18O_H2O_VSMOW, length(temps)),
lwd = 3, lty = 1, col = cols[5])
text(10, -23, expression("OH"^"–"), col = cols[6], adj = c(0, 0))
lines(temps, B_from_a(a18_H2O_OH(temps, eq = "Z20-X3LYP"), d18O_H2O_VSMOW),
lwd = 2, lty = 1, col = cols[6])
```
## Utility functions
``` {r}
# Convert between the VSMOW and VPDB scales:
to_VPDB(32)
to_VSMOW(1)
# Convert between classical delta and delta prime values:
prime(10)
unprime(9.95)
# Calculate isotope fractionation factors:
a_A_B(A = 30.40, B = 0.15)
epsilon(a_A_B(A = 30.40, B = 0.15))
```
## Datasets
Within `isogeochem` you have quick access to important datasets.
| Name | Description | Reference |
|-------------|---------------------------------------------------------------- |------------------------|
| `devilshole`| The original Devils Hole carbonate *δ*<sup>18</sup>O time series| Winograd et al. (2006) |
| `LR04` | A benthic foraminifera *δ*<sup>18</sup>O stack | Lisiecki & Raymo (2005)|
| `GTS2020` | An abridged version of the GTS2020 oxygen isotope stack | Grossman & Joachimski (2020)|
| `meteoric_water`| A compilation of meteoric water *δ*<sup>18</sup>O and *δ*<sup>17</sup>O values | Barkan & Luz (2010), Aron et al. (2021)|
For more information on the datasets please have a look at the corresponding documentation, e.g., `?devilshole`
___
## License and citation
Copyright (C) 2023 David Bajnai
This program is free software: you can redistribute it and/or modify it under the terms of the [GNU General Public License version 3](LICENSE.md), or (at your option) any later version. This program is distributed in the hope that it will be useful, but without any warranty.
Follow the citation format provided in [CITATION](CITATION.cff) when referencing `isogeochem`.
___
## See also
There are several other R packages that complement `isogeochem` and are worth checking out:
[`viridisLite`](https://github.com/sjmgarnier/viridisLite) and [`viridis`](https://github.com/sjmgarnier/viridis) produce color-blind and black-and-white printer friendly color scales.
[`clumpedr`](https://github.com/isoverse/clumpedr/) works with [`isoreader`](https://github.com/isoverse/isoreader) to read in raw measurement data and reproducibly process the results to clumped isotope values.
[`seasonalclumped`](https://github.com/nielsjdewinter/seasonalclumped) can be used to reconstruct temperature and salinity variations from seasonal oxygen and clumped isotope records.
[`deeptime`](https://github.com/willgearty/deeptime) adds geological timescales to ggplots.