Hydrothermal Carbonization of Organic Waste and Biomass: A Review on Process, Reactor, and Plant Modeling

Hydrothermal carbonization (HTC) is an emerging path to give a new life to organic waste and residual biomass. Fulfilling the principles of the circular economy, through HTC “unpleasant” organics can be transformed into useful materials and possibly energy carriers. The potential applications of HTC are tremendous and the recent literature is full of investigations. In this context, models capable to predict, simulate and optimize the HTC process, reactors, and plants are engineering tools that can significantly shift HTC research towards innovation by boosting the development of novel enterprises based on HTC technology. This review paper addresses such key-issue: where do we stand regarding the development of these tools? The literature presents many and simplified models to describe the reaction kinetics, some dealing with the process simulation, while few focused on the heart of an HTC system, the reactor. Statistical investigations and some life cycle assessment analyses also appear in the current state of the art. This work examines and analyzes these predicting tools, highlighting their potentialities and limits. Overall, the current models suffer from many aspects, from the lack of data to the intrinsic complexity of HTC reactions and HTC systems. Therefore, the emphasis is given to what is still necessary to make the HTC process duly simulated and therefore implementable on an industrial scale with sufficient predictive margins.


Introduction
The conversion of biomass into bio-based products and bioenergy is an important innovative aspect to pursue the circular economy principles, aimed to increase the amount of renewable sources and to reduce the consumption of raw material and energy [1]. The usage of second-generation biomass as starting material for the valorization stage is a linkage between circular economy and waste management. Indeed, this category of biomass does not interfere with food production or land usage and, in addition, it requires a particular disposal strategy. Second-generation biomasses, which comprise agricultural and forest residues, and industrial and organic waste, are commonly heterogeneous and characterized by a high moisture content [2]. Besides, they require significant pretreatment stages for improving their handling and storage. For example, energy-consuming drying stages are often needed to make the feedstock suitable for traditional thermochemical conversion processes.
In this context, hydrothermal carbonization (HTC) could play a significant role in the valorization of organic waste and biomass residues. Even if the first concept of HTC dates back to the early twentieth century to understand the mechanism of natural coalification, only in the past 10 years HTC has acquired significant attention inside the scientific community [3]. HTC is a thermochemical process occurring in subcritical water under relatively mild conditions, with a reaction temperature commonly maintained in a range of about 180-250 °C and pressure between 10 and 50 bar. Under these conditions, liquid water behaves both as a reactant and a solvent/reaction medium for the biomass, leading to its degradation through heterogeneous chemical reactions, which include dehydration, hydrolysis, aromatization, and polymerization. The result is a solid phase rich in carbon (the so-called hydrochar), a liquid by-product, and a gas phase rich in CO 2 . The hydrochar can be used in a variety of applications like agriculture (e.g. as a solid amendment), energy production (e.g. in co-combustion with traditional fossil fuels [4]), as adsorbent [5,6], and as a base for advanced carbon materials [7]. Occurring in water, HTC is particularly suitable for the upgrading of biomass with high moisture content and a wide variety of "wet" biomasses have been investigated, like lignocellulosic biomasses, organic wastes, sewage sludge, as well as basic compounds like pure cellulose and other carbohydrates [8][9][10]. Meanwhile, HTC has been successfully applied in the context of the third and fourth generation biomass, mainly represented by algae [11][12][13]. Even if algae are particularly attractive for the production of liquid biofuel through hydrothermal liquefaction (HTL), the adoption of HTC has been also investigated. Indeed, HTC has been suggested as an intermediate step for the production of biodiesel. Actually, the hydrochar adsorbs the fatty acids formed from the hydrolysis of the lipids naturally present in the algae and, if extracted, these fatty acids can be used as precursors for the production of biodiesel [14,15].
The hydrochar, the carbonaceous solid phase resulting after the process, presents interesting properties that make it suitable for various applications. Due to dehydration and decarboxylation reactions, it presents a higher carbon content than the biomass, with a contextual improved heating content (HHV) and lower O/C and H/C ratios [16]. This makes it similar to conventional solid fuels. It also presents an improved hydrophobicity [17]. Regarding surface properties, these are similar to those of the starting material: during the process, products resulting from polymerization deposit on the surface, blocking the pores [18]. Nevertheless, modifying the hydrochar through activation can provide properties suitable for contaminants remediation [19].
A field of rising interest is the integration of HTC with other valorization techniques. Among these, combining HTC with anaerobic digestion (AD) for the production of biogas represents a promising strategy for valorization and usage of the HTC process water (the liquor), otherwise a byproduct of the process [10,13]. Several integration strategies have been investigated for different feedstock at lab scale [13,[20][21][22][23][24], showing that the combination is effective both to energetically sustain HTC and to reduce the environmental impact of the liquor disposal [10]. Moreover, some studies investigated the combination of HTC with conventional thermochemical processes, like gasification and pyrolysis, in which HTC represents the pretreatment stage to enhance the chemical and physical properties of the raw biomass [16]. In particular, the usage of the hydrochar for gasification can improve the quality of the syngas and reduce the tar formation [25][26][27], while pyro-hydrochar exhibits improved surface area and durability than the normal biochar [25,28]. Recently, renewable and alternative technological configurations have been also proposed, like the coupling between an HTC reactor and a solar concentrator [29,30].
Until now, significant experimental activity has been done to investigate HTC reaction mechanisms and the effects of process parameters (like reaction temperature, residence time, and solid load) on the final products, whose physical and chemical properties have been extensively characterized according to the desired application. Even if some pilot-scale and commercial continuous reactors are present worldwide, most of the research has been carried out on 1 3 lab-or bench-scale reactors [3], typically batch reactors: certainly, more efforts are required to develop continuous reactors so to move towards the industrialization of HTC. Challenges must be faced to spread this technology on a commercial scale, and the development of reliable predictive models is considered by the research community as one of these challenges.
Despite the massive experimental activity performed to characterize the whole process, the modeling and simulation of HTC are lagging. Indeed, the development of validated models and tools is crucial for the optimization and design of the process, as well as for the understanding of the HTC mechanism dynamics. Through modeling, virtual scenario can be settled and the variations of process parameters investigated, so to provide useful information to obtain a product with the required characteristics. This is essential to increase the HTC competitiveness and attractiveness with respect to traditional techniques [31]. Commercial parties need to easily explore the possibility of HTC, without investing significant resources [32]. Therefore, the scarcity of reliable models in the current state of the art for HTC makes modeling a challenge that should be faced in our immediate future for the diffusion of the HTC technology [33]. Moreover, the current literature mainly focuses on the modeling of properties related to biofuels applications (like carbon content and HHV), limiting the spread of the other applications.
Until now, modeling and simulation of HTC have been applied on different system levels. These pass from the description of the process itself through reaction kinetics models (which often use pseudo-first-order reactions) and statistical correlations, to the description of the reactor and the overall plant. Overall, the current reaction kinetics and statistical approaches are highly affected by the experimental data through which they are calibrated. Indeed, they propose relations among the properties of the resulting phases and the HTC operating conditions. Meanwhile, techno-economic analyses are often associated with the large-scale models and some life cycle assessment studies have been also performed. While many reviews present the general knowledge of HTC with focus on the chemistry of the process, physical and chemical properties of the products, and applications [3,16,34,35], few works address the current state of the art on HTC modeling and simulation, without however presenting a complete overview of all aspects related to modeling in the HTC field.
This work provides a review of the main aspects related to HTC modeling, highlighting strengths and aspects that should be improved to reach an appropriate knowledge level. The different models applied to the different HTC subsystems (process, reactor, and plant) are presented. We hope that this work could guide researchers to the development of new and more reliable models that could drive the shift of HTC from a laboratory to an industrial scale.

Overview on HTC Process Modeling
HTC models and simulations that have been developed so far can be categorized according to the subsystem of analysis. Results and data of the previous subsystem are necessary for the definition of the subsequent subsystem (see below). The following categories (or subsystems), characterized by progressive broader boundaries, have been identified: (1) Process i.e. reaction kinetics and statistical models.
They focus on the reaction mechanisms and the correlations between process parameters (i.e. temperature, reaction time, and biomass load) and certain observed output variables. They are usually used to predict properties (like hydrochar yield, carbon content, and energy content) and are valid for a narrow range of experimental conditions and feedstock. Kinetics models are often based on a simplified reaction scheme and assume lumped first-order reactions following the Arrhenius law. Linear and non-linear correlations are widely adopted to statistically correlate HTC variables. (2) Reactor These models are intended to predict variations in space and time considering a multitude of aspects that include reaction kinetics, mass and heat transport phenomena, as well as the geometry of the reactor. Validated computational models developed in proper software (like Comsol Multiphysics ® ) belong to this model category. Until now, few studies cover this category and the lack of precise kinetics and thermodynamic (e.g. the HTC heat of reaction) data negatively affects their implementation and diffusion. (3) Plant They model the entire HTC process, considering the reactor and the various plant equipment (like pumps, heat exchangers, boiler, filters, and dryer). They are often coupled with an economic evaluation and use commercial software (like Aspen Plus ® and VMGSim ® ) or software codes ad hoc developed in proper computer programming languages. The current plant simulation models lack of precise data on the reactor and presume that lab-scale data can be representative at a much larger plant-scale. (4) Life Cycle Assessment (LCA) They evaluate the environmental impacts accounting for all the stages of the process, providing a holistic view of the process by applying a "cradle to grave" approach. Some examples of LCA applied to HTC are present in the literature, but more large-scale data are necessary to make them more reliable.

Reaction Kinetics Models
The modeling of the reaction kinetics has a key role in the industrialization of the HTC technology. Indeed, a model capable to predict the output on the base of input parameters is a fundamental requirement to design and optimize a process that could be reliable, competitive, and attractive when compared with other technologies. Mainly due to the complexity of the reactions involved and the difficulty of developing a general reaction mechanism, current kinetics models suffer from strong simplifying hypothesis and dependence on the experimental conditions under which they are calibrated and/or validated. Nevertheless, significant improvements have been done in the last years. In the current state of the art, many reaction kinetics models adopt simplified lumped schemes based on the Arrhenius equation, which describe the process through pseudo-first-or higher-order reactions. In addition, simplified "all-inclusive" models, like the severity and coalification models, were developed to predict hydrochar properties using only one parameter condensing the more affecting factors. Recently, new approaches propose the description and corresponding modeling of HTC through probabilistic laws.

General Concept
Until now, the kinetics of the HTC of biomass has been primarily modeled adopting simplified Arrhenius kinetics schemes, in which the overall mechanism is simplified in a network of discrete reactions comprised of few elementary steps. Few attempts consider other approaches, like the stochastic or computational one in a three-dimensional space. Authors developed models on a wide variety of feedstock, ranging from pure components to complex biomasses, hypothesizing possible reaction pathways. In particular, the main components of lignocellulosic biomass (cellulose, hemicellulose, lignin, and extractives), have been investigated both individually and when enclosed into biomass. Generally, even on the pure compounds, it is observed a great variability in the results.
A standard procedure to develop a kinetics model for HTC consists in: (1) collecting data through HTC experimental campaigns at different key variables (i.e. temperature, time, and initial biomass load); (2) hypothesizing a reaction mechanism; (3) writing one or more mass balance equations with Arrhenius-based kinetics relations with some "free parameters" (e.g. activation energies); (4) solving the equations by applying numerical methods and fitting the experimental data in order to obtain the best values of the free parameters.
The kinetics modeling of HTC is a demanding task. HTC occurs through a complex network of reactions, whose mechanisms are still not fully understood. The intrinsic complexity and variety of biomass do not facilitate the modeling. Indeed, each biomass component has its own structural and chemical properties, and interacts with the other components/chemical species during the process involving heterogeneous reactions. Therefore, a detailed and comprehensive reaction mechanism, comprised of reaction kinetics, thermodynamics, and transport phenomena, is still lacking [36].
To date, several reaction mechanisms have been identified and experimentally observed. During HTC, biomass undergoes several reactions that include hydrolysis, dehydration, decarboxylation, polymerization, and aromatization [37]. In particular, the solid residue resulting from HTC (the hydrochar) derives from two different reaction pathways, which lead to two distinct solid fractions. One is called "primary char" and results from a solid-solid conversion. Having a similar morphology and structure of the parent biomass, it represents somehow its unconverted fraction [38]. The second is known as "secondary char" and forms in the aqueous phase. Secondary char results from the condensation and re-polymerization of dissolved intermediates in the liquid phase [39]. Therefore, to model this complex mechanism, a sequence of steps of reaction is usually established. Schematic sketches of the simplified and detailed reaction mechanism are shown in Fig. 1a and b, respectively.
After having assumed a certain reaction pathway, mass balances have to be written for the various species involved in the HTC reactions: the biomass, the hydrochar(s), the gaseous molecules and the molecules present in the liquid phase. Some species can be grouped in "lumped-components" in order to simplify such a complex reacting system.
The approach adopted to model the HTC kinetics derives from typical homogenous phase reacting systems. This is a huge approximation of the reality: HTC occurs in a complex heterogeneous system, characterized by mass transfer phenomena and where reactions take place in the solid, liquid, and gaseous phases. The approximation of a homogenous system is to some extent justified by the much larger amount of water than biomass (common biomass to water ratio are around 0.10-0.15). Moreover, a constant volume system is considered. Despite the approximations, this approach is particularly useful for its simplicity and well suits an easy collection of experimental data. Therefore, a relation often found in the literature because of its simplicity and its direct link to the experimental practice is represented by Eq. (1).

3
where r represents the reaction rate (mol L −1 s −1 ), in particular the degradation rate of a certain feedstock F whose concentration is C F (mol L −1 ), and t is time (s). Through experimental activity, the variation in the concentration of the selected feedstock vs. time provides a direct indication of the reaction rate through which the same feedstock reacts. Actually, Eq. (1) is valid only when a series of conditions are fulfilled: (1) the reactor is a closed system (no matter in, no matter out during the reaction process); (2) reactor perfectly mixed; (3) reacting mixture occupying a constant volume. These three conditions can be actually achieved in a lab HTC reactor. The mass balance of species F then turns into Eq.
where the transient term on the right of the equal sign equals the generation rate (actually, the degradation rate) of F. Saying that the degradation rate of F coincides with the reaction rate contains the implicit simplification that species F reacts through a single "degradation reaction" where its stoichiometric coefficient is equal to 1.
In the case of a degradation reaction, the reaction rate can be expressed through Eq. (2).
where k is the reaction rate constant (s −1 if n = 1) and n is the reaction order. Equation (2) is strictly valid for an irreversible reaction occurring in a homogeneous phase. Moreover, according to Eq. (2), the reaction rate depends only on the concentration of species F, while in the general case the concentration of other reactants can affect r. In the case of HTC, it is reasonable to consider that the concentration of water C H20 is really high and constant, so that C H20 can be avoided in the reaction rate expression. Then, the reaction rate constant k is normally expressed through the Arrhenius Eq. (3).
where k 0 is the pre-exponential factor (s −1 ), E a the activation energy (J mol −1 ), R the universal gas constant (8.314 J mol −1 K −1 ), and T the temperature (K). The Arrhenius equation follows from the collision theory developed for reacting fluid systems: its application to reacting solid system, in particular to solids undergoing degradation reactions, is not so straightforward.
Through modeling, values for the Arrhenius kinetic parameters are estimated. The common approach consists of: (1) estimation of the rate constants k i of each reaction by applying specific numerical methods starting from experimental data of concentration variations; (2) construction of the Arrhenius plot (ln k i versus 1/T) for each i-reaction; (3) computation of E a and k 0 from the Arrhenius plot.
As previously shortly mentioned, a step forward in respect to the very basic description offered by Eqs. (1) and (2) is considering the whole set of mass balances of all the species (or lumped-components) involved in the HTC reactions: reactants, intermediates and products. This would result in a system of differential equations where time t is the only independent variable. The generation term of the various species appears in the mass balances as a function of the various reaction rates. For example, assuming that all the reactions are first order apart from the secondary char formation, Fig. 1a can be described by the set of Eqs. 4-8: where C B , C L , C HC1 , C HC2 , and C G are the molar concentrations in mol/m 3 of the i-species (like atomic carbon) respectively in the biomass (B), liquor (L), primary char (HC1), secondary char (HC2), and gas phase (G) derived from the biomass and the HTC liquor [41]. In addition to the description of HTC through an Arrhenius behavior and a lumped approach, which are however functional and hardly imaginable to be removable, the current modeling presents other limits like the following.
(1) A common approach consists in considering pseudofirst-order reactions (n = 1) which, although its constructive simplicity, hides the effective behavior of the system and omits some important aspects. Indeed, secondary char formation is highly favored at high concentrations of the intermediate products, meaning that the kinetics model should be likely described by a higher reaction order (n > 1). This issue is pointed out also by Jung et al. [43], who demonstrated that a first-order reaction is unsatisfactory for describing the growth of hydrochar from fructose and that computational results are affected by the initial fructose concentration. The formation of secondary char is thought to occur through typical second-order polymerization reactions. (2) Due to the complexity of the mathematical handling, authors often consider the temperature independent of time, so that Eqs. (1), (2) and (3) can be solved by separation of variables. This is a strong simplification. Indeed, a distinct heat-up phase and a constant temperature period characterize batch reactors, which have been the object of most kinetics studies [44,45]. During the heating period, a significant biomass conversion occurs and the transient phase could be neglected only in case of rapid heating times, which could occur only in batch micro-reactors [41,46]. (3) Common kinetics models do not consider mass and heat transport phenomena, which instead can affect the overall reaction rate. Indeed, during HTC, water and products have to enter and leave the biomass pores. Therefore, mass transfer through pores can represent a limiting factor for the HTC process kinetics. (4) Models are commonly developed on the base of experimental activity performed on a certain feedstock, which affects the chemistry of the reactions involved. Therefore, their application is often limited to the type of biomass and operating conditions used.

Pseudo-first-order reaction kinetics models
As abovementioned, most of the developed models are pseudo-first-order (n = 1), which implies that the reaction rate depends linearly on the reactant (e.g. the biomass) concentration. Even if these models neglect higher-order reactions (like polymerization), they have the advantage of being easily numerically implemented. Table 1 sums up the Table 1 Sum of first-order (n = 1) reaction models presented in this work Feedstock Mechanism and details of the process Process conditions Arrhenius kinetic parameters R ef.

Table 1
(continued) 1 3 process conditions and the results of the main literature concerning these models.
It is important to underline that the hypothesized reaction pathway highly affects the values of the kinetic parameters: each reaction stage is characterized by its activation energy and pre-exponential factor. Therefore, the same feedstock (better: feedstock degradation reaction) can present different values of kinetic parameters in dependence of the reaction scheme adopted. Moreover, the models currently present in the literature consider ashes as inert during the process, neglecting their slight decomposition/evolution during HTC and their possible catalytic effect.
Cellulose has been a widely investigated feedstock applying a pseudo-first-order approach. During HTC, due to temperature and pressure conditions, glycosidic bonds that connect glucose monomers forming cellulose are broken up to form glucose monomers [47]. However, competing reactions can also be present, which complicates the overall mechanism. One of the first kinetic studies on cellulose was performed by Schwald and Bobleter [48] who, assuming a pseudo-first-order reaction into glucose and degradation products, obtained an activation energy of cotton cellulose of 129 kJ mol −1 . Besides, Yousefifar et al. [49] considered a peculiar reaction pathway, in which the gaseous products derive from both the direct conversion of soluble intermediates and volatilization of intermediates, finding out three different kinetic parameters on the base of the chemical oxygen demand (COD). Dos Santos Rocha et al. [50] developed a more complex pathway to describe the decomposition of cellulose derived from sugarcane straw, considering cellobiose, glucose, formic acid, and hydroxymethylfurfural (HMF) in the reaction mechanism. The model fits well experimental data in a temperature range between 180 and 195 °C, while it overestimates the degradation at 210 °C since it does not consider the difference between crystalline and amorphous cellulose. In the same study, the degradation of hemicellulose into xylose, arabinose, acetic acid, glucuronic acid, and furfural was also investigated.
Other pure compounds have been investigated. Rogalisnki et al. [51] characterized the degradation of different biopolymers (starch, cellulose, and protein), focusing on the formation of glucose and its subsequent degradation. In the proposed reaction scheme, the substrate undergoes two possible parallel first-order reactions to form glucose and by-products (mainly maltose, oligomers, and oligopeptides). Then, glucose further decomposes into degradation products. Results showed that the degradation kinetics of protein was disadvantaged compared to that of cellulose and starch due to the presence of peptide bonds, which are more stable under hydrothermal conditions [52].
Moreover, several authors studied the kinetics of lipids degradation [47,53,54]. These, during HTC, due to the high temperature and pressure, become more miscible in water and degrade into fatty acids [47]. Fujii et al. [53] demonstrated that monoacylglycerides degrade into fatty acids following first-order reactions, with an activation energy of 73.9 kJ mol −1 and a pre-exponential factor of 5.65 × 10 4 s −1 when subjected to an isothermal program, while the value of the two kinetic parameters become equal to, respectively, 77.5 kJ mol −1 and 1.01 × 10 5 s −1 in a constant heating program. In the same research, they tested several fatty acids (caprylic, capric, and lauric acid) under different constant heating rates, demonstrating that their degradation can be also described through a firstorder reaction mechanism [53].
First-order reactions were also used to model other pure compounds that are representative of biomass. Bicker et al. [55] studied the degradation of glucose, fructose, and pyruvaldehyde between 200 and 300 °C using zinc as catalyst. Jing and Lu [56] characterized the hydrothermal decomposition between 180 and 220 °C of glucose, which degrades through a series of first-order reactions into 5-HMF and levulinic acid. Braghiroli et al. [57] investigated the hydrothermal decomposition of tannin. Sheehan et al. [58] carried out a study on the degradation of tetra-alanine between 170 and 230 °C, considering the effects of the pH on the kinetics.
Although the individual characterization of the biomass fractions is useful to understand their separate behavior, when enclosed inside the biomass each compound may kinetically behave differently, as demonstrated for instance for cellulose [59]. Single reactions can affect the others, causing a variation of the kinetics parameters. Therefore, research on the biomass itself is crucial for the development of accurate models predictive of the biomass behavior during HTC. Generally, however, for an understandable model simplification, the biomass is assumed to be composed of certain fractions of model compounds that decompose independently of each other and can lead to the formation of common species.
In this regard, Reza et al. [60] validated a simple kinetics model for cellulose and hemicellulose through HTC tests on loblolly pine (200-260 °C). In the proposed reaction mechanism, cellulose and hemicellulose degrade through two parallel first-order reactions, while lignin is considered as inert. In particular, they assumed that cellulose degrades into a solid phase, aqueous chemicals, and gases, while hemicellulose into aqueous chemicals and gases. The reduction of the solid concentration was used to validate the model, resulting in kinetic parameters (activation energy: 73 kJ mol −1 for cellulose and 29 kJ mol −1 for hemicellulose) lower than those computed for the separate compounds. Indeed, the degradation of aqueous extractives enhances the degradation of cellulose and hemicellulose by reducing the pH, while the degradation of hemicellulose into acids catalyzes the cellulose cleavage [60,61]. Besides, Basso et al. [62] and Baratieri et al. [63] modeled the HTC reaction kinetics of grape 1 3 seeds, a constituent of grape marc residual from the winemaking process, in a two-step first-order reaction scheme, between 180 and 250 °C. The scheme assumes that biomass degrades into intermediate products that further degrade into hydrochar, while volatiles form from both intermediates and hydrochar.
The HTC kinetics of sugarcane bagasse was studied by Iryani et al. [64], who adopted a shrinking core model according to which the reaction zone progressively moves from the external to the internal part of the reacting biomass particle. Then, they assumed that the degradation can be described by two parallel first-order reactions (hydrolysis and dehydration), characterized by their proper kinetics parameters.
First-order reactions have been also applied for the modeling of sewage sludge. Danso-Boateng et al. [65] modeled the decomposition of primary sewage sludge and synthetic feces with a simple first-order model based on solid yield conversion. They obtained values in the ranges of 70.4-77.8 kJ mol −1 and 4.0 × 10 6 -1.5 × 10 7 min −1 for the activation energy and pre-exponential factor, respectively. Still regarding sewage sludge, Yin et al. [66] developed a firstorder model based on the decomposition of macromolecules and low molecular weight constituents. The reaction scheme includes two stages: a first decomposition of the biomass into its primary constituents (proteins and saccharides, present inside microorganism and the extracellular polymeric substances of the sludge) and their degradation in water into ammonia, acetic acid, and gaseous products. Besides, they included correction terms to fit the data. Results show that the decomposition into proteins and saccharides occurs with an activation energy of 51 kJ mol −1 , while into acetic acid of 29.7 kJ mol −1 .
Adopting a liquid-phase thermogravimetric technique, Mochidzuki et al. [67] investigated the HTC kinetics of biomass residues like old newspaper, rice husk, and spent malt from brewery.
A different approach was proposed by Luo et al. [68], who developed a first-order model for the degradation of water hyacinth by considering kinetics expressions commonly used for pyrolysis and combustion of solid feedstock. They adopted a non-isothermal method and, by introducing the heating rate ( = dT/dt) constant during the tests and the conversion factor , obtained Eq. (9).
By applying the Coats-Redfern differential method, they computed the kinetic parameters in two temperature ranges (150-210 and 200-280 °C). Notably, there is a considerable different behavior between the two ranges, which is attributed to the subsequent decomposition of hemicellulose (lower temperature range, E a = 147.11 kJ mol −1 and k 0 = 1.74 × 10 16 min −1 ) and cellulose (upper temperature range, E a = 90.79 kJ mol −1 and k 0 = 5.08 × 10 8 min −1 ). Using the same approach, Liu and Balasubramanian [44] performed a kinetic model study on coconut fiber and eucalyptus leaves, treated between 150 and 300 °C.

Higher-Order Reaction Kinetics Models
Higher-order reaction kinetics models (n > 1) allow considering important aspects commonly absent in first-order models. Indeed, the formation of secondary char, the waterinsoluble solid fraction favored under harsh HTC conditions [38], occurs through polymerization and polycondensation reactions. Polymerization reactions are often characterized by a second-order kinetics or, to be more generalist, their rate increases more than proportionally with respect to the increase in concentration of the reagents. Therefore, from the experimental activity and the intrinsic kinetics of the polymerization reactions, severe HTC conditions and high biomass to water ratios (i.e. high biomass concentrations) seem to increase the probability of having higher-order reactions. Overall higher-order reaction models are more sophisticated and reliable than the first-order models. The price for this precision is a more complicated numerical handling and the need of a very well-defined reaction pathway consistent with the experimental measures. Thus, several authors developed HTC kinetics models considering higher reaction orders, as summarized in Table 2. It is worth noting that the handling of higher-order reaction models requires particular attention. Indeed, the values of the concentration of the starting feedstock can affect the reaction kinetics parameters. This is particularly true for biomasses, where the concept of molar concentration is not straightforward. Authors often try to adapt the concepts valid for homogenous systems (where the molar concentration is referred to as mol/m 3 ) also to heterogeneous systems (insoluble biomass in water), and this introduces a certain subjectivity and dependence of the results on the choices made. For example, when the initial concentration is lower than one, the reaction rate is lower for models with n > 1 than for those with n = 1, which may be counterintuitive. Intuitively, a higher order reaction (n > 1) is faster than a first order reaction, which actually occurs if the initial concentration of the reagents is greater than one-therefore it becomes decisive as the concentration in biomass-water systems is defined. Generally, when n is different than 1, the pre-exponential factor k 0 will assume the unit of measure s −1 (m 3 ) n−1 mol 1−n .
Knežević et al. [69] are the first authors who stressed the need for developing a kinetic model for HTC with a reactionorder higher than one. They investigated the decomposition 1 3 rate of glucose between 250 and 350 °C. According to their proposed reaction model, glucose degrades into water-soluble products (WSS) that can further degrade into watersolvent insoluble products (WSIS) through a second-order reaction. By measuring the yield of the products and glucose degradation at different residence times, temperatures, and initial concentrations, they obtained the reaction rate for each modeled reaction. In addition, after observing that lower initial glucose concentrations increase the overall glucose degradation rate, they considered that glucose degrades into WSS according to a reaction order smaller than one. In this regard, they fitted the data using a reaction order equal to 0.9, and they computed an activation energy of glucose degradation equal to 114 kJ mol −1 , and a pre-exponential factor of 1.0 × 10 9 g 1−n s −1 L n−1 .
Later, Jatzwauck and Schumpe [70] modeled the HTC of soft rush using a simplified scheme, where the solid substrate undergoes degradation into dissolved intermediates (n = 1), which can form gaseous molecules (n = 1) or polymerize into hydrochar according to a higher-order reaction (n = 1.53). Therefore, solid-solid reactions were not considered. The concentration was defined in terms of atomic carbon. Results show that the activation energy of the first degradation step (hydrolysis, 141 kJ mol −1 ) is much higher than that of the other subsequent reactions (75.0 and 74.3 kJ mol −1 ).
A kinetics scheme that considers the pathways of primary and secondary hydrochar formation was developed by Lucian et al. [41,71]. The model predicts the carbon distribution among reaction intermediates and the products, and was validated on several biomasses: olive trimmings, grape marc, and Opuntia Ficus Indica. Even if do not performed by the authors, after a proper calibration, the model could be used to predict the elemental composition of the considered products. The reaction scheme is reported in Fig. 2, where the secondary char is produced according to a n-order reaction and derives from the polymerization and condensation of the water soluble-species in the liquid phase. Results show that the formation of primary char is the most favored reaction pathway, while that of secondary char occurs with a reaction order between 1 and 2 (the higher the HTC temperature, the higher the order). The model comprises also the heating-up phase and both modeling and experimental results point out the importance of considering into the computations such a temperature transient phase.
Another higher-order model was presented by Jung et al. [43] for the decomposition of fructose. The authors point out the necessity of introducing a higher reaction order to account for the effect of the initial fructose concentration on the hydrochar yield. In the proposed scheme, fructose decomposes into degradation products and HMF, which can polycondensate into hydrochar (in a second-order reaction) or other degradation products. The model was calibrated using gravimetric data and predicts kinetics parameters at different feedstock concentrations.
Recently, Keiller et al. [40] proposed a model based on the relative amount variation of cellulose, hemicellulose, and lignin during HTC of Australian saltbush. Through HPLC analysis of key-compounds, kinetic parameters were computed assuming several reaction orders (0.5, 1, and 2) for the degradation of cellulose, hemicellulose, and lignin. By comparing the predicted parameters with experimental results, they showed that the first-order reaction models behave the best for hemicellulose and lignin, with activation energies of 61 and 66 kJ mol −1 , respectively. Meanwhile, a 0.5-order model fits better data for cellulose degradation, leading to an activation energy of 127 kJ mol −1 and suggesting that other mechanisms could be involved during the degradation of cellulose.
A novel methodological approach was recently proposed by Pecchi et al. [72], who used a high-pressure differential scanning calorimetry to assess the kinetic parameters of sludge and digestate degradation during HTC. The method adopts the concept of the extent of reaction α, defined as the ratio between the instantaneous heat of reaction and the maximum heat released at the highest temperature. Experimental curves were used to calibrate an nth order and an autocatalytic kinetic model. Results show reactions orders of 2.68 (sludge) and 2.46 (digestate).
Overall, what emerges from the above analysis is that, even if modeling the degradation of the same biomasses, it is difficult to compare the values of the kinetics parameters among the different models and with the first-order models. Indeed, these values are highly dependent on the hypothesized reaction mechanism. Each reaction stage is characterized by its proper parameters and reflects to some extent the overall reaction pathway.

Severity and Coalification Models
Another approach of kinetic modeling is that adopted in severity and coalification models. These models assess and predict the effects of HTC parameters on hydrochar properties (like yields, chemical composition, and energy content) in a compact way. Indeed, a unique parameter with a particular form condenses the effect of the main influencing parameters (i.e. temperature and reaction time). Overall, these "all-inclusive" models have the advantage of having a relatively easy implementation and show the effect of the HTC process parameters at first glance. However, in addition to a certain dependence on the starting feedstock, they do not reflect HTC reaction mechanisms and do not include important parameters like the biomass load (often expressed as biomass to water ratio) and the initial feedstock carbon content. Among these, we can distinguish between the severity models and the coalification models.
The severity models use the so called "severity factor". This concept was originally introduced by Geniesse and Reuter [73] to model the cracking of oils, and was developed later to characterize the hydrolytic depolymerization of lignocellulosic biomass [74,75]. The severity factor approach was applied successfully also to HTC. The model proposed by Abatzoglou et al. [75] assumes a first-order reaction with an Arrhenius dependence on temperature and defines the reaction coordinate (R 0 ) as expressed in Eq. (10).
where t R is the reaction time (in minutes) and T ref the reference temperature. is a parameter depending on the gas constant (R), the activation energy (E a ), and the reference temperature, as expressed in Eq. (12).
Montané et al. [76] demonstrated that the reference conditions do not affect significantly the model and can be chosen at the middle of the data sets and T ref is usually set at 100 °C and at 14.75 [74,77]. Being an empirical parameter, value has been sometimes varied when using other feedstocks. For example, Ko et al. [78] used an equal to 4.6 when modeling the decomposition of glucose.
To graphically represent the process severity, the reaction coordinate is expressed through the logarithmic function (log R 0 ). This is referred to as severity factor (SF) and is expressed in Eq. (13): Hoekman et al. [79] used the severity factor to compare products from fast HTC (reaction time of 20-30 s) of loblolly pine with those of typical HTC batch reactors. Starting from experimental data (relevant to 39 HTC tests), Heidari et al. [80] made an attempt to predict the effect of several severity factor values (3.83, 5.01, and 6.19) on the hydrochar mass yield, higher heating value (HHV), and carbon content. With the same aim, in a work dealing with the water recycling effect during the HTC of sawdust, Heidari et al. [81] used severity factors between 4.3 and 5.3. Similarly, Basso et al. [82] reported the solid, liquid, and gaseous phase yields resulting from HTC of grape marc constituents as a function of the severity factor. For the solid (hydrochar) mass yields, the maximum absolute error between predicted and experimental values was of 7.1%, while for the liquid and gas phases yields was higher than 10% [82]. In addition, Borrero-López et al. [83] used the severity factor to model the kinetics of furfural and 5-HMF derived from the degradation of hemicellulose, cellulose, and lignin, using two alternative models (one with and one without side reactions).
When fitting the parameters as a function of R 0 , a proper function has to be chosen. Several functions have been reported in literature. For example, Suwelack et al. [84] demonstrated through statistical analysis that the solid, gaseous, and liquid yields follow the pattern expressed by Eq. (14).
Guo et al. [85] implemented the correlation between HTC output properties (hydrochar yield, carbon content, and HHV) and the severity factor by using a dose-response function as expressed in Eq. (15). By investigating the trends of the derivative of Y (representing the output property) as a function of ln(R 0 ), they observed that the maximum variation rate (correspondent to the first derivative peak) of the solid yield, carbon content, and HHV, occurs in a severity factor range of 5.8-6.4.
The coalification model was firstly developed by Ruyter [86], who performed a series of hydrothermal treatments on a variety of low-rank materials, from biomass waste to sub-bituminous coal. The model derives from the Arrhenius relation and describes the oxygen content variation of the feedstock during HTC. One of the main assumptions is that the oxygen of the feedstock is converted to CO 2 and H 2 O at a fixed ratio, regardless of pressure, pH and feedstock, while the complete conversion has been conventionally assumed to occur when the oxygen content is 6 wt%. The model in based on the relation reported in Eq. (16).
where t is the reaction time, T the reaction temperature, O feed the feedstock oxygen content, O t the oxygen content of the feedstock undergoing carbonization at reaction time t, and A, B, and C are calibrating coefficients (A = 50, B = 0.2, and C = 3500). The coalification model was used by Kruse et al. [87] to characterize the kinetics of brewer grains during HTC, in which A, B, and C were adjusted on the basis of the experimental campaign performed.
Later, Jung and Kruse [88] evaluated the ability of the coalification model to predict carbon content, oxygen content, and hydrochar yield for a variety of feedstock (29 typologies). For doing this, they linearized Eq. (16) in Eq. (17) and determined the coefficients A, B, and C by fitting experimental data relevant to every feedstock.
Thus, to establish the accuracy of the model, results were compared with those obtained applying the severity model (linearized and with the dose-response curve). Results show that the coalification model can be applied to predict carbon content, oxygen content, and hydrochar yield with reasonable low mean absolute errors. In particular, even if the best fit was achieved when every parameter was fitted individually to every feedstock, little accuracy loss arose when all the parameters were adjusted considering all the feedstock. Authors also proposed and validated a method to predict the hydrochar properties based on the lignin content of biomasses.

Other Kinetics Models
Although the majority of the research around HTC kinetics involves pseudo-first or higher order models, some new interesting approaches have been recently proposed.
In this regard, Gallifuoco [33] demonstrated that HTC can be modeled stochastically. Indeed, it is observed that important HTC properties (like hydrochar HHV and yield) follow a sigmoidal behavior with respect to the reaction time, whatever biomass is adopted [31]. In particular, they demonstrated that an HTC generic property y(t) can be described with the empirical law as for Eq. (18).
where y 0 is the initial value of the property, y ∞ is its final value, t is the process time, t h the time corresponding at half of the increase of the studied property, and p an empirical parameter. Y(t) may describe any property function of time, like H/C and O/C ratios, HHV, solid yield, and the carbon content. Equation (18) is connected to the Hill's equation, which was originally developed for biochemistry processes. In addition, Gallifuoco and Di Giacomo [31] suggested that the Hill's equation could describe HTC properties as a logical consequence of assuming that HTC is governed by probabilistic laws, which describe the average behavior of thousands of random detachments of particles from the parent biomass that simultaneously react between each other.
In another study, Gallifuoco [33] suggested that the formation of hydrochar and the evolution of its properties over time are connected to a Markov process. The study performed was supported by experimental data obtained from several biomasses and a proper mathematical handling. The Markov approach is well established in the literature for the modeling of stochastic processes in several technological areas (like polymers science, biology, astrophysics, computer networks, etc.) and is the base for the general method of the Markov chain Monte Carlo [89]. According to the Markov process, the HTC can be described in a stochastic way, in which the dynamic of the reactions depends only on the current state, and does not depend on the previous history. Clearly, further research is needed to confirm this attractive approach to modeling, which could avoid the resolution of complicated systems of non-linear equations [33].

Statistical Models
Statistical models are often used to correlate key process variables on the base of the experimental activity. Many research articles dealing with HTC contain a proper section dedicated to a statistical analysis in which the process conditions are correlated with output data (like temperature and time with the hydrochar yield, carbon content, and energy properties). A statistically model can be used both to provide a better knowledge highlighting correlations between process variables and to optimize the system [2]. These models are not directly related to the kinetics models: an in-depth modeling investigation around HTC certainly requires the integration of these tools.
In the statistical modeling process, a desired observed quantity (the dependent variable) is correlated with a process factor of interest (the independent variable) through proper statistical correlations, like linear and non-linear regressions, or best fittings derived from a design of experiment (DOE) approach [2]. To validate the model, additional experimental activity is commonly performed and the errors (e.g. standard mean deviation or root square mean error) computed.
Most statistical models highly depend on the specific experimental activity. Thus, they are a function of parameters like the type of feedstock, the HTC conditions (i.e. temperature, residence time, biomass load), the heating rate, the type of reactor used (continuous or batch), and the geometry of the reactor. This dependence limits the prediction capability of these models to a narrow range of experimental conditions and restricts their universal applicability. Nevertheless, statistical models are of great importance for several reasons. Indeed, they are useful: (1) to make predictions of specific properties (like mass yields, hydrochar carbon content, energy properties, etc.) within the original experimental range and set-up; (2) to understand the magnitude of the correlation between the different variables; (3) to individuate the type and optimal combination of process parameters that maximize the information derived from the experimental activity.

Linear and Non-linear Regression Models
A widely adopted and easily implementable approach is the multiple linear regression. This method adopts linear functions to relate HTC outputs with the input variables of interest, in which the coefficients of the linear function are estimated from the experimental data. A generic observed variable y is written as a function of the i-th independent variable x as expressed in Eq. (19).
where a i are the coefficients estimated through the linear regression and n the number of independent variables. The advantage of this method is that it bases on a simple function easy to manage and understand. However, linear modeling is limited to specific conditions and can be not sufficient to highlight important relations [90]. For example, Mumme et al. [91] linearly correlated data derived from the HTC of anaerobically digested maize silage. They used a linear function to correlate the carbon content with temperature, residence time, and pH. Even if the mean deviations computed in the validation phase are small (1.54%), the authors stress the limitations of this kind of model. Indeed, important factors like biomass concentration were not included and the correlation results from-and thus is valid in-a narrow investigated experimental range.
Multiple linear regression was adopted also by Volpe et al. [92] to characterize hydrochars from Opuntia ficusindica cladodes obtained at different carbonization conditions. In particular, they individuated that the hydrochar yield is influenced by three process parameters (temperature, reaction time, and biomass to water ratio), while the energy yield is driven by the biomass to water ratio. However, they highlight that a multivariate model is not able to correlate (19) y = a 0 + a 1 x 1 + a 2 x 2 + ⋯ + a n x n the fixed carbon and volatile matter contents and a non-linear model should better represent these data.
An extensive study on the data available in the literature was performed by Li et al. [93], who continued a previous work carried out by the same research group [90]. These papers are of particular interest because they aggregate data derived from a multitude of papers, for a total of 313 in [90] and 138 in [93], thus providing a model based on a variety of process conditions, feedstock, and hydrochar properties. In the first work, Li et al. [90] analyzed data using two statistical approaches, which are the multiple linear regression and the regression tree models. The regression tree approach consists of producing binary trees for each observed variable by splitting them into nodes by applying a recursive portioning [94]. With respect to the linear regression model, it has the advantage of not assuming any a priori relation between the variables. Results show that both the linear and the tree regression approaches can be used to fit the data. Anyway, the authors point out that, as suggested by the highly branching of the tree model, a significantly high level of interactions between the independent variables occurs and a nonlinear model is more suitable for the fitting. Moreover, the models highlight that process conditions (temperature and time) mainly affect the solid yield, while the type of feedstock mainly influences the hydrochar carbon and energy content. The regression tree model was expanded to a random forest model by Li et al. [93,95]. Indeed, it consists of a large collection of trees, and its performance is the average between the tree models considered. The studies show that the random forest models have an improved predictive capacity with respect to the linear regression and tree models. In particular, these models were used to represent all the combinations of model parameters to describe the hydrochar yield. Results show that the relation between hydrochar yield and HTC parameters (feedstock and process conditions) is non-linear.
The multiple linear regression model and regression tree developed by Li et al. [90] were adopted by Ro et al. [96] to predict the properties of hydrochars derived from animal manure.

Design of Experiments: Response Surface Methodology
The design of experiment-response surface methodology (DoE-RSM) is a statistical approach used to investigate the interaction between independent and response variables of a process. The concept implies the usage of a set of designed experiments to obtain an optimal response (the output) [80]. This approach aims to identify the optimal conditions that lead to the optimal response so to optimize the experimental activity. The variations of the responses are often represented graphically, thus improving the understanding of the connection among variables [35]. Several DoE-RSM approaches have been adopted in the HTC literature.
For example, Heilmann et al. [11] performed two-levelfactorial experiments with replicated center points to investigate the correlation between three variables (temperature, time, and solid load) for the HTC of microalgae. After developing an orthogonal factorial design, a linear regression function was used to predict the percentage of carbon recovered. The same model was presented by Heilmann et al. [97] to observe the effects of temperature, time, and solid load on the carbon content and mass yield of the hydrochar from distiller's grains.
Mäkelä et al. [98] investigated the statistical effect of temperature, residence time, and solid load for the HTC of paper sludge residue by constructing individual response models. The responses are the solid yield, the carbon content, the O/C ratio, the energy densification, and the energy yield. These were computed by applying the multiple regression equation as represented by Eq. (20).
where y is the vector of measured or calculated responses, X the design matrix of design factors, b a vector of coefficients, and e a vector of residuals. b was solved separately for each response, while the confidence of intervals (i.e. the range of certainity of the unkown parameter) by applying equations as proposed by Montgomery [99]. Results show that hydrochar properties are mainly affected by temperature, while poorly by the solid load.
Recently, Afolabi et al. [100] presented a DoE-RSM method to investigate and optimize the combined effect of time and temperature (the "factors") on the hydrochar yield and HHV (the "responses") of spent coffee grounds. They adopted a faced-center central composite factorial design to fit the two responses with the two factors. Results show that in the investigated range (180-220 °C and 1-5 h), the best combination of hydrochar yield and HHV (63.9% and 31.6 MJ kg −1 , respectively) is achieved at a temperature of 216.4 °C and residence time of 1 hour. This optimum conditions corresponds to a severity factor equal to 5.2, which belongs to the interval of values reported in Sect. 3.4. Some results of the analysis are graphically reported in Fig. 3, where the simultaneous effect of residence time and temperature on hydrochar yield is reported. The hydrochar yield at different times and temperatures forms a curve-shaped surface, in which the effect of the temperature is more significant than the residence time.
Similarly, Heidari et al. [80] used RSM to find relations between biomass constituents and HTC severity factor, considering as responses the mass yield, HHV, carbon content, and energy recovery factor. To achieve this, they performed 39 HTC experiments on a mixture of pure hemicellulose, (20) y = Xb + e cellulose, and alkali lignin. The resulting relations consist of polynomial equations in which the output response is function of the factors. For each model, the statistical significance was determined by performing an analysis of variance (ANOVA). Results show that the optimal point for the mixture in terms of solid yield and HHV (corresponding to 67% and 25 MJ kg −1 , respectively) is achieved at a severity factor of 4.5. The authors tested the model also on 10 different biomasses, demonstrating that it is reliable at low ash and extractives content.
The Central Composite Design (CCD) approach was used by several authors to implement the DoE in the RSM to determine the optimum process variables and the number of experimental runs. The CCD builds a second-order polynomial for the RSM without using a full factorial design of experiments, thus avoiding the usage of all the possible combinations levels of the experimental factors [101]. It bases on a two-level factorial design, where the total number of experiments (N) using k number of factors (like HTC temperature and residence time) is computed as the sum of 2 k , 2k axial runs, and n centre runs, as shown in Eq. (21).
The CCD method was adopted to observe the effects and optimize the variables of the HTC of different biomasses (e.g. olive stone, tomato peels, and water hyacinth) by several authors, like Álvarez-Murillo et al. [102], Sabio et al.  [100]. Copyright © 2020, Elsevier 1 3 [103], and Román et al. [104]. Authors observed the influence of temperature, residence time and biomass to water ratio (which implied a number of random experiments equal to 18). Then, the investigated target functions Y (the solid yield and HHV) were modeled using a second-order function as Eq. (22).
where A i are coefficients (computed after the fitting of the experimental data), R is the biomass to water ratio, T the temperature, and t the residence time. Results show this statistical approach fits the data well and enables identifying important correlations which would not be highlighted with classical methods.
Moreover, the CCD approach was successfully adopted in other studies to observe the combined effect of each factor during the HTC of several biomasses (primary sludge, fish waste, poultry litter, sunflower stalk, and algae) [105][106][107].

Other Statistical Models
Other statistical models were also developed to understand the correlation between process parameters and product properties.
An alternative approach was proposed by Ismail et al. [108], who used the artificial neural network (ANN) coupled with a Kriging interpolation approach. The model was used to generate more data points from the collected data to understand the effect of temperature and residence time on the capacity of recovering carbon and inorganic phosphorus. The ANN model is useful for complex systems and is inspired by biological neural works, while the Kriging is an interpolation method based on stochastic relations. Results show that temperature has the highest impact on carbon recovery.
Besides, Xu et al. [109] applied the Taguchi statistical method, in which, adopting orthogonal arrays to perform experiments, control factors are identified to obtain the optimum process results. In particular, they investigated the HTC of microalgae, identifying as control factors reaction temperature, time, initial solid load, particle size, and catalyst amount (citric, acrylic, and sulfuric acid).

Reactor Models
The modeling of the HTC reactor is an essential step towards the design of a process that can be energetically efficient, economically feasible, and competitive with the current technologies. The reactor modeling should include several aspects, like heat transfer phenomena and the fluid-dynamics of the system as a function of the selected geometry and process conditions. Until now, little research has been carried out around this topic. Indeed, while many efforts have been done to understand the HTC reaction mechanism itself, few studies investigate the design of the reactor and its modeling. This shortage can be attributed to both the complexity of the process itself, still not fully understood and highly dependent on process conditions, and to the few data available concerning the larger-scale systems (like pilot and industrial).

Heat Transfer Models
Developing a heat transfer model is of great importance to evaluate and minimize energy consumptions. Currently, as highlighted by Heidari et al. [16], these models have been poorly investigated mainly due to a certain lack of research around: (1) the thermodynamics of HTC, i.e. the enthalpy of reaction (endothermic or exothermic process); (2) the change of physical properties (like porosity, density, and permeability) of the mixture of biomass and water during the process.
A simple model was developed by Baratieri et al. [63] to simulate the thermal behavior (also the transient phase) of an HTC lab apparatus previously designed by Fiori et al. [110]. To do this, the HTC system was reduced to a resistancecapacitance network considering thermophysical properties of the system like the reactor shell, the gaseous and liquid phases. A sketch of the system is shown in Fig. 4, where the HTC system at the set-point temperature T exchanges heat with the external environment (T U and T D are the upper and lower temperatures of the system, while T H is the temperature of the external electrical band heater). C 0 is the sum of heat capacities of pure water, stainless steel, and marble support, R 0 is the overall thermal resistance, R U and R D are the thermal resistances of the upper and lower areas, respectively. Then, heat transfer equations were written and used to simulate the thermal behavior of HTC starting from the reactor temperature profile experimentally measured. Although its easy implementation, this model is limited since only water is considered inside the HTC reactor and the reactor is dimensionless.

Computational models
Computational models can provide a detailed description of the process and a multitude of phenomena (e.g. heat transfer, transport phenomena, and reaction kinetics) considering variations in both time and space. Indeed, they consider important parameters like the specific geometry of the reactor, the porous structure of the biomass, and possibly the fluid-dynamics of the reacting mixture inside the HTC reactor. Despite their potential, they require heavy numerical computations as well the usage of proper software and properties databases, limiting their diffusion. Actually, the properties of a reactive heterogeneous system comprising biomasses are quite difficult to calculate or correlate. As also explained by Román et al. [35], a common procedure to develop computational models consists of: (1) specifying the reactor geometry; (2) establishing a proper mesh to divide the domain into cells; (3) setting of the initial and boundary conditions; (4) defining the properties of the different phases (e.g. state equations, kinetic equations); (5) defining the proper solvers and models to compute the system behavior in each cell.
Until now, according to author's knowledge, the only computational model has been developed by Álvarez-Murillo et al. [45]. They used COMSOL Multiphysics ® to characterize the decomposition of cellulose during HTC, using as inputs of the model the reactor geometry, materials, and reaction medium. Figure 5 shows the system geometry (reactor placed in an oven) as used in the program. Only water was considered to affect thermal properties of the reaction medium. The total heat transfer Q wall (in W) through the HTC walls was computed as for equation (23).
Where h conv and h rad are the convective and radiative heat transfer coefficients (W m −2 , computed by the program through traditional heat transfer relations), A S the surface of the external wall, T ∞ the bulk temperature, T S the temperature of the external wall. The kinetics of the process was evaluated by applying a first-order Arrhenius mechanism. Cellulose with a biomass to water ratio equal to 0.1 was chosen as the feedstock. Temperatures were investigated between 170 and 245 °C, with residence times between 2 and 11.8 hours. The simulated reactor was a stainless-steel cylinder with the following geometry: 18 cm height, external diameter of 9 cm, and wall thickness of 0.5 cm. The internal volume of the oven was defined as a rectangular prism of 25 × 25 × 30 cm.
Through the simulation, the variations of the investigated parameters (like relative mass loss of cellulose, HHV, H/C ratio) with time were computed at the different process conditions. The good fitting with experimental data confirms the reliability of the model. This presents also several limitations. Indeed, it considers only water affecting the thermal properties of the reaction medium and neglects the exothermic behavior of the HTC reactions. Moreover, the porosity of the cellulose is neglected and the heating method is uncommon for industrial and lab scale reactors. Anyway, despite the limits, this work represents the first computational model developed for HTC and then, it can provide the basis for further works. Reprinted with permission from Álvarez-Murillo et al. [45]. Copy-right© 2016, Elsevier Therefore, significant efforts should be done to develop this kind of models, which could lead to a better understanding and optimization of the process. A future step could be considering higher-order reactions and, as suggested by Heidari et al. [16], the investigation of the thermal effect of the HTC reactions on the heat fluxes.

Large scale modeling
The development of a validated plant simulation can provide important insight on the process dynamics [32]. Through a process simulator, operating parameters can be varied so to optimize the process considering both energy and economic aspects. Commercial suites (like Aspen Plus ® or UniSim Design ® ) and software environments (like Matlab) allow the development of such type of models. The current process simulators solve the process equations using proper databases and solving methods, which require accurate equilibrium and kinetics input data [32]. However, few data are currently available in literature regarding HTC thermodynamic and kinetics, limiting the reliability of the implemented models. Therefore, the inclusion of these aspects requires first an extensive availability and consistency of these data. Then, they can be integrated inside the specific software or suite using specific subroutines where kinetics and thermodynamics data (or equations) are specified.

HTC plant simulations
The available models on HTC are mainly techno-economic analyses, where both the energy and economic aspects are investigated. Even if important steps have been done, there are still several limitations. Indeed: (1) no reaction kinetics effects are considered in the HTC reactor, which is often modeled using equilibrium relations or yield data experimentally measured; (2) the thermodynamics of the HTC reactions (i.e. the quantification of the exothermic behavior of the process) is often not considered or over simplified using the few data available in the literature but relevant to lab-scale experiments; (3) mass transport and fluid-dynamics phenomena are always neglected; (4) large-scale and continuous processes are often modeled with experimental data obtained with batch reactors through laboratory experiments.
A steady-state model was developed by McGaughy and Reza [111] using Aspen Plus ® to simulate the HTC of food waste. The plant is thought to treat continuously 1 ton per day of biomass. The model bases on experimental data derived from an experimental campaign performed on a batch reactor and equations are forced to satisfy carbon content and solid yield according to experimental results. No reaction kinetics is considered and HTC reactions are modeled in a "Gibbs reactor", i.e. a continuous stirred tank reactor (CSTR) working at thermodynamic equilibrium which is calculated by minimizing the system Gibbs energy. In the model, food waste converts into hydrochar (that does not participate in the equilibrium), a liquid phase consisting of organic acids, and carbon dioxide. Utilities (pump, heater, filter, and dryer) are also considered, while the reaction is considered exothermic (1 kJ kg −1 ). The outputs of the model are the energy consumptions and the carbon dioxide produced by HTC. Even if results confirm the effect of temperature and moisture content of the starting feedstock on the process output parameters, the model suffers from over simplification and the optimistic assumption of considering all the heating requirements covered by the exothermic behavior of the HTC reaction.
A more complicated HTC plant model was developed by Erlach et al. [26]. The HTC flow-sheet, shown in Fig. 6, consists of several components (heat exchangers, flash tanks, a filter, coolers, pumps, pellet press) and includes a process water recirculation loop. Before entering the plant, the slurry is mixed and pre-heated.
The HTC reactor is modeled as a black box and characterized on the base of an experimental activity performed in a batch system (using poplar wood). Heat losses were assumed for all the auxiliaries. The energy consumption of the HTC reactor was evaluated based on the HHV of the various reactants and products. An exergy analysis was conducted, demonstrating that the reactor is one of the main sources of exergy loss. In addition, the coupling of the HTC plant with a gasification plant was also simulated.
The same plant scheme proposed by Erlach et al. [26] was used by Stemann et al. [112] to simulate an HTC plant for the treatment of empty palm oil bunches. The simulation was used to evaluate the auxiliary energy demand and investment costs. Results demonstrate that about 25% of the process water can be recycled, while the pumps and the pellet press are the most energy demanding equipment.
Aspen Plus ® has been recently used also by Akbari et al. [113] to perform an economic analysis on the equipment cost of two possible HTC plant configurations, which differ for their heat scheme and equipment. Both plants work with 160 tons per day of yard waste. The reactor was modeled using a "RYield reactor", in which the yield of the different phases is specified as input data on the base of experimental data. Such data were obtained in previous works. Results show that the optimal configuration leads to a cost of production of the hydrochar of 13.1 €/GJ. A techno-economic analysis was performed by Lucian and Fiori [114] to model the continuous treatment of organic wastes (20,000 tons per year) to pelletized hydrochar. All the plant equipment were considered working at steady state, and mass yield data were assumed as those obtained experimentally in a batch HTC reactor. A sketch of the adopted scheme is shown in Fig. 7. The model computes the plant efficiency (equal to 78%) and the minimum selling price of the pelletized hydrochar (equal to 157 €/ton) which makes the whole plant economically sustainable.
Saba et al. [115] performed a techno-economic analysis for the co-processing of coal and miscanthus for producing a solid fuel. A blend of coal and miscanthus (50:50, about 65,000 kg h −1 ) was treated in an HTC reactor. The analysis was limited to an economic evaluation, which leads to a breakeven selling price of 117.7 $/ton hydrocar.
Recently, Gómez et al. [116] modeled a continuous HTC plant adopting an alternative approach to previous works. Indeed, they assumed the HTC to occur is separated reaction stages. The model was developed with the UniSim Design ® process simulator and calibrated using operational data derived from a commercial company (Grenol GmbH). Thus, after the pre-heating and pumping stages, the stream of biomass (10-50 ton/day) enters several reactors, placed in sequence, and dedicated to: 1. hydrolysis (170 °C): biomass degrades into sugars; 2. intermediate compounds degradation (220 °C): sugars react to form aqueous and gaseous products (acetic acid and carbon dioxide); 3. aromatics formation (220 °C): the aqueous phase reacts to form aromatic compounds (HMF and furfurals); 4. polymerization process (220 °C): organic compounds re-polymerize while the unreacted biomass undergoes a solid-solid reaction.
Given the scarcity of reaction kinetics data, each reactor was characterized according to a hypothesized stoichiometry of each theoretical reaction. Once developed, the model was adjusted considering operational parameters from Grenol GmbH and fractional conversions. The model predicts the hydrochar composition and yield.

New Horizons: Integrated Plants and Continuous Reactors
Even if integrated plants and continuous reactors have been little or nothing investigated through modeling, we believe that these are fundamental topics for guiding authors towards future innovative modeling schemes.
Coupling HTC with other processes is of particular interest to maximize the efficiency of the whole plant. In particular, significant efforts have been done around the integration of HTC with anaerobic digestion (AD), a consolidated technology for the valorization of wet substrates into biogas. The process exploits the HTC process water (the "liquor"), which is an inevitable HTC by-product, to generate biogas. During the operation, the process water resulting from HTC feeds an AD plant to generate biogas through a series of bacterial and metabolic pathways [10,13]. Then, the biogas can be directly used for the production of electricity and heat through combustion in internal combustion engines or turbines or, alternatively, it can be injected in the natural gas grid after a purification stage to remove CO 2 and other minor constituents such as H 2 S. Several feedstock have been investigated, like digestate [21,117], sewage sludge [20,118], microalgae [12,13], orange pomace [22], spent coffee grounds [119], biowaste [120], and the organic fraction of municipal solid waste [24].
As explained by Merzari et al. [10], the chemical composition of the process water is fundamental to maximize the production of biomethane and avoid any inhibition of the AD process. Regarding the HTC process water to be recycled to AD, important parameters to control are those typical of the wastewater treatment field, like the chemical oxygen demand (COD), the total and suspended solids (TS and SS), the volatile fatty acids (VFAs), the ammonia, and the pH [10]. For AD, control parameters are the organic loading rate (OLR) and the hydraulic retention time (HRT).
Until now, promising research has been performed on the lab-scale to individuate the optimal parameters to avoid the formation of toxic or inhibiting compounds [121]. For example, the process water from sewage sludge HTC seems particularly suitable for AD if HTC temperatures and residence times are kept lower than 180-200 °C and 15-30 min, respectively, leading to methane yields between 0.022 and 0.277 L CH4 gCOD −1 [10]. Besides, Paul et al. [120], dealing with the HTC process water from biowaste, found that residence times longer than 30 min (at 240 °C) cause the formation of inhibitory chemicals for AD. Recently, Lucian et al. [24] tested the biomethane potential from the coupling of AD with HTC of organic wastes in two different scenarios: AD of the HTC liquor and AD of the HTC slurry (containing both the hydrochar and the liquor). Besides, an energy evaluation of the two scenarios was reported. On the one hand, results show that the highest biomethane potential is obtained by the slurry rather than the liquor (+ 363% and + 37% compared to the raw biomass, respectively). On the other hand, in the hypothesized plants, the use of the sole HTC liquor represents the most energy efficient pathway as the energy content of the hydrochar can be efficiently exploited through combustion.
Interestingly, Aragon-Briceno et al. [117] modeled different scenarios of an HTC-AD plant coupled with a cogeneration stage by the use of Aspen Plus ® . The analysis focuses on the energy duties of the plant, in which both hydrochar and biogas were considered as combustion resources. RYield reactors were used to model both the HTC and the AD reactors on the base of experimental data. Results show that scenarios with high solid loads lead to a lower heat production from the combustion of the biogas, but more electricity and heat are obtained from the combustion of the hydrochar. At the net, large solid load scenarios (30% the highest) lead to the highest electrical and thermal efficiencies, equal to 25.8 and 38.9%, respectively. Furthermore, Heidari et al. [23] compared an HTC-AD scenario with a direct combustion scenario for power production. The plants comprise a Rankine and a Bryton cycle to generate electricity and the thermodynamic modeling was performed through EES (Engineering Equation Solver). Results show the HTC-AD case has higher efficiencies only under high moisture contents of the starting biomass (sawdust in this case).
Regarding the integration of HTC with other technologies, our group recently developed and investigated a renewable solution consisting of a hybrid HTC-solar apparatus. The system consists of an HTC reactor directly coupled with a solar concentrator, which was experimentally demonstrated to completely cover the HTC energy needs [29]. Through this system, a "zero-energy technology" can be realized. Besides, our group also performed a techno-economic analysis of a large-scale plant consisting of a continuous HTC reactor coupled with a solar concentrator complete of thermal energy storage tanks: results show that this apparatus ensures a minimum fuel selling price of the hydrochar equal to 37.4 €/ton [30].
The development of a continuous reactor is a mandatory challenge for the industrialization of HTC and moving towards large scale plants. In addition to the intrinsic benefits of a continuous configuration (greater automation and better control, less use of manpower), continuous reactors enable to work with a wide range of HTC residence times when compared to batch reactors, instead characterized by useless long times to charge and discharge the reactors. Consistent research around this topic is still lacking and significant efforts are required. In particular, the biomass conveying system requires particular attention in the design phase [12]. Indeed, an effective and non-trivial pumping system is required to continuously feed the biomass-water stream to the high pressure HTC reactor [122]. Among the few authors that faced this topic, Hoekman et al. [79] adopted a modified two-screw extruder to realize a fully continuous reactor. The reactor consists of an initial rotating portion (with increasing pitch) for increasing the pressure of the biomass, a central part at a constant pitch, and a final part with another rotating screw. Water is injected at the end of the first section. Experimental activity was successfully performed on loblolly pine with very short residence times (20-30 s) and high temperature (290 °C). In another work, Hoekman et al. [123] designed a semi-continuous reactor for the conversion of lignocellulosic biomass. Operating with 2-3 kg of biomass per run and with short residence times, the system consists of a V-shaped reactor. The reactor consists of a feeding inclined auger and a second auger to convey the material out of the system. Moreover, Stemman and Ziegler [124] simulated a continuous HTC system through EES (Engineering Equation Solver), accounting for the exothermic behavior of HTC.

Life Cycle Assessment
The development of a Life Cycle Assessment (LCA) is an important step for the assessment of the environmental impacts of a technology. By investigating the whole life cycle, from the extraction to the end of life, LCA provides a holistic evaluation of the whole process. Until now, mainly due to the lack of large-scale data, few studies were performed on the LCA of HTC.
Recently, Gievers et al. [125] carried out a LCA of the HTC of sewage sludge derived from anaerobic digestion for four possible applications (agriculture, horticulture, co-firing in waste incineration, and co-firing in a lignite power plant). These were compared with the results of traditional processing, which implies incineration and landfilling of the ash. The global warming potential (in kg CO 2 equivalent) was assessed. According to the results obtained by the authors, total CO 2 emissions are due for the 63% to the electricity required to run the HTC plant and auxiliary process, for the 28% to the natural gas used for the heating, and for about the 9% to the treatment of the process water. Besides, the analysis demonstrates that using hydrochar for substituting fossil fuel has the highest potential to reduce the global warming potential of the sewage sludge treatment. This LCA could be improved considering more life cycle impact categories and more data for the HTC plant.
Liu et al. [126] performed a LCA to assess the impacts of the production of briquettes of hydrochar with coal (with a hydrochar percentage from 10 to 100%) to be used as a solid fuel for electricity production in a power plant. The hydrochar is produced starting from loblolly pine, and a fast HTC (30 s of residence time and temperature between 260 and 290 °C) was considered. The HHV of the hydrochar resulted equal to 24.4 MJ/kg. The stages considered in the LCA are the harvesting of the loblolly pine, the pre-processing, the HTC apparatus (a twin-screw extruder reactor), and the production of electricity in a power plant. The investigated parameters are the GHG emissions and life cycle energy use. Results show that the GHG impact of the hydrochar-coal briquettes is much lower than that of coal (from 1.033 to 0.149 kgCO 2,e / kWh from pure coal to 100% of hydrochar), while the difference in the energy use is relatively small (ranging between 10.64 and 12.63 MJ/kWh).
Benavente et al. [127] assessed the environmental impacts of the HTC of olive mill wastes compared with traditional management strategies (composting, anaerobic digestion, oil extraction, and incineration). HTC is combined with an incineration stage for heat and power production. Results show that this solution is less environmentally impacting than composting and anaerobic digestion, while no significant improvements are observed compared with incineration and oil extraction. This gap is mainly due to the water disposal stage required for HTC. However, as pointed out by Heidari et al. [16], important aspects like the impact of ash in the incineration phase were neglected in this LCA.
Other interesting LCA studies were performed by Ahamed et al. [128] to assess a combined plant HTC-oil system, and by Fornes et al. [129] regarding hydrochar derived from forest waste for horticultural applications.
Overall, the current LCA studies applied to HTC suffer from the lack of large-scale data. Moreover, future studies should include the valorization stages (like the coupling with anaerobic digestion) of the liquid phase.

Conclusions
Hydrothermal carbonization is acquiring a rising interest in the research community and industry for the valorization of organic waste and biomass residues. Through hydrothermal carbonization, biomasses with high moisture content often considered unpleasant can be upgraded into useful bio-based products, usable for energy applications, in agriculture or as a basis for advanced carbon materials. Several challenges have to be tackled to shift HTC from the laboratory activity to industry. The development of modeling and simulation tools is one of these. Through modeling, virtual scenarios can be implemented and the effects of process variations investigated. In this way, the process can be optimized so to improve its attractiveness to commercial parties.
This work presents the different modeling strategies present in the current state of art. These have been divided according to the system of analysis. Three main system are proposed: (1) HTC process itself, involving reaction kinetics and statistical models; (2) the reactor, with heat transfer and computational models that account for a variety of physical phenomena (like mass transfer and biomass porosity); (3) the HTC plant, with techno-economic analyses considering the equipment and auxiliaries. Besides, life cycle assessment models applied to HTC are also presented. Even if significant improvements have been done in the last years, more efforts are required to develop appropriate and reliable models able to drive a shift to the industrialization of HTC. Indeed, the current models suffer from the general complexity of the reaction mechanism (HTC involves heterogeneous reactions in a three-phase system), still not fully understood, and the lack of some important kinetics and thermodynamic data. Moreover, models are often limited to a narrow range of experimental conditions (i.e. type of feedstock, reaction temperature, residence time, and starting biomass load) and pilot plant data are often not available to develop reliable large-scale models.
Therefore, further research is required to develop more credible and reliable models. In this regard, the production of more reaction kinetics and thermodynamic data could be an important starting point. Then, the development of reactor simulations by using computational tools, accounting for mass and heat transfer phenomena, could drive to important improvements. Moreover, it is worth mentioning that, ideally, models should be generalized and user-friendly. In this way, they could increase their attractiveness for commercial parties, who could easily explore the possibility of adopting hydrothermal carbonization as a valuable and competitive technology for waste and biomass valorization.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.