a, ,b b T T. .A A. . G Go or rd Ev Va
email: vladimir.solozhenko@univ-paris13.fr

TEM study of boron phosphide: Discovery of rhombohedral BP

Keywords: boron phosphide, transmission electron microscopy, polymorphism, twinning

Microstructure of sphalerite (3C) boron phosphide, BP, produced by self-propagated hightemperature synthesis has been studied by high-resolution transmission electron microscopy.

Along with numerous twins on the {111} 3C plane, layers of wurtzite (2H) polymorphic modification and previously unknown for BP rhombohedral (9R) structure were found which indicates trimorphism of BP.

Introduction

Boron phosphide BP is a hard (Vickers hardness H V  30 GPa [START_REF] Solozhenko | Mechanical properties of boron phosphides[END_REF]) refractory (melting temperature at ambient pressure is 2840 K [START_REF] Solozhenko | On melting of boron phosphide under pressure[END_REF]) and low-compressible (300-K bulk modulus is 174 GPa [START_REF] Solozhenko | Boron phosphide under pressure: in situ study by Raman scattering and X-ray diffraction[END_REF]) wide bandgap (E g = 2.1 eV [START_REF] Agui | Electronic structure of BP studied by resonant soft X-ray emission spectroscopy[END_REF]) semiconductor with outstanding chemical and high-temperature stability that makes it a promising material for a wide range of applications [START_REF] Woo | BP: synthesis and properties of boron phosphide[END_REF]. Under ambient conditions, BP crystallizes in cubic sphalerite (F-43m) structure with boron and phosphorus tetrahedrally coordinated to each other. Another tetrahedral structure of boron phosphide, hexagonal wurtzite (P6 3 mc) polymorph, has also been reported [START_REF] Villars | Pearson's Handbook of Crystallographic Data for Intermetallic Phases[END_REF], however, in the literature there is no data on the structure and properties of this phase.

The main limitation for the use of boron phosphide is the lack of relatively simple and economical methods of its production, especially, of BP single crystals. The disadvantages of the existing methods are the use of toxic and aggressive reagents, complicated technical implementation and high time consumption. Recently, two new methods of boron phosphide synthesis have been developed i.e. self-propagated high-temperature synthesis [START_REF] Mukhanov | Self-propagating hightemperature synthesis of boron phosphide[END_REF] and mechanochemical synthesis [START_REF] Mukhanov | Ultra-fast mechanochemical synthesis of boron phosphides, BP and B 12 P 2[END_REF] that are characterized by simplicity, high efficiency, low cost and good perspectives for large-scale production. Here we report the results of transmission electron microscopy studies of BP produced by self-propagated high-temperature synthesis.

Experimental

Microcrystalline powder of boron phosphide has been synthesized by self-propagating hightemperature reaction of boron phosphate and metallic magnesium: BPO 4 + 4Mg = BP + 4MgO using the method described earlier [START_REF] Mukhanov | Self-propagating hightemperature synthesis of boron phosphide[END_REF]. According to X-ray diffraction study (TEXT 3000 Inel, CuKα1 radiation) the sample is well-crystallized sphalerite BP with lattice parameter a = 0.45356(9) nm, which was close to the literature value 0.4538(2) nm [START_REF] Rundqvist | Crystal structure of boron phosphide BP. Congres International de Chimie pure et applique[END_REF]. The amount of B 12 P 2 boron subphosphide impurity caused by a partial decomposition of as-forming BP due to a high local temperatures in the reaction mixture did not exceed 1 vol.%. Microstructure of boron phosphide has been studied by high-resolution transmission electron microscopy (HRTEM) using JEM-2010 microscope with energy dispersive x-ray spectroscopy (EDS) attachment.

Results and discussion

According to TEM data (Fig. 1), the boron phosphide powder consists of flat polyhedral particles having dimensions from 50 to several hundred nanometers (Fig. 1a); significant part of them contains twins. Fig. 1b shows two particles containing twins. Traces of twinning planes are shown by arrows. Such particles with a sphalerite structure predominate in the sample. At the same time, other phases were detected in small quantities. are caused by diffuse scattering from {111} type stacking faults. Similar structural transformations were previously observed in carbon (diamond → lonsdaleite) [START_REF] Kulnitskiy | Polytypes and twins in the diamondlonsdaleite system formed by high-pressure and high-temperature treatment of graphite[END_REF] and in silicon (SiI → SiIV) [START_REF] Kulnitskiy | Mutual transformation between crystalline phases in silicon after treatment in a planetary mill: HRTEM studies[END_REF].

According to [START_REF] Momida | Effects of lattice parameters on piezoelectric constants in wurtzite materials: A theoretical study using first-principles and statistical-learning methods[END_REF], the lattice parameters of wurtzite BP are а = 0.320 nm and с = 0.531 nm, which are noticeably different from the lattice parameters of its closest structural analog, wurtzite BN (а = 0.25505 nm; с = 0.4210 nm [START_REF] Kurdyumov | Lattice parameters of boron nitride polymorphous modifications as a function of their crystal-structure perfection[END_REF]). The crystal lattice in the upper part corresponds to the sphalerite (3C) structure with lattice parameter a = 0.4538 nm, whereas the structure in the lower part is a rhombohedral (9R) lattice with trippled value of interplanar distance for (111) sphalerite planes i.e. 0.2623 = 0.786 nm. The tripling is visible both on HRTEM and FFT imags. As follows from the comparison of FFT images of two neighboring fragments (Fig. 3b) the (11-1) 3C ǁ (001) 9R and [1-10] 3C ǁ [11-20] 9R orientation relation is observed. Red circles indicate reflexes from 3C crystalline lattice. According to our assessment, the lattice parameters of the rhombohedral phase are: a = 0.321 nm, c = 2.358 nm. Thus, Fig. 3 shows a twinned BP particle which also contains two fragments of rhombohedral phase. The rhombohedral phase in the right part is formed as a result of twinning. Fig. 4 shows the schematic representation of different structures of the BP composition: a) ABAB packing (2H); b) ABCABC packing (3C) and c) ABABCBCAC packing (9R). Drawings containing models of crystal structures are made in a software package for visualization and analysis Ovito [START_REF] Stukowski | Visualization and analysis of atomistic simulation data with OVITO -the Open Visualization Tool[END_REF]. We could not find any information on the existence of rhombohedral BP in the literature, however, there is the experimental evidence for formation of rhombohedral polymorphic modification of zinc sulfide (ZnS) [START_REF] Buck | Trimorphism in zinc sulfide[END_REF], the structural analogue of BP. This rhombohedral structure of ZnS was found to occur in the 600-1020C temperature range as an intermediate phase between sphalerite (3C) and wurtzite (2H) structures. In our case, the appearance of rhombohedral BP can be associated with deformation of as-formed sphalerite lattice in the course of ultrafast self-propagated high-temperature reaction, which leads to a change of interplanar distances for {111} 3С planes and, consequently, results in the change of the lattice symmetry. Fig. 5a shows a two-phase particle. The left part of the particle corresponds to B 12 P 2 composition and contains polytypes, as evidenced by the elongated streaks passing through the spots in the corresponding FFT-image (Fig. 4b). The right part contains a twice-twinned fragment of BP composition (the corresponding FFT-image is shown in Fig. 4c). The mutual orientation of these two phases corresponds to the (003) B 12 P 2 ǁ (111) BP orientation relation. The arrows denote the {111} BP planes. One of the planes separates B 12 P 2 and BP crystals, while the other two are twinning planes in BP crystal. Such a junction of twinning planes is generally denoted as Σ = 9 [START_REF] Narayan | Dislocations, twins, and grain boundaries in CVD diamond thin films: Atomic structure and properties[END_REF]. In other words, Σ = 9 (second order twin) boundary is created when two Σ = 3 (first order twin) boundaries intersect. EDS spectra of B 12 P 2 is shown in Fig. 4d.

Previously, in boron subphosphide B 12 P 2 synthesized by the same self-propagating hightemperature method we have discovered two systems of twins i.e. conventional twins on the (0003) h plane and nanotwins resulting from duplication of the rhombohedral unit cell of B 12 P 2 along one of the basic vectors [START_REF] Kulnitskiy | Nanotwinning in boron subphosphide B 12 P 2[END_REF]. In the present study, neither these nor other twins were observed in the only few B 12 P 2 particles found. The latter can be explained by the fact that B 12 P 2 is a secondary phase formed due to the thermal decomposition of as-synthesized BP at points of the local overheating in the reaction front.

Conclusions

In boron phosphide produced by self-propagated high-temperature synthesis, along with traditional sphalerite (3C) polymorphic modification, fragments of wurtzite (2H) polymorph and previously unknown for BP rhombohedral (9R) structure were found which shed new light on polymorphism of boron phosphide.

Fig. 2

 2 Fig. 2 presents TEM and HRTEM images of sphalerite BP particle containing two systems of twins. White rectangles depict stacking faults which change the layer stacking sequence from ABCABC for cubic (3C) sphalerite structure to ABA for hexagonal (2H) wurtzite structure. The corresponding fast Fourier transform (FFT) image is shown in the inset. Strong streakings parallel to (111) and (1-1-1)

Fig. 3

 3 Fig.3shows a BP particle with a complex structure. In the upper part of HRTEM image twins in sphalerite 3C structure along the (111) plane are observed. The trace of the twinning plane is indicated by a vertical line at the top of the figure. In the lower part on the right, twins along the

Fig. 1 .

 1 Fig. 1. A characteristic TEM image of boron phosphide particles(a); two particles contain twins (b). Traces of twinning planes are shown by arrows.

Fig. 2 .

 2 Fig. 2. Two systems of twins in BP structure: (a) overview of the particle and (b) high resolution image of a portion of the particle where white rectangles depict the fragments of wurtzite (2H) structure with BAB and ABA layer stacking sequence; the corresponding FFT image is shown in the inset.

Fig. 3 .Fig. 4 .

 34 Fig. 3.

Fig. 5 .

 5 Fig. 5. (a) HRTEM image of two-phase particle: the left part is of B 12 P 2 composition and contains polytypes; the right part is twice-twinned fragment of BP composition. The arrows indicate {111} BP planes. (b) FFT image of the left part of (a); (c) FFT image of the right part of (a); d) EDS-spectra from B 12 P 2 .


Sindi Baco

Marcel Klinksiek

¥ , Rashid Ismail Bedawi Zakaria

Elizabeth Antonia Garcia-Hernandez

Mélanie Mignot

Julien Legros

Christoph Held

Valeria Casson Moreno

Sébastien Leveneur
email: sebastien.leveneur@insa-rouen.fr

Solvent effect investigation on the acid-catalyzed esterification of levulinic acid by ethanol aided by a Linear Solvation Energy Relationship

Keywords: Esterification, Kinetic modeling, Bayesian statistics, ethyl levulinate, ePC-SAFT, Kamlet-Abboud-Taft equation

When processing lignocellulosic biomass materials to obtain platform molecules such as levulinic acid (LA), alkyl levulinates or -valerolactone (GVL), the choice of solvent is of prime importance for kinetics. The knowledge of relationships between reaction kinetics and solvent serves as a decision tool for process design. To determine such relationships, esterification reactions was chosen because such reaction steps are present in several biomass conversion processes. In this work, kinetic models of LA esterification by ethanol over sulfuric acid in polar aprotic solvent (GVL) and polar protic solvents (water or ethanol) were developed and evaluated by Bayesian statistics. The apparent dissociation constants in solvents were estimated by ePC-SAFT approach to distinguish the proton concentration from the rate constants. The developed models can fit the experimental concentrations of ethyl levulinate and predict the proton concentration. Using the Kalmet-Abboud-Taft equation, linear relationships between estimated rate constants and solvent properties were established at different temperatures. We observed that solvents with low polarizability and high bond acceptor capacity should be favored for this reaction. Hence, the reaction of esterification is faster in ethanol solvent than in GVL solvent than in water solvent.

Introduction

The development of lignocellulosic biomass (LCB) valorization processes is growing because this substrate is renewable, abundant, and not competing with the alimentary sector. The shift from fossil-based refineries to biorefineries is challenging due to the diversity of biomass structures [START_REF] Singhvi | Lignocellulosic biomass: Hurdles and challenges in its valorization[END_REF], the biomass gathering and seasonability [START_REF] Pyrgakis | An Engineering Tool to Screen and Integrate Biomass Valorization Paths in Multiple-Feedstock Biorefineries[END_REF], the pretreatment complexity for LCB fractionation [START_REF] Zhao | Advances in pretreatment of lignocellulosic biomass for bioenergy production: Challenges and perspectives[END_REF], the identification of platform molecules [START_REF] Shinde | Lignocellulose-derived platform molecules: An introduction[END_REF], the need to discover new catalysts [START_REF] Singhvi | Lignocellulosic biomass: Hurdles and challenges in its valorization[END_REF], etc.

Alkyl levulinate (AL) and levulinic acid (LA) are promising platform molecules issued from cellulose and hemicellulose hydrolysis or alcoholysis. The hydrogenation of AL or LA leads to the production of γ-valerolactone (GVL), which is also a platform molecule [START_REF] Alonso | Gamma-valerolactone, a sustainable platform molecule derived from lignocellulosic biomass[END_REF]Capecci et al., 2021aCapecci et al., , 2021b;;[START_REF] Piskun | Kinetic modeling of levulinic acid hydrogenation to γ-valerolactone in water using a carbon supported Ru catalyst[END_REF][START_REF] Wang | Application of the concept of Linear Free Energy Relationships to the hydrogenation of levulinic acid and its corresponding esters[END_REF][START_REF] Yan | Catalytic reactions of gamma-valerolactone: A 59 platform to fuels and value-added chemicals[END_REF]. In the fuel sector, alkyl levulinates are used as fuel oxygenate additives or as blending components for diesel [START_REF] Peixoto | Production of ethyl levulinate fuel bioadditive from 5-hydroxymethylfurfural over sulfonic acid functionalized biochar catalysts[END_REF][START_REF] Peng | Catalytic upgrading of renewable furfuryl alcohol to alkyl levulinates using AlCl3 as a facile, efficient, and reusable catalyst[END_REF], improving the fuel quality and reducing pollutant emission [START_REF] Joshi | Ethyl levulinate: A potential bio-based diluent for biodiesel which improves cold flow properties[END_REF][START_REF] Windom | Advanced Distillation Curve Analysis on Ethyl Levulinate as a Diesel Fuel Oxygenate and a Hybrid Biodiesel Fuel[END_REF].

In 2017, Tian et al. showed that methyl and ethyl levulinates have higher anti-knock quality than Euro95 gasoline [START_REF] Tian | Anti-knock quality of sugar derived levulinic esters and cyclic ethers[END_REF].

Production of alkyl levulinate is commonly done through the esterification of LA [START_REF] Badgujar | A review on catalytic synthesis of energy rich fuel additive levulinate compounds from biomass derived levulinic acid[END_REF]. The other production route is via the alcoholysis of monosaccharide or biomass-derived furanic compounds reducing the side production of humins compared to water solvolysis [START_REF] Bhat | Recent advances in the preparation of levulinic esters from biomass-derived furanic and levulinic chemical platforms using heteropoly acid (HPA) catalysts[END_REF][START_REF] Menno | Production of levulinic acid and alkyl levulinates: A process insight[END_REF][START_REF] Li | Stable and efficient CuCr catalyst for the solvent-free hydrogenation of biomass derived ethyl levulinate to γ-valerolactone as potential biofuel candidate[END_REF][START_REF] Van Zandvoort | Formation, Molecular Structure, and Morphology of Humins in Biomass Conversion: Influence of Feedstock and Processing Conditions[END_REF].

From a chemistry and chemical engineering viewpoint, we need to address the transformation challenges of LCB into valuable chemicals. The tryptic kinetics-catalystthermodynamics is fundamental to find the optimal operating conditions, i.e., using less energy, being safe, and providing the maximum of target chemicals quickly.

In this perspective, the role of solvent or reaction mixture is fundamental if we want to succeed in the biomass-industry shift. A solvent can interfere with the kinetics and thermodynamics of a reaction system. In 2014, Mellmer et al. [START_REF] Mellmer | Solvent Effects in Acid-Catalyzed Biomass Conversion Reactions[END_REF] showed that aprotic polar solvent could improve the catalytic activity of protons. They suggested that the polar aprotic solvent can change the stabilization of the acidic proton relative to the protonated transition state during the acid-catalyzed conversion of xylose into furfural. They proposed to express the rate constant by using the Gibbs energy for proton solvation. Nevertheless, such thermodynamic data can be cumbersome to estimate.

In 1981, [START_REF] Kamlet | An Examination of Linear Solvation Energy Relationships[END_REF] developed a linear relationship between the solvent effect and the reaction kinetics. This equation, also known as Kamlet-Abboud-Taft equation (KAT), considers three solvent properties: its dielectric behavior, solvent ability to donate protons and electrons. This approach requires the knowledge of solvatochromic parameters [START_REF] Islam | Kamlet Taft Parameters: A Tool to Alternate the Usage of Hazardous Solvent in Pharmaceutical and Chemical Manufacturing/Synthesis -A Gateway towards Green Technology[END_REF][START_REF] Jessop | Solvatochromic parameters for solvents of interest in green chemistry[END_REF][START_REF] Laurence | A Database of dispersion-induction DI, electrostatic ES, and hydrogen bonding α1 and δ1 solvent parameters and some applications to the multiparameter correlation analysis of solvent effects[END_REF]Lu et al., 2002a;[START_REF] Nikolć | A linear solvation energy relationship study for the reactivity of 2-(4-substituted phenyl)-cyclohe-1-enecarboxylic, 2-(4substituted phenyl)-benzoic, and 2-(4-substituted phenyl)-acrylic acids with diazodiphenylmethane in various solvents[END_REF][START_REF] Sherwood | A Method of Calculating the Kamlet-Abboud-Taft Solvatochromic Parameters Using COSMO-RS[END_REF].

KAT was used to evaluate the solvent effect on spectral properties of compounds [START_REF] Aaron | A solvatochromic study of new benzo[a]phenothiazines for the determination of dipole moments and specific solute-solvent interactions in the first excited singlet state[END_REF][START_REF] Aaron | Quantitative treatment of the solvent effects on the electronic absorption and fluorescence spectra of acridines and phenazines. The ground and first excited singlet-state dipole moments[END_REF][START_REF] Párkányi | A quantitative study of the effect of solvent on the electronic absorption and fluorescence spectra of substituted phenothiazines: evaluation of their ground and excited singlet-state dipole moments[END_REF][START_REF] Rauf | Solvent effects on the absorption and fluorescence spectra of Cu(II)-phthalocyanine and DFT calculations[END_REF], on molecule solubility or solvent-solute interactions [START_REF] Holzweber | Mutual lewis acid-base interactions of cations and anions in ionic liquids[END_REF][START_REF] Marcus | Are solubility parameters relevant to supercritical fluids?[END_REF][START_REF] Paul | Optimization in solvent selection for chlorin e6 in photodynamic therapy[END_REF][START_REF] Thadathil | An Insight into Photophysical Investigation of (E)-2-Fluoro-N'-(1-(4-Nitrophenyl)Ethylidene)Benzohydrazide through Solvatochromism Approaches and Computational Studies[END_REF][START_REF] Wu | Excited state intramolecular proton transfer in 2-(2′arylsulfonamidophenyl)benzimidazole derivatives: insights into the origin of donor substituent-induced emission energy shifts[END_REF], on radical scavenging activity [START_REF] Jabbari | Solvent dependence on antioxidant activity of some waterinsoluble flavonoids and their cerium(IV) complexes[END_REF] or on reaction kinetics [START_REF] Auxenfans | Simultaneous pretreatment and enzymatic saccharification of (ligno) celluloses in aqueous-ionic liquid media: A compromise[END_REF][START_REF] García-Río | Nitrosation of Amines in Nonaqueous Solvents. 2. Solvent-Induced Mechanistic Changes[END_REF][START_REF] Raducan | Kinetics of hydrogen peroxide decomposition by catalase: Hydroxylic solvent effects[END_REF][START_REF] Ušcumlić | The study of linear solvation energy relationship for the reactivity of carboxylic acids with diazodiphenylmethane in protic and aprotic solvents[END_REF]. Nevertheless, the use of this approach in kinetics at different reaction temperatures is seldom. The study of Lu et al. (Lu et al., 2002b) on the hydrolysis of nitroaromatic compounds in different solvents and temperatures can be cited.

Catalyst studies for ALs synthesis via esterification or transesterification are vast [START_REF] Da Silva | Esterification of levulinic acid over Sn(II) exchanged Keggin heteropolyacid salts: An efficient route to obtain bioaditives[END_REF][START_REF] Melchiorre | Solvent-free transesterification of methyl levulinate and esterification of levulinic acid catalyzed by a homogeneous iron(III) dimer complex[END_REF][START_REF] Russo | Kinetic study of Amberlite IR120 catalyzed acid esterification of levulinic acid with ethanol: From batch to continuous operation[END_REF][START_REF] Zainol | Kinetics and thermodynamic analysis of levulinic acid esterification using lignin-furfural carbon cryogel catalyst[END_REF]. While the number of research papers on the solvent effect in esterification is relatively significant [START_REF] Janssen | Kinetics of lipasecatalyzed esterification in organic media: Correct model and solvent effects on parameters[END_REF][START_REF] Lemberg | Predicting the Solvent Effect on Esterification Kinetics[END_REF][START_REF] Liu | Effect of water on sulfuric acid catalyzed esterification[END_REF][START_REF] Lotti | Effects of methanol on lipases: Molecular, kinetic and process issues in the production of biodiesel[END_REF][START_REF] Riechert | Solvent effects on esterification equilibria[END_REF], there are just a few of them concerning the solvent effect for the esterification of LA [START_REF] Di | Biocatalytic upgrading of levulinic acid to methyl levulinate in green solvents[END_REF]. The study of [START_REF] Di | Biocatalytic upgrading of levulinic acid to methyl levulinate in green solvents[END_REF] was on the esterification of methyl levulinate over immobilized enzyme, and they showed that enzymatic activity is lower when the hydrogen-bond forces of the solvent are strong. Thus, hydrophobic solvents (i.e., ionic liquids) are preferred. However, they did not quantify the contribution of polarizability, hydrogen donor and acceptor on the kinetics.

The present work aims at developing three kinetic models for the synthesis of ethyl levulinate (EL) from LA in GVL solvent (polar aprotic), ethanol solvent (polar protic) and in water solvent (polar protic). These kinetic models, assessed by Bayesian inference, allowed a fair comparison between solvents. KAT approach cannot predict the acid dissociation constant in different solvents. Thus, the acid dissociation constants were estimated by ePC-SAFT equation of state [START_REF] Voges | Thermodynamics of the 58 alanine aminotransferase reaction[END_REF]. The knowledge of acid dissociation constants allows quantifying the proton concentration contribution to the rate constant. Esterification kinetic and thermodynamic constants in the three solvents were estimated via non-linear regression [START_REF] Delgado | Kinetic model assessment for the synthesis of γ-valerolactone from n-butyl levulinate and levulinic acid hydrogenation over the synergy effect of dual catalysts Ru/C and Amberlite IR-120[END_REF][START_REF] Stewart | Computer-Aided Modeling of Reactive Systems[END_REF].

The last stage was using kinetic constants in different solvents and at different temperatures in KAT equations to understand the solvent effect on kinetics.

Material and methods

Chemicals

The following chemicals were used without further purification: γ-valerolactone (wt% >= 99%) and LA (wt% >= 97%), both purchased from Sigma Aldrich. Ethanol absolute (wt% ≥ 99,7%) and acetone (wt% ≥ 99,8 vol%) were supplied by VWR Chemicals. Finally, ethyl levulinate (wt% ≥ 98%) was purchased from Acros Organics and sulfuric acid (wt% ≥ 98%) was obtained from Honeywell Fukla Chemicals. Distilled water was used.

Analytical method

GC-FID was used to measure the concentration of EL. The apparatus is BrukerScion GC436 equipped with a flame ionization detector, an autosampler, and a capillary column (ZB5, 30 m x 0.32 mm x 0.25 µm). Helium (99.99%) was the carrier gas used at a constant flow rate of 1.2 mL min -1 to transfer the sample from the injector, through the column, and into the FID-detector. The temperature of the injector and the detector were set at 250°C, while the program temperature of the oven was as follows: from 50°C (1 min) to 200°C (1 min) at 20°C min -1 . The injection volume was 1 μL and the split ratio was 1:20.

During the reaction, different samples were collected and diluted with acetone up to a dilution factor of 10000. In this way, a very low concentration of H2SO4 was present in the GC column, avoiding any damage. To evaluate the EL concentration in samples, calibration curves were done between 0 mg L -1 to 100 mg L -1 of EL.

Experimental procedure

Esterification reaction was carried out in a 300 mL glass reactor vessel, operating in batch mode. It was equipped with a jacket, in which water circulated and heated using a thermostat. A condenser was present to avoid mass loss due to the evaporation, and water was the cooling fluid. The pressure inside the reactor was equal to the atmospheric pressure. The presence of a mechanical stirrer was required in order to obtain a homogeneous reaction mixture. The stirrer was fixed at 450 rpm. Samples at different time intervals were taken by a plastic syringe through a reactor neck during a reaction.

For the performed kinetic experiments, starting solutions with only one reactant were poured into the reactor, and the temperature was brought to the desired reaction temperature prior to the addition of the preheated second reactant. The starting solutions were an aqueous solution of LA (or EL) and H2SO4 in water solvent, non-aqueous GVL solution of LA (or EL) and H2SO4 in GVL solvent reactions, and non-aqueous H2SO4 in ethanol solvent. The preheated second reactant added was ethanol in water solvent and GVL solvent reactions, while LA was in ethanol solvent. A typical reaction mass was ca. 230 grams, and the total sampling weight percentage was lower than 15%.

Tables 1, 2 and 3 show the experimental matrix for the kinetic experiments performed in the GVL, ethanol and water solvents. Besides probing solvent influences on the reaction, the operating conditions were varied for reaction temperature, catalyst, reactant, and product initial concentrations to develop robust kinetic models. For experiments in ethanol solvent, the maximum reaction temperature was fixed to 70°C because of the ethanol boiling point, which is ca. 78°C at 1 bar.

Reaction equilibria for any chemical reaction, also for dissociation reactions, are characterized by the activity-based equilibrium constant K th . The equilibrium constant is calculated with the mole fractions and activity coefficients of the reactants and products according to Eq. (3).

K th (T, p) = K x (T, p, x) ⋅ K γ (T, p, x) = ∏(x i ⋅ γ i ) ν i i ( 3
)
13 K th depends neither on the type nor on the concentration of the solvent. The dependence on reaction conditions is taken into account by the activity coefficients 𝛾 𝑖 . Thus, once the equilibrium constant is known, it can be used to predict the equilibrium concentrations in any other solvent or solvent mixture. K x is also known as the mole-fraction ratio at equilibrium, or sometimes also denoted apparent equilibrium constant; any concentration scale can be used instead of mole fractions. In the literature, the concentration-based value K c is typically used to determine the apparent dissociation constant 𝑝𝐾 𝑎 𝑐 by Eq. ( 5). The latter is not exactly the same as the standard value for 𝑝𝐾 𝑎 𝑎 , which is true only assuming that activity coefficients of products and reactants cancel out.

𝑝𝐾 𝑎 𝑎 = -log 10 K th (4)

𝑝𝐾 𝑎 𝑐 = -log 10 K c (5)
However, activity coefficients might crucially depend on the kind of solvent. The related value K γ requires knowledge about rational activity coefficients γ i * . For all species, the infinite dilution in pure water was applied as a reference state to calculate the activity coefficients γ i * . The equilibrium constant K th,1 and K th,2 for each dissociation step is defined according to Eq. (3) as follows:

K th,1 = x H +⋅x HS𝑂 4 - x H 2 S𝑂 4 ⋅ γ H + * ⋅γ HS𝑂 4 - * γ H 2 S𝑂 4 * → K 𝑥,1 = x H + ⋅x HS𝑂 4 - x H 2 S𝑂 4 = K th,1 ⋅ ( γ H + * ⋅γ HS𝑂 4 - * γ H 2 S𝑂 4 * ) -1 (6) K th,2 = x H +⋅x S𝑂 4 2- x HS𝑂 4 - ⋅ γ H + * ⋅γ S𝑂 4 2- * γ HS𝑂 4 - * → K x,2 = x H +⋅x S𝑂 4 2- x HS𝑂 4 - = K th,2 ⋅ ( γ H + * ⋅γ S𝑂 4 2- * γ HS𝑂 4 - * ) -1 (7)
The 𝑝𝐾 𝑎 𝑥 values are then the negative logarithm (to the base of 10) of the 𝐾 𝑥 value, which allows easy conversion to 𝐾 𝑐 and 𝑝𝐾 𝑎 𝑐 , respectively. Thus, the activity coefficients of the ions are required. These were modeled in this work with the thermodynamic equation of state ePC-SAFT advanced [START_REF] Bülow | ePC-SAFT advanced -Part I: Physical meaning of including a concentration-dependent dielectric constant in the born term and in the Debye-Hückel theory[END_REF]) (S1).

Athena Visual Studio [START_REF] Stewart | Computer-Aided Modeling of Reactive Systems[END_REF] was used to solve the ordinary differential equations-algebraic equations and estimate kinetic, equilibrium and Kamlet-Abboud-Taft parameters. This software uses Bayesian inference to analyze the estimated parameters. The GREGLUS package, implemented in Athena, can provide optimal parameter estimates with the 95% confidence intervals, expressed by the highest probability density (HPD). GREGPLUS provides the normalized parameter covariance matrix. Ordinary differential equations ODEs and algebraic equations were solved simultaneously by the solver DDALPUS algorithm, via a damped Newton method, implemented in Athena [START_REF] Caracotsios | Sensitivity analysis of initial value problems with mixed odes and algebraic equations[END_REF].

Kamlet-Abboud-Taft (KAT) equation

There is a distinction between non-specific and specific solute-solvent interactions in the KAT approach. The latter interaction, i.e., the solute-solvent specific interactions, comprises solvent Lewis-acidity and solvent Lewis-basicity interactions. The relationship between the rate constant in a solvent 𝑘 𝐶 𝑆𝑜𝑙𝑣𝑒𝑛𝑡 and the solvent properties is expressed in Eq. 8 [START_REF] Kamlet | An Examination of Linear Solvation Energy Relationships[END_REF][START_REF] Nikolć | A linear solvation energy relationship study for the reactivity of 2-(4-substituted phenyl)-cyclohe-1-enecarboxylic, 2-(4substituted phenyl)-benzoic, and 2-(4-substituted phenyl)-acrylic acids with diazodiphenylmethane in various solvents[END_REF][START_REF] Weiß | The physical significance of the Kamlet-Taft: π * parameter of ionic liquids[END_REF].

ln(𝑘 𝐶 𝑆𝑜𝑙𝑣𝑒𝑛𝑡 ) = 𝐴 0 + 𝑠 • 𝜋 * + 𝑎 • 𝛼 + 𝑏 • 𝛽 (8)
where, the terms 𝜋 * , 𝛼 and 𝛽 are solvatochromic solvent parameters that were supposed to be temperature independent in this study. The solvatochromic parameter α measures the solvent hydrogen-bond donor (HBD) acidity and evaluates the ability of a solvent to release a proton in a solvent to solute hydrogen bond. The values of α vary from zero (e.g., non-HBD solvents) to 1.0 (e.g., methanol). The solvatochromic parameter π* is an index of solvent dipolarity/polarizability, representing the ability of the solvent to stabilize a charge or a dipole via its dielectric effect. The values of π* vary from 0.00 for cyclohexane to 1.00 for dimethylsulfoxide. The solvatochromic parameter β measures the solvent hydrogenbond acceptor (HBA) basicity, and evaluates the solvent ability to accept a proton in a solute-to-solvent hydrogen bond. The β-scale values are from zero for non-HBD solvents to 1 for hexamethylphosphoric acid triamide (HMPT).

The term 𝐴 0 is the regression value. The regression coefficients (a, b, and s) present the relative susceptibilities of the solvent dependence of ln(𝑘 𝐶 𝑆𝑜𝑙𝑣𝑒𝑛𝑡 ) to the solvent parameters. Table 4 displays the solvatochromic solvent parameters for each solvent/environmental condition used in this study [START_REF] Islam | Kamlet Taft Parameters: A Tool to Alternate the Usage of Hazardous Solvent in Pharmaceutical and Chemical Manufacturing/Synthesis -A Gateway towards Green Technology[END_REF]. Athena visual studio was used to estimate a, b and s via DDAPLUS solver [START_REF] Stewart | Computer-Aided Modeling of Reactive Systems[END_REF]. A priori, water could affect sulfuric acid dissociation when GVL is used as a solvent. Fig. 1 shows that the increase of water concentration decreases the kinetics of EL synthesis, meaning that water is not needed to catalyze this reaction over sulfuric acid. Based on this experimental observation, proton in GVL solvent was expressed in (H + ) state and not in the hydroxonium ion (H3O + ). In ethanol solvent, proton was also expressed in (H + ) state.

This observation was also done by Liu et al. [START_REF] Liu | Effect of water on sulfuric acid catalyzed esterification[END_REF], who observed an inhibition effect of water on sulfuric acid for the esterification of acetic acid by methanol.

In Fig. 1 and the other, the error bar represents the standard deviation, and each sample was analyzed three times. was used as an input value to calculate K th at the same temperature and pressure using the ePC-SAFT advanced predicted activity coefficients referred to an infinite dilution state in pure water according to Eq.(S1.8). This K th value was then kept constant (at the considered temperature and pressure), but new activity coefficients for the species of H 2 SO 4 in the organic solvent were required to obtain 𝑝𝐾 𝑎 𝑐 in organic solvents. These were predicted by means of ePC-SAFT advanced using Eq. (S1.8), i.e., relating the ion's fugacity coefficient at infinite dilution in organic solvent to the infinite-dilution reference state in pure water.

Based on this, a K γ value was obtained, which finally yielded K c and 𝑝𝐾 𝑎 𝑐 using the aforedetermined constant K th (Eq. ( 3)). The same procedure was done for both dissociation steps. The calculated K th,2 for the second dissociation by ePC-SAFT using 𝑝𝐾 𝑎,2 𝑐 = 1.96 [START_REF] Izutsu | Electrochemistry in Nonaqueous Solutions: Second Edition[END_REF] interaction parameter between protons and organic solvents. For DMSO and DMF the prediction is in agreement with the literature with 𝑘 𝑖𝑗 = 0.17 and 𝑘 𝑖𝑗 = 0.177, respectively.

The dissociation of H 2 SO 4 is mainly determined by the differences of the dielectric constants of organic solvent to water. This procedure allowed reliable modeling results of 𝑝𝐾 𝑎 𝑐 values of H 2 SO 4 in DMSO, DMF and ACN (Fig. 2). Thus, the same strategy was applied to predict 𝑝𝐾 𝑎 𝑐 values of H 2 SO 4 in GVL and ethanol. The first dissociation constant of sulfuric acid based on concentration in GVL results to 9.13 • 10 -3 mol.L -1 , whereas the one in ethanol is even lower with 1.12 • 10 -7 mol.L -1 (Table S. 1.4). This clarify, that the dissociation is significantly reduced in organic solvents.

Additionally, Fig. 2 reveals the influence of the dielectric constant (𝜀 𝑟,𝐷𝑀𝑆𝑂 = 46.45, 𝜀 𝑟,𝐷𝑀𝐹 = 36.7, cf. Table S3.1) on the dissociation constants in DMSO (𝑝𝐾 𝑎 𝑐 = 1.99)

and DMF (𝑝𝐾 𝑎 𝑐 = 3.1).

In the following, we described our approach to determine the dissociation constant value for the second dissociation of sulfuric acid in water and in GVL solvent, the dissociation constant of LA in GVL and in water. In this study, the first dissociation of sulfuric acid in water was considered to be total, the second dissociation constant of sulfuric acid and the one of LA in ethanol was neglected based on the ePC-SAFT analysis.

Shen et al. have observed that DMF (N,N-dimethylformamide) and GVL have similar dielectric constant values [START_REF] Shen | γ-valerolactone (GVL) as a biobased green solvent and ligand for iron-mediated AGET ATRP[END_REF]. Hence, dissociation constants in DMF solvent are similar than those in GVL solvent [START_REF] Izutsu | Acid-Base Reactions in Non-Aqueous Solvents[END_REF].

Table S1.4 shows the values of dissociation constants for sulfuric acid, its second dissociation, acetic acid, benzoic acid and LA in different solvents [START_REF] Izutsu | Acid-Base Reactions in Non-Aqueous Solvents[END_REF]. The first dissociation constant of sulfuric acid in DMF is lower than in water. The second sulfuric acid dissociation constant is higher in water solvent, whereas it is very low in DMF solvent.

The dissociation constant of LA in water was found from the work of Kopetzki et al. [START_REF] Kopetzki | Transfer hydrogenation of levulinic acid under hydrothermal conditions catalyzed by sulfate as a temperature-switchable base[END_REF] and Piskun et al. [START_REF] Piskun | Kinetic modeling of levulinic acid hydrogenation to γ-valerolactone in water using a carbon supported Ru catalyst[END_REF]. To the best of our knowledge, there are no data concerning the dissociation of LA in non-aqueous solvents like GVL, DMF. In this study, the effect of temperature on the dissociation constant was supposed to be negligible.

The following notations are introduced: KI, KII and KIII, represent the apparent constants of the first dissociation of sulfuric acid, the second dissociation of sulfuric acid, and the dissociation of LA, respectively

In GVL solvent, the following values for the dissociation constants were used: KI: 9.13 . 10 - 03 mol.L -1 from ePC-SAFT approach; KII: 6.31 . 10 -18 mol.L -1 as the one found in DMF solvent (Table S1.4) and KIII: 3.00 . 10 -13 mol.L -1 . The latter value KIII, the dissociation constant of LA in GVL solvent, was evaluated as an average one between the dissociation constant of acetic and benzoic acids in DMF. In ethanol solvent, KI was estimated to be 1.12 . 10 -07 mol.L -1 from ePC-SAFT approach, and KII and KIII were found to be negligible.

Knopf et al. have developed an expression for the dissociation reaction of the bisulfate ion [START_REF] Knopf | Thermodynamic dissociation constant of the bisulfate ion from Raman and ion interaction modeling studies of aqueous sulfuric acid at low temperatures[END_REF] (𝐻𝑆𝑂 4 -↔ 𝑆𝑂 4 2-+ 𝐻 + ) in aqueous H2SO4 solutions. In the experimental conditions of Knopf et al., the sulfuric acid concentration is low; therefore, the activity coefficient of the solutions is close to unity. By dividing the apparent second dissociation from [START_REF] Knopf | Thermodynamic dissociation constant of the bisulfate ion from Raman and ion interaction modeling studies of aqueous sulfuric acid at low temperatures[END_REF] by the water concentration in the case of high dilution, i.e., 55.5 mol.L -1 , the apparent second dissociation constant for sulfuric acid can be obtained.

In water solvent, the following values for the dissociation constants were used: 𝐾 𝐼𝐼 = 0.011 55.5

and 𝐾 𝐼𝐼𝐼 = 2.57•10 -5 55.5

. The value of KII was calculated as an average value between 50 and 80°C, and the value KIII was obtained from Piskun et al. [START_REF] Piskun | Kinetic modeling of levulinic acid hydrogenation to γ-valerolactone in water using a carbon supported Ru catalyst[END_REF].

In sum, Table 5 displays the values of KI, KII and KIII used in this study in the different solvents.

GVL solvent EtOH solvent

Water solvent KI (first dissociation of H2SO4) 9.13 . 10 -03 mol.L -1 1.12 . 10 -07 mol.L -1 Infinity KII (second dissociation of H2SO4) 6.31 . 10 -18 mol.L -1 0 mol.L -1 0.011 55.5 KIII (dissociation of LA) 3.00 . 10 -13 mol.L -1 0 mol.L -1 2.57 • 10 -5 55.5 Our group used this approach to model the hydrolysis of acetic anhydride [START_REF] Garcia-Hernandez | Kinetic modeling using temperature as an on-line measurement: Application to the hydrolysis of acetic anhydride, a revisited kinetic model[END_REF]. A quasi-equilibrium hypothesis was used to determine the concentration of protons from sulfuric acid and LA dissociation and determine the concentrations of some intermediates.

Mechanism

Kinetic expressions

In Supplementary Material, the reader can find the derivation of the kinetic expressions in GVL, ethanol, and water solvents (S5).

In GVL solvent and ethanol solvent

For the catalyzed route, the quasi-equilibrium assumption was applied to the reversible proton donator steps 6 and 8 (Fig. S4.1). The rate-determining step is reversible step 7 (Fig.

S4

.1). The rate of esterification in GVL solvent can be derived as

𝑟 𝑒𝑠𝑡𝑒𝑟𝑖𝑓𝑖𝑐𝑎𝑡𝑖𝑜𝑛 𝐺𝑉𝐿 = 𝑘 𝑐 𝐺𝑉𝐿 • [𝐻 + ] • ([𝐿𝐴] • [𝐸𝑡𝑂𝐻] - 1 𝐾 𝑐 𝐺𝑉𝐿 • [𝐸𝐿] • [𝐻 2 𝑂]) (9)
The rate of esterification in ethanol solvent can be expressed as

𝑟 𝑒𝑠𝑡𝑒𝑟𝑖𝑓𝑖𝑐𝑎𝑡𝑖𝑜𝑛 𝐸𝑡ℎ𝑎𝑛𝑜𝑙 = 𝑘 𝑐 𝐸𝑡ℎ𝑎𝑛𝑜𝑙 • [𝐻 + ] • ([𝐿𝐴] • [𝐸𝑡𝑂𝐻] - 1 𝐾 𝑐 𝐸𝑡ℎ𝑎𝑛𝑜𝑙 • [𝐸𝐿] • [𝐻 2 𝑂]) (10)
In the modeling stage, 𝑘 𝑐 𝐺𝑉𝐿 (𝑇 𝑟𝑒𝑓 ), 𝐸 𝑎 𝐺𝑉𝐿 , 𝑘 𝑐 𝐸𝑡ℎ𝑎𝑛𝑜𝑙 (𝑇 𝑟𝑒𝑓 ) and 𝐸 𝑎 𝐸𝑡ℎ𝑎𝑛𝑜𝑙 were estimated.

In water solvent

For the catalyzed route, the quasi-equilibrium assumption was applied to the reversible proton donator steps 6 and 8 (Fig. S4.2). The rate-determining step is reversible step 7 (Fig.

S4

.2), and the rate of esterification can be expressed as

𝑟 𝑒𝑠𝑡𝑒𝑟𝑖𝑓𝑖𝑐𝑎𝑡𝑖𝑜𝑛 𝑊 = 𝑘 𝑐 𝑊 • [𝐻 3 𝑂 + ] [𝐻 2 𝑂] • ([𝐿𝐴] • [𝐸𝑡𝑂𝐻] - 1 𝐾 𝑐 𝑤 • [𝐸𝐿] • [𝐻 2 𝑂]) (11)
In the modeling stage, we estimated 𝑘 𝑐 𝑊 (𝑇 𝑟𝑒𝑓 ) and 𝐸 𝑎 𝑊 .

Expression of the apparent equilibrium constant

In the low sulfuric acid concentration range used in this study, the apparent equilibrium constant (calculated based on concentrations) was supposed not to be affected by the sulfuric acid concentration itself. Hence, apparent equilibrium constants calculated based on the concentration can be assumed to be the true thermodynamic equilibrium constants.

The effect of temperature on the true thermodynamic equilibrium constant is described by the law of Van't Hoff equation:

𝑑ln (𝐾(𝑇)) 𝑑𝑇 = ∆𝐻 𝑟 0 𝑅𝑇 2 (12)
where, ∆𝐻 𝑟 0 stands for the standard reaction enthalpy change. Assuming that ∆𝐻 𝑟 0 is independent from T, the integration of Eq.( 12) from a reference temperature Tref to an arbitrary temperature T leads to:

𝑙𝑛 𝐾(𝑇) 𝐾(𝑇 𝑟𝑒𝑓 ) = -𝛥𝐻 𝑟 𝑜 𝑅 ( 1 𝑇 - 1 𝑇 𝑟𝑒𝑓 ) (13)
In the modelling stage, 𝐾 𝑐 𝑆𝑜𝑙𝑣𝑒𝑛𝑡 (𝑇 𝑟𝑒𝑓 ) and 𝛥𝐻 𝑟 𝑆𝑜𝑙𝑣𝑒𝑛𝑡 were estimated for each solvent.

Calculation of protons concentration

In order to facilitate reading, the derivation of concentration expression for the three solvents was described in detail in Supplementary Materials (S6). Autoprotolysis of water and ethanol was neglected in the three solvents. To determine the proton concentration, molar balances on sulfate and organic species were performed, and the electroneutrality principle was used.

Calculation of proton concentrations in GVL solvent and ethanol solvent

In these solvents, the expression for the proton concentration is:

[𝐻 + ] = 𝐾 𝐼 • [𝐻 2 𝑆𝑂 4 ] 𝑂 [𝐻 + ] +2•𝐾 𝐼 •𝐾 I𝐼 •[𝐻 2 𝑆𝑂 4 ] 𝑂 •( 1 [𝐻 + ] ) 2 1+𝐾 𝐼 • 1 [𝐻 + ] +𝐾 𝐼 •𝐾 𝐼𝐼 • 1 [𝐻 + ] 2 + 𝐾 𝐼𝐼𝐼 • 1 [𝐻 + ] •([𝐿𝐴] 𝑂 -[𝐸𝐿]) 1+𝐾 𝐼𝐼𝐼 • 1 [𝐻 + ] (14)
The initial proton concentration was calculated by neglecting KII and KIII.

Calculation of proton concentrations in water solvent

The first dissociation of sulfuric acid (KI) in water solvent tends to infinity. Thus, the implicit expression for the hydroxonium concentration is:

[𝐻 2 𝑂] [𝐻 3 𝑂 + ] [𝐻 2 𝑆𝑂 4 ] 𝑂 [𝐻 2 𝑂] [𝐻 3 𝑂 + ] + * 𝐾 𝐼𝐼 * [𝐻 2 𝑂] 2 [𝐻 3 𝑂 + ] 2 + 2 * 𝐾 𝐼𝐼 [𝐻 2 𝑂] 2 [𝐻 3 𝑂 + ] 2 [𝐻 2 𝑆𝑂 4 ] 𝑂 [𝐻 2 𝑂] [𝐻 3 𝑂 + ] +𝐾 𝐼𝐼 * [𝐻 2 𝑂] 2 [𝐻 + ] 2 + 𝐾 𝐼𝐼𝐼 [𝐻 2 𝑂] [𝐻 3 𝑂 + ] * [𝐿𝐴] 𝑂 -[𝐸𝐿] 1+𝐾 𝐼𝐼𝐼 [𝐻 2 𝑂] [𝐻 3 𝑂 + ] = [𝐻 3 𝑂 + ] (15)
The initial hydroxonium concentration was calculated using the same approach as in the reference [START_REF] Leveneur | Kinetic study and modeling of peroxypropionic acid synthesis from propionic acid and hydrogen peroxide using homogeneous catalysts[END_REF]. This leads to:

[𝐻 3 𝑂 + ] = 1 2 • [𝐻 2 𝑆𝑂 4 ] 0 + √ [𝐻 2 𝑆𝑂 4 ] 0 2 4 + 2 • 𝐾 𝐼𝐼 • [𝐻 2 𝑆𝑂 4 ] 0 • [𝐻 2 𝑂] + 𝐾 𝐼𝐼𝐼 • [𝐻 2 𝑂] • [𝐿𝐴] (16)
The solver DDAPLUS algorithm was used during the modeling stage to determine the proton concentration by solving Eqs ( 14) and (15).

Material balances and modeling results

The rate equation is expressed by Eqs ( 9), ( 10) or ( 11). Material balances in solvents are:

dC 𝐿𝐴 𝑑𝑡 = -𝑟 𝑒𝑠𝑡𝑒𝑟𝑖𝑓𝑖𝑐𝑎𝑡𝑖𝑜𝑛 ( 17
)
dC 𝐸𝑡ℎ𝑎𝑛𝑜𝑙 𝑑𝑡 = -𝑟 𝑒𝑠𝑡𝑒𝑟𝑖𝑓𝑖𝑐𝑎𝑡𝑖𝑜𝑛 (18) dC 𝐸𝐿 𝑑𝑡 = 𝑟 𝑒𝑠𝑡𝑒𝑟𝑖𝑓𝑖𝑐𝑎𝑡𝑖𝑜𝑛 (19) dC 𝑊 𝑑𝑡 = 𝑟 𝑒𝑠𝑡𝑒𝑟𝑖𝑓𝑖𝑐𝑎𝑡𝑖𝑜𝑛 (20)
Ordinary differential equations ODEs ( 17)-( 20) and algebraic equations ( 14) or ( 15) were solved simultaneously by the solver DDALPUS algorithm [START_REF] Caracotsios | Sensitivity analysis of initial value problems with mixed odes and algebraic equations[END_REF].

During the parameter estimation stage, the GREGPLUS package minimizes the objective function 𝑆(𝜃), and can calculate the maximum posterior probability density of the different estimated parameters 𝜃 and the values of the posterior distribution of the tested models.

𝑆(𝜃) = ∑ ([𝑦 𝑢 -𝑓 𝑢 (𝜃)]) 2 𝑛 𝑢=1 (21)
where, 𝑦 𝑢 and 𝑓 𝑢 (𝜃) are experimental and simulated concentrations of EL, respectively.

A modified Arrhenius equation was used to express the rate constants during the modeling stage to ease the parameter identification and decrease the correlation between the preexponential factor and activation energy [START_REF] Buzzi-Ferraris | Planning of experiments and kinetic analysis[END_REF].

𝑘 𝑐 (𝑇) = 𝑒𝑥𝑝 [𝑙𝑛 (𝑘 𝑐 (𝑇 𝑟𝑒𝑓 )) + 𝐸 𝑎 𝑅•𝑇 𝑟𝑒𝑓 • (1 - 𝑇 𝑟𝑒𝑓 𝑇 )] (22)
where, Tref is a reference temperature which is the average temperature of the experimental matrix.

The apparent equilibrium constant was linearized as

𝐾 𝑐 (𝑇) = 𝑒𝑥𝑝 [𝑙𝑛 (𝐾 𝑐 (𝑇 𝑟𝑒𝑓 )) + ∆𝐻 𝑟 𝑅•𝑇 𝑟𝑒𝑓 • (1 - 𝑇 𝑟𝑒𝑓 𝑇 )] (23)
The following constants were estimated: 𝑙𝑛 (𝑘 𝑐 (𝑇 𝑟𝑒𝑓 )), Table 6 shows the values of the estimated apparent constants with a reference temperature of 64 °C for the kinetic constant. One can notice that the credible intervals for each estimated constant are short, meaning that the experimental matrix was well designed (Table 1). Esterification of LA by ethanol in GVL solvent is an endothermic system. The values of 𝑘 𝑐 𝐺𝑉𝐿 (𝑇 𝑟𝑒𝑓 ), 𝐸 𝑎 𝐺𝑉𝐿 , 𝐾 𝑐 𝐺𝑉𝐿 (𝑇 𝑟𝑒𝑓 ) and 𝛥𝐻 𝑟 𝐺𝑉𝐿 are displayed in Table S7.1. Table S7.2

shows that correlations between the estimated constants are very low, meaning that the estimated constants were well identified. Globally, the parity plot (Fig. S7. [START_REF] Values | [END_REF] shows that the developed model fits well the experimental concentration of EL.

The residual analysis (Figs S7.2-S7.5) shows that the residuals are normally distributed versus the events, time, simulated and experimental EL concentrations. These normal distributions confirm that the developed model is correct.

Kinetic modeling results in ethanol solvent

Table 7 shows the values of the estimated constants with a reference temperature of 57 °C.

As in the previous case, it can be noticed that the credible interval for each estimated constant is short (except for 𝛥𝐻 𝑟 𝐸𝑡ℎ𝑎𝑛𝑜𝑙 ) showing the good design of the experimental matrix (Table 3). The values of 𝑘 𝑐 𝐸𝑡ℎ𝑎𝑛𝑜𝑙 (𝑇 𝑟𝑒𝑓 ), 𝐸 𝑎 𝐸𝑡ℎ𝑎𝑛𝑜𝑙 , 𝐾 𝑐 𝐸𝑡ℎ𝑎𝑛𝑜𝑙 (𝑇 𝑟𝑒𝑓 ) and 𝛥𝐻 𝑟 𝐸𝑡ℎ𝑎𝑛𝑜𝑙 are displayed in Table S7.3. Table S7.4 shows the correlations between the estimated constants. In general, the correlation between estimated parameters is very low, except for, 𝐾 𝑐 𝐸𝑡ℎ𝑎𝑛𝑜𝑙 (𝑇 𝑟𝑒𝑓 ) and 𝛥𝐻 𝑟 𝐸𝑡ℎ𝑎𝑛𝑜𝑙 , which is 0.857. Nevertheless, according to Toch et al. [START_REF] Toch | A systematic methodology for kinetic modeling of chemical reactions applied to n-hexane hydroisomerization[END_REF], two parameters are correlated if the binary correlation coefficient is Fig. 4 shows the fit of the model to some experimental concentrations of EL with the 95% prediction intervals and the mean estimated values. As in the system with GVL solvent, it is possible to confirm that the model fits well the experimental concentrations, and most of the experimental concentrations lie between the intervals. The initial kinetics and the equilibrium phase are well predicted with this model. Globally, the parity plot (Fig. S7.11) shows that the developed model fits well the experimental concentration of EL.

Kinetic modeling results in water solvent

The residual analysis (Figs S7.12-S7.15) shows that the residuals are normally distributed versus the events, time, simulated and experimental EL concentrations. These normal distributions confirm that the developed model is correct.

Explanation of solvent effect on kinetics

Using the kinetics and equilibrium constants obtained previously (Tables 678), we simulated reactant and product concentrations and the reaction rate of three experiments carried out in similar operating conditions in each solvent (Table S8.1). In these experiments, sulfuric acid and LA concentrations were similar.

Figs. 6 A&B show that the esterification rate is faster in ethanol solvent. At 50°C, the initial reaction rate in ethanol is ca. 8 times faster than in water solvent. The same thing can be observed at 70°C, where the esterification rate of water is still the lowest: 4 times slower than in GVL and 6 times in ethanol solvent.

Both at 50°C and 70°C (Fig. 6.A&B), one can observe that the initial kinetics in ethanol is faster than in GVL solvents, even if the difference is not that remarkable since at 50°C, the esterification rate is 1.5 times higher than in ethanol.

Figs. 6 C&D show that maximum EL yields are obtained in ethanol solvent. At 70°C, the final concentration of EL represents 99% of the initial concentration of LA in ethanol compared to 71% for the water solvent.

However, the benefits of using GVL solvent are its highest boiling point, allowing experiments at higher temperatures to increase the reaction rate.

Fig. 6. Temperature and solvent effect on esterification reaction rate and ratio

[𝐸𝐿]

[𝐿𝐴] 0 .

The Kamlet-Abboud-Taft equation was used to find the relationships between solvent properties and kinetics. The solvatochromic solvent parameters (Table 4) were used in Eq.( 8). The parameters (s, a and b) of Eq.( 8) were estimated by using the solver DDAPLUS, implemented in Athena Visual Studio. The responses were the natural logarithm of the estimated constants at different temperatures from Tables 678. Fig. 7 shows the parity plot between the estimated values of ln(𝑘 𝐶 𝑆𝑜𝑙𝑣𝑒𝑛𝑡 ) 𝑆𝑖𝑚 and the ones obtained from the regression using Kamlet-Abboud-Taft equation ln(𝑘 𝐶 𝑆𝑜𝑙𝑣𝑒𝑛𝑡 ) 𝐾𝐴𝑇 . One can notice that KAT equation can predict the kinetic constant in different solvents. Table 9 shows the values of A0, s, a and b at different temperatures (Eq. ( 8)). From Table 9, one can notice that solvent effect on the rate constant is mainly due to its ability to accept hydrogen-bond (measured by b value) and its ability to stabilize a charge (measured by s value). These two properties have antagonist effect on the rate constant value

Conclusions

This manuscript proposes a kinetic investigation to synthesize ethyl levulinate via the acidcatalyzed esterification of levulinic by sulfuric acid in three solvents: ethanol solvent (polar protic), GVL solvent (polar aprotic) and water solvent (polar protic). In order to perform a fair comparison, kinetic models were developed for each solvent at different temperatures, reactant, product and catalyst concentrations, and by taking into account the proton concentration. ePC-SAFT approach was used to estimate the sulfuric acid first dissociation in GVL solvent and ethanol solvent. Bayesian inference was used to evaluate the credible intervals of the estimated constants for each model. The developed models were found to be reliable.

In the absence of sulfuric acid, the reaction rate was found to be negligible in solvents. The esterification rate follows the order: 𝑟 𝑒𝑠𝑡𝑒𝑟𝑖𝑓𝑖𝑐𝑎𝑡𝑖𝑜𝑛 𝐸𝑡ℎ𝑎𝑛𝑜𝑙 >𝑟 𝑒𝑠𝑡𝑒𝑟𝑖𝑓𝑖𝑐𝑎𝑡𝑖𝑜𝑛 𝐺𝑉𝐿 >𝑟 𝑒𝑠𝑡𝑒𝑟𝑖𝑓𝑖𝑐𝑎𝑡𝑖𝑜𝑛

𝑊𝑎𝑡𝑒𝑟

. The conversion at the equilibrium follows the same order. For instance, at 70°C, the initial rate in ethanol solvent, GVL solvent and water solvent was calculated to be 0.15, 0.08 and 0.02 mol.L -1 .min -1 , respectively; the conversion of LA at equilibrium in ethanol and GVL solvents was calculated to be respectively 99% and 87% compared to 71% in water solvent at 70°C.

Based on the developed kinetic models, the Kamlet-Abboud-Taft linear solvation energy was used to establish a relationship between kinetics and solvent properties in different solvents and reaction temperatures. Solvents able to stabilize a charge could stabilize the cation transition state and proton, decreasing the reactivity of these species and thus decreasing the kinetics.

by ethanol by evaluating the contribution of polarizability, H donor and acceptor capacity.

By knowing the KAT properties of a solvent, one can predict the kinetics in this solvent.

One can also use this information to choose the most appropriate solvent to get faster kinetics. Solvent with low polarizability capacity and high hydrogen-bond donor capacity should be used for LA esterification over acid catalyst. The continuation of this work could be the kinetic study over heterogeneous and the temperature dependence of solvatochromic parameters.

Notation

Fig. 1 .

 1 Fig. 1. Effect of water concentration on the kinetics of esterification by comparing Runs 8, 9, 10 and 14 at 60 °C and [H2SO4]0: 0.09 mol.L -1 .

Fig. 2 .

 2 Fig. 2. 𝑝𝐾 𝑎 𝑐 values of H2SO4 in various solvents at 298.15 K and 1 bar. Literature values

Figs S4. 1

 1 Figs S4.1 and S4.2 illustrate reaction mechanisms in GVL, water and ethanol solvents. The

Fig. 3 Fig. 3 .

 33 Fig.3shows the fit of the model to some experimental concentrations of EL with the 95%

  of experimental [EL] with standard deviation Mean value of [EL] simulated 95% prediction interval case.

Fig

  Fig. S7.6 shows the parity plot of the ethanol solvent. With an R 2 = 0.94, it is possible to

Fig. 5

 5 Fig.5shows the fit of the model to some experimental concentrations of EL with the 95%

  . The solvent ability to donate hydrogen-bond has less influenced (measured by a) on the rate constant value. The negative values of s means that solvents with a high ability to stabilize a charge decrease the kinetic constants, whereas the positive value of b means that hydrogen-bond acceptor solvents increase the kinetic constants. Thus, solvents with low values for π* and high values for β should be favored for the acid-catalyzed esterification. Solvents that can stabilize a charge could stabilize the carbocation intermediate Int1 and/or proton (Figs S4.1 and S4.2), decreasing their reactivities.

  on concentration of LA esterification in GVL solvent 𝐾 𝐶 𝐸𝑡ℎ𝑎𝑛𝑜𝑙 Equilibrium constant based on concentration of LA esterification in excess of ethanol 𝐾 𝐶 𝑊 Equilibrium constant based on concentration on LA esterification in excess of water Kx Mole-fraction ratio at equilibrium or 'apparent' equilibrium constant 𝑘 𝐶 𝐺𝑉𝐿 Rate constant based on concentration for LA esterification in GVL solvent 𝑘 𝐶 𝐸𝑡ℎ𝑎𝑛𝑜𝑙 Rate constant based on concentration for LA esterification in excess of ethanol 𝑘 𝐶 𝑊 Rate constant based on concentration for LA esterification in excess of water in Kamlet-Abboud-Taft equation for the solvent hydrogenbond donor acidity 𝛽 Solvatochromic solvent parameter in Kamlet-Abboud-Taft equation for solvent hydrogen-𝐷𝑀𝐹 Relative permittivity of H2SO4 in DMF 𝜀 𝑟,𝐷𝑀𝑆𝑂 Relative permittivity of H2SO4 in DMSO

Table 1 .

 1 Experimental matrix for kinetic experiments in GVL solvent.

	Run	Temp	[H2SO4]0	[EL]0	[LA]0	[Ethanol]0	[GVL]0	[Water]0
		°C	mol.L -1	mol.L -1	mol.L -1	mol.L -1	mol.L -1	mol.L -1
		80.0	0.04	0.00	2.27	9.01	1.89	0.00
		70.0	0.05	0.00	2.30	9.11	1.94	0.00
		60.0	0.09	0.00	2.30	9.12	2.03	0.00
		50.0	0.09	0.00	2.33	9.19	2.02	0.00
		70.0	0.04	0.00	3.39	6.88	2.13	0.00
		80.0	0.04	0.00	3.34	6.74	2.12	0.00
		60.0	0.04	0.00	2.33	9.21	1.93	0.00
		60.0	0.09	0.00	2.32	9.19	1.90	0.00
		60.0	0.09	0.00	2.30	9.21	1.83	0.54
	10	60.0	0.09	0.00	2.34	9.18	1.63	1.60
	11	60.0	0.09	0.07	3.32	9.20	0.85	0.00
	12	60.0	0.09	0.40	2.32	9.20	1.34	0.00
	13	60.0	0.09	0.00	2.32	9.19	1.94	0.00
	14	60.0	0.09	0.00	2.32	9.17	1.19	4.05
	15	80.0	0.07	0.00	2.72	8.18	1.91	0.00
	16	50.0	0.07	0.00	2.84	8.48	1.93	0.00

Table 2 .

 2 Experimental matrix for kinetic experiments in water solvent.

		Temp	[H2SO4]0	[EL]0	[LA]0	[Ethanol]0	[Water]0
	Run	°C	mol.L -1	mol.L -1	mol.L -1	mol.L -1	mol.L -1
	1W	60.0	0.09	0.00	2.33	9.21	10.28
	2W	60.0	0.04	0.00	2.33	9.23	10.30
	3W	50.0	0.09	0.00	2.33	9.23	10.72
	4W	80.0	0.04	0.00	2.33	9.23	9.22
	5W	70.0	0.04	0.00	3.33	6.83	11.58
	6W	50.0	0.07	0.00	2.69	8.06	11.05
	7W	80.0	0.04	0.00	3.37	6.72	11.79
	8W	80.0	0.07	0.00	2.23	8.88	10.92
	9W	70.0	0.00	0.00	3.41	6.82	11.84
	10W	60.0	0.09	0.07	2.25	8.99	10.85
	11W	60.0	0.09	0.00	2.28	9.09	10.93
	12W	70.0	0.09	0.00	3.39	6.78	11.80
	13W	50.0	0.04	0.00	2.31	9.23	10.93
	14W	80.0	0.09	0.00	2.67	8.02	11.23
	15W	50.0	0.09	0.00	2.76	8.29	11.35

Table 3 .

 3 Experimental matrix for kinetic experiments in ethanol solvent.

		Temp	[H2SO4]0	[EL]0	[LA]0	[Ethanol]0	[Water]0
	Run	°C	mol.L -1	mol.L -1	mol.L -1	mol.L -1	mol.L -1
	1E	50.0	0.09	0.00	2.88	11.51	0.00
	2E	60.0	0.09	0.00	2.85	11.38	0.00
	3E	70.0	0.09	0.00	2.82	11.25	0.00
	4E	60.0	0.04	0.00	2.86	11.41	0.00
	5E	60.0	0.07	0.00	2.76	11.40	0.00
	6E	50.0	0.07	0.00	2.88	11.53	0.00
	7E	50.0	0.07	0.00	3.39	10.48	0.00
	8E	50.0	0.07	0.00	4.45	8.83	0.00
	9E	60.0	0.090	0.00	2.85	11.38	0.00

Table 4 .

 4 Solvatochromic solvent parameters[START_REF] Islam | Kamlet Taft Parameters: A Tool to Alternate the Usage of Hazardous Solvent in Pharmaceutical and Chemical Manufacturing/Synthesis -A Gateway towards Green Technology[END_REF].

		π* [-]	α [-]	β [-]
	GVL	0.9	0	0.6
	Water	1.2	1.2	0.5
	Ethanol	0.6	0.9	0.8

performed for water solvent, GVL solvent and ethanol solvent. It was observed that the used procedure is repeatable, that the increase of temperature and catalyst leads to accelerate the reaction rate. In the absence of sulfuric acid, the kinetics was observed to be very slow and negligible in any solvents at 70°C. For these reasons, such experiments were discarded in the modeling stage.

  Solvent effect on the dissociation constant 𝑝𝐾 𝑎 𝑐The key to predict the dissociation constant in non-aqueous solutions is the equilibrium constant K th . The availability of the thermodynamically consistent equilibrium constant K th allows predicting solvent effects on K c values and thus on the dissociation equilibrium.The acidity of an acid is expressed as dissociation constant 𝑝𝐾 𝑎 𝑐 . Because of the complications of measuring the acidity in non-aqueous solvents or mixtures with organic solvents and because of the fact that such data is missing for acids in many organic solvents, it is necessary to reliably predict the experimentally accessible 𝑝𝐾 𝑎 𝑐 values. In order to apply this procedure to H 2 SO 4 , the 𝑝𝐾 𝑎 𝑐 values of H 2 SO 4 in water, 𝑝𝐾 𝑎

	3.2
		0.70
		0.60
		0.50
	[EL]/[LA]0	0.30 0.40
		0.20					Run 8: 0.00 mol/L of W
							Run 9: 0.54 mol/L of W
		0.10					Run 10: 1.6 mol/L of W
							Run 14: 4.05 mol/L of W
		0.00
		0.00	10.00	20.00	30.00	40.00	50.00	60.00	70.00
					Time/min

𝑐 (𝑠𝑜𝑙𝑣𝑒𝑛𝑡 = 𝑤𝑎𝑡𝑒𝑟, 𝑇 = 298.15 𝐾, 𝑝 = 1 𝑏𝑎𝑟) = -3

[START_REF] Kolthoff | Treatise on analytical chemistry[END_REF][START_REF] Williams | pKa Data compiled by R.Williams[END_REF]

,

  The dissociation of the second step is very low therefore, the focus is on the first dissociation step in organic media via ePC-SAFT approach. The resulting 𝑝𝐾 𝑎 𝑐 values in five different solvents at 298.15 K and 1 bar are shown in Fig.2.

	different solvents are in good agreement with literature data (DMSO 𝑝𝐾 𝑎,2,𝑒𝑃𝐶-𝑆𝐴𝐹𝑇 𝑐	= 16.3,
	𝑝𝐾 𝑎,2,𝐿𝑖𝑡. 𝑐	= 14.7 (Izutsu, 2009)) and DMF 𝑝𝐾 𝑎,2,𝑒𝑃𝐶-𝑆𝐴𝐹𝑇 𝑐	= 18.3, 𝑝𝐾 𝑎,2,𝐿𝑖𝑡. 𝑐	= 17.2 (Izutsu,
	2009)).
					2 𝑐 in

results to 0.0174, and is in proper agreement with the literature value of 0.0119 published by Sippola and Sippola and Taskinen

[START_REF] Sippola | Critical evaluation of the second dissociation constants for aqueous sulfuric acid over a wide temperature range[END_REF][START_REF] Sippola | Thermodynamic properties of aqueous sulfuric acid[END_REF]

. Further, the calculated dissociation constants for the second step 𝑝𝐾 𝑎,

Table 6 .

 6 Estimated values at Tref =64 °C and statistical data for esterification in GVL

			solvent.
		Units	Estimated value 95% marginal HPD intervals HDP in %
	𝑙𝑛(𝑘 𝑐 𝐺𝑉𝐿 (𝑇 𝑟𝑒𝑓 ))	L 2 .mol -2 .min -1	-2.22	0.05	2.40
	𝐺𝑉𝐿 𝐸 𝑎 𝑅 • 𝑇 𝑟𝑒𝑓	-	18.06	1.84	10.21
	𝑙𝑛(𝐾 𝑐 𝐺𝑉𝐿 (𝑇 𝑟𝑒𝑓 ))	-	0.39	0.11	28.17
	𝐺𝑉𝐿 𝛥𝐻 𝑟 𝑅 • 𝑇 𝑟𝑒𝑓	-	15.49	4.34	28.00

Table 7 .

 7 Estimated values at Tref =57 °C and statistical data for esterification in ethanol

			solvent.
		Units	Estimated value 95% marginal HPD intervals HDP in %
	𝑙𝑛(𝑘 𝑐 𝐸𝑡ℎ𝑎𝑛𝑜𝑙 (𝑇 𝑟𝑒𝑓 )) L 2 .mol -2 .min -1	2.91	0.06	1.95
	𝐸𝑡ℎ𝑎𝑛𝑜𝑙 𝐸 𝑎 𝑅 • 𝑇 𝑟𝑒𝑓	-	24.26	2.84	11.72
	𝑙𝑛(𝐾 𝑐 𝐸𝑡ℎ𝑎𝑛𝑜𝑙 (𝑇 𝑟𝑒𝑓 ))	-	1.49	0.43	28.81
	𝐸𝑡ℎ𝑎𝑛𝑜𝑙 𝛥𝐻 𝑟 𝑅 • 𝑇 𝑟𝑒𝑓	-	42.32	23.14	54.68

Table 8 .

 8 Table 8 shows the values of the estimated constants with a reference temperature of 64 °C for the kinetic constants. One can notice that the credible intervals for each estimated constant are narrow. The values of 𝑘 𝑐 𝑊 (𝑇 𝑟𝑒𝑓 ), 𝐸 𝑎 𝑊 , 𝐾 𝑐 𝑊 (𝑇 𝑟𝑒𝑓 ) and 𝛥𝐻 𝑟 𝑊 are displayed in Table S7.5. Esterification of LA by ethanol in water solvent is an endothermic system. Estimated values at Tref =64 °C and statistical data for esterification in water solvent.

	Table S7.6 shows that correlations between the estimated constants are low in general.
	Nevertheless, the correlations between 𝑙𝑛(𝑘 𝑐 𝑊 (𝑇 𝑟𝑒𝑓 )) and 𝑙𝑛(𝐾 𝑐 𝑊 (𝑇 𝑟𝑒𝑓 ));	𝛥𝐻 𝑟 𝑊 𝑅•𝑇 𝑟𝑒𝑓	and	𝑊 𝐸 𝑎 𝑅•𝑇 𝑟𝑒𝑓
	can be qualified as medium, but lower than 0.95. Hence, the medium correlation observed
	for these estimated parameters did not result in wrong identification.

Table 9 .

 9 Values of A0, s, a and b of KAT equation (Eq. (8)).

	Temperature/°C	A0	s	a	b
	50	-0.252	-7.413	2.062	6.536
	52	-0.319	-7.343	2.082	6.736
	54	-0.387	-7.274	2.101	6.936
	56	-0.032	-7.411	2.138	6.737
	58	-0.050	-7.366	2.159	6.887
	60	-0.069	-7.322	2.179	7.036
	62	-0.087	-7.278	2.200	7.183
	64	-0.105	-7.235	2.219	7.327
	66	-0.123	-7.192	2.239	7.471
	68	-0.140	-7.150	2.258	7.612
	70	-0.158	-7.108	2.278	7.752

* + 𝑎 • 𝛼 + 𝑏 • 𝛽 for each solvent and different temperatures (Table

9

).

  𝐴 0Regression value for Kamlet-Abboud-Taft equation 𝑎 Regression parameter for Kamlet-Abboud-Taft equation

	𝑏	Regression parameter for Kamlet-Abboud-Taft equation
	𝐶 𝐿𝐴	Concentration of LA
	𝐶 𝐸𝑡ℎ𝑎𝑛𝑜𝑙	Concentration of ethanol
	𝐶 𝐸𝐿	Concentration of ethyl levulinate
	𝐶 𝑊	Concentration of water
	𝑐 𝑝

𝑂

Specific standard heat capacity used in Eq.17

Acknowledgments

The authors thank INSA Rouen Normandy, University of Rouen Normandy, the Centre National de la Recherche Scientifique (CNRS), European Regional Development Fund (ERDF) N• HN0001343, Labex SynOrg (ANR-11-LABX-0029), Carnot Institute I2C, the graduate school for reasearch XL-Chem (ANR-18-EURE-0020 XL CHEM) and Region Normandie for their support. GC/FID was financed by ERDF RIN Green Chem 2019NU01FOBC08 N• 17P04374. The authors thank the Maîtrise des Risques et Environementaux department, and the Erasmus programme to make the research projects of Sindi Baco and Rashid Ismail Bedawi Zakaria possible. This research was funded, in whole or in part, by the ANR (French National Research Agency) and the DFG (German Research Foundation) through the project MUST (MicroflUidics for Structure-reactivity relationships aided by Thermodynamics & kinetics) [ANR-20-CE92-0002-01 -Project number 446436621].

Abbreviations

Declaration of Competing Interest

The authors report no declarations of interest.