The increasing worldwide concerns for sustainability push chemical engineers to perfect the industrial processes according to stringent paradigms. Circular economy, green chemistry, intensification, clean production, and integration are ubiquitous keywords of the current process studies (Clark et al., 2016; Avraamidou et al., 2020; Tula et al., 2020). In this scenario, waste biomasses and biorefinery play a central role to meet both the demands of the economy of scale and the increasing environmental solicitudes (Larragoiti-Kuri et al., 2017; Sherwood, 2020; Ubando et al., 2020). Bio-waste feedstocks could integrate efficiently into the chemical supply chain at the level of medium-scale chemical plants (Guo et al., 2019). A further benefit is that, unlike other renewable energy sources, biomass conversion into heat, electricity, and fuels could be an on-demand process (Murele et al., 2020). In many cases, however, the drying of wet bio-residues and waste is an energy sink that negatively affects the overall efficiency.

Hydrothermal conversions, i.e., biomass processing in hot compressed water, bypass this limit, and generate fuels and chemicals with lesser energy consumption than other thermal conversions (Antero et al., 2020). Integrated processes for the biomass-to-energy chain currently include hydrothermal reactors (Lee et al., 2019). Among the hydrothermal treatments, carbonization (HTC) occurs at the mildest operating conditions, subcritical temperature range, 450–620 K, and autogenous pressure. HTC is advisable for non-energy conversions of mixed wastes (Antero et al., 2020; Zhan et al., 2020). The literature debates on specialized applications of the main product, the solid hydrochar (Kruse and Dahmen, 2018). HTC has now reached a level of maturity that allows researchers to develop process considerations and classify plants (Ischia and Fiori, 2020).

To comply with the state-of-art of research activity, industrial HTC reactors should treat various materials, integrate them into other biorefinery processes, and maximize the yield of valuable products (Usman et al., 2019).

A prerequisite for the design of optimal reaction conditions is the availability of numerous experimental data and reliable HTC kinetic models. Researchers face the challenge of developing a model valid for different feedstocks, having relatively few kinetic data (Heidari et al., 2018). The scarcity and heterogeneity of data motivate to increase the range of investigative tools. Clear examples of this course are prediction techniques such as non-linear random forest models (Li et al., 2020), the design of experiments using surface response techniques (Román et al., 2020), and the assessment of models using high-pressure differential scanning calorimetry (Pecchi et al., 2020). The studies are expanding, aiming to bring the HTC modeling to maturity, as occurs for other biomass thermochemical conversions, such as pyrolysis and gasification (Weber et al., 2017; Safarian et al., 2019). In the authors' opinion, stochastic techniques could contribute effectively to perfect kinetic models and analyze experimental data. This claim is undoubtedly valid for other biomass conversion processes (Dhaundiyal et al., 2019; Terrell et al., 2020), and the HTC kinetic studies should benefit from a stochastic view inside the reaction as well. Bayesian and Markov-chain methods applied to chemical engineering show a mature state-of-art, as demonstrated by textbooks and specialized papers (Beers, 2006; Shields et al., 2021). This introductory paper aims to bring Bayesian specialists' attention to the HTC modeling and stimulate researchers working on the hydrothermal conversion of biomass to consider stochastic techniques as an additional tool.

A previous study introduced probability as the possible underlying law that steers the time-course of HTC reactions network (Gallifuoco and Di Giacomo, 2018). That paper showed how to use proper cumulative frequency distributions (CFD) and probability density functions (PDF) for describing the dynamics of solid and liquid phase transformations. A more in-depth investigation proved that several HTC kinetic mechanisms, widely used in the literature, could be modeled as Markov-chain processes (Gallifuoco, 2019). Another study enlarged the adoption of CFDs and proposed their use as a general tool for supporting HTC modeling (Gallifuoco et al., 2020). The successful accordance between statistic calculations and experimental data from different residual biomasses warrants to persevere using the stochastic approach. In this way, HTC studies could take advantage of the previous knowledge gained in the statistical analysis of other chemical engineering systems, particularly chemical reaction engineering (Erban and Chapman, 2019). The present paper introduces the novelty of the Bayesian approach and Markov-chain Monte Carlo techniques (MCMC) in the HTC studies. The aim is to enlarge the panoply of methods commonly used for studying the HTC process. The paper shows the practicality of stochastic techniques analyzing both literature data and experimental results obtained on purpose.

Table 1 summarizes the datasets over which the programs run and Table 2 details the results of the experiments of this study.

The datasets encompass a wide range of operational conditions and serve as a stress-test for assessing the routines' reliability. Substrates are representative of herbaceous biomass [miscanthus (Smith and Ross, 2019)], model carbohydrates [cellulose (Paksung et al., 2020)], lignocellulosic materials (fir) and agro-food industry scraps (potato starch). Temperature ranges from 200 to 250°C, reaction duration from 2 to 24 h, the solid-liquid ratio from 1:3.5 to 1:14. The programs analyzed the time course of significant properties of the hydrochar, e.g., mass yield, volatile matter, total carbon content, and fixed carbon. The liquid phase electrical conductivity was also studied as it previously proved to be a convenient lumped parameter for monitoring the reaction progress (Gallifuoco et al., 2018).

The misfit functions used two different model equations. For the data referring to the solid phase, the following logistic equation mostly proved to give the best fitting performances (Gallifuoco and Di Giacomo, 2018; Gallifuoco, 2019):

The time-course of liquid phase conductivity follows a series mechanism given by two first-order steps, formation plus depletion (Gallifuoco et al., 2018). Accordingly, the proper equation is:

In both Equations (11, 12), t is the reaction time and P1…P4 the estimand parameters. In Equation (11), P3 is the time scale, P4 a shape factor, and P1 and P2 the final and the initial values, respectively. Equation (12) has two different time scales, P2 and P4. The initial value is P3, the final one P1. According to data of Table 1, the software performed 22 fittings of 4 parameters each, which sum up to 88 different estimates.

The use of model equations is an entry-level problem for illustrating the features of program (A). More advanced applications are possible, such as the inverse problem of fitting parameters to the system of ordinary differential equations coming from a hypothesized reaction network.

The misfit functions associated with Equations (9, 10) lie in the 5-dimensional space, and consequently, complete visualization of the random walks is not possible. Nevertheless, plots of any couple of parameters could illustrate the essential features. Figure 4 shows typical examples of the results and serves well in the discussion.

The diagrams trace the 10 random walks from the start-point (□) to the end (♢). By way of comparison, the figure also locates the non-linear least square estimate (⋆). By inspection of the left-hand diagram, which refers to dataset 12 of Table 1, one could observe that the walkers, after the initial wandering in far, not significant regions, converge on the expected target. The right-hand diagram (dataset 4) shows more straight tracks that flow into a specific portion of the parameter plane, crowded enough to hide the endpoints. The numerical outputs show that all the walkers hit the same point, coincident with the least square estimate. The intersection of the two dashed segments locates the coordinates of the final estimate (0.647, 0.933). The trajectories in the right-hand diagram signal that the valley surrounding the minimum has steep walls. Figure 5 shows examples of the satisfactory accordance of predictions with experimental data. Figure 6 reports a selection of data vs. time plots and the fitting lines connected to the parameters' endpoint estimates.

Other plots, not reported here to avoid crowding the paper, give similar results, so the discussion of Figure 6 assumes the character of generality. The adjusted R²s, shown next to the respective correlation curves, prove that the software works satisfactorily with both the model equation, whatever the HTC conditions were. The left diagram refers to Miscanthus at two reaction temperature and reports the time courses of hydrochar fixed carbon and volatile matter content. The right diagram refers to starch and records the evolution of liquid phase electrical conductivity at three different liquid/solid weight ratios.

Comparing the obtained parameter values with those coming from the traditional non-linear least-squares method (Levemberg-Marquardt algorithm) allows testing the reliability of the MCMC technique. Figure 7 summarizes the results obtained.

The diagram shows the relative percentage deviation between the value estimated by the program (MCMC) and that obtained with least squares (NLS) for each of the 88 parameters. Most of the points align with unity (full line), indicating substantial equality between the two estimates. Six values out of the 88 are within a ±2% deviation (dashed lines). The remaining 8 points (full symbol) deserve further discussion. The numbers to the right of symbols label the dataset to which the parameter belongs. So, three out of four parameters of dataset 1 differ significantly from the respective values obtained with the least-squares method.

Interestingly, in this case, the MCMC method found a global minimum, whose misfit function assumes a value lower than that of NLS. This circumstance also occurs for dataset 3, although to a lesser extent. Datasets 2 and 20 represent opposite cases, in which NLS performs better than MCMC. The increase of iterations above 30,000 did not improve the result, and one could conclude that these two cases deserve further investigation. A more precise diagnosis would require increasing the number of walkers to reduce the effect of individual deviations. For this study, however, one could conclude that the software passed the reliability test entirely. It can detect global minima that elude the NLS method.

NLS methods are fast algorithms for identifying the misfit function's global minimum with models of low dimension. With the increase in model dimensions, these techniques become inefficient and vulnerable to finding a local minimum closest to the starting point rather than the desired global minimum. To overcome this limit, one could repeat the optimization with various starting points and monitor if they all converge to the same solution. This stratagem is prone to become computationally expensive and troublesome as one complies with many parameters. Although this study analyzes few-parameters models, in some cases, NLS failed to find the global minimum. Changing the routine's conditioning and exploring up to 100 different starting points did not improve the performance. One could expect that the advantage of MCMC over NLS becomes even more evident with multi-parameter models. Conversely, the increase of the model complexity could make the computational cost of MCMC higher than that of traditional least-square methods. Hence, to deal with more detailed reaction schemes, one should resort to more sophisticated MCMC algorithms, well-established in the literature. It is worth overstepping the scope of this introductory study and prosecute the research in this way.

Finally, dataset 19 shows the value of one of its parameters which differs sensibly (+8.71%) from the corresponding one estimated by NLS, although the misfit function's value is substantially the same with the two methods. This last result warns not to accept any regression method ipso facto, especially when the estimated parameters are critical variables for the subsequent process design. The Bayesian paradigm considers parameters as random variables whose distributions are updated by the knowledge of experimental data. The steps after the burn-in period sample repeatedly from these distributions, and this furnishes valuable information on each parameter. Table 3 sums up all the results.

One can observe that all the estimands deviate from the average value within a very narrow range. The satisfactory regressions of experimental data (adjusted R-square factors appear in the last column) reinforce the precision of estimates. Figure 7 visualizes the near-convergence samplings of four selected parameters and serves as an example of the general results.

For reference, each histogram shows the normal probability distribution. The probit analysis, which uses the cumulative probabilities to test the normality of the distribution, is superimposed to the histograms. The more the points align on a straight line, the more the observed distribution approximates the normal one. In Figure 7, the most significant deviations appear in C and D, where the distributions are slightly left-skewed. More in-depth elaborations are possible, although beyond the scope of this introductory study.

Program (B) run on reaction networks reported in Table 4.

All the schemes make use of mass-action kinetics between compartments and respect essential literature findings. As the amount of gaseous phase generated by the low-temperature hydrothermal reactions negligible, no compartment for gas species appear. Hydrochar formation is a two time-scale process, and accordingly, the networks involve two distinct pseudo-kinetic constants, k₁ and k₂. Moieties released in the liquid by primary hydrochar formation contribute to build-up secondary hydrochar, and the schemes entail liquid-phase compartments. As the reaction proceeds, the solid yield (recovered solid/biomass) should decrease, and hydrochar energy density should increase.

In the following, B stands for biomass, HC1 for primary hydrochar, HC2 for secondary hydrochar, and L1 an L2, respectively, for the corresponding liquid-phases substances. Scheme 1 is the simplest one. This naïve model, a test bench for assessing the reliability of program (B), conceives the HTC as a first-order two-step process and disregards the dynamics of liquid phases species. The related system of three ordinary differential equation has a straightforward solution, reported in Table 4. Figure 8 illustrates the results of the tests.

Diagram A traces 10 simulations of the dynamics of the secondary hydrochar, each originating from the same initial condition (t = 0, HC1 = HC2 = 0, B = 100). The random trajectories follow independent paths, relatively different from each other, but matching on average the exact solution (dashed line) reported in Table 4. Expectably, the more the realizations performed, the higher the precision of the results, as demonstrated by the other diagrams of the Figure. Each of them reports the averages across all realizations for the three species and the corresponding exact solutions. The number of realizations performed to simulate the system enlarges moving clockwise from the upper right diagram to the lower left one. One can easily verify that the precision of the simulations increases, as proved by the reported adjusted R-squares of the worst fitting among the three of each diagram. One could observe that 10 realizations suffice for obtaining a satisfactory average. The remaining diagrams demonstrate how the fluctuations reduce if one performs a higher number of simulations. Diagram D displays the simulation coming from 10,000 realizations, with just five seconds of computing time. It appears that program (B) correctly behaves, since the stochastic simulation tends to coincide with the exact solution in the limit of infinite realizations, and the computing task for attaining a satisfactory approximation is affordable.

The remaining schemes are worth considering as tentative, more detailed descriptions of the HTC reactions. The literature remarks on the solid-phase autocatalytic behavior of HTC (Brown, 1997; Paksung et al., 2020), and the first reaction of schemes 2, 3, and 6 accounts for this. For activating the process, a certain amount of HC1 should be present as the reactions start (time zero, reactor heated up to the setpoint temperature). Experimental evidence confirms that a partial biomass transformation occurs during the reactor warmup. The more the transient lasts, the higher the extent of modification. For example, in datasets 1–4, the per cent fixed carbon of native biomass is 12.3, those of time zero are 16.0 and 27.9, respectively, for 200–250°C (Smith and Ross, 2019). The finding is generally evident for lignocellulosic materials, as confirmed by the experiments of this study. Non-lignocellulosic biomass, such as starchy materials, could undergo the entire first step in the warmup. In eventuality, solid yields could increase with the reaction time. A cautious choice of the initial conditions could allow the simulation run in agreement with the experimental observations. The second reactions of Table 4 account for the formation of secondary hydrochar via interactions between liquid and solid phases. The schemes envisage the progressive reduction of solid-phase species, except for network 2 (solid balance of zero), which could adapt to cases of increasing solid yield. The improvement of these illustrative networks is straightforward, e.g., by considering different reactions. The present basic form serves the scope of this paper and gives remarkable results, easily comparable with the experimental data to steer the selection of the proper scheme. Figure 9 illustrates typical results.

The simulations run on equal values of the parameters (k₁ = 20, k₂ = 1, initial HC1 = 10, initial L1 = 10). Diagram A reports the fraction of hydrochars in the solid phase as a function of the normalized reaction time (end-time equal to one). Each network displays recognizable dynamics, and this helps to link the proper model to the experimental data. Diagram B gives the corresponding distribution of the liquid to solid components. The differences between schemes persist, except for models 3 and 4, intrinsically structured to produce a constant ratio of liquid to solid products. Liquid-phase analyses, both of key components and lumped properties (Gallifuoco et al., 2018), could allow evidencing the best-fitting model. Diagram C and D show the relationships between liquid and solid phases retrievable at any time from the reactor. Definite patterns appear for all schemes in that the solid recovery decreases monotonously with the accumulation of liquid product (C), and this last gives clear trace of the conversion of biomass to hydrochar (D).

Overall, the results of Figure 9 stimulate further elaborations for comparing previsions with experiments. A detailed linking of stochastic simulations with experiments allow selecting the proper reaction network. Although the schemes of this study consider a restrict number of compartments, there are sufficient to illustrate the procedure. Figure 10 shows some examples of how to bring valuable information out of the results. Diagram A reports the fixed carbon of dataset 2 vs. the corresponding simulated property for four out of the six reaction schemes. The authors (Smith and Ross, 2019) reported the fixed carbon of the native biomass (12.3%) of the time-zero solid (16.0%), and that measured after 24 h of reaction (31.1%), reasonably due to the complete conversion. One could speculate that the fixed carbon measured at intermediate times should be due to the weighted contributes of the biomass not yet reacted and the hydrochar already produced.

Similarly, the x-axis values weigh the amounts of B and (HC1 + HC2) with the experimental data. The diagram shows that the simulations match the experiments unambiguously. By the way, data fit very well third-degree polynomials (lines). Further investigation could disclose the essence of this finding and assess the possible theoretical foundations. Diagram B is analougus and refers to schemes 2 and 4. Although in these last cases, the polynomial fitting failed, ordered patterns appear connecting simulations and data unambiguously, as shown by the lines connecting points. Part C reports the direct relationship between datasets 20 and 21 and the simulated total hydrochar of reaction scheme 3. The left y-axis reports the liquid phase electrical conductivity, the right y-axis the solid yield. The solid line connecting the yields is an exponential which fits satisfactory the data (R² = 0.99919). The dashed line is a linear correlation (R² = 0.99876) of the conductivities, excluding the point on the upper right corner. Diagram C shows satisfactory matching trends and demonstrates that scheme 3 is predictive of the datasets analyzed. Finally, diagram D shows a good coupling of the amount of liquid recovered after the reaction and the corresponding compartments. Data refer to datasets 11, 12, and 13 as matched to model 6.

Figure 10 shows some of the many ways of linking experiments with simulations. One could envisage further fruitful implementations for steering the process of model assessments, obtaining at the same time valuable feedbacks on which part of the experimental investigation to strengthen. A more in-depth analysis requires identifying and quantifying liquid-phase key chemicals to match with the compartments of models. Experimental data on the dynamics of such compounds appear in the HTC literature seldom. As available, this evidence could contribute to more accurate validations.

The state-of-art on HTC kinetic modeling shows the success of both lumping and detailed descriptions of reaction networks. Models that describe the overall biomass conversion catch the autocatalytic progression of the condensed phase transformation satisfactorily (Pecchi et al., 2020). Mechanistic descriptions, involving up to a 10th of reactions and chemically defined liquid-phase products, result in a broader knowledge of the kinetic constants (Jung et al., 2020). This paper demonstrates that compartmental modeling, a sort of intermediate approach, could find room in the HTC kinetic studies. The research is proceeding this way.

The study demonstrates the usefulness of Bayesian statistics and Monte Carlo methods for studying biomass hydrothermal carbonization kinetics. Stochastic simulation of HTC reactions is a flexible tool for testing hypothesized networks and could improve the knowledge of the mechanism of biomass conversion. The approach could face the limit of computational expensiveness when extended to a thorough description of the reaction kinetics. The matter deserves future researches. Despite the basicity of the routines, the results are satisfactory. The estimation of the parameters furnishes regression coefficients as high as 0.9999 and detects the global minimum of the space parameters for all the datasets. The test of possible reaction schemes is straightforward. The upgrading to more sophisticated and efficient algorithms is clear-cut. It could exploit the vast library of software, readily available and long used in other areas of kinetics applied to chemical engineering. The proposed method has the potential to guide the selection of the correct kinetic model. It can flexibly simulate the dynamics of any experimentally measured property in both the solid and liquid phases.

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

AG conceived the work, developed the modelling, and the computers routines. AP and LT contributed equally in carrying out the experiments. All authors contributed to the article and approved the submitted version.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.