Kinetic Modeling for the Gas-Phase Hydrogenation of the LOHC γ -Butyrolactone–1,4-Butanediol on a Copper-Zinc Catalyst

: Liquid organic hydrogen carriers (LOHCs) are an interesting alternative for hydrogen storage as the method is based on the reversibility of hydrogenation and dehydrogenation reactions to produce liquid and safe components at room temperature. As hydrogen storage involves a large amount of hydrogen and pure compounds, the design of a three-phase reactor requires the study of gas and liquid-phase kinetics. The gas-phase hydrogenation kinetics of LOHC γ -butyrolactone/1,4-butanediol on a copper-zinc catalyst are investigated here. The experiments were performed with data, taken from the literature, in the temperature and pressure ranges 200–240 ◦ C and 25–35 bar, respectively, for a H 2 / γ -butyrolactone molar ratio at the reactor inlet of about 90. The best kinetic law takes into account the thermodynamic chemical equilibrium, is based on the associative hydrogen adsorption and is able to simulate temperature and pressure effects. For this model, the conﬁdence intervals are at most 28% for the pre-exponential factors and 4% for the activation energies. Finally, this model will be included in a larger reactor model in order to evaluate the selectivity of the reactions, which may differ depending on whether the reaction takes place in the liquid or gas phase.


Introduction
The threat of climate change as well as geopolitical tensions are gradually reducing or even prohibiting the use of fossil fuels.Hydrogen, produced from renewable energy sources, is widely studied as a replacement [1].Among the hydrogen storage technologies, the most mature are compression or liquefaction and one of the most studied is physisorption.However, they have drawbacks regarding safety, due to high compression pressures and the risk of leakage and evaporation of flammable hydrogen [2][3][4].With regard to chemical storage of hydrogen in carrier molecules, this technology can involve a wide variety of chemical reactions and storage solutions.The most common ways involve powerto-X systems, by CO 2 and/or N 2 hydrogenation to fuels capable of being stored under reasonable conditions over extended periods [1].However, these systems often lack in reversibility because hydrogen is hardly retrieved [1,5].Another solution consists in binding hydrogen to metals under a solid state [6].Anyway, the most promising solutions in terms of hydrogen gravimetric density are currently not reversible and the reversible ones are limited in terms of hydrogen content.Another way is the use of a liquid organic hydrogen carrier (LOHC), which involves reversible reactions, is easy to transport and stable over long periods [1,4].Indeed, a LOHC consists in a couple of hydrogen-lean|hydrogen-rich molecules, which are respectively catalytically hydrogenated and dehydrogenated into the other molecule in order to store and release hydrogen.These organic liquids are fluids sharing similar properties with hydrocarbon fuels or common chemicals.Their stability at ambient temperature and atmospheric pressure in the liquid phase makes them compatible with the existing fuel infrastructures and thus LOHCs can be easily transported in tankers, trucks and/or pipelines [1,4,7].
Multiple LOHC systems have already been identified.The main ones highlighted in the literature are aromatics|cycloalkanes, such as toluene|methylcyclohexane (TOL|MCH), dibenzyltoluene|perhydro-dibenzyltoluene (H0-DBT|H18-DBT), benzyltoluene|perhydrobenzyltoluene (H0-BT|H12-BT) and couples based on N-heterocycles, such as N-ethylcarba-zole|perhydro-N-ethylcarbazole (H0-NEC|H12-NEC).The first cycloalkane-based LOHCs were developed because of their good hydrogen capacities (5.8-6.2%),their wide availabilities and the amount of knowledge on the processes of hydrogenation and dehydrogenation.However, reaction enthalpies, which dictate how much energy is lost during the storage process, are high for these LOHCs at around 65-69 kJ•mol H2 −1 , although NEC LOHC has a slightly lower reaction enthalpy (52 kJ•mol H2 −1 ).Safety and toxicity profiles have also been improved from the early light aromatics LOHCs to NEC, BT and DBT LOHCs, due to their higher flash points and lesser carcinogenic properties.The liquid state criterion, which is a key point for the process designs, is hardly met for some LOHCs at some conditions.For example, H12-NEC is solid at ambient conditions.Moreover, DBT LOHCs have reaction intermediates with high viscosities.NEC has also reversibility issues with high degradation when subject to harsh dehydrogenation temperatures.Finally, the use of expensive noble metal catalysts and relatively severe reaction conditions hinder the economical or environmental viability of the aromatic-based LOHCs [1,3,[7][8][9][10][11][12].CEA has recently identified a bioavailable LOHC: γ-butyrolactone|1,4-butanediol (GBL|BDO) [13].The latter, described for the first time by Onoda et al. [14], has low reaction enthalpies of 42 kJ•mol H2 −1 in liquid phase and 31 kJ•mol H2 −1 in gas phase, thus making it the most energy efficient.Its lower cycle operating temperature range (160-220 • C), compared to that of cycloalkane-based LOHCs _(300-350 • C), makes its process less energy consuming.BDO is partly produced from biomass, by hydrogenating bio-based succinic acid into GBL or by producing BDO directly from a fermentation of sugars over microorganisms [14][15][16][17], although most of it is currently produced from acetylene and formaldehyde [15][16][17].It has a H 2 capacity of 4.5 wt%, which is lower than that of the other LOHCs.The melting point of pure BDO at 20 • C could limit its use in colder weather or regions [16], but as for the stoichiometric ratio in feed, it remains liquid at reaction operating temperature and pressure, so reaction product separation is easier.Although GBL and BDO have good safety and toxicity profiles, GBL is known as a drug precursor because of its psychotropic properties and chemical resemblance to γ-hydroxybutyric acid (GHB) [18].They are commonly used as solvents and monomers and their uses are allowed in industrial processes in most countries.See Table S1 in Supplementary Materials for LOHC comparisons [1,16,[19][20][21][22].
For all LOHCs, except for TOL|MCH, hydrogenation and dehydrogenation are performed in three-phase reactors.Although LOHCs are usually compared with respect to their mass H 2 capacities, the volume H 2 capacity is therefore also useful for the design of suitable reactors as it indicates the volume of hydrogen involved during the reactions.For all LOHCs, the ratio H 2 /LOHC is significant (613-856 L H2 •L −1 LOHC ) and could act as a barrier to liquid-solid transfer by forming dry zones on the surface of the catalyst pellets.Thus, the knowledge of the reaction kinetics in both the liquid and gas phases is necessary to accurately design a three-phase reactor.As the main quality of LOHCs is the perfect reversibility of their reactions, the formation of co-products in both phases and reactions must also be assessed.This study is therefore a part of a larger study on GBL|BDO.Here, the paper is focused on the kinetics of the hydrogenation in the gas phase.The kinetics of the hydrogenation in the liquid phase will be the subject of another article.
The reaction pathway describing the equilibrium between GBL and BDO has been debated within the literature.Most have studied the reaction on copper-based catalysts in the liquid phase with temperature and pressure ranging from 200 to 250 • C and from 30 to 90 bar [14,15,23,24], respectively.Two parallel reactions, resulting in a dehydration towards tetrahydrofuran (THF) and n-Butanol (BuOH), occur at high temperature and are catalyzed by acid sites present in the catalyst support [23][24][25][26][27]; however, their mechanisms are unclear.Reaction intermediates, such as 2-Hydroxytetrahydrofuran, have been proposed as they could be involved in the formation of secondary products [19,24,27].[27].Chaudhari et al. proposed a mechanism where both GBL and BDO can lead up to the formation of THF and BuOH [23].However, no article seems to indicate that the latter are formed from intermediates, despites hydrogen imbalances between the reactants and the secondary products (except for the dehydration of BDO into THF).For instance, the experimental data in Schlander's Ph.D. thesis [25], performed in the gas phase, seem to indicate that GBL concentrations decrease faster than those of BDO as the secondary products build up in the mixture.As for Chaudhari et al., they contend THF is likely to be formed from the hydrogenation of GBL, rather than the dehydration of BDO [23]; we have also used Schlander's reaction pathway, shown in Figure 1, for the kinetic modeling.Kinetic models for the hydrogenation of GBL usually involve a simple power law model such as that proposed by Schlander [15,23,25].However, it (Model 1 in §2.2.1) proved unable to correctly predict the experimentally described pressure effect.Hermann et al. [24] proposed a Langmuir-Hinshelwood type model because an intermediate is thought to be formed during the hydrogenation of GBL, but this approach was first ruled out because of the number of parameters needed to be estimated, giving the model too many degrees of liberty and thus making it hard to extrapolate.This is why we decided to revise the kinetic study of the GBL gas-phase hydrogenation.The aim was to build a kinetic model, able to predict the effect of temperature and pressure, to be included in a larger reactor model including all phases involved.As the data from Schlander's works [25,26] were obtained under operating conditions that can mimic the composition of the gas phase in a trickle bed, the published data are used here to determine the appropriate kinetic law.

Ichikawa et al. proposed a model taking into account its formation without measuring it
The GBL hydrogenation was done in gas phase on a CuZnO catalyst in an isothermal tubular reactor.The fixed bed was loaded on a reactor length of 80 mm with crushed catalyst and glass of 400 µm-diameter particles.The catalyst loading was about 8 g.The molar ratio H 2 /GBL at the reactor entrance was about 90.Temperature and pressure ranges from 200 to 240 • C and from 25 to 35 bar, respectively, were screened.The experimental data presented in Figure 2 are plotted as a function of the residence time t V , which is the ratio between the reactor volume V R and the total volume flow rate Q.The bimetallic catalyst was prepared by co-precipitation and contained 15 mol% Cu and 85 mol% ZnO [28].See Section S2 in Supplementary Materials, and references [25,26,28] for further details.Regarding the potential transfer limitation, Schlander showed that the apparent kinetics could not be explained either by internal or external mass transfer (Section S3 in Supplementary Materials).
ns 2022, 3, FOR PEER REVIEW data presented in Figure 2 are plotted as a function of the residence time tV, which ratio between the reactor volume VR and the total volume flow rate Q.The bimetalli alyst was prepared by co-precipitation and contained 15 mol% Cu and 85 mol% ZnO See Section S2 in Supplementary Materials, and references [25,26,28] for further de Regarding the potential transfer limitation, Schlander showed that the apparent kin could not be explained either by internal or external mass transfer (Section S3 in Su mentary Materials).
In this paper, based on these previous bibliographic results, a reactor model wa described and then several kinetic laws were tested and kinetic parameters estimat order to predict the effect of temperature and pressure on the gas-phase hydrogenat GBL into BDO.The kinetic laws were established in interaction with the thermodyn to take into account the reversibility of the reactions.

Reactor Model
To establish mass balance equations suitable for simulation of the experiments, criteria were checked.GBL was highly diluted by hydrogen in the feed.As the H2 molar ratio at the reactor inlet was about 90, the GBL conversion induced at most a v tion of about 1.8% of the volumetric flow rate along the catalytic bed and a negl change in the gas phase composition for all experiments considered in the kinetics.In this paper, based on these previous bibliographic results, a reactor model was first described and then several kinetic laws were tested and kinetic parameters estimated in order to predict the effect of temperature and pressure on the gas-phase hydrogenation of GBL into BDO.The kinetic laws were established in interaction with the thermodynamics to take into account the reversibility of the reactions.

Reactor Model
To establish mass balance equations suitable for simulation of the experiments, some criteria were checked.GBL was highly diluted by hydrogen in the feed.As the H 2 /GBL molar ratio at the reactor inlet was about 90, the GBL conversion induced at most a variation of about 1.8% of the volumetric flow rate along the catalytic bed and a negligible change in the gas phase composition for all experiments considered in the kinetics.Consequently, the variation of the gas velocity along the reactor could be neglected as well as the possible temperature gradients induced by the reaction.So, the catalytic bed was considered as isothermal.Moreover, as the reactor diameter (d t ) to particle diameter (d p ) ratio was greater than 10 (Equation ( 1)), then wall effects causing a radial velocity profile resulting in deviations from plug flow could be neglected.In addition, as the reactor length (L B ) to particle diameter ratio was greater than 50 (Equation ( 2)), the axial dispersion resulting in consequences on the residence time distribution could also be neglected.
Consequently, an isothermal plug flow model was applied.Steady state was supposed to be reached.The material balance followed the evolution of a compound j within an elementary volume of catalyst V cata : where i is the reaction index and j is the compound index.The apparent kinetic rate is: As the diffusion has been estimated as a non-limiting step, the efficiency factor η i is considered equal to 1 (See Section S3 in Supplementary Materials).Capillary condensation has also been assessed in Section S3 in Supplementary Materials.
The experimental data used the scale t V , which is the residence time.Thus with x j is the molar fraction of condensable compound j equal to the ratio between the compound partial pressure P j and the partial pressure of GBL at the reactor entrance P GBL,0 and r j the production rates of compound j.

Kinetic Modeling
Model 1 and Model 2 were built upon power laws with Model 2 including equilibrium constants.Model 3 and Model 4 included in addition an inhibition parameter on the hydrogen pressure.All the models were based on the reaction scheme described in Figure 1.All the models presented in this article (whether kinetic or thermodynamic) were constructed on the basis of the ideal gas hypothesis (see Section S4 in Supplementary Materials).

Model 1
Model 1 is a simple power law model first investigated in modeling the hydrogenation of GBL into BDO.It is the model proposed by J.H. Schlander in his thesis [25].To avoid the interpretation bias due to the numerical methods, this model has been assessed here with the same method as used for the other models.That explains the difference between the values of the kinetic parameters estimated in this paper and given in the thesis.
Based on the reactions given in Table 1, the production rates of all active compounds are as follows: where k j are the kinetic parameters.

Reactions Kinetic Parameters Hypothesis
In order to study the effect of pressure into the model, the terms of the production rates were modified to introduce the hydrogen pressure and the chemical equilibrium.k THF and k BuOH were grouped under a single parameter k SP to better fit the single kinetic following of both side reactions.The reaction mechanism for this power law model is given in Table 2.Moreover, because of the large excess of hydrogen, the total pressure within the reactor was assumed to be the partial pressure of hydrogen.The activity of hydrogen for side reactions was neglected because of its low effect on kinetics.
Table 2. Reaction mechanism used for the construction of the second power law model: Model 2.

Reactions Kinetic or Thermodynamic Parameters Hypothesis
Model 3 is based on power law kinetics combined with an inhibition parameter added to the hydrogen partial pressure, which would describe slower kinetics due to a harder accessibility of the reactants to the active copper sites.Here, the adsorption mechanism of hydrogen is associative, meaning that each hydrogen molecule takes one site on the catalyst.The activity of hydrogen for side reactions was neglected because of its low effect on kinetics.The following reaction steps were considered for the construction of Model 3 (Table 3):

Reactions
Kinetic or Thermodynamic Parameters Hypothesis In this model, the production rates also depend on the surface coverage of hydrogen on the catalyst θ H 2 as well as the surface coverage of free sites θ s , as follows: The production rate r SP for the side products (SP) is the same as Equation ( 12) of Model 2.
A steady state is assumed for the adsorption of H 2 onto the catalyst and leads to expressions of the surface coverage of hydrogen on the catalyst θ H 2 (Equation ( 15)) as well as the surface coverage of free sites θ s (Equation ( 16)).
Moreover, because of the large excess of hydrogen, the total pressure within the reactor is assumed to be the partial pressure of hydrogen.The production rates used for Model 3 are:

Model 4
Model 4 is based on Model 3, assuming a dissociative H 2 sorption instead of an associative H 2 sorption.For the dissociative H 2 sorption, each hydrogen atom takes up one site on the catalyst.Model 4 was built similarly to Model 3. The activity of hydrogen for side reactions was neglected because of its low effect on kinetics.
The proposed reaction mechanism is given in Table 4.The induced production rates for the dissociative H 2 sorption are written in Equations ( 19) and (20).As for Model 3, the production rate r SP for the side products (SP) is the same as Equation ( 12) of Model 2 and a steady state is assumed for the adsorption of H 2 onto the catalyst, and leads to the expressions of the surface coverage of hydrogen on the catalyst θ H (Equation ( 21)) as well as the surface coverage of free sites θ s (Equation ( 22)).21) Moreover, because of the large excess of hydrogen, the total pressure within the reactor is assumed to be the partial pressure of hydrogen.The terms of production rates used for the simulation are:

Thermodynamic Equilibrium Modeling
The overall reaction equilibrium constant K is calculated through the standard Gibbs free energy change of reaction, ∆rG • : The standard Gibbs free energy change of reaction, ∆rG • , is calculated by Hess's law, using formation enthalpies ∆fH • j and entropies ∆fH • j for GBL and BDO.Thermodynamics properties were extracted from the NIST database and completed by the ProPhyPlus database and can be seen in Table 5.
Table 5. Thermodynamic data used for the modeling of the equilibrium constant.Data were extracted from the databases of NIST and ProPhyPlus.Formation enthalpies and entropies are calculated at a given temperature T according to Kirchhoff's laws: Finally, the heat capacities of the component C pj are calculated according to Hess's law, between the component and its hypothetical decomposition into elements:

Numerical Methods
All numerical operations were done using the software MATLAB.Differential equation systems were solved using a modified Rosenbrock formula of order 2 (function ode23s in MATLAB).
Parameters were estimated using a trust-region-reflective least squares algorithm (function lsqnonlin in MATLAB).The minimized objective function J(a) was the sum of gaps between the experimental data and the modeled data squared, represented through Equation (31).
For Model 1 and 2, parameters were estimated for each set of reaction conditions independently.For Model 3 and 4, parameters were estimated for the experiments at 200 and 240 • C independently and for the two experiments at 220 • C simultaneously in order to obtain kinetic constants independent of pressure conditions.For the last parameter estimation for Model 3, wherein both pre-exponential and activation energy parameters were simultaneously estimated, all the experimental data were used.
Because the uncertainties on experimental data were not available, the precision of estimated parameters only takes into account the uncertainties due to the modeling.Uncertainties on estimated parameters were calculated as follows: The covariance matrix Cov k was calculated through the Jacobian matrix given by the function lsqnonlin of MATLAB and t N-Mp is the Student number.

Kinetic Model Comparison
For Model 1, the model proposed by J.H. Schlander [25], the estimated kinetic constants and their uncertainties for a 97.5% confidence interval are given in Table 6.The confidence intervals for k THF and k BuOH show that these parameter values obviously are not significant.Table 6.Estimated kinetic parameters and their uncertainty for a 97.5% confidence interval for Model 1.

Experiment
This is due to the very large number of parameters to be estimated compared to the amount of experimental data available.For example, k THF and k BuOH should not be evaluated separately because all the secondary products produced by the two parallel dehydrations are drawn in a single dataset from the experiments.Moreover, k GBL and k BDO can easily be made independent of one another by a proper modeling of the reaction equilibrium, which is usually reached in the experimental section.The redundancy of kinetic constants makes the model not sensitive towards the estimated parameters, since a default in one kinetic constant can be balanced by another.This gives very large uncertainties on the estimated parameters, especially for the side reactions, thus making this model impossible to extrapolate.Moreover, the model does not take into account the effect of pressure.These unreliable results motivated the development of Model 2 (Table 7).Although we took into account the hydrogen pressure in a new model (Model 2), this model displayed a significant gap for the k kinetic parameter at 220 • C between 25 and 35 bar.Reaction rates are slower at higher pressure, which is not taken into account in that approach.Moreover, the confidence intervals for the k constant did not overlap for the experiments at 220 • C, meaning that the effect of pressure was not correctly modeled.So, it was impossible to fit an Arrhenius law for this model.Additionally, there was a nonsignificant k kinetic parameter at 240 • C. Indeed, uncertainties were low for every constant except for the k constant at 240 • C.This is explainable by the fact that the equilibrium was immediately reached right after the first point, and thus made the regression of this k constant happen between just the two first points.Therefore, simple power law models (Model 1 and 2) are inefficient when accounting for the drop of reaction rates as the pressure rises.
Model 3 (Table 8) was able to fit both experiments at 220 • C with a single set of constants in accordance to Arrhenius's law where kinetic parameters were only dependent on temperatures.Parameter estimation showed that the adsorption constant K H 2 was almost constant in the studied range of temperatures.This result is not surprising, considering that activation energies for adsorption constants are usually low, thus making the dependency of adsorption constants towards temperature weak.Table 8.Estimated kinetic parameters and their uncertainty for a 97.5% confidence interval for Model 3. Regarding Model 4, the estimated kinetic constants and their uncertainties for a 97.5% confidence interval are given in Table 9.The estimated constants for k SP and K H 2 kinetic parameters were quite similar to those of Model 3.For the k parameter, some differences in terms of value appeared for all the temperatures, while as for Model 3, the k parameter remained non-significant at 240 • C. Comparing Model 3 and 4 shows that, while kinetic constants were slightly different, the accuracy of the models was exactly the same, whether the adsorption was dissociative or associative.The model needed to slow down the reaction rates between the two pressure points of 25 and 35 bar regardless of the mathematical expression of this phenomenon.To evaluate which phenomenon, hydrogen dissociative or associative adsorption, was the limiting one taking place, further experiments with a wider range of pressures need to be performed.

Temperature Extrapolation
In order to design a LOHC process, a kinetic model able to predict the kinetics of the GBL hydrogenation into BDO at a given temperature and pressure, between 200 and 240 • C and between 25 and 35 bar, is needed.This paragraph is devoted to the extrapolation of Model 3, which was one of the only models successful in describing the slower kinetics at higher pressures.Extrapolation of Model 4 is seen to be unnecessary due to its similarity with Model 3. Extrapolation of pressures is built in the mathematical expression of the different models.For the extrapolation of temperatures, kinetic constants are modeled by Arrhenius's law: Pre-exponential factors A i and activation energies Ea i were re-estimated by replacing the kinetic parameters by their respective Arrhenius's laws and using all the experimental data at once.Because of the added parameters, the problem is highly dependent on initial conditions.To find a good set of initial parameters, exponential fitting of k and k SP constants into an Arrhenius's law were made.However, the pre-exponential factor for the k parameter given by the exponential fitting proved to be several orders of magnitude off a good initialization point.In order to find a better set of initial parameters, we noted that the uncertainty calculation on the k constant at 240 • C displayed a lot of play in the constant k, because the equilibrium was obtained right after the first point (Figure 3).Therefore, the k constant of experiment (b) was manually adjusted in order to fit perfectly into an Arrhenius's law.This adjustment was made through an exponential fit made with the estimated k constants for 200 and 220 • C, in furtherance of calculating the k constant of experiment (b) upon the exponential fit.Initial constants found with the retrofitted k constant showed decent results in the model fitting.Exponential fitting of k SP constants into an Arrhenius's law proved to be good enough and no adjustment was made here.
Model testing showed that the adsorption constant K H 2 was almost constant in the studied range of temperatures.Therefore, the K H 2 constant was manually adjusted to 1.03 (which fitted all confidence intervals) and treated as a constant across the studied range of temperatures.Table 10 displays the initial parameters used.Estimation of parameters using the aforementioned method led to a good fit of the extrapolated model to the experimental data (Figure 4).
The estimated pre-exponential factors and activation energies are given in Table 11 along with parity diagrams in Figure 5.For all the parameters, the confidence intervals were narrow and showed good reliability of the parameter values.The estimated pre-exponential factors and activation energies are given in Table 11 along with parity diagrams in Figure 5.For all the parameters, the confidence intervals were narrow and showed good reliability of the parameter values.

Conclusions
In this paper, a kinetic model has been developed for the gas-phase hydrogenation of γ-butyrolactone in the context of a hydrogen storage process.The data were gathered from experiments conducted by J.H. Schlander on the valorization of succinic acid.Indeed, several works addressing the hydrogenation of both maleic anhydride and succinic Surface coverage of free sites on the catalyst (-) Stoichiometric coefficient for reaction i and for compound j (j: GBL, BDO, THF, BuOH, side products) (-)

Figure 1 .
Figure 1.Reaction scheme used for the construction of the kinetic model.

Figure 3 .
Figure 3. Exponential fitting of k constants of Model 3 between 200 and 240 • C before adjustment.

Table 1 .
Reaction mechanism for the construction of the first power law model: Model 1.

Table 3 .
Reaction mechanism used for the construction of the kinetic model with inhibition of hydrogen partial pressure: Model 3.

Table 4 .
Reaction mechanism used to build the kinetic model based on a dissociative adsorption of H 2 : Model 4.

Table 7 .
Estimated kinetic parameters and their uncertainty for a 97.5% confidence interval for Model 2.

Table 9 .
Estimated kinetic parameters and their uncertainty for a 97.5% confidence interval for Model 4.

Table 10 .
Initial parameters used for the parameters estimation of the extrapolated Model 3.

Table 11 .
Estimated parameters for Model 3 with 95% confidence intervals. A

Table 11 .
Estimated parameters for Model 3 with 95% confidence intervals.