## Generate a synthetic dataset
# Losely based on the example of the Cayambe volcano
# Samaniego P, Martin H, Robin C, Monzier M (2002) Transition from calc-alkalic to adakitic magmatism at Cayambe volcano, Ecuador: Insights into slab melts and mantle wedge interactions. Geology 30 (11):967-970
##########################################################
# Working dir, need to be adjusted to your system !!
setwd("C:/Users/moyen/Documents/Recherche/Springer Book")
#Definitions
############
mjrs<-c("SiO2","Al2O3","FeOt","MgO","MnO","CaO","Na2O","K2O","TiO2")
trc<-c("Rb","Sr","Ba","Zr","V","Ni","Y","La","Nb","Ta","X")
mins<-c("Pl","Hb","Mt","Zrc")
# Molecular weights, needed for M
# MW does exist in GCDkit, but it is (re)defined here to allow the code to run without GCDkit
MWt<-c(60.0848,101.9613,71.8464,40.3044,70.9374,56.0774,61.9789,94.1954,79.8658)
names(MWt)<-mjrs
mw.coefs<-c(1,2,1,1,1,1,2,2,1)
names(mw.coefs)<-mjrs
# C0 -- the most mafic lava of the suite (CAY56 from Samaniego et al.)
c0.maj<-c(56.5,17.08,6.90,4.24,0.11,7.08,3.54,1.31,0.88)
names(c0.maj)<-mjrs
c0.trc<-c(30,570,620,140,200,50,18,18,5.5,0.4,1)
names(c0.trc)<-trc
c0<-c(c0.maj,c0.trc)
# Mineral compositions
# We use only 4 minerals here (plag, amphibole, magnetite, zircon)
# As the "real" modelling shows that this the likely fractionnating assemblage is amp+pg+mt
min.comp<-matrix(c(53.11,28.98,0.57,0.08,0,11.71,4.8,0.32,0.09,
42.24,12.15,12.5,14.04,0.2,11.39,2.63,0.61,2.82,
0.6,1.29,94,0.39,0.14,0.01,0,0.02,3,
50,0,0,0,0,0,0,0,0),
nrow=4,byrow=T)
colnames(min.comp)<-mjrs
rownames(min.comp)<-mins
# Partition coefs
Kd<-matrix(c(
0.04,4.4,0.5,0.01,0.01,0.38,0.055,0.4,0.025,0.05,0,
0.014,0.022,0.044,0.02,10,12,5,0.74,0.9,0.5,0,
0.00001,0.00001,0.00001,0.00001,8.6,8.6,0.00001,0.22,0.00001,0.04,0,
0.00001,0.00001,0.00001,3800,0.00001,0.00001,200,2,25,50,1000),
nrow=4,byrow=T)
colnames(Kd)<-trc
rownames(Kd)<-mins
# Mineral proportions (from actual model)
min.props<-c(0.45,0.5,0.05,0)
names(min.props)<-mins
# Range of F values (from actual model, tweaked to make sure no major gets < 0)
fmin<-0.46
fmax<-1
fstep<--0.01
ff<-seq(fmax,fmin,fstep)
# Temperature defined as a function of F
# Function adjusted so that we have
# ~1020 °C at F = 1 (SiO2 = 56%)
# ~720 °C at F = 0.45 (SiO2 = 72%)
# Probably realistic based on experimental literature
t.c<-function(f.n){
return(460+560*f.n)
}
# Water content
# Also tweaked as a function of F
# It's main use is to have totals < 100 % (LOI)
tot.anh<-function(f.n){
return(99.6-1.5*(1-f.n))
}
# Calculate the cumulate and D
cs<-min.props%*%min.comp
D<-min.props%*%Kd
# Theoretical evolution (fwd model)
###################################
# For each value of F, we calculate major, traces, and Zr saturation
qq<-sapply(ff,function(f.n){
# Major elements
maj<-(c0[mjrs]-(1-f.n)*cs[1,mjrs])/f.n
maj<-maj/sum(maj)*tot.anh(f.n) # Majors renormalized to match the expected anydrous total
# M parameter (Watson & Harrison)
ee<-maj[mjrs]/MWt[mjrs]*mw.coefs[mjrs]
ee<-ee/sum(ee)
M<-(ee["Na2O"]+ee["K2O"]+2*ee["CaO"])/(ee["SiO2"]*ee["Al2O3"])
# Traces (first pass, without Zr sat)
tr<-c0[trc]*f.n^(D[1,trc]-1)
# Zr sat
DZr<-exp(-3.8-0.85*(M-1)+12900/(t.c(f.n)+273))
Zr.sat<-497644/DZr
# Correction for Zrc
# Technically wrong, zrc amount treated as batch, probably close enough...
if(tr["Zr"]