# Reactive Distillation of Glycolic Acid Using Heterogeneous Catalysts: Experimental Studies and Process Simulation.

Carole Mutschler, Juliana Aparicio, Ilham Mokbel, Mickael Capron, Pascal Fongarland, Marcia-Carolina Araque Marin, Clémence Nikitine

<figure><img src="image_1.png" /></figure>

## To cite this version:

Carole Mutschler, Juliana Aparicio, Ilham Mokbel, Mickael Capron, Pascal Fongarland, et al. "Reactive Distillation of Glycolic Acid Using Heterogeneous Catalysts: Experimental Studies and Process Simulation." Frontiers in Chemistry, 2022, Frontiers in Chemistry, 10, pp.909380. 10.3389/fchem.2022.909380. [hal-04044528]

HAL Id: hal-04044528

https://hal.univ-lille.fr/hal-04044528

Submitted on 24 Mar 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques 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.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

Check for
updates

# Reactive Distillation of Glycolic Acid Using Heterogeneous Catalysts: Experimental Studies and Process Simulation

Carole Mutschler¹, Juliana Aparicio², Ilham Mokbel³, Mickaël Capron², Pascal Fongarland¹, Marcia Araque² and Clémence Nikitine¹*

¹CP2M, UMR CNRS 5128, University Lyon 1, CPE Lyon, Villeurbanne, France, ²University Lille, CNRS, Centrale Lille, ENSCL, University Artois, UMR 8181 - UCCS - Unité de Catalyse et Chimie du Solide, Lille, France, ³LMI, UMR CNRS 5615, University Lyon 1, Villeurbanne, France

The glycerol oxidation reaction was developed leading to selective catalysts and optimum conditions for the production of carboxylic acids such as glycolic acid. However, carboxylic acids are produced in highly diluted mixtures, challenging the recovery and purification, and resulting in high production costs, polymerization, and thermal degradation of some of the products. The protection of the acid function by esterification reaction is one of the most promising alternatives through reactive distillation (RD); this technique allows simultaneously the recovery of carboxylic acids and the elimination of most part of the water. The reactive distillation, experimental and simulation, of glycolic acid was performed, based on kinetic and thermodynamic models developed. For the thermodynamic model, binary parameters of the missing couples were determined experimentally, and the non-random two-liquid (NRTL) model was selected as the most suitable to represent the binary behavior. The kinetic study of the esterification in the presence of homogeneous and heterogeneous catalysis concluded that the heterogeneous reaction can be accurately described either by a pseudo-homogeneous model or the Langmuir-Hinshelwood (L-H) adsorption model. Reactive distillation was conducted in a distillation column filled with random packing sulfonated ion-exchange resin, Nafion NR50®, or with extruded TiO₂-Woₓ. The conversion rate of glycolic acid in reactive distillation increases from 14% without catalyst to 30% and 36% using Nafion NR50® and TiO₂-Woₓ, respectively. As opposed to the batch reactor study, the conversion rate of glycolic acid was better with TiO₂-Woₓ than with sulfonated ion-exchange resin. The better performance was related to an increase in the hydrodynamics inside the column. Tests using water in the feed confirm the hypothesis by increasing the conversion rate because of the decrease in the mass transfer resistance by reducing the average diffusion coefficient. The simulation of the reactive distillation column with ProSim® Plus showed that the yield of the ester increased operating at a low feed rate with reactive stripping. In the presence of water in the feed, nonreactive stages are required, including an enrichment region to separate water vapor.

**Keywords:** vapor-liquid equilibrium, glycolic acid, butyl glycolate, reactive distillation, reaction kinetics

## OPEN ACCESS

**Edited by:**
Patrick Cognet,
National Polytechnic Institute of
Toulouse, France

**Reviewed by:**
César Ramírez-Márquez,
Universidad de Guanajuato, Mexico
Aida Luz Villa,
Universidad de Antioquia, Colombia

*Correspondence:
Clémence Nikitine
clemence.nikitine@cpe.fr

**Specialty section:**
This article was submitted to
Green and Sustainable Chemistry,
a section of the journal
Frontiers in Chemistry.

**Received:** 31 March 2022

**Accepted:** 05 May 2022

**Published:** 15 June 2022

**Citation:**
Mutschler C, Aparicio J, Mokbel I, Capron M, Fongarland P, Araque M and Nikitine C (2022) *Reactive Distillation of Glycolic Acid Using Heterogeneous Catalysts: Experimental Studies and Process Simulation*. *Front. Chem.* 10:909380.
doi: 10.3389/fchem.2022.909380
# 1 INTRODUCTION

Glycerol is generally used as an additive or as a raw material in a wide variety of processes, including the production of food additives, tobacco, and pharmaceuticals, and synthesis of trinitroglycerin, alkyd resins, and polyurethanes. It is also used in the manufacture of lacquers, varnishes, inks, adhesives, synthetic plastics, regenerated cellulose, explosives, and other industrial uses (Veluturla et al., 2018). Despite its versatility, new ways of using glycerol have been proposed. Among these valorizations, the oxidation of glycerol in aqueous solution is an interesting alternative, leading to a mixture of not only carboxylic acids including oxalic and formic acids, but also α-hydroxy acids such as lactic, glycolic, glyceric, and tartronic acids (Skrzyńska et al., 2014). Using classical noble metals such as Au and Pt, the main products are glyceric and tartronic acids. The selectivity of the reaction turns to glycolic acid when Ag-based catalysts are used (Skrzyńska et al., 2016). In recent years, the applications (hence the economic potential) of glycolic acid have increased. The molecule has two functionalities: an alcohol group and a moderately strong acid group. These qualities make the glycolic acid perfect for a wide range of applications, for example, in the pharmaceutical industry, skin care products, and the food industry—as a flavoring agent and preservative—in adhesives and plastics, and in the textile industry as dyeing and in organic synthesis (Bianchi et al., 2005; Skrzyńska et al., 2016). However, the product obtained after the catalytic reaction is a highly diluted mixture of at least three acids: glyceric, glycolic, and formic acids. In the case of Ag-based catalyst, 69% of selectivity toward glycolic acid were reached up, with the formation of 20% of formic acid and 10% of glyceric acid. An increase in the initial glycerol concentration from 0.3 to 1.5 M, a full transformation of glycerol, was achieved with the production of 42% of glycolic acid, glyceric acid, and formic acid (15%). CO₂ production was also reported (Tavera Ruiz et al., 2021).

The production of these acids from renewable resources is of growing interest; however, new issues regarding separation emerge with the use of these resources. Unlike the petrochemical industry, the mixtures involved are highly diluted. The efficient elimination of water is a major issue in this type of process. In addition, the nonselectivity of the reactions and impurities present in the raw materials lead to the formation of coproducts and the production of highly complex mixtures. The separation processes from the petrochemical sector, such as distillation and extraction, are difficult to use and expensive. The separation by extraction or distillation is limited by the phase separation and the distribution of the components involved in the system. Therefore, the costs associated with the recovery, concentration, and purification of these carboxylic acids can represent 60%–70% of the cost of the product, making these technologies unviable.

In order to provide an integrated process, the use of reactive distillation (RD) is envisaged. It combines both reaction and separation within the same equipment. This process is clearly part of the current process evolution toward the design of multifunctional, compact, and performance-enhancing devices. RD is applied specifically to reversible chemical reactions in the

liquid phase, in which reaction equilibrium limits the conversion rate of the reactants (Talnikar and Mahajan, 2014). Reactive distillation (RD) has been proposed as a promising technique for the recovery of a short-chain carboxylic acid with high purity and high yield (Kumar et al., 2006). RD improves selectivity, increases conversion rate, allows a better heat control and an effective utilization of reaction heat, increases the scope for difficult separations, and helps to void azeotropes. As the products in RD are continuously separated from the reaction zone, no limiting chemical equilibrium can be established, and thus, the reaction is maintained at a high rate, resulting in higher yields. Other benefits of RD can include the minimization of side reactions and the utilization of the reaction heat for the mass transfer within the same column. Therefore, by acting simultaneously on the distillation and reaction (RD), both investment and operating costs can be reduced compared with conventional processes and can yield benefits such as reduction in recycling, separation optimization, and lower requirements of pumps, instrumentation, and piping (Talnikar and Mahajan, 2014) (Komesu et al., 2015).

The design of the process is based on the knowledge and understanding of the chemical reaction, the phase equilibrium, and the feasibility analysis of reactive distillation. Then, the reactive distillation column is designed and synthesized to specify the configuration of the column (number of theoretical stages, position of the reactive zone, number and position of the feed, etc.) and the operating parameters of the process (reflux rate, heating power, etc.).

Esterification reactions have been repeatedly carried out in reactive distillation processes. However, the esterification of alpha-hydroxy acids in catalytic distillation has been little investigated except for the recovery of lactic acid from fermentation by esterification (Rao et al., 2014; Kumar and Mahajani, 2007). Reactive distillation (RD) has also been used for the recovery of lactic acid from aqueous solution, using different alcohols for the esterification (e.g., ethanol, butanol, methanol, and 2-propanol) (Komesu et al., 2017; Liu et al., 2006). This work presents the experimental and simulation study of the reactive distillation of glycolic acid by esterification with butanol (Rxn. 1).

$$C_2H_4O_3 + C_4H_{10}O \leftrightarrow C_6H_{12}O_3 + H_2O \quad (Rxn1)$$

Both kinetic and thermodynamic models were developed. The reaction kinetics were determined using heterogeneous catalysts, and the catalytic results were compared with those obtained using homogeneous sulfuric acid (H₂SO₄). The thermodynamic model was fitted using vapor–liquid and liquid–liquid equilibrium data, and the coefficients were also used for the estimation of the kinetic model via the fitting of the activities. This study also presents experimental results obtained on a continuous distillation pilot, where the influence of the flow rate, feed composition, reflux rate, and the catalyst mass was studied. Finally, the experimental results were compared with the simulation results of equilibrium stagewise model using ProSim* Plus, and an optimization of operating parameters was proposed.
<figure><img src="image_4.png" /><figcaption>**FIGURE 1** | Apparatus FISCHER® LABODEST® VLE 602 used in this work. 1: Reservoir with immersed rod heater. 2: Condenser. 3: Vapor-phase temperature sensor. 4: Liquid-phase temperature sensor. 5: Mixer chamber. 6: Coolant water connection. 7: Pressure control line. 8: Liquid sampling port. 9: Vapor sampling port. 10: Condenser. 11: Liquid-phase septum. 12: Vapor-phase septum. 13: Side reservoir.</figcaption></figure>

This work develops—from bench to pilot scale—the efficient implementation of a separation technology in a more cost-effective and environmental-friendly manner, which reduces the water content before the implementation of more complex technologies for the separation of a highly diluted mixture of carboxylic acids. The implementation of separation and reaction in a single unit has high potentials, and the development of this type of technologies will accelerate the shift between petrochemical reagents and the biomass-producing pathways.

## 2 MATERIALS AND METHODS

The design and simulation of the reactive distillation column for the production of butyl glycolate was performed based on four different studies: 1) thermodynamic equilibria, in order to determine the activity coefficients used for the separation and reaction calculations; 2) $K_{Eq}$, and kinetic study, in order to determine the feasibility of reactive distillation and to integrate these studies in simulation; 3) reactive distillation experiments with a laboratory column; and 4) column simulation, in order to investigate the use of a conventional equilibrium stage model to simulate reactive distillation columns and determine the best

configuration to optimize the RD process with this chemical system. The materials used for the experimental part in the first three studies are summarized in **Supplementary Table S1**.

### 2.1 Thermodynamic Measurements

The esterification reaction of glycolic acid with butanol (Rxn. 1) represents six different binaries: 1-butanol-butyl glycolate (binary 1); 1-butanol-water (binary 2); butyl glycolate-water (binary 3); butyl glycolic acid (GA) (binary 4); 1-butanol-glycolic acid (binary 5); and water-glycolic acid (binary 6). Binary interaction parameters are only reported for binary 2, and the other equilibriums have not yet been assessed in the open literature. Besides, few experimental data concerning the saturated vapor pressure of butyl glycolate are available. Thus, vapor pressure of BG was determined and vapor-liquid equilibrium (VLE) was proposed for the determination of the binary interaction parameter of binary 1. Liquid-liquid equilibrium (LLE) was used for the determination of interaction parameters of binary 3, and solid-liquid equilibrium (SLE) was used for the determination of the binary interaction parameters of binaries 4, 5, and 6. For this study, the purity of the chemicals was verified by gas chromatography and a vacuum distillation was performed for the butyl glycolate (BG) in order to increase the purity to 99 wt.%.

#### 2.1.1 Vapor Pressure

The vapor pressure of BG was measured in two different equipment: From 1 to 79 mbar, a static device was used and described in another publication (Mokbel et al., 1995); and from 100 to 1013.25 mbar, the measurements were performed in the FISCHER® LABODEST® VLE 602 unit (Figure 1).

The VLE 602 unit is equipped with a COTTRELL pump, and it is based on the principle of the circulation method, assuring the contact between phases and a quickly reach of the equilibrium. For the vapor pressure measurements, the mixed chamber (5, Figure 1) and the side reservoir (13, Figure 1) were filled with 80–90 ml of the pure component. The determination of the equilibrium points was performed by fixing the pressure to 10 kPa and increasing the heating power in order to obtain one or two drops per second. The equilibrium conditions were considered to be attained when the variation in temperature was lower than 0.1 K during 30 min. Once the equilibrium condition was attained, a sample of liquid was taken (40 µl) for HPLC analysis. Before the test, the unit was degassed and dried under vacuum (25 mbar) during 1 h. Then, a flow of nitrogen was introduced and the starting working pressure was set.

#### 2.1.2 Vapor-Liquid Equilibrium

For BuOH-BG binary, VLE measurements at three pressures were carried out: 300.00, 700.00, and 1013.25 mbar. The measurements were performed in the FISCHER® LABODEST® VLE 602 unit as previously introduced. The mixed chamber (5, Figure 1) and the side reservoir (13, Figure 1) were filled (80–90 ml) with the less volatile component, and the second pure component is added *via* liquid-phase septum (11, Figure 1), according to previously fixed quantities. Before test, the unit was degassed and dried under vacuum (25 mbar) during 1 h. Then, a
flow of nitrogen was introduced, the starting working pressure was fixed, and the heating and stirring systems were turned on. Once the equilibrium condition was attained, samples of liquid and condensed vapor were taken (40 µl) for HPLC analysis. The accuracy of the temperature and the pressure measurements were ± 0.01 K and ± 0.1 mbar, respectively, as indicated by the supplier.

### 2.1.3 Liquid-Liquid Equilibrium

LLE data for BG-W binary were determined using magnetically agitated vials of 30 ml, placed in an isothermal oil bath equipped with a temperature control system. Mixtures of different compositions were prepared and stirred at a constant temperature for 1 h. Then, the temperature of the mixture was measured, using an electronic thermocouple (± 0.1°C), and the samples were centrifuged for 2 min at 4000 rpm. The two liquid phases were recovered separately and weighed. The composition of the samples was measured by HPLC. Samples of 400 µl were taken and diluted in 10-ml volumetric balloons with a solution of 42%-58% volume W-CAN.

### 2.1.4 Solid-Liquid Equilibrium

For solid-liquid equilibrium data, solubility measurements were performed. The solubility was measured by preparing a GA/solvent saturated mixtures in a 30-ml glass vial. The mixture was then immersed in an oil bath equipped with a temperature control system and agitated for 1 h. After agitation, the temperature was measured inside the vial, with an electronic thermocouple (± 0.1°C), and the liquid phase was recovered and centrifuged at 4000 rpm for 2 min. A sample of 400 µl was taken and weighed for HPLC analysis.

## 2.2 Batch Measurements: K<sub>Eq.</sub> and Kinetic

In general, GA esterification was carried out in a 250-ml three-necked flask, open to the atmosphere and equipped with a cooling system. The flask was placed in an oil bath under magnetic stirring, with a temperature control of within ±0.5°C. First, GA was solubilized in n-butanol at 50°C and the time zero of the reaction “t<sub>0</sub>” was considered when the catalyst was added, just before the reaction temperature was reached (T = T<sub>reaction</sub>-2°C). For K<sub>Eq.</sub> tests, the reaction was performed during 20 h using H<sub>2</sub>SO<sub>4</sub> (0.257 g). The molar ratio of butanol to glycolic acid was 3, and the reaction was performed at three temperatures: 50°C, 60°C, and 70°C. For the effect of the catalyst, the esterification was carried out at 70°C with a molar ratio of butanol to glycolic acid of 3. The reaction was performed during 4 h under reflux at constant temperature and stirring. The samples were taken at 10, 20, 30, and 60 min after t<sub>0</sub> and then every hour for 4 h. H<sub>2</sub>SO<sub>4</sub> was used as a homogeneous catalyst (0.257 g), and cation-exchange acid resins (Amberlyst 15, Amberlyst 16, Amberlyst 36, Nafion, and Dowex) were used for the heterogeneous catalytic tests (1.32 wt.%). The kinetic study was performed at different conditions: 50°C-70°C; molar ratio of butanol to GA of 1:1; 1:3; 1:6, and 1:10; and catalyst quantity of 0.3-1.5 wt.%.

For diffusional studies, the stirring rate ranged from 300 to 700 rpm, keeping the other operating parameters constant. The reaction was performed during 4 h under reflux at 70°C with

Amberlyst 36 and Nafion NR50 and a molar ratio of butanol to GA of 1:10. The effect of intra-particle diffusion in the reaction was studied for the Amberlyst 36. Two different particle sizes were screened, between 250 and 500 µm and greater than 500 µm. The same reaction conditions as the one previously mentioned were used for both tests.

## 2.3 Reactive Distillation

The experimental setup for reactive distillation studies was performed in a pilot-scale column illustrated in Figure 2. The column consists of three glass sections with an internal diameter of 32 mm and a packing height of up to 25 cm for each segment. The setup includes a total condenser, an electronic reflux splitter for reflux ratio control, and a reboiler (1 L) with an overflow outlet and a maximum heating capacity of 450 W. A heat carrier fluid (ethylene glycol) flowing in an outer wall aims to reduce heat losses in the length of the column. Several ports in the entire column, from boiler to condenser, allow internal temperature measurement. Additionally, the assembly has a pump allowing to feed the column in stages 2 and 3 after passing through a preheating system. All experiments were performed in reduced pressure at 370 mbar, and other conditions are indicated in Table 1.

## 2.4 Analysis

The samples obtained in the kinetic, VLE, LLE, and solubility measurements were analyzed by a high-performance liquid chromatograph (HPLC, SHIMADZU) with a refractive index detector (RID). The HPLC was equipped with a column Luna Omega C18 (octadecyl, inverse phase), 250 mm in length, 4.6 mm as an internal diameter, and a particle size of 5 µm. The mobile phase was acetonitrile/water 58-42 (%v/v) and acidified 0.06 g L<sup>-1</sup> with respect to water, at a constant flow rate of 0.5 ml.min<sup>-1</sup>. The column oven temperature was kept constant at 30°C for 15 min. The volume injected for each analysis was 10 µl. The samples were diluted in a solution of equal concentration of the mobile phase, without acidification. The calibration of the products was carried out in triplicates to obtain the repeatability within 0.5% in moles. Acetonitrile (HPLC grade, EMD) was used in the HPLC analysis.

The concentrations obtained by HPLC analysis allowed for calculating the evolution of the molar concentration. From the data obtained, conversion (X) and yield (Y) were calculated by the following equations:

$$X = \frac{n_{GA}^{ini} - n_{GA}}{n_{GA}^{ini}} \qquad (1)$$

$$Y = \frac{n_{BG}}{n_{GA}^{ini}} \qquad (2)$$

where n<sub>GA</sub><sup>ini</sup> and n<sub>GA</sub> represent the moles of acid at t = 0 and at the corresponding sampling time; n<sub>BG</sub> corresponds to moles of butyl glycolate produced. Due to the lack of an instrument capable to determine the water amount such as Karl Fischer, the mass balance could not be determined. However, to corroborate the reliability of the analysis method, a carbon balance was completed as illustrated by Eq. 3.
<figure><img src="image_5.png" /><figcaption>FIGURE 2 | Reactive distillation pilot design by PIGNAT®.</figcaption></figure>

$$CB(\%) = \frac{(2n_{GA}^{ini} + 4n_{BuOH}^{ini}) - (2n_{GA} + 4n_{BuOH} + 6n_{BG})}{(2n_{GA}^{ini} + 4n_{BuOH}^{ini})} \quad (3)$$

Determination of acidic sites in water for two resins was carried out using an acid-base feedback method described by Vilcocq et al. (2015). With an acid site number of 0.88 mol/g, the Nafion NR50* ion-exchange

resin has a catalytic potential of interest for the reaction studied.

## 2.5 Column Simulation

Two simulations studies have been performed. The first one is a feasibility study based on reactive residue curves. The reactive residue curve map (rCRM) is a very useful tool to obtain
<table><caption>TABLE 1 | Operating conditions for reactive distillation and results obtained with m<sub>AG</sub>/m<sub>BuOH</sub> 1/10 g/g.</caption><thead><tr><th colspan="3">Feed</th><th rowspan="2">R</th><th rowspan="2">m<sub>cata</sub> (g)</th><th rowspan="2">Flow (kg.h<sup>-1</sup>)</th><th colspan="3">Boiler</th><th rowspan="2">X<sub>exp</sub> (%)</th></tr><tr><th>F (kg h<sup>-1</sup>)</th><th colspan="2">Mass fraction</th><th>GA</th><th>BG</th><th>BuOH</th></tr></thead><tbody><tr><td>0.6</td><td>0.1</td><td>0</td><td>1</td><td>8</td><td>0.549</td><td>0.081</td><td>0.046</td><td>0.873</td><td>22</td></tr><tr><td>0.6</td><td>0.1</td><td>0</td><td>1</td><td>8</td><td>0.552</td><td>0.079</td><td>0.044</td><td>0.877</td><td>23</td></tr><tr><td>0.6</td><td>0.1</td><td>0</td><td>1</td><td>17</td><td>0.451</td><td>0.075</td><td>0.057</td><td>0.868</td><td>29</td></tr><tr><td>0.5</td><td>0.1</td><td>0</td><td>0.5</td><td>17</td><td>0.416</td><td>0.076</td><td>0.063</td><td>0.861</td><td>30</td></tr><tr><td>0.5</td><td>0.1</td><td>0</td><td>0</td><td>17</td><td>0.395</td><td>0.078</td><td>0.077</td><td>0.845</td><td>34</td></tr><tr><td>0.5</td><td>0.1</td><td>0</td><td>5</td><td>17</td><td>0.499</td><td>0.069</td><td>0.056</td><td>0.875</td><td>28</td></tr><tr><td>0.6</td><td>0.1</td><td>0</td><td>1</td><td>35</td><td>0.486</td><td>0.054</td><td>0.086</td><td>0.86</td><td>43</td></tr><tr><td>0.6</td><td>0.1</td><td>0</td><td>1</td><td>22</td><td>0.522</td><td>0.063</td><td>0.070</td><td>0.867</td><td>30</td></tr><tr><td>0.2</td><td>0.1</td><td>0</td><td>1</td><td>22</td><td>0.212</td><td>0.053</td><td>0.075</td><td>0.872</td><td>42</td></tr><tr><td>0.2</td><td>0.084</td><td>0.044</td><td>1</td><td>22</td><td>0.158</td><td>0.056</td><td>0.1</td><td>0.844</td><td>47</td></tr><tr><td>0.2</td><td>0.084</td><td>0.100</td><td>1</td><td>22</td><td>0.099</td><td>0.074</td><td>0.159</td><td>0.767</td><td>54</td></tr><tr><td>0.2</td><td>0.073</td><td>0.147</td><td>1</td><td>22</td><td>0.07</td><td>0.080</td><td>0.217</td><td>0.703</td><td>65</td></tr></tbody></table>

conceptual designs of reactive distillation columns. The methodology for determining them has already been described (Niang and Mikitenko, 1998; Barbosa and Doherty, 1988). Such maps allow us to verify the existence of stable nodes and to predict the numbers of feed and their positions. The kinetic model is incorporated into the differential equations for the computation of the reactive residue curve as:

$$ \frac{dX_i}{dt} = \frac{L}{V} (X_i - Y_i) \quad (4) $$

where X<sub>i</sub> and Y<sub>i</sub> are the transformed molar fraction in the liquid and vapor phase, respectively.

For this reactive system, BG was chosen as a reference compound. Thus, the three other transformed fractions (GA, n-butanol, and water) are calculated as below:

$$ X_i = x_i - \nu_i x_{ref} \quad (5) $$

$$ Y_i = y_i - \nu_i y_{ref} \quad (6) $$

The pressure is fixed at 380 mbar, and thus, the degree of freedom of the system is equal to 2. Therefore, by choosing two independent variables X<sub>BuOH</sub> and X<sub>GA</sub>, it is possible to solve the differential equations (Nc-1) for k initial coordinate points using the MATLAB software interfaced with the Simulink ProSim software for liquid-vapor equilibrium resolution.

The second one is a simulation of the RD process. To represent the continuous reactive distillation system, a model was developed with the ProSim® Plus software. The thermodynamic data obtained in this work were implemented in the NRTL model of the software. An equilibrium stage model was considered, and the developed kinetic law was implemented via a Cape open source file. It is worth mentioning that nonequilibrium models normally provide more details and more precise information to the simulation than the equilibrium models in the case of conventional packed distillation columns. However, the availability of reliable mass transfer correlations for the catalytic packing would be a prerequisite for the use of a nonequilibrium stage model. Such models and correlations are not available in the software used.

<figure><img src="image_6.png" /><figcaption>FIGURE 3 | BG vapor pressure reported by Crosby and Berthold (1960) (♦), reported by Forman et al. (1941) (♦), measured in FISCHER® LABODEST® VLE 602 (♦), measured in LMI (♦), and calculated (---).</figcaption></figure>

addition, the diffusional limitations of L/S mass transfer are not taken into account in the software. To overcome this problem, it was required to reduce the amount of catalyst in the reaction stage, relative to the amount of catalyst used in the experiments. Finally, in order to validate the simulation parameters, the number of equilibrium stages was set to 4 as obtained experimentally, and the conversion rate, the ester recovery rate, and this purity were compared for the same feed, distillate, and residue flow rates.

# 3 RESULTS

## 3.1 Thermodynamic Measurements

### 3.1.1 Vapor Pressure of Butyl Glycolate

The experimental and so far reported values for the BG vapor pressure are summarized in Supplementary Table S2. A regression analysis of Antoine equation parameters (Eq. 7)
<figure><img src="image_7.png" /><figcaption>FIGURE 4 | T-x-y diagram for the system butanol-butyl glycolate measured experimentally at 300 mbar (▲), 700 mbar (■), and 1013.25 mbar (●); and calculated by NRTL adjusted (---).</figcaption></figure>

was performed by minimization of the relative least sum of squares defined in Eq. 8.

$$ \ln P^o(T) = A - \frac{B}{T + C} \qquad (7) $$

$$ \sigma = \sum \left( \frac{P_{exp} - P_{cal}}{P_{exp}} \right)^2 \qquad (8) $$

In Eq. 7, $P^o(T)$ represents the vapor pressure reported in bar, $T$ is the corresponding equilibrium temperature in Kelvin, and $A, B$, and $C$ are the Antoine equation parameters. The parameters regressed are $11.94 \pm 0.3, 5140 \pm 247$, and $-30.52 \pm 8$, respectively. The standard deviation associated is 0.837. The predicted and measured/reported vapor pressure values were compared (Figure 3).

As observed, it is a good agreement between the reported and the measured data, and between the measured data and the calculated values using the regressed parameters. The maximum deviation between the measured and calculated vapor pressure was of 26 Pa between 30°C and 190°C.

### 3.1.2 Vapor-Liquid Equilibrium of Butanol-Butyl Glycolate System

The activity coefficients ($\gamma_i$) of the mixture BuOH_BG were determined according to the simplified equilibrium equation at moderate pressures (Eq. 9).

$$ \gamma_i = \frac{\varphi_i \gamma_i P}{x_i P^o(T)} \qquad (9) $$

$P$ is the total pressure of the system, $x_i$ and $\gamma_i$ are the liquid and vapor molar fractions, respectively, and $\varphi_i$ is the fugacity coefficient of component i in the mixture.

Ideal behavior was assumed for the gas phase, given the low working pressures, and the fugacity coefficient was set to 1. The VLE measured data and the calculated $\gamma_i$ are summarized in Supplementary Table S3. In all cases, the consistency of the VLE data was verified using the tool for thermodynamic consistency

available in Aspen Plus V10. For the three pressures, the consistency test—similar to the Redlich-Kister total area test (Duran et al., 2013)—was passed with tolerance intervals lower than 10%.

The VLE data were correlated by the NRTL thermodynamic model (Renon and Prausnitz, 1968) using the Aspen Plus V10 regression tool. The interaction parameters were adjusted by minimizing the maximum-likelihood objective function (Eq. 10), where $n$ and $\sigma$ are the number of data points and the standard deviation, respectively.

$$ OF = \sum_{i=1}^{n} \left[ \left( \frac{T_i^{exp} - T_i^{cal}}{\sigma^{T'}} \right)^2 + \left( \frac{x_i^{exp} - x_i^{cal}}{\sigma^x} \right)^2 + \left( \frac{y_i^{exp} - y_i^{cal}}{\sigma^y} \right)^2 \right] \qquad (10) $$

The NRTL binary interaction parameters were calculated at 300.00, 700.00, and 1013.25 mbar, and are summarized in Supplementary Table S4. The VLE predicted values and the VLE measured experimentally are compared and shown in Figure 4. It is observed that for the three pressures, it is a good agreement between the predicted and the measured equilibrium values. The total differences between the values predicted by the NRTL model and the experimental data are 1.1%.

### 3.1.3 Liquid-Liquid Equilibrium of Butyl Glycolate-Water

Figure 5 shows the LLE data of the BG-W system within a temperature range of 32.05°C-61.15°C, at atmospheric pressure. From these data and using Aspen Plus V10 regression tool, the binary parameters for the NRTL and UNIQUAC models were determined. The predicted values using the two models are also presented in Figure 5. It is observed both models can accurately predict the LL equilibria. For both models, the root mean square error (RMSE) percentage for temperature prediction was of 0.4%. The equilibria were more accurately predicted at high water (W) compositions, rather than at high butyl glycolate (BG) compositions.

The binary interaction parameters for both models are summarized in Supplementary Table S4. For NRTL model, $\alpha_{ij}$ was set to 0.2 (Toikka et al., 2018).

<figure><img src="image_8.png" /><figcaption>FIGURE 5 | T-x-y diagram for the system BG + W at 1013.25 mbar. Exp (●), NRTL (-).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**FIGURE 6** | Glycolic acid solubility in water (---), butyl glycolate (--- --), and butanol (-).</figcaption></figure>

### 3.1.4 Solid-Liquid Equilibria

The solid-liquid equilibrium is widely used in the literature to determine the activity coefficients of a solute (Negadi et al., 2006; Gutiérrez et al., 2019). For GA as solute, the equilibrium equation can be expressed as follows:

$$ 
\begin{aligned}
\ln(a_{GA}) &= \ln(x_{GA}\gamma_{GA}) \\
&= \frac{\Delta H_{fus}}{R} \left[ \frac{1}{T_{fus}} - \frac{1}{T} \right] + \left( \frac{1}{RT} \right) \int_{T}^{T_{fus}} (C_{pL} - C_{pS})dT \\
&\quad - \left( \frac{1}{R} \right) \int_{T}^{T_{fus}} \frac{(C_{pL} - C_{pS})}{T}dT - \lambda_{PT}
\end{aligned}
\quad (11) $$

$a_{GA}$, $x_{GA}$, and $\gamma_{GA}$ are respectively the activity, molar fraction, and activity coefficient of glycolic acid in saturated solution. $\Delta H_{fus}$ is the fusion enthalpy, $C_{pL}$ and $C_{pS}$ are the heat capacities of GA in the liquid and solid state, respectively; and T and $T_{fus}$ are the temperatures of the system and melting point of GA, respectively.

The activity coefficients for GA ($\gamma_{GA}$) were calculated from the experimental solubility curves (T, $x_{GA}$), obtained in water, butyl glycolate, and butanol (Figure 6). For the calculation; the melting point and enthalpy of fusion of GA were taken from the literature (Emel'yanenko et al., 2010), and the heat capacities of liquid and solid were calculated using the Aspen Plus V10 estimation tool. The values and expressions are summarized in Supplementary Table S4.

The binary parameters were fitted for the NRTL model by minimizing the difference between the calculated activity coefficients using the binary NRTL parameters and the experimental values determined using Eq. 8. The minimization was performed using the MATLAB 9.4 tool, and the results are presented in Supplementary Table S4.

## 3.2 Kinetic Studies

### 3.2.1 Equilibrium Measurement

The equilibrium constants ($K_{eq}$) were determined using Eq. 12, where $K_x$ represents the ratio of the molar fractions of the mixture at the equilibrium, and $K_y$ represents the relation of the activity coefficients calculated at the same conditions. The activity coefficients of the components of the reaction mixture

were calculated by the NRTL method using the binary parameter determined in the thermodynamic study (Supplementary Table S4). The equilibrium constants are presented in Table 2.

$$ K_{Eq} = K_x K_y = \frac{(x_{EqC} * x_{EqD})}{(x_{EqA} * x_{EqB})} \frac{(\gamma_{EqC} * \gamma_{EqD})}{(\gamma_{EqA} * \gamma_{EqB})} \quad (12) $$

The comparison between the values of $K_x$ and $K_{eq}$ concludes that at low temperatures, the system tends to follow an ideal behavior. Indeed, by comparing Eqs 13, 14, which represent the evolution of $K_{eq}$ and $K_x$ as a function of temperature, we notice that $K_{eq}$ is equal to $K_x$ for a temperature lower than or equal to 45°C. However, when the temperature increases, there is a significant difference between the ideal and nonideal system. By including the effect of nonideality in the calculation of the equilibrium constant, an increase in the total value is observed. These results are in agreement with those presented by Orjuela et al. (2012). They found, for the succinic acid and ethanol system, an increase in a factor of four for the equilibrium constant for a temperature range of 70°C-120°C; in this study, the increase in factors was found to be between 1 and 1.5. However, they consider that for their reaction, the effect of temperature was negligible. In the present study, the effect of temperature was considered in order to reduce possible sources of error as much as possible and to achieve the best possible adjustment. Considering the equilibrium constant, as a function of the activities, the curve corresponded to an endothermic reaction and the standard enthalpy of reaction can be calculated from the Van't Hoff equation (see Eq. 13). In the literature, there are no works similar to our case study; the closest result is the value reporter by Xu et al. (2012) who worked on the hydrolysis of methyl glycolate (reverse reaction). They found a standard reaction enthalpy of -15.52 kJ mol⁻¹, which is in agreement with our calculated value.

$$ \ln K_{eq} = -\frac{2858.5}{T} + 9.41 \quad (13) $$

$$ \ln K_x = -\frac{1125.7}{T} + 3.96 \quad (14) $$

### 3.2.2 Catalyzed Reaction

#### 3.2.2.1 Diffusional Limitations

The study to determine diffusion limitations was performed for the esterification reaction of GA and butanol. Different stirring rates were evaluated to determine the effect of external mass transfer resistance. At the same reaction time, the conversion rate remained constant despite the change in the stirring rate, suggesting that the external mass transfer resistance is not the speed-controlling step (Supplementary Figure S1). This result

<table><caption>**TABLE 2** | Equilibrium constants of glycolic acid esterification with butanol.</caption><thead><tr><td>Temp (°C)</td><td>Kx</td><td>Keq</td></tr></thead><tbody><tr><td>50</td><td>1.60</td><td>1.73</td></tr><tr><td>60</td><td>1.80</td><td>2.39</td></tr><tr><td>70</td><td>1.96</td><td>2.90</td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>FIGURE 7 | Effect of catalyst type on glycolic acid conversion rate with butanol.</figcaption></figure>

agrees with that of different works carried out with NAFION NR50® (Lopez et al., 2007) or Amberlyst 36 (Akyalçın, 2017).

The effect of intra-particle diffusion in the reaction was studied for the Amberlyst 36. Two different particle sizes were screened, between 250 and 500 µm and greater than 500 µm. Through the experiments, it was observed that there were no evident differences in reaction rates with change in particle sizes, which shows that internal resistance to mass transfer can be neglected. This work could not be done with Nafion® NR50, which has a larger nominal size than most ion-exchange resins and cannot be ground. For each experiment, the Weisz modulus was calculated (Supplementary Table S6), taking into account the effective diffusivity of GA since it is the limiting reagent.

$$ \phi'_{st} = \frac{r_p L^2}{D_{eAG} C_{AG}^{st}} \qquad (15) $$

where $\bar{r}_p$ represents the apparent kinetic rate; L, the characteristic length of catalyst equal to volume and surface ratio; $D_{eAG}$, the effective diffusion coefficient of glycolic acid; and $C_{AG}^{st}$, the concentration of glycolic acid on the surface of the catalyst.

For each experiment, the Weisz modulus was less than 0.1. These results are in agreement with those of the extensive literature, which state that external diffusion and intra-particle resistances are usually negligible for most reactions catalyzed by Amberlyst-type resins (Orjuela et al., 2012; Nguyen et al., 2018; Toor et al., 2011; Pappu et al., 2013). This allows the kinetic study to be carried out with the certainty of working in a kinetic regime.

### 3.2.2.2 Effect of the Catalyst

Figure 7 compares the evolution of the glycolic acid conversion rate using different commercial ion-exchange resins as catalysts. The test performed using H₂SO₄ under the same conditions is also presented in Figure 7 as a benchmark. All resins showed similar catalytic behavior, with high activities and slower kinetics compared with H₂SO₄. After 240 min, only the Amberlyst 36 achieved a conversion rate equal to the one obtained with H₂SO₄

(i.e., 86%). Amberlyst 16 also showed high activity with a conversion rate of only 5% lower compared to that obtained with the H₂SO₄. In this case, it is noticed that the activity does not depend on the morphology of the catalyst, but it can be related to its ion-exchange capacity. Dowex and Amberlyst15 despite having a capacity superior or equal to 4.7 eq.kg⁻¹ have similar conversions, both having 10 activity points less than Amberlyst 36. Amberlyst 36 has been widely used for esterification reactions due to its advantages in terms of catalytic activity and thermal stability, as already mentioned (Pappu et al., 2013; Akyalçın, 2017; Tsai et al., 2011; de Aguiar et al., 2017). Liu et al. (2006) have shown that in order for Nafion resin to have a similar activity to H₂SO₄, its polymeric structure must be modified; the authors also found that the activity of this catalyst is strongly inhibited by the presence of water. This does not seem to be the case here, which is why this resin will be studied in more detail.

### 3.2.2.3 Kinetic Model

For the kinetic study, the data were adjusted using a pseudo-homogeneous (PH) model and adsorption-based models as presented in the literature (see Table 3). Among the authors who used the PH model, Steinigeweg and Gmehling (2002) demonstrated that this model is sufficient to describe the profiles of the reactive distillation columns if there are small or medium amounts of water in the system. Orjuela et al. (2012) also demonstrated that the PH model can describe the esterification reaction between succinic acid and ethanol; they also considered the dehydration reaction of ethanol, a secondary reaction in the studied conditions. The model most used in the literature to describe the kinetic behavior of the esterification reaction in the presence of heterogeneous catalysts is the Langmuir–Hinshelwood (LH) model; this model considers that all compounds are adsorbed on the surface of the catalyst (Gangadwala et al., 2003). Finally, another model also mentioned and studied in the literature is the Eley–Rideal (ER) model. This model considers that the esterification reaction occurs between the adsorbed alcohol and non-adsorbed acid to obtain the non-absorbed ester and adsorbed water molecules (Ju et al., 2011).

The mass balance is a system of ordinary differential equations and is solved using an ordinary differential equation solver ode45 in MATLAB 9.4*. The fourth-order Runge–Kutta method was used to numerically integrate the kinetic model. The optimization of the kinetic parameters was performed by minimizing the residual sum of squares (SRS) between experimental ($x_{exp}$) and calculated ($x_{cal}$) species mole fractions using the following objective function, represented in a general way in the equation:

$$ SRS = \frac{1}{n} \sum_{Allsamples}^{Nc} (x_{exp} - x_{cal})^2 \qquad (16) $$

Where n is the number of experimental samples taken from the batch reactors in all experiments performed, and NC is the number of components considered in each sample.

To confirm that the parameters found correspond to the global optimum, the “Multistart” function was used, generating different combinations of initial points. For the
<table><caption>**TABLE 3** | Kinetic models for esterification reaction.</caption><thead><tr><th>Model</th><th>Kinetic law</th></tr></thead><tbody><tr><td>pseudo-homogeneous (PH)</td><td>The reaction takes place as if it were a homogeneous system</td></tr><tr><td>Eley–Rideal (ER)</td><td>Only one reactant and one product are found as adsorbed species</td></tr><tr><td>Langmuir–Hinshelwood (LH)</td><td>All reactants and products are present as adsorbed species</td></tr></tbody></table><table><tbody><tr><td>k<sup>+</sup> (a<sub>GA</sub>a<sub>BuOH</sub> - <sup>a<sub>BG</sub>a<sub>H<sub>2</sub>O</sub></sup>/<sub>k<sub>eq</sub></sub>)</td></tr><tr><td>k<sup>+</sup> (a<sub>GA</sub>a<sub>BuOH</sub> - <sup>a<sub>BG</sub>a<sup>2</sup>H<sub>2</sub>O</sup>/<sub>k<sub>eq</sub></sub>)</td></tr><tr><td>(1 + K<sub>GA</sub>a<sub>GA</sub> + K<sub>BG</sub>a<sub>BG</sub>)</td></tr><tr><td>k<sup>+</sup> (a<sub>GA</sub>a<sub>BuOH</sub> - <sup>a<sub>BG</sub>a<sup>2</sup>H<sub>2</sub>O</sup>/<sub>k<sub>eq</sub></sub>)</td></tr><tr><td>(1 + K<sub>GA</sub>a<sub>GA</sub> + K<sub>BuOH</sub> a<sub>BuOH</sub> + K<sub>BG</sub>a<sub>BG</sub> + K<sub>H<sub>2</sub>O</sub>a<sub>H<sub>2</sub>O</sub>)<sup>2</sup></td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>**FIGURE 8** | Parity diagram for esterification of GA with butanol and Nafion NR® as catalyst.</figcaption></figure>

adjustment of the variables, the thermodynamic consistency was considered, forcing positive values of the kinetic and adsorption constants.

As the last stage of verification, the absolute and relative error or mean relative deviation of the optimized constants was determined as presented in the following equations:

$$E_{Abs} = \frac{1}{n} \sum_{Allsamples}^{Nc} |x_{exp} - x_{cal}| \quad (17)$$

$$E_{Rel} = \frac{100*1}{n} \sum_{Allsamples}^{Nc} \frac{|x_{exp} - x_{cal}|}{x_{exp}} \quad (18)$$

Kinetic parameters were determined for two catalysts (Amberlyst 36 and Nafion NR50®) with 400 experimental data in each case (10 experiments * 10 samples * 4 compounds). The obtained parameters for each reaction with the three tested models are listed in **Supplementary Table S5. Figure 8** presents an example of the parity diagram between the experimental molar fraction and the modeled ones using the parameters given in **Supplementary Table S5**. A very good fit is obtained between the experimental compositions and the calculated ones with values within an error of 10%.

The results of the adjustment show that the error decreases when the adsorption-based models are considered, although the difference between the results obtained using the PH model and adsorption-based models varies by a maximum of 1.4%. This coincides with the conclusion of studies carried out by

Steinigeweg and Gmehling (2002) and Orjuela et al. (2012) in which it is shown that the PH model can describe with great precision this type of reaction in the presence of ion-exchange resins such as Amberlyst type.

Considering the models that represent the adsorption phenomenon, there is a slight increase in the activation energy from 53 to 56 kJ.mol⁻¹. The activation energy of the esterification reaction using Amberlyst 36 is 56 kJ.mol⁻¹; the value is similar regardless of the model used for the adjustment. However, when considering the adsorption-based models, the activation energy increases, although on this particular reaction, no information is available in the literature, the values obtained are within the range of values using the same type of catalyst, type of alcohol, and primary and carboxylic acids such as carbon numbers between C2 and C4. Steinigeweg and Gmehling (2002) reported the activation energy value of 56.65 kJ mol⁻¹ for the system between acetic acid and butanol in the presence of Amberlyst 15 as catalyst. This value could be compared with the values obtained in this study since in both cases it is a carboxylic acid with two carbons. Additionally, the glycolic acid contains an -OH group. No secondary reaction associated with the presence of this -OH group was observed under the studied conditions. It is worth noting that the activation energy increases by approximately 20% when switching from a homogeneous catalyst to a heterogeneous catalyst.

The best fit of the experimental data was achieved using the ER model; the adjustment of the experimental data demonstrates that it is only the butanol, the only molecule that strongly adsorbed from the reagents, and the water molecule in the products. However, the results obtained with the pH model and the LH model are also accurate.

## 3.3 Reactive Distillation Experiments and Simulation

### 3.3.1 Reactive Residue Map

The reactive residue curve map obtained shows that there not exist reactive azeotropes and the heteroazeotrope between butanol and water and the azeotrope between GA and BG remains despite the reaction (Supplementary Figure S2). Pure n-butanol, pure water, and azeotrope between GA and ester are saddle nodes, so they cannot be obtained by reactive distillation unless a second feed is added. On the contrary, GA and BG are stable nodes and the water/n-butanol azeotrope is an unstable node. Finally, there is a distillation boundary between the unstable node (the water/n-butanol heteroazeotrope) and the saddle point (GA/BG azeotrope). Thus, it can be concluded that
<figure><img src="image_12.png" /><figcaption>**FIGURE 9** | Steady-state establishment with feed of 0.5 kg h⁻¹, p = 380 mbar, and R = 1. (A) Temperature in column T₁ (boiler) (—), T₂ (column bottom) (—), T₃ (mid-column) (—), T₄ (column top) (—), and T₅ (condenser) (—) versus time; (B) boiler composition in BA, GB, and butanol versus time.</figcaption></figure>

with a rich n-butanol reacting mixture, a reactive separation process in order to produce pure BG is feasible. Finally, it is mainly the position of the distillation frontier that is impacted by the increase in pressure (presented in **Supplementary Figure S2**). Thus, by decreasing the pressure, the operating range for ester production is also reduced.

### 3.3.2 Repeatability, Steady State, and Mass Balance

**Figure 9** shows the composition of the residue and the temperatures along the column versus the time. Steady state is achieved after 4 h. To limit transient regimes, the start-up strategy was adjusted using the composition of the previous experiment in the boiler. For each experiment carried out, the carbon balance is verified to within 0.1%, and it can be stated that there are no carbon by-products detected during the esterification between GA and n-butanol. Among other things, polycondensation reaction does not take place. Experiment with 8 g of Nafion® N50 catalyst, a feed rate of 0.66 kg h⁻¹ with 1:10 acid/butanol mass ratio, and a reflux of 1 was performed twice in order to verify the repeatability of the experiments and the catalytic performance. As shown in **Table 1**, the conversion rates are similar, as well as the flow rates and residue compositions. The flow rate in distillate differs, possibly due to heat losses as the column head is not jacketed. Experiments can be considered repeatable. The catalyst used in these experiments will have been stable over more than 12 h of operation.

### 3.3.3 Effect of Operating Parameters

To study the effect of reflux ratio in the process, this was changed keeping constant the rest of parameters. The ratios studied were 0, 0.5, 1, and 5. There is no significant change in the molar composition in the distillate. Moreover, the recovery rate of the ester is 100% at the residue in all the experiments, which is particularly interesting from the point of view of recovery of the acid, for example, by hydrolysis of the ester. At high reflux ratios, the reactants are separated too effectively from each other, which reduces the reaction rate and causes the less conversion rate of

glycolic acid. If a high reflux ratio is used, ethanol leaves the column with the distillate. For optimum reaction conditions, the concentrations of both reactants must be high in the reaction zone. The highest conversion rate of glycolic acid is 33%, obtained at zero reflux. The column behaves like a reactive stripping column, similar to the result reported by Fields et al. (2008).

**Figure 10** shows the conversion rate of GA versus total feed flow rate. The reboiler temperature is 90°C, the feed mass ratio is 1:10, and the reflux ratio is 1. The conversion rate decreases when the total feed rate is increased from 0.2 to 0.6 kg h⁻¹, mostly due to less residence time of reactants in the reactive zone. At the highest feed flow rate, it is possible to increase the conversion rate by increasing the amount of catalyst and therefore the reactive zone. However, by having twice the contact time between liquid and solid, that is, by increasing the catalyst mass from 10 to 20 g, the conversion rate only increases by 26%. This result suggests that there is external and/or internal diffusional limitation in the column. Moreover, inadequate column hydrodynamics can lead to partial wetting of the catalyst. With the maximum amount of catalyst (i.e., a whole column section), a conversion rate of 65% is achieved, which is still below the conversion rate at thermodynamic equilibrium (85%).

**Figure 10** shows the increase in the conversion rate of GA and the distillate flow rate with the mass ratio of water in feed flow rate. The GA conversion rate increases with an increase of water in feed. The excess water causes more vapor generation in the reboiler, which leads to high vapor flow rate in the reactive zone and an increase in the distillate flow rate. On the contrary, when the vapor flow rate is higher, the contact time of liquid on solid catalyst is lower. Hence, the conversion rate should decrease. However, the estimation of the external resistance fraction (Eq. 19) and the Weisz modulus (Eq. 15) showed that by increasing the amount of water in the liquid phase, the diffusional limitations decrease (Table 4). Although the chemical regime is not yet achieved, under the conditions of our experiments the presence of water in the feed is beneficial for the recovery of GA.
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>**FIGURE 10** | Conversion rate of glycolic acid versus total feed flow rate (A) versus water mass fraction (B) (R = 1, m<sub>cat</sub> = 22 g, p = 380 mbar and total feed rate 0.2 kg h<sup>-1</sup>).</figcaption></figure>

<table><caption>**TABLE 4** | Mass transfer resistance data vs. %.mt H₂O.</caption><thead><tr><th>x<sub>H2O</sub></th><th>f<sub>ex</sub></th><th>φ'<sub>su</sub></th><th>η = η<sub>su</sub> (1 - f<sub>ex</sub>) (%)</th></tr></thead><tbody><tr><td>0</td><td>0.447</td><td>21.1</td><td>3</td></tr><tr><td>0.14</td><td>0.257</td><td>4.7</td><td>26</td></tr><tr><td>0.28</td><td>0.137</td><td>3.1</td><td>28</td></tr><tr><td>0.36</td><td>0.107</td><td>1.9</td><td>47</td></tr></tbody></table>

$$f_{ex} = \frac{\overline{r_p}L}{k_{DAG}C_{AGex}} \qquad (19)$$

### 3.3.4 Simulation

The parameters of the simulations were confirmed with the experimental results for GA conversion (Supplementary Table S7). The relative difference between simulations and experimental results is less than 6% for the conversion rate (e.g., 27% by simulation versus 25%) and recovery rate of the ester to residue (e.g., 99% versus 100%). Nevertheless, the ester purity is lower in simulations (by 15% on average, e.g., 4.9% versus 5.7%), which can be explained by a lower reboiler heat in simulations. Indeed, this results in a lower vapor rate and a less efficient separation. Moreover, the ester is highly diluted so the purity is always under 10%. A small deviation on the moles recovered in the boiler obtained by simulation can lead to a large deviation between simulation and experimental data on this parameter.

Figure 11 shows the conversion rate of GA, purity, and recovery rate of BG in residue versus the boiler heat. This heat duty has a low impact on the conversion rate. However, by increasing the heating, the BG purity increases while keeping a recovery rate above 98%. However, it is important to note that increasing the heat duty leads to a high flow of steam and increases the risk of clogging the column. This parameter must therefore be controlled. To avoid clogging, to favor the gas/liquid flow, and to keep the catalytic efficiency, it would be

<figure><img src="image_15.png" /><figcaption>**FIGURE 11** | Conversion rate of glycolic acid (---), purity (---), and recovery rate of butyl glycolate in residue (---) versus the boiler heat.</figcaption></figure>

possible to use a reactive divide wall column. The addition of a partial wall inside the column or and internal tube where the steam flows inside the column allows the purity of the products recovered at the outlet of the unit (Weinfeld et al., 2018; Von Harbou et al., 2011). Finally, a fully reactive column was compared with the initial configuration (i.e., 2 separation stages and just one reactive stage). A conversion rate of 80% and a recovery rate of 99.5% of the ester to the residue were obtained versus 68% and 99.4%, respectively, for the initial simulation. Therefore, there is no impact on the liquid phase composition of the presence of a pure separation zone. This result confirms that the column behaves as a reactive stripping column. The influence of the position of the feed on the conversion rate was studied by simulation. The reactive section is number 3
(Section 4 is the boiler). The conversion rate increases from 39% to 60% depending on whether the feed is positioned at stage 1 or 3, respectively. On the contrary, the recovery rate of ester to residue is 100% regardless of configuration. Finally, a feed in the boiler leads to a decrease in the ester recovery rate. Finally, the purity of the ester to residue is the highest (17%) when the feed is located at the top of the distillation column. Indeed, most of the butanol is recovered at the distillate in this configuration. The concentration of ester decreases strongly with a boiler feed, due to the decrease in the conversion rate. The associated simulations are provided in **Supplementary Figure S3**.

## 4 CONCLUSION

In this study, a thermodynamic model based on the NRTL model was developed and validated by fitted experimental data. A screening of catalysts was carried out. It was shown that Amberlyst 36 was the most efficient. Nafion resins were also chosen for their ease of use in the distillation columns. The kinetic study of the esterification in the presence of homogeneous and two heterogeneous catalysis was also carried out. It was concluded that the heterogeneous reaction can be accurately described either by a pseudo-homogeneous model or the Langmuir-Hinshelwood (L-H) adsorption model. A parametric study of this esterification in a reactive distillation pilot showed no significant effect of reflux ratios, but the conversion rate of GA increases with the residence time in the column. The results obtained with a different mass of catalyst suggest that there is external and/or internal diffusional limitation in the column. With the maximum amount of catalyst (i.e., a whole column section), a conversion rate of 65% is achieved, which is still below the conversion rate at thermodynamic equilibrium (85%). Based on the kinetic and thermodynamic models developed, the simulation of the reactive distillation column with ProSim Plus showed that to increase the ester yield, operating at a low feed rate with reactive stripping was sufficient.

## REFERENCES

Akyalcın, S. (2017). Kinetics of Esterification of Acetic Acid with 2-ethylhexanol in the Presence of Amberlyst 36. *Anadolu Univ. J. Sci. Technol. A- Appl. Sci. Eng.* 18, 1098-1106. doi:10.18038/aubtda.310676

Barbosa, D., and Doherty, M. F. (1988). Design and Minimum-Reflux Calculations for Single-Feed Multicomponent Reactive Distillation Columns. *Chem. Eng. Sci.* 43, 1523-1537. doi:10.1016/0009-2509(88)85144-3

Bianchi, C. L., Canton, P., Dimitratos, N., Porta, F., and Prati, L. (2005). Selective Oxidation of Glycerol with Oxygen Using Mono and Bimetallic Catalysts Based on Au, Pd and Pt Metals. *Catal. Today* 102-103 (103), 203-212. doi:10.1016/j.cattod.2005.02.003

Crosby, D. G., and Berthold, R. V. (1960). n-Butyl 5-Chloro-2-Pyrimidoxyacetate-A Plant Growth Regulator Analog. *J. Org. Chem.* 25, 1916-1919. doi:10.1021/jo01081a023

de Aguiar, V. M., de Souza, A. L. F., Galdino, F. S., da Silva, M. M. C., Teixeira, V. G., and Lachter, E. R. (2017). Sulfonated Poly(divinylbenzene) and Poly(styrene-Divinylbenzene) as Catalysts for Esterification of Fatty Acids. *Renew. Energy* 114, 725-732. doi:10.1016/j.renene.2017.07.084

Duran, J. A., Córdoba, F. P., Gil, I. D., Rodríguez, G., and Orjuela, A. (2013). Vapor-liquid Equilibrium of the Ethanol+3-Methyl-1-Butanol System at 50.66, 101.33 and 151.99kPa. *Fluid Phase Equilibria* 338, 128-134. doi:10.1016/j.fluid.2012.11.004

## DATA AVAILABILITY STATEMENT

The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.

## AUTHOR CONTRIBUTIONS

CM and JA designed the methodology, investigated and visualized the study, and wrote, reviewed, and edited the manuscript. IM contributed to thermodynamic measurements and reviewed the manuscript. MC and PF conceptualized the design, supervised the study, provided resources, acquired funding, and wrote, reviewed, and edited the manuscript. MA and CN conceptualized the design, supervised the study, provided resources, and wrote, reviewed, and edited the manuscript.

## FUNDING

This work was performed, in partnership with the SAS PIVERT, within the frame of the French Institute for the Energy Transition Institut pour la Transition Energétique (ITE) P.I.V.E.R.T (www.institut-pivert.com) selected as an Investment for the Future ("Investissements d'Avenir"). This work was supported, as part of the Investments for the Future, by the French Government under the reference ANR-001-01.

## SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem.2022.909380/full#supplementary-material

Emel'yanenko, V. N., Verevkin, S. P., Stepurko, E. N., Roganov, G. N., and Georgieva, M. K. (2010). Thermodynamic Properties of Glycolic Acid and Glycolide. *Russ. J. Phys. Chem.* 84, 1301-1308. doi:10.1134/S0036024410080054

Fields, N., Miller, D. J., Asthana, N. S., Kolah, A. K., Vu, D., and Lira, C. T. (2008). *Reactive Distillation for Esterification of Bio-Based Organic Acids*. Golden, CO, United States: National Corn Growers Association. doi:10.2172/937553

Forman, S. E., Carr, C. J., and Krantz, J. C. (1941). Alkyl Nitrites. VII. Synthesis of Some Organic Nitrites and nitrates**Contribution from the Department of Pharmacology, School of Medicine, University of Maryland. *J. Am. Pharm. Assoc. Sci. ed.* 30, 132-133. doi:10.1002/jps.3030300504

Gangadwala, J., Mankar, S., Mahajani, S., Kienle, A., and Stein, E. (2003). Esterification of Acetic Acid with Butanol in the Presence of Ion-Exchange Resins as Catalysts. *Ind. Eng. Chem. Res.* 42, 2146-2155. doi:10.1021/ie0204989

Gutiérrez, M. F., Suaza, A., Rivera, J. L., and Orjuela, A. (2019). Solid-Liquid Equilibria and Characterization of the Reaction Mixture to Produce Sucrose Palmitate in Solvent-free Media. *J. Chem. Eng. Data* 64, 2052-2061. doi:10.1021/acs.jced.8b01026

Ju, I. B., Lim, H.-W., Jeon, W., Suh, D. J., Park, M.-J., and Suh, Y.-W. (2011). Kinetic Study of Catalytic Esterification of Butyric Acid and N-Butanol over Dowex 50Wx8-400. *Chem. Eng. J.* 168, 293-302. doi:10.1016/j.cej.2010.12.086

Komesu, A., de Oliveira, J. A. R., Martins, L. H. d. S., Wolf Maciel, M. R., and Maciel Filho, R. (2017). Lactic Acid Production to Purification: A Review. *BioRes* 12, 4364-4383. doi:10.15376/biores.12.2.4364-4383
Komesu, A., Martinez, P. F. M., Lunelli, B. H., Filho, R. M., and Maciel, M. R. W. (2015). Lactic Acid Purification by Reactive Distillation System Using Design of Experiments. *Chem. Eng. Process. Process Intensif.* 95, 26–30. doi:10.1016/j.cep.2015.05.005

Kumar, R., and Mahajani, S. M. (2007). Esterification of Lactic Acid with N-Butanol by Reactive Distillation. *Ind. Eng. Chem. Res.* 46, 6873–6882. doi:10.1021/ie061274j

Kumar, R., Mahajani, S. M., Nanavati, H., and Noronha, S. B. (2006). Recovery of Lactic Acid by Batch Reactive Distillation. *J. Chem. Technol. Biotechnol.* 81, 1141–1150. doi:10.1002/jctb.1444

Liu, Y., Lotero, E., and Goodwinjr, J. (2006). A Comparison of the Esterification of Acetic Acid with Methanol Using Heterogeneous versus Homogeneous Acid Catalysis. *J. Catal.* 242, 278–286. doi:10.1016/j.jcat.2006.05.026

Lopez, D., Goodwinjr, J., and Bruce, D. (2007). Transesterification of Triacetin with Methanol on Nafion Acid Resins. *J. Catal.* 245, 381–391. doi:10.1016/j.jcat.2006.10.027

Mokbel, I., Rauzy, E., Loiseleur, H., Berro, C., and Jose, J. (1995). Vapor Pressures of 12 Alkylcyclohexanes, Cyclopentane, Butylcyclopentane and Trans-decahydronaphthalene Down to 0.5 Pa. Experimental Results, Correlation and Prediction by an Equation of State. *Fluid Phase Equilibria* 108, 103–120. doi:10.1016/0378-3812(95)02707-L

Negadi, L., Wilken, M., and Gmehling, J. (2006). Solid–Liquid Equilibria for Binary Organic Systems Containing 1-Methoxy-2-Propanol and 2-Butoxy Ethanol. *J. Chem. Eng. Data* 51, 1873–1876. doi:10.1021/je060229f

Nguyen, V. C., Bui, N. Q., Mascunan, P., Vu, T. T. H., Fongarland, P., and Essayem, N. (2018). Esterification of Aqueous Lactic Acid Solutions with Ethanol Using Carbon Solid Acid Catalysts: Amberlyst 15, Sulfonated Pyrolyzed Wood and Graphene Oxide. *Appl. Catal. A General* 552, 184–191. doi:10.1016/j.apcata.2017.12.024

Niang, M., and Mikitenko, P. (1998). Analyse de faisabilité d'une distillation réactive par la méthode des courbes de résidus. *Rev. Inst. Fr. Pét.* 53, 439–462. doi:10.2516/ogst.1998037

Orjuela, A., Yanez, A. J., Santhanakrishnan, A., Lira, C. T., and Miller, D. J. (2012). Kinetics of Mixed Succinic Acid/acetic Acid Esterification with Amberlyst 70 Ion Exchange Resin as Catalyst. *Chem. Eng. J.* 188, 98–107. doi:10.1016/j.cej.2012.01.103

Pappu, V. K. S., Kanyi, V., Santhanakrishnan, A., Lira, C. T., and Miller, D. J. (2013). Butyric Acid Esterification Kinetics over Amberlyst Solid Acid Catalysts: The Effect of Alcohol Carbon Chain Length. *Bioresour. Technol.* 130, 793–797. doi:10.1016/j.biortech.2012.12.087

Rao, V. V. B., Kumar, P. S., Sailu, C., and Rao, S. R. M. (2014). Recovery of Lactic Acid by Reactive Distillation. *J. Appl. Sci.* 14, 1289–1293. doi:10.3923/jas.2014.1289.1293

Renon, H., and Prausnitz, J. M. (1968). Local Compositions in Thermodynamic Excess Functions for Liquid Mixtures. *AIChE J.* 14, 135–144. doi:10.1002/aic.690140124

Skrzyńska, E., Wondolowska-Grabowska, A., Capron, M., and Dumeignil, F. (2014). Crude Glycerol as a Raw Material for the Liquid Phase Oxidation Reaction. *Appl. Catal. A General* 482, 245–257. doi:10.1016/j.apcata.2014.06.005

Skrzyńska, E., Zaid, S., Addad, A., Girardon, J.-S., Capron, M., and Dumeignil, F. (2016). Performance of Ag/Al2O3 Catalysts in the Liquid Phase Oxidation of Glycerol - Effect of Preparation Method and Reaction Conditions. *Catal. Sci. Technol.* 6, 3182–3196. doi:10.1039/C5CY01581B

Steinigeweg, S., and Gmehling, J. (2002). n-Butyl Acetate Synthesis via Reactive Distillation: Thermodynamic Aspects, Reaction Kinetics, Pilot-Plant

Experiments, and Simulation Studies. *Ind. Eng. Chem. Res.* 41, 5483–5490. doi:10.1021/ie020179h

Talnikar, V. D., and Mahajan, Y. S. (2014). Recovery of Acids from Dilute Streams: A Review of Process Technologies. *Korean J. Chem. Eng.* 31, 1720–1731. doi:10.1007/s11814-014-0202-4

Tavera Ruiz, C. P., Dumeignil, F., and Capron, M. (2021). Catalytic Production of Glycolic Acid from Glycerol Oxidation: An Optimization Using Response Surface Methodology. *Catalysts* 11, 257. doi:10.3390/catal11020257

Toikka, M., Trofimova, D., and Samarov, A. (2018). Liquid-liquid Equilibrium and Critical States for the Quaternary System Propionic Acid-N-Butanol-N-Butyl Propionate-Water at 303.15 K. *Fluid Phase Equilibria* 460, 17–22. doi:10.1016/j.fluid.2017.12.023

Toor, A. P., Sharma, M., Kumar, G., and Wanchoo, R. K. (2011). Kinetic Study of Esterification of Acetic Acid with N-Butanol and Isobutanol Catalyzed by Ion Exchange Resin. *Bull. Chem. React. Eng. Catal.* 6, 23. doi:10.9767/bcrec.6.1.665. 23-30

Tsai, Y.-T., Lin, H.-m., and Lee, M.-J. (2011). Kinetics Behavior of Esterification of Acetic Acid with Methanol over Amberlyst 36. *Chem. Eng. J.* 171, 1367–1372. doi:10.1016/j.cej.2011.05.049

Veluturla, S., Archna, N., Subba Rao, D., Hezil, N., Indraja, I. S., and Spoorthi, S. (2018). Catalytic Valorization of Raw Glycerol Derived from Biodiesel: a Review. *Biofuels* 9, 305–314. doi:10.1080/17597269.2016.1266234

Vilcocq, L., Spinola, V., Moniz, P., Duarte, L. C., Carvalheiro, F., Fernandes, C., et al. (2015). Acid-modified Clays as Green Catalysts for the Hydrolysis of Hemicellulosic Oligosaccharides. *Catal. Sci. Technol.* 5, 4072–4080. doi:10.1039/C5CY00195A

Von Harbou, E., Schmitt, M., Parada, S., Grossmann, C., and Hasse, H. (2011). Study of Heterogeneously Catalysed Reactive Distillation Using the D+R Tray-A Novel Type of Laboratory Equipment. *Chem. Eng. Res. Des.* 89, 1271–1280. doi:10.1016/j.cherd.2011.01.011

Weinfeld, J. A., Owens, S. A., and Eldridge, R. B. (2018). Reactive Dividing Wall Columns: A Comprehensive Review. *Chem. Eng. Process - Process Intensif.* 123, 20–33. doi:10.1016/j.cep.2017.10.019

Xu, Y., Dou, W., Zhao, Y., Huang, G., and Ma, X. (2012). Kinetics Study for Ion-Exchange-Resin Catalyzed Hydrolysis of Methyl Glycolate. *Ind. Eng. Chem. Res.* 51, 11653–11658. doi:10.1021/ie3004489

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

**Publisher's Note:** All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Copyright © 2022 Mutschler, Aparicio, Mokbel, Capron, Fongarland, Araque and Nikitine. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
GLOSSARY

a activity

c concentration, mol l⁻¹

cₚ heat capacity, J.mol⁻¹ K⁻¹

Dₐ effective diffusion coefficient, m² s⁻¹

Eₐbs absolute error, %

Eₐrel relative error, %

fₐex external fraction resistance

ΔHₐfus fusion enthalpy, J.mol⁻¹

KₐDAG convection coefficient of glycolic acid, m s⁻¹

Kₐeq equilibrium constant

Kₚ ratio of the molar fractions of the mixture

L volume-to-surface ratio, m

n number of experimental samples

P total pressure of the system, bar

P° vapor pressure, bar

R ideal gas constant, J.mol⁻¹ K⁻¹

r̄ₚ apparent kinetic rate, mol.s⁻¹ kg⁻¹

SRS residual sum of squares

T temperature, °C or K

V vapor flow, mol.s⁻¹

Xᵢ transformed molar fraction of compound i in liquid phase

xᵢ molar fraction of compound I in liquid phase

Yᵢ transformed molar fraction of compound I in vapor phase

yᵢ molar fraction of compound I in vapor phase

vᵢ algebraic stoichiometric coefficient

φ' Weisz modulus

φ fugacity coefficient

γ activity coefficient

σⱼ standard deviation of j parameter

BuOH butanol

BG butyl glycolate

calc calculated

ex external to the catalyst

exp experimental

fus melting point

GA glycolic acid

i component i

l liquid

s solid

su surface of the catalyst