The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

Introduction to biocharkit

library(biocharkit)

All data in this vignette is synthetically generated with a fixed random seed, purely to demonstrate the package’s functions end to end. It is not real experimental data.

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)

Langmuir isotherm fit with fitted curve overlaid on observed points

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")

Pseudo-second-order kinetics fit with fitted curve overlaid on observed points

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)

Synthetic FTIR spectrum with five characteristic biochar functional-group peaks

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

9. Proximate and ultimate analysis

proximate_analysis(moisture_pct = 3.2, volatile_matter_pct = 28.5, ash_pct = 12.1)
#>   moisture_pct VM_pct ash_pct fixed_carbon_pct
#> 1          3.2   28.5    12.1             56.2
ultimate_ratios(C_pct = 65.2, H_pct = 2.8, O_pct = 18.4, N_pct = 1.1)
#>    HC_ratio  OC_ratio   NC_ratio
#> 1 0.5117161 0.2118637 0.01446702

10. TGA: DTG, decomposition stages, and Kissinger kinetics

A single synthetic TGA curve, with a moisture step and one main devolatilization step:

Temp <- seq(25, 800, by = 5)
Weight <- 100 - 5 * (Temp > 110) -
  40 / (1 + exp(-(Temp - 340) / 20)) + rnorm(length(Temp), 0, 0.15)
tga <- data.frame(temperature_C = Temp, weight_pct = pmax(Weight, 20))

dtg <- tga_dtg(tga)
peaks <- find_dtg_peaks(dtg, min_prominence = 0.05)
peaks
#>   peak_temperature_C    dtg_max prominence
#> 4                340 0.50836054 0.53096217
#> 2                115 0.35255407 0.37938075
#> 3                220 0.03104519 0.05364682
assign_dtg_peaks(peaks$peak_temperature_C)
#>   peak_temperature_C         stage
#> 1                340     Cellulose
#> 2                115          <NA>
#> 3                220 Hemicellulose
#>                                                        description
#> 1 Cellulose decomposition (usually the sharpest, largest DTG peak)
#> 2                                                             <NA>
#> 3      Hemicellulose decomposition (shoulder or distinct DTG peak)

tga_stages() reads moisture/volatile-matter/ash off the curve at three temperature breakpoints and hands them to proximate_analysis():

tga_stages(tga)
#>   moisture_pct   VM_pct  ash_pct fixed_carbon_pct moisture_end_C vm_end_C
#> 1     0.231178 44.69406 54.79319        0.2815708            110      650
plot_tga(tga)

TG and DTG curves against temperature

Kissinger kinetics needs the same decomposition step observed at several heating rates, not a single curve; here from DTG peak temperatures recorded across four heating rates:

beta <- c(5, 10, 15, 20)          # deg C / min
Tp_C <- c(320, 335, 344, 351)     # DTG peak temperature at each beta
kis <- tga_kinetics_kissinger(beta, Tp_C)
kis[c("Ea_kJmol", "A_min1", "R2")]
#> $Ea_kJmol
#> [1] 127.9724
#> 
#> $A_min1
#> [1] 40827506685
#> 
#> $R2
#> [1] 0.9998787
fit_confint(kis)
#>     parameter     estimate        lower        upper
#> 1 (Intercept)     14.79099     13.94572     15.63627
#> 2           x -15392.39420 -15908.28127 -14876.50713

11. 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.9500000
#> EC_dSm            0.9486833 0.9500000 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 8.762524e-05
#> EC_dSm         9.584591e-05 8.762524e-05 0.000000e+00
#> 
#> $method
#> [1] "spearman"

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.