Skip to content

Commit 891cefb

Browse files
ActiveMarginsActiveMargins
authored
Add files via upload
1 parent a73196a commit 891cefb

2 files changed

Lines changed: 203 additions & 5 deletions

File tree

‎R/LithoGas-package.R‎

Lines changed: 10 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,7 @@
55
#' @importFrom ggplot2 ggplot aes geom_line geom_ribbon theme_bw xlab ylab
66
#' scale_x_log10 scale_y_continuous sec_axis labs scale_y_log10
77
#' @importFrom rlang sym
8-
#' @importFrom magrittr %>%
8+
#' @importFrom dplyr %>% mutate bind_rows rename_with ends_with
99
"_PACKAGE"
1010

1111
utils::globalVariables(c(
@@ -38,8 +38,7 @@ utils::globalVariables(c(
3838
"Fe3FeTInitalRatMean", "Fe3FeTInitalRatSD",
3939
"Fe3FeTRatCurMin", "Fe3FeTRatCurMax",
4040
"Fe3FeTRatCurMean", "Fe3FeTRatCurSD",
41-
"mass",
42-
"SerpH2_molm3", "SerpH2Rate_molm3yr",
41+
"mass", "SerpH2_molm3", "SerpH2Rate_molm3yr",
4342

4443
#monteSum()
4544
"sampleField",
@@ -52,6 +51,12 @@ utils::globalVariables(c(
5251
#plotting - monteH2Plot() and monteHePlot()
5352
"vol_km", "colorField",
5453
"H2RateMin_molyr", "H2RateMean_molyr", "H2RateMax_molyr",
55-
"HeRateMin_molyr", "HeRateMean_molyr", "HeRateMax_molyr"
56-
54+
"HeRateMin_molyr", "HeRateMean_molyr", "HeRateMax_molyr",
55+
56+
#deepTime - deepTimeProd()
57+
"U238ppm","U238ppm_T","U235ppm","U235ppm_T", "Uppm_T", "Thppm_T", "Kpct_T","EKA_T",
58+
"EThA_T", "EUA_T", "W_T", "EThB_T", "EUB_T","EKB_T", "EKG_T", "EThG_T", "EUG_T",
59+
"ENetA_T", "sysDen_T", "ENetB_T", "ENetG_T", "YH2A_T", "YH2B_T", "YH2G_T",
60+
"Uppm", "Thppm", "Kpct", "fluDen", "Fe2O3T","Fe2O3","FeO","FeT","Fe3FeTRatCur",
61+
"Fe3FeTInitalRat","Fe3FeTRatDiff","Fe3O4Diff_wt","molFe3O4","SerpModel", "RadModel"
5762
))

‎R/deepTimeProd.R‎

Lines changed: 193 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,193 @@
1+
#' Deep-time production of H2 and He following radioactive decay and average
2+
#' rate of serpentinization
3+
#'
4+
#' For radiolysis models, present-day measured concentrations of uranium, thorium,
5+
#' and potassium, are back calucualted to deep-time concentrations at specified ages
6+
#' by reversing the radioactive decay equations for each relevant isotope (²³⁸U, ²³⁵U, ²³²Th, ⁴⁰K).
7+
#' These deep-time concentrations are then used to calculate H2 and He production rates
8+
#' via radiolysis.
9+
#'
10+
#' For serpentinization models the average rate of H2 generation by serpentinization
11+
#' that is established by \code{\link{monteProd}} is used.
12+
#'
13+
#' @param monteProdDF Dataframe. A data frame returned by \code{\link{monteProd}}.
14+
#' Long-dataframe including all trial of input Monte Carlo.
15+
#' @param startAge_Ma Numeric. Geologic age in millions of years (Ma) for the
16+
#' model to start at (youngest age). startAge_Ma must be less than endAge_Ma.
17+
#' @param endAge_Ma Numeric. Geologic age in millions of years (Ma) for the
18+
#' model to end at (oldest age). endAge_Ma must be greater than start.
19+
#' @param stepAge_Ma Numeric. Geologic duration of each step in the age model
20+
#' in millions of years (My). Must be positive.
21+
#' @param rad Logic. If true, back calculation of radiolysis using radioactive
22+
#' decay will be computed. This is explained below in.
23+
#' @param serp Logic. If true, back calculation of serpentinization rate
24+
#' will be computed.
25+
#'
26+
#' @details
27+
#' The back-projection uses the standard radioactive decay law in reverse:
28+
#'
29+
#' \deqn{N_{\text{past}} = N_{\text{now}} \cdot e^{+\lambda t}}
30+
#'
31+
#' where \eqn{\lambda} is the decay constant for a given isotope and \eqn{t}
32+
#' is the elapsed time in years.
33+
#'
34+
#' **Decay constants and half-lives used:**
35+
#' \tabular{lll}{
36+
#' Isotope \tab Half-life \tab Assumed abundance \cr
37+
#' \eqn{{}^{238}}U \tab 4.468 Ga \tab 99.2745\% \cr
38+
#' \eqn{{}^{235}}U \tab 703.8 Ma \tab 0.7204\% \cr
39+
#' \eqn{{}^{232}}Th \tab 14.05 Ga \tab 100\% \cr
40+
#' \eqn{{}^{40}}K \tab 1248 Ma \tab 0.01167\% \cr
41+
#' }
42+
#'
43+
#' @return A the input dataframe repeated for each time step. All Monte Carlo trials from monteProd() back
44+
#' calculated at each time step. If \code{rad==True}, additional columns will be added:
45+
#' \describe{
46+
#' \item{\code{timeStep_Ma}}{The time step of the model in millions of years (Ma).}
47+
#' \item{\code{stepDur_Ma}}{The duration of the time step, millions of years (Ma). Same as input \code{stepAge_Ma}}
48+
#' \item{\code{Uppm}}{Total uranium concentration at time \code{timeStep_Ma} (ppm).}
49+
#' \item{\code{Thppm}}{Total thorium concentration at time \code{timeStep_Ma} (ppm).}
50+
#' \item{\code{Kpct}}{Total potassium concentration at time \code{timeStep_Ma} (wt\%).}
51+
#' \item{\code{Kppm}}{Total potassium concentration at time \code{timeStep_Ma} (wt\%).}
52+
#' \item{\code{U238ppm}}{²³⁸U concentration at time \code{timeStep_Ma} (ppm).}
53+
#' \item{\code{U235ppm}}{²³⁵U concentration at time \code{timeStep_Ma} (ppm).}
54+
#' \item{\code{RadH2Rate_molm3yr}}{H2 generation rate (mol/m3/year), back calculated at time \code{timeStep_Ma}}
55+
#' \item{\code{RadHeRate_molm3yr}}{He generation rate (mol/m3/year), back calculated at time \code{timeStep_Ma}}
56+
#' \item{\code{RadH2_molm3}}{Cumulative volume of H2 generated over the \code{timeStep_Ma} in mols H2 per m3. Calculated by \code{RadH2Rate_molm3yr} * \code{timeStep_Ma}}
57+
#' \item{\code{RadHe_molm3}}{Cumulative volume of He generated over the \code{timeStep_Ma} in mols He per m3. Calculated by \code{RadHeRate_molm3yr} * \code{timeStep_Ma} }
58+
#' }
59+
#'
60+
#' If \code{serp==True}, additional columns will be added:
61+
#' \describe{
62+
#' \item{\code{timeStep_Ma}}{The input time in Ma (returned for reference).}
63+
#' \item{\code{SerpH2Rate_molm3yr}}{H2 generation rate (mol/m3/year), back calculated at time \code{timeStep_Ma}}
64+
#' }
65+
#'
66+
#' @examples
67+
#' data("structuredDF")
68+
#' monteProdDF <- monteProd(structuredDF,50, rad=TRUE, serp=TRUE, allowTotalFeSerp =TRUE)
69+
#' deepDF <- deepTimeProd(monteProdDF,0,2000,100,rad=TRUE)
70+
#'
71+
#' library(ggplot2)
72+
#' ggplot(deepDF) + geom_boxplot(mapping=aes(x=as.factor(timeStep_Ma),y=Uppm, fill=Sample)) +labs(x="Age (Ma)",y="U (ppm)", title="Decrease in U by decay law")
73+
#'
74+
#' ggplot(deepDF) + geom_boxplot(mapping=aes(x=as.factor(timeStep_Ma),y=RadHeRate_molm3yr, fill=Sample)) +labs(x="Age (Ma)",y="He generation rate (mol/m3/yr)", title="Decrease in He generation rate by decay law")
75+
#'
76+
#' @importFrom dplyr %>% mutate bind_rows rename_with ends_with
77+
#'
78+
#' @export
79+
deepTimeProd <- function(monteProdDF, startAge_Ma, endAge_Ma, stepAge_Ma, rad=FALSE, serp=FALSE){
80+
#Add any missing columns to the monte carlo results
81+
desired_cols <- c("Uppm", "Thppm", "Kpct", "rockDen", "porosity", "fluDen",
82+
"sysDen", "W", "EKA", "EKB", "EKG", "EThA", "EThB", "EThG",
83+
"EUA", "EUB", "EUG", "ENetA", "ENetB", "ENetG", "YH2A", "YH2B",
84+
"YH2G", "RadH2Rate_molm3yr", "RadHeRate_molm3y", "Fe2O3T",
85+
"Fe2O3", "FeO", "FeT", "Fe3FeTRatCur", "Fe3FeTInitalRat", "Fe3FeTRatDiff",
86+
"Fe3O4Diff_wt", "molFe3O4", "SerpH2_molm3", "SerpH2Rate_molm3yr",
87+
"SerpModel") #these are the columns we need for the code to run
88+
89+
missing_cols <- setdiff(desired_cols, names(monteProdDF)) # find cols in desired but not in A
90+
monteProdDF[missing_cols] <- NA #set any of the missing columns to NA
91+
92+
#Constants that are used explicitly (not as parts of distributions) just set here for easier programming
93+
#Energies used for H2 radiolysis (see Warr et al., 2023)
94+
SA <- 1.5
95+
SB <- 1.25
96+
SG <- 1.14
97+
GH2A <- 1.32
98+
GH2B <- 0.6
99+
GH2G <- 0.25
100+
AConstant <- 6.023E23
101+
102+
#Decay constants (per year)
103+
lambda_U238 <- log(2) / 4.468e9 # 238U half-life: 4.468 Ga
104+
lambda_U235 <- log(2) / 703.8e6 # 235U half-life: 703.8 Ma
105+
lambda_Th232 <- log(2) / 14.05e9 # 232Th half-life: 14.05 Ga
106+
lambda_K40 <- log(2) / 1248e6 # 40K half-life: 1.248 Ga
107+
108+
#Natural isotopic abundances
109+
abund_U238 <- 0.992745 # fraction of total U
110+
abund_U235 <- 0.007204 # fraction of total U
111+
abund_K40 <- 1.167e-4 # fraction of total K
112+
113+
#Set up model objects
114+
deepTimeProdDF <- NULL #final dataframe that will write to
115+
modelSteps <- seq(from=startAge_Ma,to=endAge_Ma,by=stepAge_Ma)
116+
117+
#Go into for loop for each model step
118+
for (t_Ma in modelSteps){
119+
print(t_Ma)
120+
stepDF <- monteProdDF #sub-dataframe that each time step will write to
121+
122+
if(rad==TRUE){
123+
stepDF <- stepDF %>%
124+
mutate(
125+
U238ppm_T = (Uppm*abund_U238) * exp(lambda_U238 * (t_Ma*1000000)),
126+
U235ppm_T = (Uppm*abund_U235) * exp(lambda_U235 * (t_Ma*1000000)),
127+
Uppm_T = U238ppm_T + U235ppm_T,
128+
Thppm_T = Thppm * exp(lambda_Th232 * (t_Ma*1000000)),
129+
Kpct_T = ((Kpct* 1e4* abund_K40) * exp(lambda_K40 * (t_Ma*1000000))) / 1e4,
130+
sysDen_T = (rockDen*(1-(porosity/100)))+(fluDen*(porosity/100)),
131+
W_T = ((porosity/100)*fluDen)/((1-(porosity/100))*rockDen),
132+
EKA_T = Kpct_T*0 ,
133+
EKB_T = Kpct_T*4.88085E15,
134+
EKG_T = Kpct_T*1.51668E15,
135+
EThA_T = Thppm_T*3.80732E14,
136+
EThB_T = Thppm_T*1.7039E14,
137+
EThG_T = Thppm_T*2.93351E14,
138+
EUA_T = Uppm_T*1.3065E15,
139+
EUB_T = Uppm_T*9.11259E14,
140+
EUG_T = Uppm_T*7.0529E14,
141+
ENetA_T = ((EKA_T+EThA_T+EUA_T)*W_T*SA)/(1+W_T*SA),
142+
ENetB_T = ((EKB_T+EThB_T+EUB_T)*W_T*SB)/(1+W_T*SB),
143+
ENetG_T = ((EKG_T+EThG_T+EUG_T)*W_T*SG)/(1+W_T*SG),
144+
YH2A_T = ((ENetA_T*GH2A)/AConstant)*sysDen_T*10,
145+
YH2B_T = ((ENetB_T*GH2B)/AConstant)*sysDen_T*10,
146+
YH2G_T = ((ENetG_T*GH2G)/AConstant)*sysDen_T*10,
147+
RadH2Rate_molm3yr_T = (YH2A_T + YH2B_T + YH2G_T), #H2 generation rate at time t_Ma - because U, Th, and K concentrations were back calculated to time t_Ma, this is accurate for time t_Ma
148+
RadHeRate_molm3yr_T = ((((3.115E6+1.272E5)*Uppm_T+7.71E5*Thppm_T)/AConstant)*sysDen_T*1E6), #He generation rate at time t_Ma - because U, Th, and K concentrations were back calculated to time t_Ma, this is accurate for time t_Ma
149+
RadH2_molm3 = RadH2Rate_molm3yr_T * (stepAge_Ma*1000000), #number of mols H2 created in that model step assuming H2 generation rate at time t_Ma for the entire step
150+
RadHe_molm3 = RadHeRate_molm3yr_T * (stepAge_Ma*1000000), #number of mols He created in that model step assuming He generation rate at time t_Ma for the entire step
151+
timeStep_Ma=t_Ma,
152+
stepDur_Ma = stepAge_Ma)
153+
stepDF <- stepDF %>%
154+
select(-sysDen,-W,-EKA,-EKB,-EKG,-EThA,-EThB,-EThG,-EUA, -EUB,-EUG,
155+
-ENetA,-ENetB,-ENetG,-YH2A,-YH2B,-YH2G, -RadH2Rate_molm3yr,
156+
-RadHeRate_molm3yr, -RadHeRate_molm3y, -Thppm, -Kpct,
157+
-Uppm ) #Drop some original columns calculation variables
158+
159+
stepDF <- stepDF %>%
160+
rename_with(~ sub("_T$", "", .x), ends_with("_T")) #rename current columns of calculation variables
161+
162+
#stepDF <- cbind(stepDF,radDF)
163+
} else {
164+
stepDF <- stepDF %>%
165+
select(-Uppm, -Thppm, -Kpct, -sysDen,
166+
-W, -EKA, -EKB, -EKG, -EThA, -EThB, -EThG, -EUA,
167+
-EUB, -EUG, -ENetA, -ENetB, -ENetG, -YH2A, -YH2B,
168+
-YH2G, -RadH2Rate_molm3yr, -RadHeRate_molm3yr, -RadModel,
169+
-RadHe_molm3, -RadH2_molm3, -fluDen)
170+
}
171+
172+
if(serp==TRUE){
173+
stepDF$SerpH2Rate_molm3yr_T <- stepDF$SerpH2Rate_molm3yr #This is the already averaged rate of generation, so every step is the
174+
stepDF$SerpH2_molm3_T <- stepDF$SerpH2Rate_molm3yr*(stepAge_Ma*1000000) #this will produce the same amount of H2 per step
175+
stepDF$timeStep_Ma <- t_Ma
176+
stepDF$stepDur_Ma <- stepAge_Ma
177+
stepDF <- stepDF %>%
178+
select(-SerpH2Rate_molm3yr,-SerpH2_molm3,-Fe2O3T, -molFe3O4)
179+
stepDF <- stepDF %>%
180+
rename_with(~ sub("_T$", "", .x), ends_with("_T")) #rename current columns of calculation variables
181+
}else {
182+
stepDF <- stepDF %>%
183+
select(-Fe2O3T, -Fe2O3, -FeO, -FeT, -Fe3FeTRatCur, -Fe3FeTInitalRat, -Fe3FeTRatDiff,
184+
-Fe3O4Diff_wt, -molFe3O4, -SerpH2_molm3, -SerpH2Rate_molm3yr,
185+
-SerpModel)
186+
}
187+
188+
deepTimeProdDF <- bind_rows(deepTimeProdDF,stepDF)
189+
190+
}#next time step
191+
192+
return(deepTimeProdDF)
193+
}

0 commit comments

Comments
 (0)