Daniele Di Menno

Di Bucchianico

Antonella Cipolla

Jean-Christophe Buvat

Mélanie Mignot

Valeria Casson Moreno

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

Daniele Di

Menno Di Bucchianico

Valeria Casson Moreno

Kinetic Study and Model Assessment for n -Butyl Levulinate Production from Alcoholysis of 5-(Hydroxymethyl)furfural over Amberlite IR-120

Keywords: Lignocellulosic biomass, Cellulose valorization, Butyl levulinate production, 5-(hydroxymethyl) furfural) alcoholysis, Amberlite IR-120, Humins production, Kinetic modeling, K-fold method

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Fossil raw materials in modern society plays a crucial role. However, using such raw materials has several drawbacks: non-renewable, difficult to extract (i.e., not full extraction), resources are not well-distributed at global scale and negative environmental impacts such as perturbation in the carbon cycle (released of greenhouse gas) or oil spills.

In order to overcome these drawbacks, the shift from fossil non-renewable raw materials to biomass renewable raw materials can aid our modern societies to continue their developments with a minimum impact on the environment. Biomass valorization faces several challenges such as biomass collection network, catalyst development, diversity of biomass content, oxygen and water presence, pretreatment, 1 etc. The use of biomass for chemical production is not new in chemical processes, one can cite the production of alcohol from cereals or sugar juice. Nevertheless, the use of such raw materials competes with food. To avoid the food versus fuel dilemma, agricultural waste or non-alimentary plants should be favored.

Lignocellulosic biomass can avoid such a dilemma and enhance the circular economy. [START_REF] Culaba | Smart Sustainable Biorefineries for Lignocellulosic Biomass[END_REF][START_REF] Ning | Recent Advances in the Valorization of Plant Biomass[END_REF][START_REF] Singhvi | Lignocellulosic Biomass: Hurdles and Challenges in Its Valorization[END_REF] Lignocellulosic biomass has a fascinating structure made of lignin (natural aromatic polymer) and cellulose/hemicellulose which are sugar polymers. The composition of these elements is different for each plant, implying that operating conditions should be tuned for each biomass.

The first stage is to separate lignin, cellulose and hemicellulose, i.e., the fractionation stage.

The valorization of lignin is still under development for industrial application. [START_REF] Lu | Reductive Catalytic Depolymerization of Semi-Industrial Wood-Based Lignin[END_REF][START_REF] Wang | From Lignin to Valuable Products-Strategies, Challenges, and Prospects[END_REF][START_REF] Sun | Bright Side of Lignin Depolymerization: Toward New Platform Chemicals[END_REF][START_REF] Schutyser | Chemicals from Lignin: An Interplay of Lignocellulose Fractionation, Depolymerisation, and Upgrading[END_REF] The valorization of cellulose and hemicellulose is more developed, [START_REF] Nzediegwu | Chemo-Catalytic Transformation of Cellulose and Cellulosic-Derived Waste Materials into Platform Chemicals[END_REF][START_REF] Chen | Recent Progress in Homogeneous Lewis Acid Catalysts for the Transformation of Hemicellulose and Cellulose into Valuable Chemicals, Fuels, and Nanocellulose[END_REF][START_REF] Šivec | Multiscale Modeling of (Hemi)Cellulose Hydrolysis and Cascade Hydrotreatment of 5-Hydroxymethylfurfural, Furfural, and Levulinic Acid[END_REF][START_REF] Dusselier | Top Chemical Opportunities from Carbohydrate Biomass: A Chemist's View of the Biorefinery[END_REF][START_REF] Delidovich | Cellulose and Hemicellulose Valorisation: An Integrated Challenge of Catalysis and Reaction Engineering[END_REF] notably for the industrial production of cellulosic ethanol.

From the cellulose or hemicellulose valorization, it is possible to develop several platform molecules such as: 1,4-diacid, 5-(hydroxymethyl)furfural (5-HMF), 2,5-furandicarboxylic acid (2,5-FDCA), aspartic acid, glutamic acid, glucaric acid, itaconic acid, sorbitol, acetone-butanolethanol (ABE process) or levulinic acid [START_REF] Takkellapati | An Overview of Biorefinery-Derived Platform Chemicals from a Cellulose and Hemicellulose Biorefinery[END_REF][START_REF] Isikgor | Lignocellulosic Biomass: A Sustainable Platform for the Production of Bio-Based Chemicals and Polymers[END_REF][START_REF] Kohli | Bio-Based Chemicals from Renewable Biomass for Integrated Biorefineries[END_REF] . Some of them are already at the industrial scale, such as levulinic acid one. [START_REF] Bozell | Cellulose and Renewable Materials Division[END_REF][START_REF] Kapanji | Techno-Economics of Lignocellulose Biorefineries at South African Sugar Mills Using the Biofine Process to Co-Produce Levulinic Acid, Furfural and Electricity along with Gamma Valeractone[END_REF] like γ-valerolactone [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] Capecci | Bayesian Statistics to Elucidate the Kinetics of γ-Valerolactone from n-Butyl Levulinate Hydrogenation over Ru/C[END_REF][START_REF] Casson Moreno | A Consequences-Based Approach for the Selection of Relevant Accident Scenarios in Emerging Technologies[END_REF] , succinic acid, [START_REF] Carnevali | Levulinic Acid Upgrade to Succinic Acid with Hydrogen Peroxide[END_REF][START_REF] Podolean | Ru-Based Magnetic Nanoparticles (MNP) for Succinic Acid Synthesis from Levulinic Acid[END_REF][START_REF] Adeleye | A Review on the Conversion of Levulinic Acid and Its Esters to Various Useful Chemicals[END_REF] 2-methylene γ-valerolactone [START_REF] Al-Naji | Sustainable Continuous Flow Valorization of γ-Valerolactone with Trioxane to α-Methylene-γ-Valerolactone over Basic Beta Zeolites[END_REF] or 2methyltetrahydrofuran (2-MTHF). [START_REF] Phanopoulos | Catalytic Transformation of Levulinic Acid to 2-Methyltetrahydrofuran Using Ruthenium -N-Triphos Complexes[END_REF][START_REF] Licursi | Cascade Strategy for the Tunable Catalytic Valorization of Levulinic Acid and γ-Valerolactone to 2-Methyltetrahydrofuran and Alcohols[END_REF][START_REF] Obregõn | One-Pot 2-Methyltetrahydrofuran Production from Levulinic Acid in Green Solvents Using Ni-Cu/Al2O3 Catalysts[END_REF] The classical production of LA is from the hydrolysis of cellulose using homogeneous catalysts. [START_REF] Tan-Soetedjo | Experimental and Kinetic Modeling Studies on the Conversion of Sucrose to Levulinic Acid and 5-Hydroxymethylfurfural Using Sulfuric Acid in Water[END_REF][START_REF] Khan | Efficient Conversion of Lignocellulosic Biomass to Levulinic Acid Using Acidic Ionic Liquids[END_REF][START_REF] Dussan | Kinetics of Levulinic Acid and Furfural Production from Miscanthus × Giganteus[END_REF][START_REF] Liang | Production of Levulinic Acid from Corn Cob Residue in a Fed-Batch Acid Hydrolysis Process[END_REF] Nevertheless, the low yields in LA production, together with high costs in separation and purification of this acid, make the direct production of alkyl levulinates more beneficial via sugar alcoholysis [START_REF] Menno | Production of Levulinic Acid and Alkyl Levulinates: A Process Insight[END_REF] .

Among alkyl levulinates, methyl, ethyl and butyl levulinates are the most studied and used in the literature. As shown in the review of Di Menno Di Bucchianico et al., [START_REF] Menno | Production of Levulinic Acid and Alkyl Levulinates: A Process Insight[END_REF] alcoholysis of sugars over heterogenous acid catalysts should be favored for the production of alkyl levulinates compared to the traditional homogeneous acids. Zhou et al., [START_REF] Zhou | Conversion of Carbohydrate Biomass to Methyl Levulinate with Al 2(SO4)3 as a Simple, Cheap and Efficient Catalyst[END_REF] reported the production of methyl levulinate from glucose, catalyzed by aluminum sulfate at 160 °C for 3 h, obtaining a production yield of 64%. Saravanamurugan and Riisager [START_REF] Saravanamurugan | Zeolite Catalyzed Transformation of Carbohydrates to Alkyl Levulinates[END_REF] studied the production of ML from fructose over zeolite (H-USY) at 160 °C for 20 h with a final ML yield of 51 %. The production of methyl levulinate was also investigated starting from 5-HMF by Hu et al., [START_REF] Hu | One-Pot Conversion of Biomass-Derived Xylose and Furfural into Levulinate Esters via Acid Catalysis[END_REF] using Amberlyst 70 at 120 °C for 2 h and resulting in a final yield of 61.7 %. Wang et al., [START_REF] Wang | Molecular Oxygen-Promoted Synthesis of Methyl Levulinate from 5-Hydroxymethylfurfural[END_REF] reported the total conversion of 5-HMF and the final yield in ML of 69.6 % over the zeolite H-beta-40. From a green-chemistry aspect, ethyl and butyl levulinates production can be considered greener because ethanol and butanol are produced from renewable raw materials. Ethyl levulinate production has been studied over a wide range of heterogeneous catalysts and starting materials. Liu et al., [START_REF] Liu | Conversion of Fructose into 5-Hydroxymethylfurfural and Alkyl Levulinates Catalyzed by Sulfonic Acid-Functionalized Carbon Materials[END_REF] have investigated the conversion of fructose over Amberlyst-15 at 120 °C for 24 h, obtaining a final yield of 73 % for ethyl levulinate and 89 % for butyl levulinate.

Considering different feedstocks, Quereshi et al., [START_REF] Quereshi | Insights into the Metal Salt Catalyzed Ethyl Levulinate Synthesis from Biorenewable Feedstocks[END_REF] showed the increase in EL yield passing from glucose (23 %), to fructose (25 %) and to 5-HMF (37.2 %) over a metal salt. 5-HMF alcoholysis in ethanol was also studied on zeolite (MZSM-5) by Chithra and Srinivas 40 with a final EL yield of 85.2 % at 150 °C for 12 h. Studies on the production of butyl levulinate are mainly from sugar alcoholysis, as the investigation of Di Menno Di Bucchianico et al., [START_REF] Menno | Role of Solvent the Production of Butyl Levulinate from Fructose[END_REF] on fructose conversion over Amberlite IR120 with a final yield 60.4 % at 110 °C for 7 h or the work of An et al., [START_REF] An | Efficient One-Pot Synthesis of n-Butyl Levulinate from Carbohydrates Catalyzed by Fe2(SO4)3[END_REF] with Fe2(SO4)3 reporting a final BL yield of 62.8 % at 190 °C. Fewer literature data are on the production of butyl levulinate, although, from an energy point of view, the use of butyl levulinate as blend in fuels has certain advantages compared to ethyl levulinate in terms of conductivity, cold flow properties, diesel lubrication, reduced vapour pressure, corrosivity and miscibility with fuels even at very low temperatures. [START_REF] Menno | Role of Solvent the Production of Butyl Levulinate from Fructose[END_REF] Some kinetic models of sugar hydrolysis to levulinic acid and esterification (levulinic acid to alkyl levulinates and formic acid to alkyl formate) can be found in the literature [START_REF] Tan-Soetedjo | Experimental and Kinetic Modeling Studies on the Conversion of Sucrose to Levulinic Acid and 5-Hydroxymethylfurfural Using Sulfuric Acid in Water[END_REF][START_REF] Kumar | 2nd Generation Biomass Derived Glucose Conversion to 5-Hydroxymethylfurfural and Levulinic Acid Catalyzed by Ionic Liquid and Transition Metal Sulfate: Elucidation of Kinetics and Mechanism[END_REF][START_REF] Ariba | Physicochemical Properties for the Reaction Systems: Levulinic Acid, Its Esters, and γ-Valerolactone[END_REF][START_REF] Konaka | A Study of Reaction Rates. Esterification of Formic Acid with Ethanol[END_REF][START_REF] Bart | Kinetics of Esterification of Succinic Anhydride with Methanol by Homogeneous Catalysis[END_REF][START_REF] Kokare | Response Surface Optimization, Kinetic Study and Process Design of n-Butyl Levulinate Synthesis[END_REF][START_REF] Van Putten | a Versatile Platform Chemical Made from Renewable Resources[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] Sacia | Biomass Conversion to Diesel via the Etherification of Furanyl Alcohols Catalyzed by Amberlyst-15[END_REF] . However, the development of kinetic models for the alcoholysis of sugars to alkyl levulinates is scarce, and the kinetic constant comparison is difficult due to the different ways to derive reaction rates. [START_REF] Wang | Experimental and Theoretical Studies on Glucose Conversion in Ethanol Solution to 5-Ethoxymethylfurfural and Ethyl Levulinate Catalyzed by a Brønsted Acid[END_REF] The kinetic modeling of this system is complex because several reactions occur simultaneously in parallel and consecutive ways (Figure 1). Besides, humins production is complex to track. It is vital to have reliable kinetic and thermodynamic models to draw process flow diagrams for alkyl levulinate production. From the knowledge of these models, it is possible to determine the optimum reaction conditions.

The kinetics of 5-HMF alcoholysis is in central position in the reaction scheme (Figure 1), because it can lead to the production of LA or BL. By considering solely the butanolysis of 5-HMF, one can evaluate the production of humins from 5-HMF at different operating conditions.

To the best of our knowledge, such model was not developed for the production of butyl levulinate over heterogeneous catalysts.

This manuscript proposes to study the kinetics of 5-HMF butanolysis over Amberlite IR-120. Indeed, Ramirez et al., [START_REF] Quereshi | Insights into the Metal Salt Catalyzed Ethyl Levulinate Synthesis from Biorenewable Feedstocks[END_REF] showed the benefit of using cation exchange resin to produce BL from fructose butanolysis. Also, Di Menno Di Bucchianico et al., [START_REF] Menno | Role of Solvent the Production of Butyl Levulinate from Fructose[END_REF] showed the benefits of using Amberlite IR-120 and γ-valerolactone (GVL) as a co-solvent for the kinetics of fructose butanolysis. The second stage of this manuscript is the development of kinetic models to explain the reaction mechanism. Special attention was paid to vary widely the 5-HMF loading, catalyst loading and reaction temperature. The knowledge of a kinetic model for concentrated 5-HMF solution is important for the alcoholysis of concentrated sugar solution. Besides, independent esterification experiments were performed to ease the estimation stage. A crossvalidation stage was included to evaluate the developed models. 5-(hydroxymethyl)furfural (5-HMF, 99% purity), 5-(ethoxymethyl)furfural (EMF, 97% purity), levulinic acid (LA, ≥ 97% purity), formic acid (FA, ≥ 95% purity) and γ-valerolactone (≥ 99% purity) were purchased from Sigma-Aldrich. 1-Butanol (BuOH, ≥ 99.5 % purity), butyl levulinate (BL, ≥ 98% purity), butyl formate (BF, ≥ 97% purity) and acetone (≥ 99.9% purity) from VWR chemicals. Amberlite IR120 H + form (ion-exchange resin, native harmonic mean size: 0.620 to 0.830 mm) commercial catalyst was provided by Acros Organics. Nitrogen gas (N2 purity > 99.999 vol%) came from Linde. All chemicals were employed without further purification.

Analytical methods

Reaction samples were diluted in acetone and analysed by gas chromatograph Bruker Scion 456-GC, equipped with a VF-1701ms Agilent column (60.0 m length x 250 µm inner diameter x 0.25 µm film thickness) and a flame ionization detector (FID). The injector was set with a temperature of 250 °C, injection volume of 1 µL and split ratio of 1:20. Helium (>99.999 vol%) is used as carrier gas at a constant flow rate of 1.0 mL/min to transfer the sample from the injector, through the column and into the FID-detector, whose temperature was set to 250 °C.

The oven temperature was programmed from 40 °C to 240 °C with 20 °C/min of ramp rate.

Reference samples with standard solutions of pure 5-HMF, LA, EMF, BL and BF were used as standards to prepare daily calibration curves. EMF was used as standard for BMF due to its commercial unavailability. A repetition of three injections for each sample was considered to estimate the uncertainty, measured by the standard deviation, in the analytical method. In this study, the standard deviation for the different concentrations was lower than 0.03 mol/L.

Kinetic experiments

Under isothermal and isobaric conditions, kinetic experiments were carried out in a 300 mL Parr stainless steel batch reactor, equipped with an electric heating system. The presence of a gas entrainment impeller (diameter 2.5 cm) with a hollow shaft provided a uniform mixing of the reacting mixture. As shown in Table 1, the reactor was loaded with HMF initial mass varying between 1.6 and 7 g, catalyst mass up to 15.3 g and a fixed solvent volume of butanol-GVL with ratio 70:30 wt%. [START_REF] Menno | Role of Solvent the Production of Butyl Levulinate from Fructose[END_REF] In runs 10-16, levulinic acid and levulinic-formic acid mixture esterification (initial concentration range 0.05 -0.5 mol/L) were performed in the same solvent system and by using the same catalyst, i.e, Amberlite IR-120. Nitrogen was used to pressurize the system at 20 bars to limit gas-liquid partition of the liquid phase. [START_REF] Menno | Role of Solvent the Production of Butyl Levulinate from Fructose[END_REF] Rotation speed was set to 800 rpm, since previous studies on fructose butanolysis by using Amberlite IR120 in the same apparatus have demonstrated the negligibility of external and internal mass transfer limitations by using the native catalyst particle size distribution at this rotation speed. [START_REF] Menno | Role of Solvent the Production of Butyl Levulinate from Fructose[END_REF] Amberlite IR120 was pre-treated as described by Di Menno Di Bucchianico et al. [START_REF] Menno | Role of Solvent the Production of Butyl Levulinate from Fructose[END_REF] The catalyst was washed with distilled water and, then with pure butanol, before being dried in an oven at 90 °C for 5 hours under atmospheric pressure. For each kinetic experiment (Table 1), we used fresh catalysts. Catalyst reusability was also tested by re-treating the used-catalyst, as described above, before reusing it in a second run. In kinetic runs, temperature was set between 80 and 115 °C, without exceeding 120 °C. Indeed, sulfonic active sites from Amberlite IR-120 can leach for temperature higher than 120°C. During the heating phase, the temperature ramp of the system has been recorded for each experiment. After a first sample collected at ambient temperature to measure the initial concentration of the loaded reagents, intermediate samples were collected to monitor the catalytic reaction, since the catalyst was loaded into the system at the beginning and therefore in contact with the reagents during the heating step. Once the set temperature was reached, a sample was taken and the next ones at 5 min, 30 min, and then every hour for up to 7 h.

Results

Repeatability

Besides analyzing each sample three times, two experiments were repeated: 1 and 6. Figures S1. 1 shows the normalized ratio of each species to the initial HMF concentration calculated at initial reaction temperature (ambient temperature) versus time (

[𝐻𝑀𝐹] [𝐻𝑀𝐹] 0 , [𝐿𝐴] [𝐻𝑀𝐹] 0 , [𝐵𝑀𝐹] [𝐻𝑀𝐹] 0 and [𝐵𝐿] [𝐻𝑀𝐹] 0
) for Experiment 1 and 1-repeated. One can notice that kinetic profiles are similar for both experiments. Also, LA concentration is lower than BL or BMF concentrations. For the ratio

[𝐵𝑀𝐹] [𝐻𝑀𝐹] 0
, one can notice that the kinetics of BMF generation is slightly affected. and one can notice that the kinetics of BL and BF are similar (Figure 3F). This similarity is higher when the temperature increases, meaning that the production of LA from HMF can be considered negligible. Mass balance on the initial concentration of HMF shows a loss of HMF to humins, which varies from 15 to 20%. F.

Nevertheless, the ratio

Effect of 5-HMF loading

To evaluate the effect of HMF loading on the kinetics, Runs 2, 5 and 17 were compared. In these experiments, the molar ratio nHMF0/nH+,apparent were kept the same. Figure 4A shows that the kinetics of 5-HMF consumption is slightly faster for high HMF loading, i.e., Run 17.

Interestingly, this higher amount of 5-HMF leads to a higher production of LA in Run 17 (Figure 4B), which is usually lower. This significant production of LA decreases the linearity/parity between BL and BF concentrations in Run 17 (Figure 4E). High HMF loading increases the BMF concentration, leading to a higher BL production in Run 17. The linearity decrease between BL and BF in Run 17 is shown in Figures 4D and4E, where the ratio

Effect of catalyst loading

To evaluate the catalyst loading effect, i.e., From Figures 5, one can observe that the generation kinetics increase from 0 to 4.9 gram of dried catalyst. Then, from 4.9 to 9.2 grams of catalyst, there is no significant effect on the rate of HMF consumption or BL production. For Run 5, the amount of protons is two times higher than the amount of HMF; in other words, more acid sites are available than chemical reactions. The blank experiment (Run 20) shows that during the first 50 minutes, HMF is consumed, but then the kinetics of decomposition is slower. These two kinetic regimes might be linked to the production of humins and BMF. One cannot say that there are no chemical reactions without acid catalysts. Nevertheless, this production could be initiated by some small amount of proton and then by LA production. Figure 5F shows that in the absence of catalysts, the linearity between BF and BL decreases, which might be due to a longer lifetime of LA and FA. Figures

𝑛

Esterification reactions

The hydrolysis of HMF leads to the production of LA and FA. These species can be esterified into BL and BF. It is very challenging to track the concentrations of LA and FA during the alcoholysis because their concentrations can be low, and FA is volatile. To overcome this issue, we performed independent experiments of esterification using the same solvent environment and LA and FA as starting materials. Figures 6 show the effect of temperature on the kinetics of esterification. One can notice that in these conditions, the esterification reactions can be considered irreversible because no water was added to the reaction mixture. The kinetics of FA esterification is faster than LA esterification. Interestingly, under these operating conditions, the temperature increase does not significantly affect the rate of esterification for FA esterification.

Discussion on kinetics

The reaction mechanism of 5-HMF hydrolysis is still under debate. This difficulty is increased by the side-reaction of humin production, where the hydroxyl group of 5-HMF plays an important role. [START_REF] Van Zandvoort | Towards the Valorization of Humin By-Products : Characterization , Solubilization and Catalysis[END_REF][START_REF] Zheng | Production of 5-Hydroxymethyl Furfural and Levulinic Acid from Lignocellulose in Aqueous Solution and Different Solvents. Biofuels, Bioprod[END_REF][START_REF] Ramírez | Assessment of Ion Exchange Resins as Catalysts for the Direct Transformation of Fructose into Butyl Levulinate[END_REF] The solvent role also interferes with the kinetics and thermodynamics of this reaction system and its capacity to dissolve 5-HMF or humins. [START_REF] Zhang | Synthesis of γ-Valerolactone from Carbohydrates and Its Applications[END_REF] Online analysis to track the different intermediates of this reaction system is cumbersome, so some research groups have used DFT calculation to unravel some elementary steps. [START_REF] Wang | Experimental and Theoretical Studies on Glucose Conversion in Ethanol Solution to 5-Ethoxymethylfurfural and Ethyl Levulinate Catalyzed by a Brønsted Acid[END_REF][START_REF] Nikbin | A First Principles-Based Microkinetic Model for the Conversion of Fructose to 5-Hydroxymethylfurfural[END_REF][START_REF] Yang | Mechanism of Brønsted Acid-Catalyzed Conversion of Carbohydrates[END_REF] From the Results section, one can consider that all esterification reactions are irreversible and humins production must be considered in these operating conditions. The fact that esterification are irreversible in this system is due to the very low concentration of water, i.e., no water was added. We assumed that the etherification of HMF to BMF is an irreversible reaction.

Figure 7 shows a simplified reaction scheme for the production of BL from 5-HMF alcoholysis.

Figure 7. Simplified reaction steps for the butanolysis of 5-HMF.

Based on the work of Wang et al., [START_REF] Wang | Experimental and Theoretical Studies on Glucose Conversion in Ethanol Solution to 5-Ethoxymethylfurfural and Ethyl Levulinate Catalyzed by a Brønsted Acid[END_REF] we have proposed another reaction scheme illustrated in Figure 8. We assumed that there is another rate-determining step before the production of LA and BL. Another possible side reaction is the etherification of butanol to dibutyl ether [START_REF] Iborra | Liquid-Phase Synthesis of Butyl Levulinate with Simultaneous Water Removal Catalyzed by Acid Ion Exchange Resins[END_REF] , but such a reaction was not observed in this study. For humins production from 5-HMF, several authors expressed its reaction rate as a first-order according to 5-HMF concentration. If humins production is a polymerization reaction system, this first-order might be wrong. In this manuscript, both orders will be evaluated.

All the reaction steps were considered to be catalyzed by the protons from Amberlite IR-120 (Figures 5). The catalytic effect due to LA was considered to be negligible due to the low concentration of LA. A pseudo-homogeneous approach was used to describe the role of proton from the resins. [START_REF] Musante | Kinetic of Liquid-Phase Reactions Catalyzed by Acidic Resins: The Formation of Peracetic Acid for Vegetable Oil Epoxidation[END_REF] In other words, protons from the sulfonic groups are considered to have a high degree of freedom. The developed models did not consider the adsorption and desorption steps on sulfonic groups. The proton concentration was calculated based on the acid capacity of Amberlite IR-120. [START_REF] Leveneur | Green Process Technology for Peroxycarboxylic Acids: Estimation of Kinetic and Dispersion Parameters Aided by RTD Measurements: Green Synthesis of Peroxycarboxylic Acids[END_REF][START_REF] Leveneur | Modeling of a Liquid-Liquid-Solid Heterogeneous Reaction System: Model System and Peroxyvaleric Acid[END_REF][START_REF] Leveneur | Synthesis of Peroxypropionic Acid from Propionic Acid and Hydrogen Peroxide over Heterogeneous Catalysts[END_REF] [𝑃𝑟𝑜𝑡. ] =

𝑚 𝑑𝑟𝑖𝑒𝑑 𝑐𝑎𝑡𝑎𝑙𝑦𝑠𝑡 •𝐶𝑎𝑝𝑎𝑐𝑖𝑡𝑦 ( 𝑚𝑜𝑙 𝑜𝑓𝑃𝑟𝑜𝑡𝑜𝑛 𝑚 𝑑𝑟𝑖𝑒𝑑 𝑐𝑎𝑡𝑎𝑙𝑦𝑠𝑡 ) 𝑉 𝑅𝑒𝑎𝑐𝑡𝑖𝑜𝑛 Eq. (1)
The volume of reaction VReaction was calculated based on the total mass ratio to the density of butanol at the corresponding temperature. [START_REF] Ariba | Physicochemical Properties for the Reaction Systems: Levulinic Acid, Its Esters, and γ-Valerolactone[END_REF] The effect of sampling was taken into account.

Models

Four models were developed based on Figures 6 and7, and by testing first and second order for humins production rate.

A. Model 1

From Figures 6, each step was considered an elementary one; thus, the reaction rates can be expressed as

𝑅 1 = 𝑘 1 • [𝐻𝑀𝐹] • [𝑃𝑟𝑜𝑡. ]
Eq. ( 2)

𝑅 2 = 𝑘 2 • [𝐻𝑀𝐹] • [𝐵𝑢𝑂𝐻] • [𝑃𝑟𝑜𝑡. ]
Eq. ( 3)

𝑅 3 = 𝑘 3 • [𝐵𝑀𝐹] • [𝐵𝑢𝑂𝐻] • [𝑃𝑟𝑜𝑡. ]
Eq. ( 4)

𝑅 4 = 𝑘 4 • [𝐿𝐴] • [𝐵𝑢𝑂𝐻] • [𝑃𝑟𝑜𝑡. ]
Eq. ( 5)

𝑅 5 = 𝑘 5 • [𝐹𝐴] • [𝐵𝑢𝑂𝐻] • [𝑃𝑟𝑜𝑡. ]
Eq.( 6)

𝑅 6 = 𝑘 6 • [𝐻𝑀𝐹] • [𝑃𝑟𝑜𝑡. ]
Eq.( 7)

Water was not added to the reaction mixture but was present in the atmosphere at an infinitesimal concentration and produced during the reaction. Thus, the quantification of water in this reaction system is not possible. For that reason, water concentration was merged into the rate constant k1. The degradation of HMF into humins was assumed to be of first order.

Ideal flow conditions were assumed for the batch reactor. Thus, material balances can be derived as

𝑑[𝐻𝑀𝐹] 𝑑𝑡 = -𝑅 1 -𝑅 2 -𝑅 6
Eq.( 8)

𝑑[𝐵𝑀𝐹] 𝑑𝑡 = 𝑅 2 -𝑅 3
Eq.( 9)

𝑑[𝐿𝐴] 𝑑𝑡 = 𝑅 1 -𝑅 4
Eq.( 10)

𝑑[𝐵𝐿] 𝑑𝑡 = 𝑅 3 +𝑅 4
Eq.( 11)

𝑑[𝐵𝑢𝑂𝐻] 𝑑𝑡 = -𝑅 2 -2 • 𝑅 3 -𝑅 4 -𝑅 5
Eq.( 12)

𝑑[𝐹𝐴] 𝑑𝑡 = 𝑅 1 -𝑅 5
Eq.( 13)

𝑑[𝐵𝐹] 𝑑𝑡 = 𝑅 3 + 𝑅 5
Eq.( 14)

𝑑[𝐻𝑢𝑚𝑖𝑛𝑠] 𝑑𝑡 = 𝑅 6
Eq.( 15)

B. Model 2

The difference with Model 1 is that reaction R6 was considered to be of second-order with respect to HMF concentration.

𝑅 6 = 𝑘 6 • [𝐻𝑀𝐹] 2 • [𝑃𝑟𝑜𝑡. ] Eq.( 16
)
Material balance is different for HMF and is expressed as

𝑑[𝐻𝑀𝐹] 𝑑𝑡 = -𝑅 1 -𝑅 2 -2 • 𝑅 6
Eq.( 17)

C. Model3

From Figure 7, each step was considered an elementary one. Thus, the reaction rates can be expressed as

𝑅 1 = 𝑘 1 • [𝐻𝑀𝐹] • [𝑃𝑟𝑜𝑡. ]
Eq. ( 18)

𝑅 2 = 𝑘 2 • [𝐼𝑁𝑇1] • [𝑃𝑟𝑜𝑡. ]
Eq. ( 19)

𝑅 3 = 𝑘 3 • [𝐻𝑀𝐹] • [𝐵𝑢𝑂𝐻] • [𝑃𝑟𝑜𝑡. ]
Eq. ( 20)

𝑅 4 = 𝑘 4 • [𝐵𝑀𝐹] • [𝐵𝑢𝑂𝐻] • [𝑃𝑟𝑜𝑡. ]
Eq. ( 21)

𝑅 5 = 𝑘 5 • [𝐼𝑁𝑇1] • [𝐵𝑢𝑂𝐻] • [𝑃𝑟𝑜𝑡. ]
Eq.( 22)

𝑅 6 = 𝑘 6 • [𝐹𝐴] • [𝐵𝑢𝑂𝐻] • [𝑃𝑟𝑜𝑡. ]
Eq.( 23)

𝑅 7 = 𝑘 7 • [𝐿𝐴] • [𝐵𝑢𝑂𝐻] • [𝑃𝑟𝑜𝑡. ]
Eq.( 24)

𝑅 8 = 𝑘 8 • [𝐻𝑀𝐹] • [𝑃𝑟𝑜𝑡. ]
Eq.( 25)

The degradation of HMF into humins was assumed to be of first-order.

Material balances can be derived as

𝑑[𝐻𝑀𝐹] 𝑑𝑡 = -𝑅 1 -𝑅 3 -𝑅 8
Eq.( 26)

𝑑[𝐵𝑀𝐹] 𝑑𝑡 = 𝑅 3 -𝑅 4
Eq.( 27)

𝑑[𝐼𝑁𝑇1] 𝑑𝑡 = 𝑅 1 -𝑅 2 + 𝑅 4 -𝑅 5
Eq.( 28)

𝑑[𝐿𝐴] 𝑑𝑡 = 𝑅 2 -𝑅 7
Eq.( 29)

𝑑[𝐵𝐿] 𝑑𝑡 = 𝑅 5 +𝑅 7
Eq.( 30)

𝑑[𝐵𝑢𝑂𝐻] 𝑑𝑡

= -𝑅 3 -𝑅 4 -𝑅 5 -𝑅 6 -𝑅 7 Eq.( 31)

𝑑[𝐹𝐴] 𝑑𝑡

= 𝑅 1 -𝑅 6 Eq.( 32)

𝑑[𝐵𝐹] 𝑑𝑡 = 𝑅 4 + 𝑅 6 Eq.( 33
𝑑[𝐻𝑢𝑚𝑖𝑛𝑠] 𝑑𝑡 = 𝑅 8
Eq.( 34)

D. Model4

The difference with Model 3 is that reaction R8 was considered to be of second order with respect to 5-HMF concentration.

𝑅 8 = 𝑘 8 • [𝐻𝑀𝐹] 2 • [𝑃𝑟𝑜𝑡. ] Eq.( 35
)
Material balance is different for HMF and can be expressed as

𝑑[𝐻𝑀𝐹] 𝑑𝑡 = -𝑅 1 -𝑅 3 -2 • 𝑅 8
Eq.( 36)

Modeling

Simulation and parameter estimation stages were carried out by commercial software Athena Visual Studio. [START_REF] Stewart | Athena Visual Studio[END_REF][START_REF] Stewart | Computer-Aided Modeling of Reactive Systems[END_REF] The experimental concentrations of HMF, LA, BMF, BF and BL were used as observables. For such multi-response parameter estimation case, the Bayesian framework is more suitable than the non-linear least squares approach. [START_REF] Kopyscinski | Comprehensive Kinetic Study for Pyridine Hydrodenitrogenation on (Ni)WP/SiO2 Catalysts[END_REF][START_REF] Stewart | Parameter Estimation from Multiresponse Data[END_REF] Indeed, the minimization of the objective function requires the determination of the determinant criterion. [START_REF] Boekel | Statistical Aspects of Kinetic Modeling for Food Science Problems[END_REF] The Ordinary Differential Equations (ODEs) issued from material balances were integrated by the DDAPLUS solver, included in Athena Visual Studio software. This solver is developed based on a modified Newton algorithm. [START_REF] Caracotsios | Sensitivity Analysis of Initial Value Problems with Mixed Odes and Algebraic Equations[END_REF] The minimization of the objective function (OF), the determination of the credible intervals for each estimated parameter and the calculation of the normalized parameter covariance were done by the subroutine GREGPLUS.

GREGPLUS uses successive quadratic programming to minimize OF. [START_REF] Stewart | Computer-Aided Modeling of Reactive Systems[END_REF][START_REF] Stewart | Parameter Estimation from Multiresponse Data[END_REF] 𝑂𝐹 = (𝑎 + 𝑏 + 1) • 𝑙𝑛|𝜐| Eq. (37) where, 𝑏 is the number of responses, 𝑎 is the number of events in response, and |𝜐| is the determinant of the covariance matrix of the responses. Each element of this matrix is

𝜐 𝑖𝑗 = ∑ [𝐶 𝑖𝑢 -𝐶 ̂𝑖𝑢 ] • [𝐶 𝑗𝑢 -𝐶 ̂𝑗𝑢 ] 𝑛 𝑢=1
Eq. (38) with Ciu the experimental concentration and 𝐶 ̂𝑖𝑢 the estimated value for response i and event u; 𝐶 𝑗𝑢 the experimental concentration and 𝐶 ̂𝑗𝑢 the estimated value for response j and event u.

The interval estimates for each estimated parameter are calculated from the final quadratic expansion of the OF. The uncertainty of the estimated parameters was evaluated by the marginal highest posterior density (HPD), which was calculated by GREGPLUS package.

The temperature dependence of rate constants was expressed by a modified Arrhenius equation. Indeed, there is a strong correlation between the pre-exponential factor and activation energy. This strong correlation can increase the confidence intervals, and thus the uncertainty. To avoid such effect, Buzzi-Ferraris recommended to linearize the Arrhenius equation, [START_REF] Buzzi-Ferraris | Planning of Experiments and Kinetic Analysis[END_REF] hence the following expression was used to express rate constants,

𝑘(𝑇) = 𝑒𝑥𝑝 [𝑙𝑛 (𝑘(𝑇 𝑟𝑒𝑓 )) + 𝐸 𝑎 𝑅•𝑇 𝑟𝑒𝑓 • (1 - 𝑇 𝑟𝑒𝑓 𝑇 )]
Eq. (39) where 𝑇 𝑟𝑒𝑓 is the reference temperature, which was the median of the different reaction temperatures.

Preliminary modeling shows that Runs 19 and 20 could not be modeled correctly. Run 20 was performed without catalyst, and 19 was performed with low catalyst loading. The reaction mechanism in the presence of a low amount of Brønsted acid sites is different (Figures 5), thus they were discarded in the modeling stage.

Our modeling strategy included some kinetic experiments with only esterification reactions (Table 1) in the same operating conditions as the butanolysis ones. This inclusion was due to the fact that LA and FA concentrations were low, and thus the estimation of their rates constants just from BL and BF could lead to some estimation bias.

Table 2 summarizes the sum of squared residuals (SSR) for each compound and model, the number of estimated parameters for each model, the Akaike Information Criterion for each model, [START_REF] Capecci | Bayesian Statistics to Elucidate the Kinetics of γ-Valerolactone from n-Butyl Levulinate Hydrogenation over Ru/C[END_REF][START_REF] Mcdonald | Kinetic Model Discrimination of Penicillin G Acylase Thermal Deactivation by Non-Isothermal Continuous Activity Assay[END_REF] , and coefficients of determination for each compound.

It is complicated to find the best model based on SSR or R 2 values (Table 2). Therefore, the criterion AIC was also used.

𝐴𝐼𝐶 = 𝑁𝑢𝑚𝑏𝑒𝑟 𝑜𝑓 𝑖𝑛𝑑𝑒𝑝𝑒𝑛𝑑𝑎𝑛𝑡 𝑒𝑣𝑒𝑛𝑡𝑠 • 𝑙𝑛 ( 𝑆𝑆𝑅 𝑁𝑢𝑚𝑏𝑒𝑟 𝑜𝑓 𝑖𝑛𝑑𝑒𝑝𝑒𝑛𝑑𝑎𝑛𝑡 𝑒𝑣𝑒𝑛𝑡𝑠 ) + 2 • 𝑁𝑢𝑚𝑏𝑒𝑟 𝑜𝑓 𝑒𝑠𝑡𝑖𝑚𝑎𝑡𝑒𝑑 𝑝𝑎𝑟𝑎𝑚𝑒𝑡𝑒𝑟𝑠 Eq.( 40
)
This criterion takes into account the number of estimated parameters. Lower is the AIC value, and higher is the model quality. Table 2, one can notice that Model 3 is more reliable than the others.

The values of the estimated constants as well as the credible intervals are given in Table S2.1

and Table S2.2 (Supporting information S2). In Table S2.3, the estimated activation energies obtained from this work and found in literature [START_REF] Konaka | A Study of Reaction Rates. Esterification of Formic Acid with Ethanol[END_REF][START_REF] Bart | Kinetics of Esterification of Succinic Anhydride with Methanol by Homogeneous Catalysis[END_REF][START_REF] Kokare | Response Surface Optimization, Kinetic Study and Process Design of n-Butyl Levulinate Synthesis[END_REF][START_REF] Van Putten | a Versatile Platform Chemical Made from Renewable Resources[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] Sacia | Biomass Conversion to Diesel via the Etherification of Furanyl Alcohols Catalyzed by Amberlyst-15[END_REF] are displayed. It is not possible to compare the activation energies, because solvent/reactant and proposed models by other authors (reaction rate expression) were not the same.

For Models 1 and 2, one can consider that the credible intervals are low, meaning that the values were well identified. The credible intervals for the estimated kinetic constants for the hydrolysis of 5-HMF (Step 1) and esterification of FA was slightly wider, which is due to the low concentration of LA and FA. For Models 3 and 4, the credible intervals are wider for the transformation of Int1 to LA and FA and to BL and BF. These wider credible intervals are linked to the fact that the concentration of Int1 was not tracked.

Parity plots for each observable and each model are displayed in Supporting information (S3).

One can notice that parity is good for HMF, LA, BF and BL concentrations because the coefficient of determination is higher than 0.95. The parity between experimental and estimated BMF concentrations is slightly lower, i.e., R2 is ca. 0.9. The mechanism could explain this lower value from HMF to BMF is more complex.

Residual plots for each compound and model are displayed in Supporting information (S4).

One can notice that residual values versus experimental and estimated concentrations are randomly distributed, meaning that the defined rate equations were well identified.

Correlation matrix for the estimated parameters were included in Supporting information (S5).

The estimated parameters were not correlated to each other because all binary correlation coefficients are lower than 0.95. 70

Cross-validation

The validation stage is important to finalize the assessment of the different models. A crossvalidation approach was used, and more specifically the K-fold method. [START_REF] Slotboom | Critical Assessment of Steady-State Kinetic Models for the Synthesis of Methanol over an Industrial Cu/ZnO/Al2O3 Catalyst[END_REF] The 18 experiments were randomly divided into 6 folds (Table 3) by paying attention that each fold contains at least one esterification experiments but not 3. The training (a.k.a regression) were made on 5 folds and testing (a.k.a validation) on the remaining fold as shown by Table 4. Eq. (41) As the 𝐶𝑉 (𝐾) value is low, then the model is predictable. Table 5 shows that 𝐶𝑉 (𝐾) numbers are lower for Models 1 and 3.

Conclusions

The butyl levulinate (BL) production from butanolysis of 5-HMF, without the addition of water, over Amberlite IR-120 was investigated in batch conditions and 20 bars of nitrogen. Butanolysis of 5-HMF is the central reaction system for the production of BL from glucose or fructose butanolysis. GVL was used as a co-solvent. It was found that experiments were repeatable and that catalyst did not undergo any significant deactivation. The analytical investigation showed that 5-(butoxymethyl)furfural (BMF), levulinic acid (LA), butyl formiate (BF) and butyl levulinate can be tracked during an experiment. We have also noticed a side-reaction of humins production from mass balance analysis. Preliminary observations showed that the increase of reaction temperature and catalyst loading increase the kinetics of all reaction steps. From our data, esterification reactions were found to be irreversible. The Etherification reaction was assumed to be irreversible.

The second part of this manuscript was to develop and assess different kinetic models to explain the kinetics of this system. The assessment was done via the determination of AIC criteria for each model, and by the use of cross-validation. To ease the parameter estimation stage, independent LA and FA esterification experiments were performed. Operating conditions such as temperature, catalyst and 5-HMF loading were widely varied to develop reliable kinetic model. Basically, two types of models were tested: one with 5 reaction steps and the other one with 8 reaction steps. We have also found that the production of humins was first-order. Based on the evaluation of the statistical output, AIC criterion and the cross-validation stage, it was found that Model 3, i.e., with 8 reaction steps and with a first-order reaction for humins production was the most probable model. Model 3 can simulate the alcoholysis of 5-HMF for experiments with a 5-HMF loadings from 9.7 to 43.7 g/L, a dried Amberlite IR-120 loading from 29.1 to 93.8 g/L and a reaction temperature from 80 to 115°C.

To improve the kinetic model of this reaction system, more investigations are needed to understand the mechanim of humins production and determine which intermediates are formed for the production of levulinates. Using quantum mechanics and advanced analytical techniques could unravel the reaction mechanism for this system and establish accurate material balances.

Notation

Figure 1 . 53 2

 153 Figure1. Reaction steps for the alcoholysis of fructose.[START_REF] Ramírez | Assessment of Ion Exchange Resins as Catalysts for the Direct Transformation of Fructose into Butyl Levulinate[END_REF]

Figures S1. 2

 2 Figures S1.2 display the kinetic profiles for the different species during the alcoholysis of 5-HMF. One can observe that Experiment 6 and 6-repeated are similar from a kinetic standpoint.

  both experiments, maybe because of the production of BL from LA esterification. Catalyst deactivation can be considered negligible.

Figures 2 .

 2 Figures 2. Effect of catalyst deactivation on the alcoholysis of HMF-5.

Figures 3 .

 3 Figures 3. Effect of temperature on the kinetics of HMF-5 alcoholysis.

  [𝐵𝐹][𝐻𝑀𝐹] 0 reached a plateau and not ratio[𝐵𝐿] [𝐻𝑀𝐹] 0. In Run 17, the fact that the ratio[𝐵𝐹] [𝐻𝑀𝐹] 0reached a plateau and the increase of the ratio[𝐿𝐴][𝐻𝑀𝐹] 0 could be linked. If the concentration of LA increases, thus the FA concentration also increases, but the volatility of the latter is higher than LA. Thus, there was less FA than LA in the liquid phase limiting the esterification of FA by butanol to produce BF.For this series of experiments, the mass balance on the initial concentration of HMF shows a side-reaction of HMF degradation into humins.

Figures 4 .

 4 Figures 4. Effect of HMF loading on the kinetics.

5Figures 5 .

 5 Figures 5. Effect of catalyst loading on the kinetics of HMF butanolysis.

Figure 8 .

 8 Figure 8. Reaction mechanism considering an intermediate Int1.

In the case of 5 -

 5 HMF-free experiments (Experiments 14 and 16; Figures 11 and 15), Model 3 slightly underestimates experimental concentrations of BL and BF for Experiment 14, but overestimates the BL experimental concentrations for Experiment 16.

Figure 9 .

 9 Figure 9. Fit of Model 3 to experimental concentrations (Experiment 8).

Figure 10 .

 10 Figure 10. Fit of Model 3 to experimental concentrations (Experiment 9).

Figure 11 .Figure 12 .Figure 13 .Figure 14 .

 11121314 Figure 11. Fit of Model 3 to experimental concentrations (Experiment 14)

Figure 15 .

 15 Figure 15. Fit of Model 3 to experimental concentrations (Experiment 16).

Table 1 .

 1 Experimental matrix

	Exp.	Temp. °C	mHMF0 g	mBuOH0 g	mGVL0 g	mdried cat. g	mLA0 g	mFA0 g	[GVL]0 mol.L -1	[BuOH]0 mol.L -1	[HMF]0 mol.L -1	[LA]0 mol.L -1	[FA]0 mol.L -1
	1		1.6	85.2	36.5	4.9	0.00	0.0	2.40	7.57	0.08	0.00	0.00
	2		1.6	85.2	36.5	4.9	0.00	0.0	2.40	7.57	0.08	0.00	0.0
	3	90	1.6	85.2	36.5	4.9	0.00	0.0	2.40	7.57	0.08	0.00	0.00
	4		7.0	85.2	36.5	8.0	0.00	0.0	2.31	7.27	0.31	0.00	0.00
	5		3.0	85.2	36.5	9.2	0.00	0.0	2.38	7.49	0.16	0.00	0.00
	6		3.0	85.2	36.5	4.9	0.00	0.0	2.38	7.49	0.15	0.00	0.00
	7	80	3.0	85.2	36.5	4.9	0.00	0.0	2.38	7.49	0.15	0.00	0.00
	8	85	7.0	85.2	36.5	8.0	0.00	0.0	2.31	7.27	0.32	0.00	0.00
	9		3.0	85.2	36.5	10.0	0.00	0.0	2.38	7.49	0.14	0.0	0.00
	10	80	0.0	80.9	34.7	4.0	8.13	0.0	2.28	7.17	0.00	0.51	0.00
	11	95	0.0	80.9	34.7	4.0	8.13	0.0	2.28	7.17	0.00	0.51	0.00
	12		0.0	80.9	34.7	3.0	8.13	0.0	2.28	7.17	0.00	0.51	0.00
	13	80	0.0	84.0	36.0	2.0	1.64	0.63	2.39	7.53	0.00	0.10	0.09
	14		0.0	84.0	36.0	2.0	1.63	0.64	2.39	7.53	0.00	0.10	0.09
	15	95	0.0	84.5	36.3	6.0	1.00	0.40	2.40	7.57	0.00	0.06	0.06
	16		0.0	84.0	36.0	8.0	1.63	0.64	2.39	7.53	0.00	0.12	0.09
	17		5.0	85.2	36.5	15.3	0.00	0.0	2.34	7.37	0.28	0.00	0.00
	18		3.0	85.2	36.5	10.0	0.00	0.0	2.37	7.49	0.17	0.00	0.00
	19		3.0	85.2	36.5	1.0	0.00	0.0	2.37	7.49	0.17	0.00	0.00
	20		3.0	85.2	36.5	0.0	0.00	0.0	2.37	7.49	0.17	0.00	0.00

Table 2 .

 2 Summary of modeling results.

		Model 1	Model 2	Model 3	Model 4
	Number of Parameters	12	12	16	16
	SSR_HMF	0.11	0.12	0.11	0.12
	SSR_BMF	0.08	0.09	0.08	0.09
	SSR_LA	0.09	0.10	0.08	0.08
	SSR_BL	0.21	0.20	0.21	0.20
	SSR_BF	0.03	0.03	0.03	0.03
	SSR_all	0.52	0.54	0.51	0.53
	AIC	-23435	-23312	-23495	-23378
	R 2 _HMF	0.97	0.96	0.97	0.96
	R 2 _BMF	0.89	0.89	0.90	0.89
	R 2 _LA	0.99	0.99	0.99	0.99
	R 2 _BL	0.97	0.97	0.97	0.97
	R 2 _BF	0.98	0.98	0.98	0.98

Table 3 .

 3 Distribution of the 18 experiments in the 6 folds.

	Fold	Experiment
		12
	Fold 1	11
		5
		6
	Fold 2	3
		10
		7
	Fold 3	14
		8
		16
	Fold 4	4
		17
		1
	Fold 5	9
		13
		18
	Fold 6	15
		2

Table 4 .

 4 Different sets for training and testing. The kinetic constants 𝑙𝑛 (𝑘(𝑇 𝑟𝑒𝑓 )) and 𝐸 𝑎 𝑅•𝑇 𝑟𝑒𝑓 for each reaction are estimated from each training stage, and these estimated constants are used for the testing stage. The prediction capacity of a model is evaluated by the 𝐶𝑉 (𝐾) number.

	Set	Train	Test
	Set 1	Folds 1-2-3-4-5	Fold 6
	Set 2	Folds 6-1-2-3-4	Fold 5
	Set 3	Folds 5-6-1-2-3	Fold 4
	Set 4	Folds 4-5-6-1-2	Fold 3
	Set 5	Folds 3-4-5-6-1	Fold 2
	Set 6	Folds 2-3-4-5-6	Fold 1

Table 5 .

 5 CV(K) and standard deviation for each model

		Model 1	Model 2	Model 3	Model 4
	𝐶𝑉 (𝐾)	0.130	0.144	0.126	0.248

Acknowledgments

This investigation was carried out in the framework of the ARBRE project. ARBRE is cofounded by European Union through the European Regional Development Fund (ERDF) and by Normandy Region, via the support of "pôle CTM (Continuum Terre-Mer) and EP2M (Énergies, Propulsion, Matière, Matériaux) de Normandie université". The Daniele Di Menno Di Bucchianico's doctoral thesis is done in the framework of "ARBRE RIN 100%" funded by Region Normandie, and this thesis is a part of the research activities of "pôle CTM (Continuum Terre-Mer) de Normandie Université". For the analytical part, this work has been partially supported by University of Rouen Normandy, INSA 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 by Region Normandie. GC was financed by FEDER RIN Green Chem 2019NU01FOBC08 Number: 17P04374.

Supporting information

This information is available free of charge via the Internet at http://pubs.acs.org/. Repeatability (S1); parameters estimated and from the literature (S2); parity plots (S3); residuals (S4); and correlation matrices (S5). His group aims to understand the reactivity of complex chemical systems: develop robust and reliable kinetic models in different thermal modes. They use different calorimeters (RC1, C80, ARSST, etc), process intensification systems (micro-reactor, microwave irradiation) and statistical methods (Bayesian inferences, cross-validation, etc).

Abbreviations


domain: sdv chim

text:
Jean-Marc Moulis
email: jean-marc.moulis@cea.fr

Zorica Bulat
email: zorica.bulat@pharmacy.bg.ac.rs

Aleksandra Buhadjordjevic
email: aleksandra.buha@pharmacy.bg.ac.rs

Aleksandra Buha Djordjevic

Threshold in the toxicology of metals: Challenges and pitfalls of the concept

Keywords:

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

General trends, definitions, and scope

Even though the chemistry of life is generally associated with the diversified carbon, i.e., organic, (bio)chemistry, with the involvement of a handful of mainly metalloid elements (C, O, S, H, P, N, etc.), the terrestrial environment is qualitatively dominated by metals, some being mandatory for life, while others are not [START_REF] Maret | The metals in the biological periodic system of the elements: Concepts and conjectures[END_REF]. The essentiality of metals in Biology is not directly related to their natural abundance.

Hence every life form has to cope with the exposure to metals, and the question addressed here is whether animals, particularly humans, can bear some exposure to metals without harmful effects. If such a safe level exists it defines a threshold for each element/metal. The biological need for metals of the first two columns of the periodic table, such as sodium, magnesium or calcium, is high. Similarly, transition metals, such as iron, copper, or zinc, are mandatory; they must be sufficiently available to avoid deficiency which implies that the onset of toxicity occurs at high concentrations relative to the required availability range. The definition of a threshold in such cases shifts to that of a toxicological level, and it does not need further elaboration here. Instead the focus will be on other, non-essential metals often with exchangeable electrons that usually occur as cations in biological systems.

It should be stated that non-essential is a provisional notion defined by default as already discussed [START_REF] Maret | The metals in the biological periodic system of the elements: Concepts and conjectures[END_REF], since the demonstration of essentiality in a given species sometimes only awaits the application of proper, often new, methods under the suitable experimental conditions. Also, a metal which is mandatory for some species may be non-essential for others; nickel is required for many microorganisms (bacteria and archaea), but it seems merely non-essential (toxic) for animals.

Metal ions are not metabolized since they are not transformed per se; they only exchange ligands and sometimes electrons in low energy reactions. In contrast to organic pollutants that experience modifications on or of their carbon skeleton, the toxicology of metals instead relies on such phenomena as speciation, traffic, or storage, and their impact on the normal functions of the exposed organism.

The issue of the use and existence of a toxicological threshold, in general, has already been deeply considered in the literature [START_REF] Crump | Use of threshold and mode of action in risk assessment[END_REF], and the present Opinion addresses the specific situation of metal ions, with a focus on mechanistic considerations that may prove crucially informative in the Toxicology of metals.

Peculiar features of chronic exposure: detoxification, adaptation, and exposure threshold for metals

Although acute exposure to toxic metals can occur with relatively short-term adverse health effects, most frequent and critical situations involve long-term, low concentrations of metal exposure with equivocal health effects. In the latter case, dose needs to be precisely described as two independent components define it, the (metal) concentration which may vary with time, and the length of exposure, up to life-long, spanning all developmental stages from gestation to aging. Although exposure is most often considered in these circumstances, a better parameter should COTox-S-19-00011 Revised p. 4 be the burden of the toxic compound [START_REF] Jacquet | Impact of maternal low-level cadmium exposure on glucose and lipid metabolism of the litter at different ages after weaning[END_REF] which is the functional parameter that implicitly take sources of variability such as the mode of administration/exposure, single nucleotide polymorphisms, or imprinting effects into account. Indeed, beyond the individual sensitivity to the metal toxicant, the rate of uptake and processing varies with external concentration and considered organs as demonstrated for Cr(VI) for instance [START_REF] Kirman | Physiologically based pharmacokinetic model for rats and mice orally exposed to chromium[END_REF].

For chronic exposure, the longer the length of exposure and the lower the concentration to minimize immediate damage, the more likely the ability of the exposed organism to develop some adaptive response to the insult. Indeed, organisms are endowed with specific or serendipitous detoxification systems, be-it export and homeostasis in general, or repair as for genotoxic compounds. These mechanisms counterbalance the potential direct toxic effect of a metal ion, and they contribute to set up an actual high limit of the dose -a genuine threshold-below which no harm occurs, as long as they can cope with the insult. In such circumstances the question of the existence of a threshold is whether non-essential metals trigger regulatory mechanisms and export systems.

Alternatively, a genuine threshold may occur if the metal needs to reach a high enough concentration to interfere with cellular or organismal function, as chronic gradual uptake of metal pollutants is unlikely to instantly saturate biochemical systems. If there is no cellular binding and uptake of the metal cation (unlikely), if the cells efficiently expel the metal, or if they can safely store it, the biochemical pathways, except the detoxifying ones and those that support them (e.g. energy provision), are not affected. About metal storage, are animal cells endowed with systems that safely and lastingly trap non-essential metal cations as plant or yeast do in vacuoles [START_REF] Nguyen | Isolated saccharomyces cerevisiae vacuoles contain low-molecular-mass transition-metal polyphosphate complexes[END_REF]? Metallothioneins are supposedly metal detoxification molecules, and ferritin and albumin are serendipitous ones. However, these and other cation-binding proteins experience dynamic turnovers in mammals, which may shift the toxicological problem to some other compartment or organ, but without solving it for the whole organism if these metal complexes are not safely excreted.

At the cellular level, the macroscopic behavior of the studied system results from discrete molecular reactions in which a given non-essential metal cation is expected to have multiple targets of similar affinities. This multiplicity of similarly sensitive targets comes from the limited set of binding anion groups in the multiple biological metal ligands. The particular property of metal ions interacting with several molecular targets with similar affinities is likely to promote the occurrence of multiple modes of action at low concentrations leading to overlapping discrete molecular responses and affected functions [START_REF] Höfer | Dose-and routedependent hormonal activity of the metalloestrogen cadmium in the rat uterus[END_REF][START_REF] Renieri | Nonlinear responses to waterborne cadmium exposure in zebrafish. An in vivo study[END_REF][START_REF] Zheng | In utero exposure to hexavalent chromium disrupts rat fetal testis development[END_REF]. Such combination may sum up to an apparent linear response at low doses, i.e., no threshold, and the existence of individual thresholds at a given level, e.g., the cell, may result in a quasi-linear, non-threshold, response at higher ones, e.g., organ or organism [START_REF] Bogen | Linear-no-threshold default assumptions for noncancer and nongenotoxic cancer risks: A mathematical and biological critique[END_REF][START_REF] Crump | Bogen's critique of linear-no-threshold default assumptions[END_REF]. But in any case managing interactions with non-essential metals mobilizes resources, which is detrimental in molecular and energetic terms.

COTox-S-19-00011 Revised p. 5

Furthermore, the sensitivity of the exposed organism depends on various characteristics such as the development stage and the mode of exposure [START_REF] Jacquet | Impact of maternal low-level cadmium exposure on glucose and lipid metabolism of the litter at different ages after weaning[END_REF]. A special case, not treated here, is associated with exposure to non-essential metals released by medical devices [START_REF] Bijukumar | Systemic and local toxicity of metal debris released from hip prostheses: A review of experimental approaches[END_REF][START_REF] Granchi | Biological effects of metal degradation in hip arthroplasties[END_REF][START_REF] Marchioni | Silver nanoparticle fate in mammals: Bridging in vitro and in vivo studies[END_REF]. Non-essential metal ions (e.g., Cd 2+ ) may avidly bind to receptors (e.g., estrogen receptor) of naïve cells [START_REF] Stoica | Activation of estrogen receptor-alpha by the heavy metal cadmium[END_REF] and act as endocrine disruptors [START_REF] Buha | Overview of cadmium thyroid disrupting effects and mechanisms[END_REF][START_REF] Iavicoli | The effects of metals as endocrine disruptors[END_REF]. An adaptive response may develop up to the point where disruption of endocrine functions strongly perturbs the organism [START_REF] Rattan | Exposure to endocrine disruptors during adulthood: Consequences for female fertility[END_REF]. The consequences vary between induction of (cell) death, sham stimulus for growth or differentiation or loss of functionality, and they change over time [START_REF] Jacquet | Impact of maternal low-level cadmium exposure on glucose and lipid metabolism of the litter at different ages after weaning[END_REF]. Mechanisms determine the chronology of events with a sharp difference between the fast receptor-mediated and hormone-like actions [START_REF] Zheng | In utero exposure to hexavalent chromium disrupts rat fetal testis development[END_REF] or epigenetics [START_REF] Bommarito | Effects of prenatal exposure to endocrine disruptors and toxic metals on the fetal epigenome[END_REF], and delayed oxidative stress [START_REF] Jacquet | Chronic exposure to low-level cadmium in diabetes: Role of oxidative stress and comparison with polychlorinated biphenyls[END_REF] that is initially controlled by the antioxidant response.

In summary, if a metal does not have any known function, a safe assumption is that it will be toxic at some point if it interacts with cells. But this relies on relevant data 1) to accurately identify the source of input and the metal traffic and storage, including full reversibility in case exposure ceases, and 2) to estimate the damage or perturbation the interacting metal can do. This scientific area is vast and leaves plenty of room to future biochemical toxicology.

Experimental dose-response curves for metals

The threshold concept implies that low doses below a defined value are harmless.

However, toxicological effects, particularly at low doses, are plagued with experimental imprecisions and various sources of uncertainty as recently listed in the case of cadmium-induced thyroid dysfunction [START_REF] Buha | Overview of cadmium thyroid disrupting effects and mechanisms[END_REF]. For example, any parameter or effect must be measured over a baseline which corresponds to the background exposure and is rarely null with environmentally relevant metals, even in laboratory settings [START_REF] Jacquet | Impact of chronic and low cadmium exposure of rats: Sex specific disruption of glucose metabolism[END_REF]. In addition, taking for reference functional concentrations of physiological metals such as zinc [START_REF] Wilson | Picomolar concentrations of free zinc(II) ions regulate receptor protein-tyrosine phosphatase beta activity[END_REF] or iron [START_REF] Pourcelot | Iron for proliferation of cell lines and hematopoietic progenitors: Nailing down the intracellular functional iron concentration[END_REF], monitoring signals at very low doses, i.e., at or below part per billion levels, may mean not more than a few metal cations per cell or per sub-cellular compartment at which level quantitative speciation is a huge technological challenge. The issue of probing departure points in doseresponse curves is sketched in Figure 1.

As a further concern, the most often assessed functional effects of non-essential metals may not be detected before the onset of structural ones as observed with cadmium in various rodent organs [START_REF] Brzóska | Changes in the structure and function of the kidney of rats chronically exposed to cadmium. I. Biochemical and histopathological studies[END_REF][START_REF] Jacquet | Mitochondrial morphology and function of the pancreatic beta-cells ins-1 model upon chronic exposure to sub-lethal cadmium doses[END_REF]. Additionally, consideration of negative results should not be neglected in estimating a threshold as long as the level of sensitivity of the associated methods is adequately assessed.

Due to the sources of uncertainty, non-essential metals have shown dose-response curves that depart from linearity [START_REF] Renieri | Nonlinear responses to waterborne cadmium exposure in zebrafish. An in vivo study[END_REF]. Comparison of markers of exposure to cadmium and effect in the general Japanese population revealed the presence of multi-linearity [START_REF] Ikeda | Bi-linear dose--response relationship in general populations with low-level cadmium exposures in nonpolluted areas in japan[END_REF] with significant geographical variations [START_REF] Sakuragi | Variation in benchmark dose (bmd) and the 95% lower confidence limit of benchmark dose (bmdl) among general japanese populations with no anthropogenic exposure to cadmium[END_REF]. As an example of a possible threshold, decreased glomerular filtration rate correlated with increased β2-COTox-S-19-00011 Revised p. 6 microglobulin in human subjects displaying more than 0.5 µg of urinary cadmium/g creatinine [START_REF] Satarug | Urinary cadmium threshold to prevent kidney disease development[END_REF]. This cut-off value, which is significantly lower than toxicity values set by international regulatory organizations, signals functional damage to the kidney, but it does not inform on the intimate status of the glomerulus. Hence the operational threshold estimated in this study cannot be converted into evidence for no (molecular and possibly cellular) effects on the nephron. A new meta-analysis examining the association between cadmium burden and diabetes mellitus in human adults [START_REF] Guo | Evaluation of the association between urinary cadmium levels below threshold limits and the risk of diabetes mellitus: A dose-response meta-analysis[END_REF] implemented an arbitrary cubic spline curve to the selected data, claimed a 2.43 µg cadmium in urine/g creatinine threshold value for the onset of correlation with diabetes, but could not conclude on the shape of the dose-response curve at lower doses (≤2.4 µg urinary cadmium/g creatinine) within the 95% confidence intervals.

A recent study designed to describe the relationship between environmental lead (Pb) exposure and maternal stress among pregnant women suggested a nonlinear association between the investigated parameters [START_REF] Li | The non-linear association between low-level lead exposure and maternal stress among pregnant women[END_REF]. Similarly, the non-linear relationship between exposure and burden of oral hexavalent chromium [Cr(VI)] [START_REF] Kirman | Physiologically based pharmacokinetic model for rats and mice orally exposed to chromium[END_REF] was associated with rodent duodenal hyperplasia and long-term cancerous lesions (reviewed in [START_REF] Thompson | Assessment of the mode of action underlying development of rodent small intestinal tumors following oral exposure to hexavalent chromium and COTox-S-19-00011 Revised p. 14 relevance to humans[END_REF]). Among the mechanisms of action of chromium are hypermethylation and down-regulation of DNA repair genes [START_REF] Hu | Cr(vi)-induced methylation and down-regulation of DNA repair genes and its association with markers of genetic damage in workers and 16hbe cells[END_REF].

All the above and other examples with metals exemplify non-monotonous doseresponse (NMDR) curves. This type of curve suggests, but does not necessarily imply, that the threshold concept should be used in interpreting such data. It may also include effects such as hormesis [START_REF] Shi | Progress in the studies on hormesis of low-dose pollutants[END_REF] which, when really present, can be seen as the result of overlapping mechanisms depending on the toxicant dose and an early step of adaptation. In such cases extrapolation from effects observed at high doses (most commonly used in standard animal models) to environmental doses for humans is highly questionable. Indeed, studies of mutagenic, carcinogenic, developmental and reproductive toxicity of metals are dominated by relatively high concentrations, irrelevant delivery routes, and short times of exposure adjusted to the life of the model species (and to the patience and means of the experimenter) [START_REF] Buha | Cadmium exposure as a putative risk factor for the development of pancreatic cancer: Three different lines of evidence[END_REF].

Thus, high-dose testing that is not quantitatively relevant to human, or other species, is not justified, as has been translated into the Kinetically-derived Maximum Dose concept for dose selection strategy [START_REF] Bus | The dose makes the poison": Key implications for mode of action (mechanistic) research in a 21st century toxicology paradigm[END_REF]. The number of critical pathways in which non-essential metals can be involved calls for alertness in the studies of their toxicology.

Consequences for risk assessment of non-essential metal exposure

When it comes to threshold importance for the regulatory implications of doseresponse curves, the societal perception now shifts the relevant question from the statistical definition of an acceptable healthcare and economic risk for a given human population to one of potentially finding an exposure value below which no individual of any population will experience toxic effects. The data to address this shifting purpose are hardly available [START_REF] Beausoleil | Review of nonmonotonic dose-responses of substances for human risk assessment[END_REF].

COTox-S-19-00011 Revised p. 7

Most toxicological models defining benchmark doses (BMD), no-(adverse)-effects levels, and other alternative parameters implicitly include the threshold concept [START_REF] Crump | Use of threshold and mode of action in risk assessment[END_REF][START_REF] Thompson | Integration of mechanistic and pharmacokinetic information to derive oral reference dose and margin-ofexposure values for hexavalent chromium[END_REF]. They all define some threshold and with the use of safety and uncertainty factors can be transferred into doses that are unlikely to cause any effects in humans. The notion of threshold of toxicological concern refers to this kind of considerations, by targeting a dose value below which worry becomes negligible even in the face of cocktails of potential toxicants [START_REF] Hartung | Thresholds of toxicological concern -setting a threshold for testing below which there is little concern[END_REF]. But metals were excluded from the list of chemicals considered with this concept, and mechanistic information is only marginally included in such risk assessment procedure [START_REF] Hartung | Thresholds of toxicological concern -setting a threshold for testing below which there is little concern[END_REF].

To provide improved risk assessment of chronic low-dose non-essential metal exposure to populations, detailed and quantitative mechanistic data obtained via more sensitive and more comprehensive (-omics) analytical methods must be implemented. The development of suitable computational tools to integrate these data in the description of such complex systems must also follow a parallel progress.

Presently, specific issues still need further scrutiny [START_REF] Saracci | The hazards of hazard identification in environmental epidemiology[END_REF]: relevance of high dose effects for exposure at low doses, appropriateness of used methodology to determine threshold doses, and accuracies in establishing NMDR are crucial items.

Historically, the linear non-threshold form of the dose-response curve has been favored for genotoxic compounds, including metals, implying that a single metal molecule may generate DNA lesions. Indeed, redox-active non-essential (for humans) transition metals such as nickel or chromium may catalyze oxidative modifications of nucleotides, and others, such as cadmium, may indirectly lead to the same outcome by blurring the antioxidant response. But this neglects the existence of DNA repair mechanisms that avoid the persistence of modifications eventually leading to tumors, at least up to some density of damage. This type of resistance mechanism should support the threshold concept in this example of cancer as long as the inflection point in the dose-response curve can be precisely established in a range of exposure conditions, and the fate of the DNA damaging metal is monitored.

Both of these propositions usually remain elusive. At best, extensive mechanistic and pharmacokinetic information with benchmark dose modeling, as in the case of chromate exposure [START_REF] Haney | Consideration of non-linear, non-threshold and threshold approaches for assessing the carcinogenicity of oral exposure to hexavalent chromium[END_REF], enables derivation of a Reference Dose, i.e., the acceptable daily lifetime intake for human subjects that is supposedly safe. But the difficulty in precisely measuring toxicant exposure and burden for humans, the uncertainties plaguing dose-response curves (Figure 1), and the degree of subjectivity associated with what is acceptable open such parameters as Reference Dose to controversy.

Observations adding to the complexity of the problem include the interaction of the molecular networks handling essential and non-essential metals. For instance, low daily doses of Cd given to rats for 4 weeks (BMDL5= 10 µg Cd/kg b.w.) induced a 5% decrease in bone Zn levels that mirrored the inverse association between bone Cd content and that of Zn or bone mineral density of human samples [START_REF] Buha | Bone mineral health is sensitively related to environmental cadmium exposure-experimental and human data[END_REF].

COTox-S-19-00011 Revised p. 8

For some metals such as lead, a consensus emerges for zero tolerable blood levels based on physiological and epidemiological arguments rejecting the occurrence of a threshold in the dose-response relative to the nervous system and the kidneys [START_REF] Shefa | Both physiology and epidemiology support zero tolerable blood lead levels[END_REF].

Building up the successive quantitative and temporal events occurring upon very low, i.e chronic and environmental, metal exposure, is slowly but regularly updated.

Unambiguously ruling on the existence or absence of thresholds for non-essential metal ions would shift risk assessment from a computationally derived level to a knowledge-based value with a thorough insight in case of changing conditions, such as upon co-exposure to other potential toxicants or in the study of different populations.

Conclusions

The demonstration of a safe threshold for any given metal ion requires very detailed experimental studies and attention to molecular mechanisms relevant to environmental conditions. Modeling of dose-response curves is still helpful for practical purposes, but it does not help when looking for intervention for corrective actions. For prevention let's state sublata causa, tollitur effectus (get rid of the cause and the effect will disappear): easier said than done, especially in present day technological societies that have recently found fancy applications of previously unused and unmined metals such as lithium or lanthanides.

What is needed is knowing whether the non-essential metal ions can interact with the organism/tissue/cell under study (generally yes) and to define their targets at infinite low concentrations. As such, safe exposure due to the existence of thresholds is likely to remain an illusion even though the human body contains seemingly nonessential metals (e.g., strontium, rubidium, titanium, barium and others [START_REF] Maret | The metals in the biological periodic system of the elements: Concepts and conjectures[END_REF]) at levels similar to some essential metals. The source, pathway(s), disposal, and biological impact of such non-essential metals are completely unknown.

Uncertainties in the mechanisms of toxicity relevant to low-level human exposures remains a challenging current issue, but it seems unlikely that any informed progress on the matter, beyond current and somewhat arbitrary risk assessment, can be reached without in-depth mechanistic knowledge of the biological action of metals.

What happens around and below current benchmark doses is terra incognita and exploring this unknown is the challenge ahead, particularly considering the increased societal scrutiny of concepts such as "acceptable risk".

COTox-S-19-00011 Revised p. 9

Acknowledgements

The work on low-level toxicity of metals in JMM's group was supported by a grant from Agence Nationale de la Recherche ANR-13-CESA-008-Cadmidia. The work on low-level toxicity of metals in AB's group was partly supported by the Ministry of Education, Science and Technological Development of Serbia (Project III 46009).

The quotations in this Opinion are voluntarily biased toward recent publications: it does not imply that previous data addressing the topic should be neglected.

COTox-S-19-00011 Revised p. 11

Conflict of Interest

Aleksandra Buha Djordjevic is one of the Editors of this special issue, but had no personal involvement in the reviewing process, or any influence in terms of adjudicating on the final decision, for this article. The other authors declare that they have no competing interests regarding this work.

COTox-S-19-00011 Revised p. 10

Figure legend

Virtual dose-response curve based on experimental data. The error bars usually increase as doses are lowered which explains the uncertainty of the suitable curve to fit the data. The two extreme situations are represented by the (1) and (2) curves, with or without threshold. A decision tree is sketched to determine whether a threshold does exist and may support any safe exposure.


Yasser Diab
email: yasser.diab@gmail.com

Pascal Venet

Hamid Gualous

Gerard Rojat

Electrical, Frequency and Thermal Measurement and Modelling of Supercapacitor Performance

Keywords: Characterisation techniques, supercapacitor characteristics, coulumbic efficiency, self-discharge modelling, pseudo-inductance modelling

This paper presents an evaluation of commercial supercapacitors performance (ESR, C, self-discharge, P max , E max , coulumbic efficiency, etc), under different conditions. Characterization of supercapacitor performances as a function of current, voltage, frequency and thermal constraints, is carried out by constant current test, electrochemical impedance spectroscopy, and voltamperometry. A comparative study of measurement techniques is conducted in detail in this work. The supercapacitor self-discharge amplitude is determined by measuring the decline of open-circuit voltage. The mechanisms of the self-discharge, leakage current and diffusion of ions at the electrode-electrolyte interfaces are modelled by a new electrical equivalent circuit. In addition, we model the electrochemical impedance of self-discharge with diffusion-controlled in frequency domain. Supercapacitor is represented by two complementary electrical equivalent circuits: two branches and non-homogeny pores. The first represents the supercapacitor time behaviour. The second represents supercapacitor dynamic behaviour. Finally, the origin and a new model of pseudoinductance are described.

INTRODUCTION

The determination of supercapacitors parameters under conditions recreated in laboratory can predict their behaviours in industrial applications. Supercapacitors generate variable and superimposed physical phenomena. Supercapacitors characterization methodologies are not standardized and it is difficult to compare the results obtained with those of manufacturers. At present time no international standards are existed. The results are depending on different parameters as ambient temperature, current intensity, frequency and voltage. Therefore, reliable and reproductive experimental results are required; we propose in this paper a comparative study on measurement methods applied on MAXWELL commercialized supercapacitors (see tab. 1). In addition, we need a whole equivalent circuit model for representation supercapacitor phenomena in power electronics simulation. So, we use and compare two complementary models: two branches and non-homogeny pores, in frequency domain. We develop, in this paper, new equivalent circuits for pseudo-inductance and self-discharge phenomena.

Two branches model

A method of determining supercapacitor parameters is to perform constant-current charging or discharging. This method consists in representing of a supercapacitor by two time constants immediate and delayed branches. The supercapacitor equivalent circuit can be represented by fig. 3.

The immediate parameters can be identified carrying out a high charge/discharge current. . The differential capacitance C 1 (u) is presented by two capacitances: a constant capacitance C 0 and a linear voltagedependent capacitance k as given in equation 1 [START_REF] Zubieta | Characterization of doublelayer capacitors (dlcs) for power electronics applications[END_REF]:

u k C u C . ) ( 0 1 + = (1)
The resistance R 1 can be determined by the ratio of voltage drop at the beginning of charging to current. The constant C 0 in the equation 1 can be determined from the definition of differential capacity I/du(0)/dt determined near to initial instant of charging. The slope k can be identified by the evolution of voltage ∆U t and Q=I. ∆t t charge of supercapacitor stored during charging as follow [START_REF] Zubieta | Characterization of doublelayer capacitors (dlcs) for power electronics applications[END_REF]:

U ) U • C - t I. ( 2 = k 2 t t 0 t ∆ ∆ ∆ (2)
We founded experimentally that the R 1 is practically constant but, the capacitance C 1 (u) elements depend on the rate of charge/discharge current (e.g. the component BCAP010 see fig. 4). This difference between the charge and discharge capacitance is probably due to Faradic current (redox) [START_REF] Kurzweil | Electrochemical and Spectroscopic Studies on Rated Capacitance and Aging Mechanisms of Supercapacitors Based on Acetonitrile[END_REF] and to charges trapped in pores with low accessibility.

Coulumbic (discharge/charge) efficiency

In order to quantify the difference between the discharge and charge capacitance, we define a factor called coulumbic efficiency. The last represents the discharge capacitance relative to that of charge for the rated voltage. Its variation with charge/discharge current is shown in fig. 5 for the component BCAP010.

Voltamperometry measurement

This method is based on the current measurement resulting from a linear potential sweep dU/dt between two chosen limits U 1 and U 2 . It allows directly determining the capacitance [START_REF] Pell | Analysis of power limitations at porous supercapacitor electrodes under cyclic voltammetry modulation and dc charge[END_REF].

The variation of coulumbic efficiency with cycle number for Maxwell supercapacitors is shown in fig. 6. The capacitance determination for the studied supercapacitors was done with the same relationship capacitance-sweep. The increase of coulumbic efficiency during subsequent cycling is attributed to irreversible phenomena [START_REF] Tripathi | Electrochemical redox supercapacitors using PVdF-HFP based gel electrolytes and polypyrrole as conducting polymer electrode[END_REF].

Ragone plot

Ragone plots provide the available energy of a supercapacitor for constant power request [START_REF] Christen | Theory of Ragone plots[END_REF][START_REF] Christen | Optimizing energy storage devices using Ragone plots[END_REF]. Fig. 7 shows the Ragone plot for two Maxwell supercapacitors. The supercapacitors are charged as described above during 1h (see § 5.1) then they are discharged with constant power from U R to U R /2. The half-order model has been called to represent the non-homogeny pores at the electrodes phenomenon [START_REF] Riu | Half-order modelling of supercapacitors[END_REF] (see fig. 9). This model allows developing the coth term of eq. 4 (left side) to secondary order as follows:

ωτ ωτ τ ω j C j Z dl p . ) coth( . ) ( = dl p C j j Z . . 1 ) ( ω ω τ ω + ≅ (4)
In order to take into account non-homogeny pores phenomenon in impedance pore model, we can use the fitting.

The resulting model is called the non-homogeny pores model. This model represents a trade-off between precision and simplicity and can be used to model supercapacitors in power electronics simulation.

Fig. 9: Comparison between the different models in Nyquist diagram

It is essentially the self-inductance of the conductive loop formed by the terminals, connectors and electrodes of the cell inductive behaviour can also result from non-homogeneous current distribution [START_REF] Karden | A frequency-domain approach to dynamical modeling of electrochemical power sources[END_REF] and electrochemical reactions at supercapacitor [START_REF] Bisquert | Inductive behaviour by charge-transfer and relaxation in solid-state electrochemistry[END_REF][START_REF] Savova-Stoynov | Analysis of the inductance influence on the measured electrochemical impedance[END_REF]. In these cases, it represents an error in the EIS measurement; this is why it should be represented in the equivalent circuit (see fig. 9).

COMPARISON OF MEASUREMENT TECHNIQUES

From fig. 12, we constant that, on the one hand, the capacitance measured by EIS is generally lower than that obtained by the two transit techniques: charge/discharge with constant current and voltamperometry. This is due to the Faradique charges generated by redox which can take place during the charge. During the discharge (Faradique effect is practically non present), the results obtained by voltamperometry approach to those obtained by EIS. One the other hand, for tensions higher than 0.5 V, we find a good agreement between the capacitance values obtained by the two last techniques (see fig. 12).

Fig. 12: Capacitance voltage dependence by different techniques

In order to compare the different techniques of measurement, the nominal capacitance is determined by the mean of the capacitance between U R and U R /2 [START_REF] Harzfeld | Capacitance and Series Resistance determination in high power ultracapacitors[END_REF].

SELF-DISCHRAGE CHARACTERISATION AND MODELING 5.1. Self-discharge measurements

It has been found that a gradual decrease in voltage across the charged supercapacitor occurs in long term open-circuit. This is due to self-discharge. In this case, the self-discharge has two mechanisms: leakage current and diffusion of ions at the electrode-electrolyte interfaces [START_REF] Ricketts | Self-discharge of carbon-based supercapacitors with organic electrolytes[END_REF][START_REF] Conway | Diagnostic analyses for mechanisms of self-discharge of electrochemical capacitors and batteries[END_REF]. The second is called the self-discharge with diffusion-controlled. The selfdischarge can be characterized by measuring the decline of open-circuit voltage of the supercapacitor.

In order to limit the relaxing phenomena effect on the selfdischarge measurement, the supercapacitor is charged with a constant voltage supply for one hour. After one hour, all the components of the supercapacitor are approximately charged to the same voltage, assuming that the relaxing phenomena have time constants of less than one hour. Consequently, when the charging current is stopped, the decrease in voltage across the supercapacitor follows only the effect of self-discharge [START_REF] Ricketts | Self-discharge of carbon-based supercapacitors with organic electrolytes[END_REF]. Fig. 13 shows the self-discharge of the aforementioned supercapacitors at an ambient temperature of 25 °C and for an initial voltage of 2.5 V. It illustrates that the self-discharge rates of components are different. It is noticed that the components MC260 has the highest self-discharge rates, whereas the components BCAP010 has the smallest rate.

Voltampero

Self-discharge modelling

The self-discharge with diffusion-controlled is due to Faradic impurity reactions and it dominates the self-discharge during several hours (3h to 45h) [START_REF] Conway | Electrochemical Supercapacitors: Scientific Fundamentals and Technological Applications[END_REF]. From the diffusion equation of charges concentration at the interfaces carbon-electrolyte, we can deduce the supercapacitor current and voltage during the self-discharge [START_REF] Ricketts | Self-discharge of carbon-based supercapacitors with organic electrolytes[END_REF] (for BCAP010 is less than 10 h) (see eq. 5 and 6).

        - = t D h C m t i diff 2 0 exp . . ) ( ( 5)                 +                 - - - = t D h erfc D h t D h t m U t u 0 2 0 0 exp 1 ) ( (6)
Where, U 0 is the initial voltage across the supercapacitor. 2h 0 is the thickness of a region of electrolyte immediately adjacent to each electrolyte-carbon interface of the supercapacitor. m is the diffusion parameter which can be calculated by equation 4 [START_REF] Ricketts | Self-discharge of carbon-based supercapacitors with organic electrolytes[END_REF].

π 12 0 C D q C m R = (7)
where D is the diffusion coefficient of the ions in the electrolyte. C 12 is the series combination of the capacitances at the two interfaces of the capacitor, per unit area of one of the interfaces. C R0 is the initial ionic species concentration at the electrode-electrolyte interfaces. Also, q is the charge carried to the carbon surface by each ion [START_REF] Ricketts | Self-discharge of carbon-based supercapacitors with organic electrolytes[END_REF]. These parameters can identified by a physico-chemical analysis of supercapacitor electrolyte.

Fig. 14 shows the evolution of current and voltage during the self-discharge with diffusion-controlled. We obverse the voltage across the supercapacitor firstly decreases during several hours, and then it comes back very slowly to its initial value U 0 during several days. The increase in a supercapacitor voltage is called recovery voltage phenomenon. From the theory of self-discharge with diffusion-controlled and by using Laplace transform, we could deduce an electrochemical impedance model (see eq. 8).

C D j Dirac h D D j h BesselI jh D j D j h BesselI h Z ) ( . . . 2 , 0 . . . 2 , 1 . ) ( 2 2 ω ω ω ω ω ω +           -         - - = (8)
Nevertheless, this model requires some physical parameters which are difficult to determine and to be established in time domain. We propose modelling the self-discharge with diffusion-controlled and leakage current by the equivalent circuit depicted in fig. 15. The equivalent circuit elements R le , R r and C r can be determined by time behaviour of selfdischarge.

Experimental approach

The circuit shown in fig. 15 demonstrates that part of the charge stored immediately in capacitance C 1 will be redistributed into capacitance C 2 and C r . C 2 has a negligible effect. The phenomenon of the redistribution is relatively low because the capacitor is charged during an hour. When the supply current is stopped the charge Q 0 stored in C 1 will be diffused to C r as the C r has an almost zero voltage. At the end of the diffusion-controlled, both capacitances, C 1 and C r , will have the same voltage, U r , with total charge, Q r . The capacitance C r can be determined by using the charge balance:

0 Q Q r = (9)
where, ) ( .

) ) ( .( 1 1 0 0 0 r r r r U C U Q C U C U Q = + = (10) As a result, r r r r U U U C U U C C ). ( ). ( 1 1 0 0 - = (11)
The circuit schematised on the figure 15 can be reduced to that shown in the figure 16 during the self-discharge with diffusion-controlled From this figure, we can write the following equations:

dt t du C dt t du C t i t u R t i t u r C r r c r ) ( ) ( ) ( ) ( ). ( ) ( 1 - = = + = (12)
Thus, we can deduce the time evolution of supercapacitor voltage during the self-discharge with diffusion-controlled:

0 ) ( ) ( 1 2 2 1 = + + dt t du C C C dt t u d R C r r r ( 13
)
The solution of the above equation with presented initial conditions can be estimated by the equation

        - = r t U t u τ exp ) ( 0 ( 14
)
τ r is the constant time of self-discharge due to diffusioncontrolled, which can be calculated with the following equation :

r r r r C C R C C + = 1 1 τ (15)
Assuming that the capacitance C 1 >>C r the time constant can be reduced as follow:

r r r R C = τ (16)
The time constant of the exponential described in equation 16 was determined experimentally by fitting (with software of non-linear least squares fit) of the experimental result for the first hours. The resistance R r ascribed to the evolution of Faradic reactions at the electrode-electrolyte interfaces can be calculated by equation 17. Equivalent parallel resistance R le will impact long term storage performance because it is basically a leakage effect. The last can be negligible compared with the self-discharge with diffusion-controlled within the first hours There is a slight difference between the circuit simulation and the experimental results (1.2 % mean relative error) (see fig. 17). The experimental results demonstrate the self-discharge can be entirely represented by an equivalent circuit with low error. This circuit is more advantageous than the analytical model described in the previous section because it is easy to be established, and allow estimating the magnitude of the electrochemical storage of energy in electrical double layer capacitor. Tab. 4 shows the self-discharge parameters normalized to those given in the manufacturer's datasheet. We remark: the capacitance C r rises with an increase in capacitance of supercapacitor because the redox increases with an increase in electrode-electrolyte interfaces. Although this, the component BCAP010 has a relative capacitance, C r /C R , lower than that of BCAP013 because the concentration of impurities per unit area of electrode-electrolyte interface is lower in highest capacitances [START_REF] Tripathi | Electrochemical redox supercapacitors using PVdF-HFP based gel electrolytes and polypyrrole as conducting polymer electrode[END_REF]. These results correspond to those founded by voltamperometry (see fig. 6). The component MC2600 has the highest capacitance, C r , and the lowest resistances, R le , which signifies this component dissipates an important stored energy compared to that of others. Tab. 4: Self-discharge parameters of studied electrochemical capacitors

EFFECT OF VOLTAGE AND TEMPERATURE ON SUPERCAPACITOR PARAMETERS

The supercapacitors can operate in a temperature range as wide as -35 °C to 65 °C. In fig. 18, we show the self-discharge rate for a temperature of -25 °C, 0 °C, 25 °C, 45 °C and 65 °C and for an initial voltage of 1 V, 1.5 V, 2 V, and 2.5 V. Consequently, because self-discharge rates are dependent on the rates of ionic transport [START_REF] Ricketts | Self-discharge of carbon-based supercapacitors with organic electrolytes[END_REF], the self-discharge rates accelerate strongly with an increase in temperature or in initial voltage.

The impedance measurements are performed by EIS at different DC potentials (0 V, 0.5 V, 1V, 1.5 V, 2 V and 2. This result can be interpreted by impact of temperature on redox. In negative temperatures with fast charge/discharge the redox is deactivated. So, the coulumbic efficiency is almost 100%. In high temperatures with low charge/discharge (as by voltamperometry), it is very low (77%).

CONCLUSION

Supercapacitor parameters (R s , R el , C, coulumbic efficiency, self-discharge) are complex functions of voltage, time, and temperature. We studied the supercapacitor characteristics as a function of current, voltage, frequency and temperature by the help of different characterisations techniques and a comparative study is performed. Self-discharge with diffusioncontrolled is modelled in both the time and frequency domain as well as the pseudo-inductance. The different presented resultants were validated for all the studied supercapacitor.

Tab. 1 : 5 . 1 .

 151 Studied supercapacitors parameters as given in manufacturer's datasheet2. DC CHARACTERISATION AND MODELLING Classical equivalent RC circuitSupercapacitor performance can be approximated by the first order response of RC circuit. This model includes an equivalent series resistance (ESR) and an ideal capacitance (C). There are several measurement methods able to determine the both parameters (C and ESR). One of which is the standard. The voltage across the supercapacitor is maintained at its rated value U R during 1h then the supercapacitor is discharged with constant current I. The capacitance determination, C, was made during DC discharge between 0.8U R and 0.4U R according the linear relation C=I.∆t c /∆U c . The ESR can be determined at the start of discharge by ESR=∆U 0 /I. The voltage drop ∆U 0 is calculated as the intersection of extrapolated curve of experimental voltage and the straight of instantaneous voltage drop, occurred when the supply current is switched on (see fig.1)[START_REF] Kurzweil | Capacitance Characterization Methods and Ageing Behaviour of Supercapacitors[END_REF][START_REF] Iec | 62391-1, Fixed electric double layer capacitors for use in electronic equipment[END_REF][START_REF] Iec | 62391-2, Fixed electric double layer capacitors for use in electronic equipment[END_REF].

Fig. 1 :

 1 Fig. 1: Determination of supercapacitor parametersTherefore, this method results depend on discharge current as shown in fig.2:The capacitance, C, varies slightly but the ESR increases considerably with a reducing of current.

Fig. 2 :

 2 Fig. 2: Characterisation of supercapacitor parameters by standard technique

Fig. 3 :

 3 Fig. 3: Two branches model [4]

Fig. 4 :

 4 Fig. 4: Capacitance dependence on rate of charge/discharge current

Fig. 5 :

 5 Fig. 5: Coulumbic efficiency vs. charge/discharge current

Fig. 6 :

 6 Fig. 6: Coulumbic efficiency vs. cycling number

Fig. 7 :

 7 Fig. 7: Ragone plot of supercapacitors 3. FREQUENCY MEASUREMENT AND MODELLING 3.1. Frequency model Electrochemical impedance spectroscopy (EIS) is used for the electrical/ electrochemical characterization of supercapacitors.Its different parameters R s high frequency resistance, R el electrolyte resistance, inductance L s , low frequency capacitance C dl , etc, can be determined independently in different frequency ranges[START_REF] Buller | Modeling the dynamic behavior of supercapacitors using impedance spectroscopy[END_REF]. An electrical equivalent model can be represented as shown in fig.8with a finite number n of RC cells. The resistance R n and C dln are calculated by eq. 3.

Fig. 8 :Fig. 9

 89 Fig.8: Impedance pore model[START_REF] Buller | Modeling the dynamic behavior of supercapacitors using impedance spectroscopy[END_REF]

Fig. 10 :Fig. 11 :

 1011 Fig. 10: Supercapacitor complex pore model tacking into account pseudoinductance In order to validate experimentally the established model, we carried out a high frequency test by EIS from 60 Hz to 10 kHz (for component BCAP010). We compare in Nyquist diagram (see fig. 11) the experimental result with the two models. The first is schematized on fig. 8 with a simple inductance and the second is schematized on fig. 10. We observe the proposed model improve very much the supercapacitor frequency response.

Fig. 13 :

 13 Fig. 13: Self-discharge rate of studied supercapacitors

Fig. 14 :

 14 Fig. 14: Evolution of supercapacitor voltage and current during selfdischarge

Fig. 15 :

 15 Fig. 15: Proposed equivalent circuit of self-discharge

Fig. 16 :

 16 Fig. 16: Equivalent circuit of capacitor during self-discharge with diffusion-controlled

  Equivalent parallel resistance R le will impact long term storage performance because it is basically a leakage effect. In this case, the voltage across the capacitor decreases as an exponential function of time as follows: where τ le is the time constant of leakage current. The time constant is determined experimentally; the decline of opencircuit voltage is approximated by an exponential function. The Faradic leakage resistance R le is calculated as follows:

Fig. 17 :

 17 Fig. 17 : Comparison of the equivalent circuit of the self-discharge with experimental

Fig. 18 :

 18 Fig. 18: Self-discharge rate vs. time for different initial voltages and temperature (for component BCAP010 ) The dependence of the self-discharge parameters, C r and R le , on both the temperature and on the initial voltage across the supercapacitor can be shown in fig.19 and 20.

Fig. 19 :Fig. 20 :

 1920 Fig. 19: Variation of Cr v. temperature and initial voltage

  5 V), and different temperatures (-25 °C, 0 °C, 25 °C, 35 °C, 45 °C, 55 °C and 65 °C). Fig. 21 (a and b) shows the variation of supercapacitor parameters versus the temperature and voltage.

Fig. 21 :

 21 Fig. 21: Supercapacitor parameters vs. temps and voltage

Fig. 22 :

 22 Fig. 22: Coulombic efficiency vs. temperature for different current charge/discharge

  We present in the tab. 3 the different resistances of supercapacitor relative to that given by manufacturer. We can observe that the resistance values are dispersed; each method gives a value different from the others, but measurements by the methods standard and EIS (low frequency) give close values.

		Resistance (µΩ)	Rs (φ=0°)	EIS ESR 55mHz	400 A	Standard	50 A	Charge/discharg (400A)
		ESR	322 620 587 792 447
	BCAP010	ESRR (DC) ESR/ESRR (DC)% ESRR (1kHz)	700 700 700 700 700 46 89 84 113 64 300 300 300 300 300
		ESR/ESRR (1kHz)% 107 201 196 264 149
		ESR	272 439 422 701 360
	MC2600	ESRR (DC) ESR/ESRR (DC)% ESRR (100Hz)	400 400 400 400 400 69 110 106 175 90 280 280 280 280 280
		ESR/ESRR (100Hz)% 97	157 150 250 129

The results of the aforementioned techniques are compared in tab. 2. The capacitance values obtained both by impedance spectroscopy and by DC characteristics (discharge of voltamperometry and Norm) are in a good agreement, while the methods charge/discharge with constant current (C 1 ) and voltamperometry give higher values.

Tab. 2: Comparing of nominal capacitance values

We can distinguish in manufacturer's documentation two resistances: DC resistance. (determined by discharge with high constant current) and AC resistance (determined for a given frequency).

Tab. 3: Comparing of resistances values

- metry (10 mV/s) Standard Capacitance (F)

			Charge/discharge (400A)	EIS	Charge	Discharge	400 A	50 A
		CR	2600 2600 2600 2600 2600 2600
	BCAP010	C	2863 2614 2792 2596 2750 2623
		C/CR 110% 101% 107% 100% 106% 101%
	MC2600	C	3130 2842 2956 2869 2857 2861
		C/CR 120	109	114	110	110	110

Table 4

 4 illustrates that leakage resistance increases with a decrease of capacitance because membrane surface and quantity of impurities are lower.

	Component	BCAP010 BCAP013 BCAP350 MC2600
	Capacitance CR (kF)	2.6	0.45	0.35	2.6
	Rr (Ω)	58.1	291	384	45.2
	Cr (F)	201	47.4	38.8	367
	Cr/CR %	8	11	11	14
	Manufacturer
	leakage resistance	0.50	0.83	2.5	0.54
	(Rlem) (kΩ)
	Rle (kΩ)	1.34	4.99	8.0	0.73
	Rle/ Rlem %	269	599	321	135
	Rle/ CR (Ω/F)	0.52	11.1	23.0	0.28