### Appendix F
### Included in Denny et al., 2019: A History of Pore Water Oxygen Isotope Evolution in the Cretaceous Travis Peak Formation.

## burmin variable input list and definitions:
# "init_detmin_d18o" - Initial d18O value of detrital quartz, in permil
# "init_wat_d18o" - Initial d18O value of water, in permil
# "depth_starting" - Initial depth (start of model)
# "depth_final" - Final depth (end of model)
# "depth_increment" - A positive number, depth_starting is the lower bound
# "INPUT_porosity" - A dataframe for inputing the formation's historical porosity curve and consisting of two columns, the first containing depth numbers, and the second containing porosity values (as % of total rock)
# "INPUT_mech" - A dataframe for inputing the current proportion of porosity loss that is occurring due to mechanical compaction and consisting of two columns, the first containing depth numbers, and the second containing proportion mechanical compaction values (a decimal number equal to [incremental porosity loss due to mechanical compaction]/[total incremental porosity loss])
# "INPUT_geotherm" - A dataframe for inputing the local historical geotherm and consisting of two columns, the first containing depth numbers, and the second containing temperatures (in Celcius)
# "cem_diss_prop" - Multiplier of cement redissolved by pressure solution, relative to the volume ratio of overgrowth quartz to total quartz (default setting is 0.5)


## EXAMPLE INPUTS (FOR INTERNAL CODE TESTING)
#porosity_input_std <- data.frame(depth_por = c(0,1000,1995,2239,2300,2393,2494,2521,2580,2646,2828,2872,2915,3063,3133,3216,3291,3386,3546,3786,5000),porosity = c(40,24,13.24,11.12,10.59,9.92,9.24,9.04,8.65,8.29,7.27,7.04,6.81,6.14,5.81,5.50,5.21,4.88,4.35,3.66,2));
#mech_input_std <- data.frame(depth_std = c(0, 1000, 1500, 2000, 3000, 4000, 5000),mechanical_contribution = c(1, 0.7, 0.5, 0.3, 0.05, 0.025, 0.01));
# 40 degree Celcius/km geotherm
#geotherm_input_std <- data.frame(depth_std = c(0, 1000, 2000, 3000, 4000, 5000),geotherm = c(20, 60, 100, 140, 180, 220));
#init_detmin_d18o <- 13
#init_wat_d18o <- 0
#depth_starting <- 0
#depth_final <- 5000
#depth_increment <- 10
#INPUT_porosity <- porosity_input_std
#INPUT_mech <- mech_input_std
#INPUT_geotherm <- geotherm_input_std
#cem_diss_prop <- 0.5;

## burmin function
burmin <- function(init_detmin_d18o,init_wat_d18o,depth_starting,depth_final,depth_increment,INPUT_porosity,INPUT_mech,INPUT_geotherm,cem_diss_prop=0.5)
  {

  # Fractionation A and B numbers
  qtzA = 3.38; qtzB = -2.90;
  calA = 2.78; calB = -2.89
  
  # Mineral densities
  DEN_wat = 1;
  DEN_qtz = 2.65;
  DEN_cal = 2.71;
  
  # Define some other constants
  AMU_O = 15.999; AMU_H = 1.008; AMU_Si = 28.085; AMU_C = 12.0107; AMU_Ca = 40.078
  avo = 6.0221409*(10^23); # avogaddro's number
  AMU_wat = (AMU_H*2)+(AMU_O); # Atomic mass of water
  AMU_qtz = (AMU_Si)+(AMU_O*2); # Atomic mass of qtz
  AMU_cal = (AMU_Ca)+(AMU_C)+(AMU_O*3); # Atomic mass of cal
  
  PRP_qtz_o = 2; # Ratio of oxygen atoms to formula molecules in quartz
  PRP_qtz_si = 1; # Ratio of silicon atoms to formula molecules in quartz
  PRP_cal_o = 3; # Ratio of oxgygen atoms to formula molecules in calcite
  PRP_cal_ca = 1; # Ratio of calcium atoms to formula molecules in calcite
  PRP_wat_o = 1; # Ratio of oxygen atoms to formula molecules in water
  
  VSMOW = 2005.2/1000000; # ratio of 18o/16o, 2005.2 parts 18o to 1 million 16o
  
  # Assign mineral-specific values to the mineral phase being buried
  minA = qtzA; minB = qtzB; DEN_min = DEN_qtz;
  AMU_min = AMU_qtz; PRP_min_cat = PRP_qtz_si; PRP_min_o = PRP_qtz_o;
  #minA = calA; minB = calB; DEN_min = DEN_cal;
  #AMU_min = AMU_cal; PRP_min_cat = PRP_cal_ca; PRP_min_o = PRP_cal_o;
  
  i = 1; # i will keep track of current row in dataframe
  
  depth_vect = seq(depth_starting,depth_final,depth_increment)
  df_l = length(depth_vect)
  
  df = data.frame(depth = depth_vect, # Initialized at start
                  in_por = as.double((approx(INPUT_porosity[[1]],INPUT_porosity[[2]],depth_vect))[[2]]),
                  in_mech = as.double((approx(INPUT_mech[[1]],INPUT_mech[[2]],depth_vect))[[2]]),
                  in_temp = as.double((approx(INPUT_geotherm[[1]],INPUT_geotherm[[2]],depth_vect))[[2]]),
                  porosity_increment = double(df_l),
                  porosity_chem = double(df_l),
                  num_min_o_diss = double(df_l),
                  num_min_o_precip = double(df_l),
                  num_ovgmin_18o = double(df_l),
                  num_ovgmin_16o = double(df_l),
                  num_detmin_18o = double(df_l),
                  num_detmin_16o = double(df_l),
                  num_wat_o = double(df_l),
                  num_wat_18o = double(df_l),
                  num_wat_16o = double(df_l),
                  current_volume = double(df_l),
                  cum_porosity_chem_loss = double(df_l),
                  cum_porosity_mech_loss = double(df_l),
                  d18o_ovgmin_new = double(df_l))
  
  rat_wat = ((init_wat_d18o/1000)+1)*VSMOW # initial ratio of 18o to 16o in water/soluble phases
  rat_min = ((init_detmin_d18o/1000)+1)*VSMOW # initial ratio of 18o to 16o in min
  num_wat_o_init = (df$in_por[1]/100)*DEN_wat*(avo/AMU_wat)*PRP_wat_o; # initial number of o atoms in water

      # Set initial row values
  df$porosity_increment[-1] = -diff(df$in_por); # Change in porosity with each depth increment
  df$porosity_chem = df$porosity_increment * (1-df$in_mech); # Porosity change with each depth increment not lost to mechanical compaction;
  df$num_min_o_diss = NA; # Number of o atoms dissolved
  df$num_min_o_precip = NA; # Number of o atoms precipitated
  df$num_min_o[1] = ((100-df$in_por[1])/100)*DEN_min*(avo/AMU_min)*PRP_min_o; # Sum of 16o and 18o atoms in mineral phase
  df$num_ovgmin_16o[1] = 0; # Number of 16o atoms in overgrowth mineral phase
  df$num_ovgmin_18o[1] = 0; # Number of 18o atoms in overgrowth mineral phase
  df$num_detmin_16o[1] = df$num_min_o[1]/(rat_min+1); # Number of 16o atoms in detrital mineral phase
  df$num_detmin_18o[1] = df$num_min_o[1]-df$num_detmin_16o[1]; # Number of 18o atoms in detrital mineral phase
  df$num_wat_o[1] = num_wat_o_init; # Total number of atoms in water, as determined by current volume and remaining porosity
  df$num_wat_16o[1] = num_wat_o_init/(rat_wat+1); # Number of 16o atoms in water
  df$num_wat_18o[1] = num_wat_o_init-df$num_wat_16o[1]; # Number of 18o atoms in water
  df$cum_porosity_chem_loss = cumsum(df$porosity_chem); # Cumulative porosity lost to dissolution
  df$cum_porosity_mech_loss = cumsum(df$porosity_increment*df$in_mech); # Cumulative porosity lost to mech compaction
  df$current_volume = (100 - df$cum_porosity_mech_loss - df$cum_porosity_chem_loss)/100; # Modeled volume relative to original volume (1); volume is reduced by mechanical compaction and chemical compaction, but chemical compaction only reduces volume if si is not imported
  df$d18o_ovgmin_new[i] = NA;

   # Bury the rock incrementally
  while(df$depth[i] < depth_final){
    i = i + 1
    df$num_min_o_diss[i] = (1-(df$in_por[i]/100))*(df$porosity_chem[i]/100)*DEN_min*(avo/AMU_min)*PRP_min_o; # Number of o atoms dissolved with current porosity change;
    df$num_min_o_precip[i] = (1-(df$in_por[i]/100))*(df$porosity_chem[i]/100)*DEN_min*(avo/AMU_min)*PRP_min_o; # Number of o atoms that will be precipitated with current porosity change;
    
    # Now to calculate how much dissolves, and from where it dissolves
    if (df$num_min_o_diss[i]*cem_diss_prop > df$num_ovgmin_18o[i-1]+df$num_ovgmin_16o[i-1]){ # If there aren't enough atoms in the overgrowth cement, sum overgrowth and detrital atoms and dissolve the combined amount
      PRPmin_diss = df$num_min_o_diss[i]/(df$num_ovgmin_18o[i-1] + df$num_ovgmin_16o[i-1] + df$num_detmin_18o[i-1] + df$num_detmin_16o[i-1]) #proportion of oxygen atoms dissolved from all min o atoms
      dslv_18o = PRPmin_diss * (df$num_ovgmin_18o[i-1] + df$num_detmin_18o[i-1])
      dslv_16o = PRPmin_diss * (df$num_ovgmin_16o[i-1] + df$num_detmin_16o[i-1])
      df$num_ovgmin_18o[i] = df$num_ovgmin_18o[i-1] - (df$num_ovgmin_18o[i-1] * PRPmin_diss)
      df$num_ovgmin_16o[i] = df$num_ovgmin_16o[i-1] - (df$num_ovgmin_16o[i-1] * PRPmin_diss)
      df$num_detmin_18o[i] = df$num_detmin_18o[i-1] - (df$num_detmin_18o[i-1] * PRPmin_diss)
      df$num_detmin_16o[i] = df$num_detmin_16o[i-1] - (df$num_detmin_16o[i-1] * PRPmin_diss)
    } else { # If there are enough atoms in the overgrowth cement to process
      PRPdetmin_diss = (df$num_min_o_diss[i]*(1-cem_diss_prop))/(df$num_detmin_18o[i-1] + df$num_detmin_16o[i-1]) #proportion of oxygen atoms dissolved from detrital quartz
      PRPovgmin_diss = (df$num_min_o_diss[i]*cem_diss_prop)/(df$num_ovgmin_18o[i-1] + df$num_ovgmin_16o[i-1]) #proportion of oxygen atoms dissolved from overgrowth quartz
      dslv_18o = (PRPovgmin_diss * df$num_ovgmin_18o[i-1]) + (PRPdetmin_diss * df$num_detmin_18o[i-1])
      dslv_16o = (PRPovgmin_diss * df$num_ovgmin_16o[i-1]) + (PRPdetmin_diss * df$num_detmin_16o[i-1])
      df$num_ovgmin_18o[i] = df$num_ovgmin_18o[i-1] - (df$num_ovgmin_18o[i-1] * PRPovgmin_diss)
      df$num_ovgmin_16o[i] = df$num_ovgmin_16o[i-1] - (df$num_ovgmin_16o[i-1] * PRPovgmin_diss)
      df$num_detmin_18o[i] = df$num_detmin_18o[i-1] - (df$num_detmin_18o[i-1] * PRPdetmin_diss)
      df$num_detmin_16o[i] = df$num_detmin_16o[i-1] - (df$num_detmin_16o[i-1] * PRPdetmin_diss)
    }
    
    # Move dissolved atoms to water
    df$num_wat_18o[i] = df$num_wat_18o[i-1] + dslv_18o;
    df$num_wat_16o[i] = df$num_wat_16o[i-1] + dslv_16o;

    # Calculate number of atoms that will precipitate as overgrowth
    ratio = (exp((((minA*1000000)/((273.15+df$in_temp[i])^2))+minB)/1000)) * (df$num_wat_18o[i]/df$num_wat_16o[i]) #calculates equilibrium ratio value of 18o/16o for mineral phase
    num_ovgmin_16o_new = df$num_min_o_precip[i]/(ratio + 1)
    num_ovgmin_18o_new = df$num_min_o_precip[i] - num_ovgmin_16o_new
    
    # Remove oxygen atoms from water pool during cement precipitation
    df$num_wat_18o[i] = df$num_wat_18o[i]-num_ovgmin_18o_new
    df$num_wat_16o[i] = df$num_wat_16o[i]-num_ovgmin_16o_new

    # Add oxygen atoms from water pool to overgrowth pool during cement precipitation
    df$num_ovgmin_16o[i] = df$num_ovgmin_16o[i] + num_ovgmin_16o_new
    df$num_ovgmin_18o[i] = df$num_ovgmin_18o[i] + num_ovgmin_18o_new

    df$num_wat_o[i] = (df$in_por[i]/100)*df$current_volume[i]*DEN_wat*(avo/AMU_wat)*PRP_wat_o; # Number of o atoms in water/soluble phases per unit interaction volume (recalculate here with new porosity, so water can be expelled)

    wat_o_tot = df$num_wat_18o[i] + df$num_wat_16o[i];
    df$num_wat_18o[i] = df$num_wat_18o[i]*(df$num_wat_o[i]/(wat_o_tot)); # correct num_wat_18o to new water volume
    df$num_wat_16o[i] = df$num_wat_16o[i]*(df$num_wat_o[i]/(wat_o_tot)); # correct num_wat_16o to new water volume

    df$d18o_ovgmin_new[i] = ((num_ovgmin_18o_new/num_ovgmin_16o_new)/(VSMOW)-1)*1000;
  }
  
  # Calculations made at end of simulation
  df$d18o_wat = ((df$num_wat_18o/df$num_wat_16o)/(VSMOW)-1)*1000;
  df$d18o_ovgmin_bulk = ((df$num_ovgmin_18o/df$num_ovgmin_16o)/(VSMOW)-1)*1000;
  df$d18o_allmin_bulk = (((df$num_ovgmin_18o+df$num_detmin_18o)/(df$num_ovgmin_16o+df$num_detmin_16o))/(VSMOW)-1)*1000;
  df$num_min_o = df$num_ovgmin_18o + df$num_ovgmin_16o + df$num_detmin_18o + df$num_detmin_16o
  df$prop_wat_of_original_vol = df$num_wat_o/(DEN_wat*(avo/AMU_wat)*PRP_wat_o)
  df$prop_det_of_original_vol = (df$num_detmin_16o+df$num_detmin_18o)/(DEN_min*(avo/AMU_min)*PRP_min_o)
  df$prop_cem_of_original_vol = (df$num_ovgmin_16o+df$num_ovgmin_18o)/(DEN_min*(avo/AMU_min)*PRP_min_o)
  df$prop_wat_of_current_vol = df$prop_wat_of_original_vol/df$current_volume
  df$prop_det_of_current_vol = df$prop_det_of_original_vol/df$current_volume
  df$prop_cem_of_current_vol = df$prop_cem_of_original_vol/df$current_volume
  df$prop_cem_of_min = (df$num_ovgmin_18o+df$num_ovgmin_16o)/(df$num_ovgmin_18o+df$num_ovgmin_16o+df$num_detmin_18o+df$num_detmin_16o) # Fraction cement of total mineral volume

  return(df)
}

## EXAMPLE RUN AND OUTPUT (FOR TESTING)
#test<-burmin(init_detmin_d18o,init_wat_d18o,depth_starting,depth_final,depth_increment,INPUT_porosity,INPUT_mech,INPUT_geotherm)
#library(reshape2)
#library(ggplot2)
#test_d <- melt(test, "depth")
#ggplot(test_d, aes(value)) + 
#  geom_histogram(bins=30) + 
#  facet_wrap(~variable, scales = "free")
#ggplot(test_d, aes(value,-depth)) + 
#  geom_path() + 
#  facet_wrap(~variable, scales = "free")