Experimentally Constrained Woody Biochar Models

• Experimentally constrained woody biochar models are detailed simulations of woody biomass-derived biochars that match experimental pyrolysis data to capture key compositional, structural, and functional traits.
• They integrate atomistic and continuum methods with spectral, porosimetric, and elemental analyses to iteratively optimize and validate structure–function relationships.
• These models enable accurate predictions of adsorption isotherms and reactive transport, advancing applications in carbon sequestration and pollutant remediation.

Experimentally constrained woody biochar models are atomistic or continuum representations of biochar materials produced from woody biomass, in which all compositional, structural, and physical parameters are quantitatively matched to experimental data from pyrolysis at defined temperatures (typically 400–800 °C). These models serve as a foundation for predicting biochar properties, guiding application-specific optimization, and elucidating structure–function relationships underlying carbon sequestration and pollutant remediation. They integrate data from elemental analysis, spectroscopic and porosimetric characterization, and surface functionality measurements, employing iterative parameterization and validation by computational methods such as molecular dynamics and reactive-transport simulation.

1. Experimental Foundations: Characterization and Target Property Definition

The construction of experimentally constrained woody biochar models begins with the compilation of a suite of experimental descriptors characterizing biochars derived from lignocellulosic woody feedstocks processed under controlled pyrolysis conditions. These descriptors include:

• Molar atomic ratios $\mathrm{H}/\mathrm{C}$ and $\mathrm{O}/\mathrm{C}$ determined by ultimate analysis.
• Aromatic carbon fraction, typically quantified by $^{13}$C NMR, FTIR, Raman spectroscopy (e.g., $I_D/I_G$ for defect density), and cross-validated with XPS surface speciation.
• True density $\rho_{\text{true}}$ measured by He pycnometry.
• Microporosity, defined by cumulative pore volume for pore widths $<2$ nm, evaluated using CO$_2$/N$_2$ isotherms and geometric probe-insertion algorithms.
• Functional group population statistics, from Boehm titration (carboxyl, phenolic, lactonic) and spectroscopic peak assignments.

Empirical confidence intervals for these descriptors (e.g., $\mathrm{H}/\mathrm{C} = 0.23 \pm 0.13$, $\mathrm{O}/\mathrm{C} = 0.07 \pm 0.08$, aromaticity $96\% \pm 7\%$ at 600°C) are enforced as hard constraints throughout model parametrization (Wood et al., 2023, Ngambia et al., 2024, Wood et al., 27 Dec 2025).

2. Molecular Model Construction: Building Block Selection and Topology

Atomistic molecular models are assembled via a modular workflow:

• Selection of basic structural units (BSUs) tuned to experimental $\mathrm{H}/\mathrm{C}$, $\mathrm{O}/\mathrm{C}$, and aromaticity. BSUs encompass polyaromatic domains of predefined domain size (ADS), functionalized with oxygenated groups (–OH, –C=O, –COOH) and aliphatic bridges. For example, BSU I/II (ADS 22–33) matches $\mathrm{H}/\mathrm{C} \approx 0.28$, BSU IV (large polycondensed, ADS $\sim$425) enforces $\mathrm{H}/\mathrm{C} \approx 0.11$ (Ngambia et al., 2024).
• Combinatorial packing of BSUs in a periodic box to produce target atom ratios. Base models such as BCMA and BCMB are constructed with distinct BSU ratios, engineered to match target elemental analysis.
• Topological control over non-hexagonal ring content (pentagons, heptagons) is used to alter packing density and curvature, contributing to intrinsic ultra-microporosity (Wood et al., 2023).

Bond connectivity, crosslink density, and block arrangements are optimized to avoid steric clashes (Packmol) and human-imposed overlap constraints. Parameter assignment is executed via OPLS-AA force field and geometry optimization routines (LigParGen, PolyParGen).

3. Simulation Protocols: Condensation, Equilibration, and Porosity Engineering

To densify and equilibrate the amorphous carbon matrix, melt–quench molecular dynamics protocols are employed:

• Initial NPT simulations at elevated temperature (1000 K, 200 bar) allow reorganization, followed by staged cooling to ambient (down to 300 K) (Ngambia et al., 2024).
• Introduction of “virtual atoms” (VAs)—massless, strongly repulsive Lennard-Jones spheres (e.g., $\sigma_V = 1.0$ or 3.0 nm, $\epsilon_V = 10^{-6}$ kJ·mol$^{-1}$)—imposes controlled microporosity. Subsequent melt–quench simulation results in the formation of pores with target volume and width distributions. Removal of VAs is followed by re-equilibration to assess permanent pore retention (Ngambia et al., 2024).
• Simulation parameters: timestep 1–2 fs, LJ cutoff 1.4–1.5 nm, PME electrostatics, thermostats (velocity-rescale, $\tau_T$ = 0.1 ps), barostats (Nosé–Hoover or Berendsen, $\tau_P$ = 10 ps), and total run lengths up to 50 ns for equilibrated sampling (Wood et al., 2023, Wood et al., 27 Dec 2025).

Surface models are generated by box elongation and surface relaxation, with solvent-accessible surface area (SASA) calculated via Monte Carlo probe rolling ($r_{probe}$ = 0.18 nm for N$_2$) (Wood et al., 2023).

4. Property Computation and Model Validation

Post-equilibration, structural and functional properties are computed directly from atomistic trajectories:

• Atomic ratios: $\mathrm{H}/\mathrm{C} = N_\mathrm{H}/N_\mathrm{C}$, $\mathrm{O}/\mathrm{C} = N_\mathrm{O}/N_\mathrm{C}$
• Aromatic carbon fraction: $$f_{\mathrm{arom}} = N_\mathrm{C}^{\text{arom}}/N_\mathrm{C}^{\text{total}}$$
• Bulk and true density: $$\rho_{\text{bulk}} = m_{\text{tot}} / V_{\text{box}}$$, $\rho_{\text{true}} = \rho_{\text{bulk}} / (1-\phi)$ with porosity $\phi$ from helium probe-insertion (Ngambia et al., 2024)
• Pore-size distributions: maximal sphere fitting (MoloVol, Zeo++), differential volume histograms, cumulative pore volume up to $D = 2$ nm
• Surface functionality population densities: $$\rho_{\text{func}} = N_{\text{func}} / A_{\text{surface}}$$ (groups/nm$^2$) (Wood et al., 27 Dec 2025)
• Validation: simulated observables are compared against experimental targets, including XRD-derived stacking parameters ($d_{002}$, $L_c$, $L_a$), FTIR and Raman spectra, BET SSA, Boehm titration values, and functional group counts. Root-mean-square deviations (RMSD) between simulation and experiment serve as quantitative validation metrics, with accepted models exhibiting deviations within $\pm5$\% of benchmarks (Wood et al., 2023, Ngambia et al., 2024).

For analytical and continuum (packed-bed) representations, isotherm (Langmuir, Freundlich) and kinetic (pseudo–first-order rate constant $k$) parameters are fitted to batch and column adsorption experiments. These are integrated into breakthrough curve models and pore-scale transport simulations using lattice-Boltzmann codes (Pettersson et al., 2024).

5. Structure–Function Relationships and Insight into Porosity Origins

Systematic variation of BSU size, topology, and porosity engineering methods unveils the determinants of biochar microporosity and functionality:

• Densely-packed BSUs with low curvature yield nearly non-porous networks; introduction of large, non-planar BSUs (high ADS, odd-membered rings) creates ultra-micropores intrinsic to the structure (Ngambia et al., 2024).
• Virtual atom-mediated porosity allows fine control of pore width distributions and total accessible volume, mimicking the volatile loss process during pyrolysis.
• Surface functional group population influences adsorptive interactions—at low HTT, biochars retain –OH and –COOH groups enabling polar and hydrogen-bond interactions, while higher HTT favors graphitic domains and $\pi$–$\pi$ stacking (Wood et al., 27 Dec 2025).
• Correlations between pore volume, SASA, and adsorption capacity are empirically established, supporting model transferability for prediction of pollutant uptake (Wood et al., 2023, Pettersson et al., 2024).

6. Applications in Adsorption, Remediation, and Macroscale Modeling

Experimentally constrained woody biochar models support a spectrum of applications:

• Prediction and optimization of adsorption isotherms for organic pollutants and anionic herbicides, with atomistic simulations corroborating experimental sorption energies and mechanisms (aromatic stacking, polar anchoring, cation bridging) (Wood et al., 27 Dec 2025).
• Packed-bed reactive transport modeling of contaminant breakthrough in soil and engineered filtration systems, with quantification of kinetic ($k$) and equilibrium (Langmuir $q_{\max}$, $K_L$) parameters (Pettersson et al., 2024).
• Impact of particle size uniformity (monodisperse vs. polydisperse beds): narrow distributions (0.5–1.0 mm) suppress channeling and maximize adsorption efficiency; heterogeneous beds display preferential flow and lower utilization (Pettersson et al., 2024).
• Transfer of model parameters to field-scale simulation by upscaling kinetic and equilibrium constants, embedding molecularly derived descriptors in continuum codes (Pettersson et al., 2024).

7. Software, Data Resources, and Reproducibility

Model construction, simulation, and analysis are implemented using:

• GROMACS (versions 2021/2022.3) with OPLS-AA force field and automated parameter assignment (LigParGen, PolyParGen), geometry packing (Packmol), and porosity calculation (MoloVol, in-house scripts) (Wood et al., 2023, Ngambia et al., 2024).
• Pore-scale and hydrodynamics: Palabos or equivalent lattice-Boltzmann solvers, with explicit boundary reaction modeling (Pettersson et al., 2024).
• Model libraries, input decks, and example files are openly accessible: https://github.com/Erastova-group/Biochar_MolecularModels (Wood et al., 2023).
• Fitted parameter tables (e.g., $q_{\max}$, $K_L$, $k$, porosity $\epsilon$, particle diameter $d$) and validation metrics (RMSE, MAE, $R^2$) are reported to enable direct comparison and replication.

By constraining each modeling decision with experimentally derived quantities, spanning bulk, surface, porosity, and adsorption metrics, these models achieve physical fidelity and predictive utility for biochars derived from woody biomass across the spectrum of pyrolysis conditions.