terrestrialProx.qmd

---
title: "Terrestrial Proxies"
author: "Gabriel Bowen"
date: "`r format(Sys.time(), '%Y-%m-%d')`"
editor: source
format: pdf
---

```{r setup, include=FALSE}
source("code/prepData.R")
load("bigout/post.rda")

pm = post$BUGSoutput$median

col.c = rgb(1, 0.2, 0.6)
col.a = rgb(0.2, 0.6, 1)

knitr::opts_knit$set(global.par = TRUE)
```

```{r setpar, include=FALSE}
par(mar = c(5, 5, 4, 1))
```

Let's review the high-latitude (terrestrial) proxy data and its interpretation. I'll refer to the sites by number for convenience, here's the key:

```{r sites, echo=FALSE}
print("Carnian = pink")
cont_sites.car

print("Anisian = blue")
cont_sites.an
```


## Al/Si

Below I'm plotting the calibration data set for the proxy (black) and the paleo-proxy values vs. the posterior median site MAT (pink = Carnian, blue = Anisian).

```{r alsi, echo=FALSE,fig.height=4}
plot(AlSi ~ MAT, data = alsi_cal, ylab = "Al/Si",
     xlim = range(c(as.data.frame(alsi_cal)$MAT, pm$t_cont.car[alsi_sites.ind.car],
                    pm$t_cont.an[alsi_sites.ind.an])),
     ylim = range(c(as.data.frame(alsi_cal)$AlSi, alsi_samples.car$AlSi,
                    alsi_samples.an$AlSi)))

points(pm$t_cont.car[alsi_sites.ind.car], alsi_samples.car$AlSi,
       pch = 20, col = col.c)

points(pm$t_cont.an[alsi_sites.ind.an], alsi_samples.an$AlSi,
       pch = 20, col = col.a)

text(pm$t_cont.car[alsi_sites.ind.car], par("usr")[4], alsi_sites.ind.car, 
     xpd = NA, pos = 3, col = col.c)

text(pm$t_cont.an[alsi_sites.ind.an], par("usr")[4], alsi_sites.ind.an, 
     xpd = NA, pos = 3, offset = 1.25, col = col.a)
```

* The data are very noisy.
* Al/Si values for all the Carnian sites are similar. Because of the influence of other proxies some have posterior MATs that fall below and others above the values you would estimate from Al/Si alone. 
* The one Anisian site has a MAT estimate ~3 $^\circ$C hotter than what would be estimated from MAT alone.

## CIA

Same thing for CIA.

```{r CIA, echo=FALSE,fig.height=4}
plot(CIA ~ Temp, data = cia_cal, xlab = "MAT",
     xlim = range(c(as.data.frame(cia_cal)$Temp, pm$t_cont.car[cia_sites.ind.car],
                    pm$t_cont.an[cia_sites.ind.an])),
     ylim = range(c(as.data.frame(cia_cal)$CIA, cia_samples.car$CIA,
                    cia_samples.an$CIA)))

points(pm$t_cont.car[cia_sites.ind.car], cia_samples.car$CIA,
       col = col.c, pch = 20)
points(pm$t_cont.an[cia_sites.ind.an], cia_samples.an$CIA,
       col = col.a, pch = 20)

text(pm$t_cont.car[cia_sites.ind.car], par("usr")[4], cia_sites.ind.car, 
     xpd = NA, pos = 3, col = col.c)
text(pm$t_cont.an[cia_sites.ind.an], par("usr")[4], cia_sites.ind.an, 
     xpd = NA, pos = 3, offset = 1.25, col = col.a)
```

* Posterior MATs are pretty consistent with the proxy values, within the range of proxy uncertainty.
* On average, posterior MATs are slightly higher at all Carnian sites than would be predicted from CIA values alone...this proxy is pulling temperatures down. 

## PWI

Now paleosol weathering index, available at one Anisian and two Carnian sites. I've added the calibration curve using the median curve fit parameters in the posterior for reference.

```{r PWI, echo=FALSE,fig.height=4}
plot(PWI ~ MAT, data = pwi_cal, 
     xlim = range(c(as.data.frame(pwi_cal)$MAT, pm$t_cont.car[pwi_sites.ind.car],
                    pm$t_cont.an[pwi_sites.ind.an])),
     ylim = range(c(as.data.frame(pwi_cal)$PWI, pwi_samples.car$PWI, 
                    pwi_samples.an$PWI)))
xs = seq(par("usr")[1], par("usr")[2], length = 100)
ys = exp((xs + c(pm$pwi_a)) / c(pm$pwi_b))
lines(xs, ys)

points(pm$t_cont.car[pwi_sites.ind.car], pwi_samples.car$PWI,
       col = col.c, pch = 20)
points(pm$t_cont.an[pwi_sites.ind.an], pwi_samples.an$PWI,
       col = col.a, pch = 20)

text(pm$t_cont.car[pwi_sites.ind.car], par("usr")[4], pwi_sites.ind.car, 
     xpd = NA, pos = 3, col = col.c)
text(pm$t_cont.an[pwi_sites.ind.an], par("usr")[4], pwi_sites.ind.an, 
     xpd = NA, pos = 3, offset = 1.25, col = col.a)
```

Here the reconstructions fit the calibration data well. Obviously the interpretation is based on an extrapolation at Carnian site 3, but if one takes the proxy at face value it's pretty tough for this site to have a MAT higher than, say, 5 or 10 $^\circ$. Note that this same site has a cool bias in the first plot (based on Al/Si).

## Clumped

Here I'm using the calibration function from Julia Kelson's paper, which is represented as a line. Only one Anisian site.

```{r clumped, echo=FALSE,fig.height=4}
plot(20:45, 6.36e4 / (20:45 + 273) ^ 2 + -4.7e-3, type = "l", 
     xlab = "MAT", ylab = expression(Delta*47),
     ylim = range(c(6.36e4 / (20:45 + 273) ^ 2 + -4.7e-3, 
                    clump_samples.an$cap_47)))

points(pm$t_cont.an[clump_sites.ind.an], clump_samples.an$cap_47,
       pch = 20, col = col.a)

points(pm$t_seas.an[clump_sites.ind.an], clump_samples.an$cap_47,
       col = col.a)

text(pm$t_cont.an[clump_sites.ind.an], par("usr")[4], clump_sites.ind.an, 
     xpd = NA, pos = 3, col = col.a)
```

* The JPI posterior is a little warm at this site for what you'd expect accounting for seasonality. The closed symbols are MAT, the open ones are summer temperature.
* Note that this same site has a bunch of CIA data shown above, and the posterior MAT is also a bit warmer than expected given those data. This is a sub-tropical site, and the fairly warm SSTs are pulling the MAT values at this and other sub-tropical Anisian sites toward higher values.

## Oxygen isotopes

I'm using carbonate $\delta^{18}O$ data from 3 Carnian and 4 Anisian sites. The lines on the plot show predicted average continental carbonate values as a function of temperature for the model prior (black) and posteriors (colors as throughout). This is derived in two steps:

* Precipitation $\delta^{18}O$ value is a function of temperature. The prior uses the classic modern empirical relationship. The slope and intercept of the relationship are model parameters which are sampled separately and optimized for each stage.
* Carbonate $\delta^{18}O$ value is a function of the local precipitation $\delta^{18}O$ value and temperature. Previously we had been using MAT to describe fractionation in this step, which is reflected in the solid lines below. Now we're using summer season temperature. Because the relationship between MAT and summer temperature is a function of both MAT and latitude we can't plot a single line representing this, but I've approximated the relationship by plotting the predicted carbonate $\delta^{18}O$ values using summer season temperature for a site at 45$^\circ$ latitude (dotted lines).

The data are plotted as above, with the posterior median MAT on the x axis and the measured proxy value on the y axis. 

```{r oxygen, echo=FALSE,fig.height=4}
mat = 0:40
c_cont = -mat + 1.4836 - 0.2738 * 45
t_seas = (-0.9189 + sqrt(0.9189 ^ 2 - 4 * -0.0015 * c_cont)) /
      (2 * -0.0015)
d18O.prior = -14 + 0.59 * mat
d18Oc.prior = (d18O.prior + (2.78e6 / (mat + 273) ^ 2 + -2.98)) * 0.97002 - 29.98
d18Oc.prior.seas = (d18O.prior + (2.78e6 / (t_seas + 273) ^ 2 + -2.98)) * 0.97002 - 29.98

d18O.car.post = c(pm$d.car) + c(pm$e.car) * mat
d18Oc.car.post = (d18O.car.post + (2.78e6 / (mat + 273) ^ 2 + -2.98)) * 0.97002 - 29.98
d18Oc.car.post.seas = (d18O.car.post + (2.78e6 / (t_seas + 273) ^ 2 + -2.98)) * 0.97002 - 29.98

d18O.an.post = c(pm$d.an) + c(pm$e.an) * mat
d18Oc.an.post = (d18O.an.post + (2.78e6 / (mat + 273) ^ 2 + -2.98)) * 0.97002 - 29.98
d18Oc.an.post.seas = (d18O.an.post + (2.78e6 / (t_seas + 273) ^ 2 + -2.98)) * 0.97002 - 29.98

plot(mat, d18Oc.prior, type = "l", 
     ylim = range(c(d18Oc.prior, d18Oc.an.post, d18Oc.car.post,
                    pc_samples.car$d18O, pc_samples.an$d18O)),
     xlab = "MAT", ylab = expression(delta^{18}*"O terrestrial carb"))
lines(mat, d18Oc.prior.seas, lty = 3)
lines(mat, d18Oc.car.post, col = col.c)
lines(mat, d18Oc.car.post.seas, col = col.c, lty = 3)
lines(mat, d18Oc.an.post, col = col.a)
lines(mat, d18Oc.an.post.seas, col = col.a, lty = 3)

points(pm$t_cont.car[pc_sites.ind.car], pc_samples.car$d18O, pch = 20, col = col.c)
points(pm$t_cont.an[pc_sites.ind.an], pc_samples.an$d18O, pch = 20, col = col.a)

text(pm$t_cont.car[pc_sites.ind.car], par("usr")[4], pc_sites.ind.car, 
     xpd = NA, pos = 3, col = col.c)

text(pm$t_cont.an[pc_sites.ind.an], par("usr")[4], pc_sites.ind.an, 
     xpd = NA, pos = 3, offset = 1.25, col = col.a)
```

Here again we have a mix. The Anisian sites, where there is limited non-isotope terrestrial proxy data (only site 1 has other data) tend to have relatively warm posterior MATs for their $\delta^{18}O$ values. For the Carnian sites, the new one from Morocco (site 10) stands out as the clear outlier...very high oxygen isotope values and a low MAT. This is driven by the CIA data from this site, which prefer a lower MAT. The CIA wins out because the oxygen isotope model has more degrees of freedom - site level temperature and precipitation $\delta^{18}O$ can both adjust to accomodate the higher carbonate isotope values, whereas the only freedom in the CIA proxy model is the uncertainty in the calibration relationship.