1. Sample IDs
Pyrolysis samples are often named to encode their production
conditions, e.g. SBC450-30 for a spirulina-derived biochar
pyrolysed at 450 degC for 30 minutes. parse_sbc_id() splits
these back out:
ids <- c("SBC300-15", "SBC300-30", "SBC300-60",
"SBC450-15", "SBC450-30", "SBC450-60",
"SBC600-15", "SBC600-30", "SBC600-60")
parse_sbc_id(ids)
#> id feedstock temperature_C residence_time_min
#> 1 SBC300-15 SBC 300 15
#> 2 SBC300-30 SBC 300 30
#> 3 SBC300-60 SBC 300 60
#> 4 SBC450-15 SBC 450 15
#> 5 SBC450-30 SBC 450 30
#> 6 SBC450-60 SBC 450 60
#> 7 SBC600-15 SBC 600 15
#> 8 SBC600-30 SBC 600 30
#> 9 SBC600-60 SBC 600 60
2. Adsorption capacity and removal efficiency
set.seed(1)
C0 <- rep(50, 9)
temps <- rep(c(300, 450, 600), each = 3)
Ce <- pmax(1, 42 - 0.045 * temps + rnorm(9, 0, 2))
qe_batch(C0, Ce, V = 0.05, m = 0.1)
#> [1] 11.37645 10.56636 11.58563 12.52972 13.79549 14.94547 17.01257 16.76168
#> [9] 16.92422
removal_efficiency(C0, Ce)
#> [1] 45.50582 42.26543 46.34251 50.11888 55.18197 59.78187 68.05028 67.04670
#> [9] 67.69687
3. Isotherms: fitting and comparing multiple models
set.seed(2)
Ce_iso <- c(2, 5, 10, 20, 35, 50, 70, 90)
qe_iso <- (5.2 * 0.12 * Ce_iso) / (1 + 0.12 * Ce_iso) + rnorm(length(Ce_iso), 0, 0.08)
langmuir <- fit_langmuir(Ce_iso, qe_iso)
freundlich <- fit_freundlich(Ce_iso, qe_iso)
temkin <- fit_temkin(Ce_iso, qe_iso)
dr <- fit_dr(Ce_iso, qe_iso)
sips <- fit_sips(Ce_iso, qe_iso)
data.frame(
model = c("Langmuir", "Freundlich", "Temkin", "Dubinin-Radushkevich", "Sips"),
R2 = round(c(langmuir$R2, freundlich$R2, temkin$R2, dr$R2, sips$R2), 4)
)
#> model R2
#> 1 Langmuir 0.9976
#> 2 Freundlich 0.9257
#> 3 Temkin 0.9858
#> 4 Dubinin-Radushkevich 0.8830
#> 5 Sips 0.9977
Comparing R² across models like this is the usual way to decide which
isotherm best describes a given system — here, unsurprisingly, Langmuir
and Sips (which nests Langmuir) come out on top, since the data was
simulated from a Langmuir relationship.
plot_isotherm(langmuir, Ce_iso, qe_iso)
Confidence intervals on the fitted parameters:
fit_confint(langmuir)
#> parameter estimate lower upper
#> 1 Qmax 5.2027170 5.0801297 5.3253043
#> 2 KL 0.1206835 0.1090092 0.1323579
Batch fitting across many samples at once
In practice you’ll often have an isotherm for every sample in a
factorial design, not just one. fit_isotherm_batch() fits a
model separately per group and returns one tidy row per sample:
set.seed(3)
batch_df <- do.call(rbind, lapply(c("SBC300-30", "SBC450-30", "SBC600-30"), function(id) {
Qmax <- runif(1, 3, 6); KL <- runif(1, 0.05, 0.2)
Ce <- c(2, 5, 10, 20, 35, 50, 70)
data.frame(id = id, Ce = Ce,
qe = (Qmax * KL * Ce) / (1 + KL * Ce) + rnorm(length(Ce), 0, 0.05))
}))
fit_isotherm_batch(batch_df, "id", "Ce", "qe", model = "langmuir")
#> id Qmax KL R2 method
#> 1 SBC300-30 3.543857 0.1634148 0.9990243 nonlinear
#> 2 SBC450-30 3.330104 0.1526459 0.9978597 nonlinear
#> 3 SBC600-30 3.421292 0.1261427 0.9986497 nonlinear
4. Kinetics
set.seed(4)
t_kin <- c(5, 15, 30, 60, 120, 180, 240, 360)
qt_kin <- (0.02 * 4.3^2 * t_kin) / (1 + 0.02 * 4.3 * t_kin) + rnorm(length(t_kin), 0, 0.05)
pfo <- fit_pfo(t_kin, qt_kin)
pso <- fit_pso(t_kin, qt_kin)
elovich <- fit_elovich(t_kin, qt_kin)
intraparticle <- fit_intraparticle(t_kin, qt_kin)
data.frame(
model = c("Pseudo-first-order", "Pseudo-second-order", "Elovich", "Intraparticle diffusion"),
R2 = round(c(pfo$R2, pso$R2, elovich$R2, intraparticle$R2), 4)
)
#> model R2
#> 1 Pseudo-first-order 0.9716
#> 2 Pseudo-second-order 0.9980
#> 3 Elovich 0.9233
#> 4 Intraparticle diffusion 0.7267
plot_kinetics(pso, t_kin, qt_kin, model = "pso")
5. Thermodynamics
Given a distribution coefficient Kc measured at several
temperatures, fit_vant_hoff() recovers standard enthalpy,
entropy, and Gibbs free energy of adsorption:
temperature_K <- c(298, 308, 318, 328)
Kc <- c(1.8, 2.3, 2.9, 3.4)
vh <- fit_vant_hoff(temperature_K, Kc)
vh[c("deltaH_kJmol", "deltaS_Jmolk", "R2")]
#> $deltaH_kJmol
#> [1] 17.42241
#>
#> $deltaS_Jmolk
#> [1] 63.44336
#>
#> $R2
#> [1] 0.9954749
vh$table
#> temperature_K Kc deltaG_kJmol
#> 1 298 1.8 -1.483710
#> 2 308 2.3 -2.118144
#> 3 318 2.9 -2.752577
#> 4 328 3.4 -3.387011
6. FTIR: baseline correction, peak picking, and functional
groups
wn <- seq(4000, 400, by = -4)
gaussian_peak <- function(x, center, height, width) height * exp(-((x - center)^2) / (2 * width^2))
baseline <- 0.08 + 0.00003 * (4000 - wn)
spectrum_raw <- data.frame(
wavenumber_cm1 = wn,
intensity = baseline +
gaussian_peak(wn, 3400, 0.35, 120) + gaussian_peak(wn, 2920, 0.15, 40) +
gaussian_peak(wn, 1710, 0.22, 25) + gaussian_peak(wn, 1600, 0.28, 20) +
gaussian_peak(wn, 1030, 0.30, 35)
)
plot_ftir_spectrum(spectrum_raw)
spectrum <- baseline_correct(spectrum_raw)
find_ftir_peaks(spectrum, window = 10, min_prominence = 0.05)
#> peak_cm1 intensity prominence
#> 1 3400 0.3499989 0.3499989
#> 7 1028 0.2995104 0.2995104
#> 5 1600 0.2800133 0.2800133
#> 4 1708 0.2192968 0.2192968
#> 2 2920 0.1501165 0.1501165
functional_group_density(spectrum)
#> group range_low_cm1 range_high_cm1 area mean_intensity
#> 1 O-H stretch 3200 3550 88.796395 0.25368785
#> 2 Aliphatic C-H 2850 2950 10.734066 0.11040199
#> 3 C=O stretch 1690 1750 9.638130 0.16860458
#> 4 Aromatic C=C 1580 1620 9.562558 0.23278514
#> 5 C-O / C-O-C stretch 1000 1300 21.162456 0.07098025
#> 6 Si-O / mineral 1000 1100 20.560622 0.20247465
7. XRD: crystallinity index via peak deconvolution
set.seed(5)
two_theta <- seq(10, 40, by = 0.1)
intensity <- 300 * exp(-((two_theta - 22)^2) / (2 * 0.8^2)) +
150 * exp(-((two_theta - 29)^2) / (2 * 0.6^2)) +
rnorm(length(two_theta), 0, 3)
intensity <- pmax(0, intensity)
xrd <- xrd_deconvolve(two_theta, intensity, peak_centers = c(22, 29))
xrd$peaks
#> center height width area
#> 1 21.99692 300.0166 0.7994207 601.1883
#> 2 29.00238 149.7794 0.6052389 227.2317
xrd$crystallinity_index
#> [1] 0.9712719
If you already have crystalline/amorphous peak areas from your own
XRD software, [xrd_crystallinity_index()] computes the index directly
without needing to fit anything in R.
8. BET surface area
P_P0 <- c(0.05, 0.10, 0.15, 0.20, 0.25, 0.30)
Q <- c(12.1, 15.8, 18.9, 21.6, 24.1, 26.8)
bet_surface_area(P_P0, Q)
#> $Qm_cm3g
#> [1] 20.63382
#>
#> $C
#> [1] 21.36484
#>
#> $SBET_m2g
#> [1] 89.81031
#>
#> $R2
#> [1] 0.9987648
#>
#> $model
#>
#> Call:
#> stats::lm(formula = y ~ P_P0)
#>
#> Coefficients:
#> (Intercept) P_P0
#> 0.002268 0.046196
10. Correlation matrix
pH <- 6.5 + 0.004 * temps + rnorm(9, 0, 0.15)
EC <- 0.8 + 0.003 * temps + rnorm(9, 0, 0.08)
df <- data.frame(temperature_C = temps, pH = pH, EC_dSm = EC)
correlation_matrix(df, method = "spearman")
#> $estimate
#> temperature_C pH EC_dSm
#> temperature_C 1.0000000 0.9486833 0.9486833
#> pH 0.9486833 1.0000000 0.9166667
#> EC_dSm 0.9486833 0.9166667 1.0000000
#>
#> $p_value
#> temperature_C pH EC_dSm
#> temperature_C 0.000000e+00 9.584591e-05 9.584591e-05
#> pH 9.584591e-05 0.000000e+00 5.066191e-04
#> EC_dSm 9.584591e-05 5.066191e-04 0.000000e+00
#>
#> $method
#> [1] "spearman"