Liquid-Liquid Equilibrium of Deep Eutectic Solvent-Aromatic-Aliphatic Ternary Systems: Experimental Study with COSMO Model Predictions

: Common solvents used for aromatic extraction from aliphatics typically degrade into toxic compounds, while green alternatives perform poorly compared to the state-of-the-art solvents. Deep eutectic solvents (DES) are a novel solvent type made of hydrogen bond donors (HBD) and hydrogen bond acceptors (HBA). DES have been applied in various applications, including advanced separa-tions. In this study, DES were studied experimentally and using the Conductor-like Screening Model (COSMO) to separate benzene from cyclohexane as model compounds for an aromatic:aliphatic system. Both equilibrium and kinetic studies were performed to determine the liquid liquid equilibrium (LLE) and mass transfer rate for the DES-based separation. Selected HBAs including tetrabutylammonium bromide (N4444Br), tetrahexylammonium bromide (N6666Br), choline chloride (ChCl), and methyltriphenylphosphonium bromide (METPB) were paired with HBDs including ethylene glycol (EG) and glycerol (Gly). COSMO was used, with adjustments to reﬂect DES speciﬁc interactions, to predict the liquid-liquid equilibrium (LLE). COSMO results showed that ChCl and N6666Br-based DES extracted too little benzene or too much cyclohexane, respectively, to be considered for experimental evaluation. Overall, the COSMO model predictions for LLE of EG-based DES were very accurate, with root-mean-square deviations (RMSD) below 1% for both N4444Br:EG and METPB:EG. The glycerol systems were less accurately modeled, with RMSD’s of 4% for N4444Br:Gly and 6% for METPB:Gly. The lower accuracy of glycerol system predictions fmay be due to limitations in COSMO for handling glycerol’s inﬂuence on polarizability in the DES that is not seen in EG-based DES. Mass transfer kinetics were determined experimentally for DES and the results were ﬁt to a ﬁrst order kinetics model. METPB:Gly had the highest mass transfer coefﬁcient at 0.180 min − 1 , followed by N4444Br:EG at 0.143 min − 1 . N4444Br:Gly and METPB:EG had the lowest mass transfer coefﬁcients at 0.096 min − 1 and 0.084 min − 1 , respectively. It was found that mass transfer rate was not directly related to maximum benzene solubility, as N4444Br:EG and METPB:Gly had the highest and lowest benzene removal, respectively, but had similar mass transfer coefﬁcients.


Introduction
Benzene, toluene, and xylene (collectively referred to as BTX) are common contaminants in hydrocarbon mixtures that must be removed before the hydrocarbons can be used in most processes. C 6 -C 8 mixtures derived from petroleum fractionalization typically have 10-20 mass% BTX. For these mixtures, benzene is particularly problematic as it has nearly the same boiling point as cyclohexane at around 80-81 • C. Thus, mixtures containing benzene and cyclohexane cannot be separated by thermal distillation and must be purified with solvent extraction. Sulfolane is the most common solvent used in aromatic extraction [1]. Though effective, with a 99% removal rate for aromatics, sulfolane is a hazardous chemical that breakdowns into toxic compounds in the environment [2,3]. Sulfolane is used at high temperatures above 120 • C and must therefore also be at high pressures (6-9 atm) to keep species such as benzene and cyclohexane in the liquid phase [1]. If a solvent with a lower melting point but similar polarity could be used, the separation process may become more energy efficient, and extraction can be performed at lower temperatures and pressures. Furthermore, if a green solvent is used then the hazardous byproducts made from e.g., sulfolane degradation may be avoided. Other chemicals, particularly polar compounds that are immiscible with nonpolar organics, can also be used in aromatic extraction processes. This includes compounds such as ethylene glycol, glycerol, and propylene carbonate [4]. However, these compounds do not typically perform well. For instance, ethylene glycol has a low selectivity for aromatics over aliphatic at 10:1 [5,6]. New solvents are needed that are selective while being able to be used at low temperatures.
Deep eutectic solvents (DES) are a type of multicomponent solvent, usually made from a salt or metal halide and acid, that combine to form a liquid at ambient conditions. There is usually one specific ratio, the eutectic ratio, that has a minimized melting point which is below the melting point of either individual component. Oftentimes DESs' physical properties vary by changing the component pairing as well as the ratio that the components are mixed at. This means that DES can be tailored to meet task-specific requirements like viscosity, density, or polarity just by adjusting the species and mixing ratio. Components that have been shown to form DES through hydrogen interactions include hydrogen bond acceptors (HBA) of quaternary phosphonium and ammonium salts along with hydrogen bond donors (HBD) of nearly any compound with an acidic hydrogen [7][8][9]. These components form Type III DES, which are the most promising candidates for selecting green solvents (Type I use metal chlorides and Type II use metal hydrates which often contain heavy or toxic metals).
DES are excellent candidates for the removal of BTX from aliphatic phases as common Type III DES have already been shown to have high selectivity for aromatic compounds. Kareem et al found that methyltriphenylphosphonium bromide (METPB) and ethylene glycol (EG) at a 1:2 molar ratio forms a Type III DES that is more selective for aromatics in hexane than sulfolane (selectivity ratio of 80 compared to 60 for DES and sulfolane respectively) [10]. Existing literature on DES used for the removal of BTX from aliphatic solvents includes a limited range of tetraalkylammonium salts, monohydroxyl ammonium salts, and phosphonium salts. Hadj-Kali et al. summarized the literature on this topic in a 2017 review [11]. In general, the review found that DES have a selectivity between 20 and 50 of benzene over cyclohexane at infinite dilution, so while some DES outperform the state of the art in BTX selectivity, not all DES do. Furthermore, the experimental methods used did not consider kinetics of mass transfer for the aromatics into the DES, which is needed to evaluate any process model for extraction.
To select and evaluate best performing DES, modeling methods are used to expedite experimental evaluation. Most methods need adjustments to account for the hydrogen bonding that occurs in the DES, but the extent and type of bonding varies even within the category of Type III DES [12][13][14]. Even models based on empirical data can still fit poorly though. Non-random two-liquid (NRTL) modeling, for example, has been used to predict sodium chloride solubility accurately in DES made of ethyltriphenylphosphonium bromide:ethylene glycol (1:3 molar ratio) with error below 4%. However, the same method had an error of 11-25% compared to experimental measurements for sodium chloride solubility in tetrabutylphosphonium bromide:ferric (III) chloride [12]. An additional complication is that the hydrogen bonding in the eutectic interact not only influences the polar portions of the DES components, but all types of electrostatic interactions [15]. While empirical modeling methods can be modified to consider these types of complex electrostatic influences caused by solvent dipolarity, the same method might not work for similar systems of DES. Predictive modeling based on molecular properties, such as electron density, would better fit the unique interactions seen in DES.
The Conductor Screening Model (COSMO) is part of the COSMOthermTR software suite from Dassault Systèmes SE (previously Cosmologic.de) that predicts electron density distribution in molecules via density functional theory (DFT) [16]. This method of first Processes 2021, 9, 1169 3 of 15 modeling electron density on the surface of a molecule and then examining electrostatic interactions better fits the complex ways that DES components interact with themselves and solutes. The maximum solubility of benzene and hexane have been accurately predicted in DES using COSMO, with an average root mean square deviation (RMSD) of 5% [17,18]. Though this is not complete LLE data, it is promising. Adjustments need to be made to the molecular modeling method to limit the shortcomings that COSMO and generate accurate LLE predictions for a wide variety of DES.
One common simplification for aliphatic-aromatic separation is assuming that the system is made entirely of cyclohexane and benzene. However, Hadj-Kali et al.'s 2017 review found that cyclohexane systems were particularly difficult to model in COSMO, with root-mean square deviations (RMSD) averaging around 30% [11,18]. Meanwhile, Salleh et al. were able to use COSMO to qualitatively rank various Type III DES for their ability to remove benzene from cyclohexane, although, their work was not able to achieve RMSD's below 30% using the COSMO modeling method [18]. Predicted selectivities of benzene compared to cyclohexane were generally below 5, with a few exceptions. These are far lower than the selectivities reported in the other studies included in Hadj-Kali et al.'s review, so while maximum solubility can be predicted with acceptable accuracy, more complicated phase equilibria require further corrections.
The COSMO method can be modified to increase accuracy and this has been demonstrated with aqueous phase extraction processes recently with furfural, where modeling error was below 5% [19]. Modeling the DES components as they behave in solution results in this increased predictive accuracy. HBA are modeled as wholly independent cations and anions with their ionic charges distributed across the molecule, which maintains charge neutrality in the solvent. Additionally, using methodology that better considers the effect of polarization increases accuracy. The tri-zeta valence-polarized (TZVP) method can be used, at the expense of additional computational time, to consider the polarizability effects in DES. Even with increased accuracy COSMO does not give information on the kinetics of mass transfer, which are vital for both evaluating DES and for designing processes to compete with the state-of-the-art. To the authors' knowledge there is currently no publication on the mass transfer kinetics of any DES used in a strictly hydrocarbon aliphatic-aromatic extraction process.
The COSMOtherm software suite, along with an understanding of electrostatic modeling, can be used to quickly determine whether a DES is compatible for aromatic removal from aliphatics. If this method can be completed with a higher accuracy than what has been reported in existing literature, this selection methodology would be far faster than simply experimentally evaluating all possible DES. This would only give the equilibrium conditions though, where mass transfer kinetics are also needed to complete a process model that could evaluate the performance of candidate DES. To this end, this study will use COSMO modeling for DES selection by first screening solvents then taking experimental measurements that confirm model accuracy as well as determine mass transfer kinetics. Exploring new solvent extraction processes with DES requires both equilibrium and kinetics data. Screening with accurate model prediction followed by kinetics measurements would aid in the development of DES based extraction processes.

Cosmo Methods
The following versions were used for each module of the COSMO modeling suite which is copyright of Dassault Systems (formerly copyrighted by COSMOLogic); TmoleX optimal geometry. SP 3 hybridized atoms had a tetrahedral configuration, SP 2 atoms had a trigonal configuration, etc. Bond lengths and dihedral angles were determined by the automatic 'clean up' operation in TmoleX that optimized for geometric spacing. Final values for bond lengths, dihedral angles and bond orders were determined by the TmoleX -based on molecular geometry, electron distribution and intramolecular electrostatic interactions [16]. For salts, only the cation was modeled in TmoleX; with a designated positive charge. All modeled components used the BP86 functionals for determining electron density.
Equation (1) shows the fundamental equation of COSMO that determines phase equilibrium. µ s (σ) refers to the potential of a particular surface segment caused by a difference in electron density difference. E pair is the energy associated with the interaction between the segment of interest and some nearby surface (σ'). ρ s (σ) is the sigma potential of the solution and T and k and the temperature and Boltzmann constant, respectively. This calculation is completed for each surface segment of a molecule and all of its neighbors iteratively [16,18]. By modeling the cation independently of the anion, the ions are able to be considered separately for LLE modeling, while also influencing how TmoleX considered molecular orbitals of the independent species. Additionally the total surface area considered by Equation (1) is increased, since the cation and anion are no longer spatially restricted to one another: The output of the TmoleX molecular modeling was also used with the COSMOConfx program for automatic conformer generation. The preset BP-TVZP-COSMO (Basis Pursuit function-Tri Zeta Valence Polarized-COnductor Screening MOdel) job available in version 4.3 was used to generate the conformers. The 'Liquid-Liquid Equilibrium (LLE)' function was used for generating ternary data and tie lines between a DES phase, cyclohexane as a model compound for the aliphatic phase, and benzene as a model for the aromatic phase. COSMO's default variables for bonding interaction and parameterization for temperature were not changed.
The specific DES ratios selected for evaluation was based on literature for the HBA:HBD ratio that gave the lowest melting point (eutectic ratio) or closest ratio that was a liquid at 25 • C. Additionally, evaluating specifically at the eutectic ratio also gives some consistency between the different families of DES that are being modeled. Table 1 shows the selected DES and their eutectic ratio.

Mass Transfer Kinetics and Equilibrium Measurement Methodologies
Mass transfer kinetics were measured experimentally to complement LLE values. Additionally, these experiments were used to verify the predicted equilibrium conditions. Mass transfer kinetics were measured for benzene transfer from cyclohexane into select DES; METPB:EG (1:3), METPB:Glycerol (1:3), N4444Br:EG (1:2), N4444Br:Glycerol (1:2). These DES were chosen due to promising equilibrium predictions (relatively low solubility of cyclohexane and high solubility of benzene). For the mass transfer kinetics experiments, initial concentrations of benzene ranged from 2-10 mass% with a solvent:solute ratio of Stirring rate was varied from 400-600 RPM. The organic phase was regularly sampled and analyzed using a DR6000 UV-Vis spectrometer (Hach, Loveland, CO, USA) with a quartz cuvette rated for 190-1100 nm. A Beer-Lambert law calibration at 238.8 nm was made using benzene in cyclohexane ranging from 3.5 × 10 −3 to 0.1 mass% benzene. Samples were then taken at regular intervals until concentration of the benzene was constant to find concentration at equilibrium. Both benzene's UV-Vis profile and the Beer Lambert calibration can be found in the Supplementary Information.
Change in the concentration of benzene in the cyclohexane phase over time can be defined by Equation (2): where, C(t) is the concentration of benzene, t is time in minutes, and k is the mass transfer coefficient in minutes −1 . Since the fraction of benzene removed was found to be constant across the concentration range examined, and this was consistent with literature, benzene concentration was normalized against the initial concentration [18]. This is shown in Equation (3): where, C n is the normalized concentration of benzene and C i is the initial concentration of benzene in the cyclohexane phase. Solving Equation (2) with the substitution of normalized concentration from Equation (3) yields Equation (4), which can be plotted according to Equation (5) to find the mass transfer coefficient: where C e is the equilibrium concentration of benzene in the cyclohexane phase. Note that Equation (4) has already been solved for the boundary conditions at times of zero and infinity: The equilibrium concentration of the DES-cyclohexane-benzene system was modeled in COSMOtherm using the same process as previously described. The predicted concentration of benzene was compared to the measured concentration to evaluate the accuracy of the COSMO Model. To ensure that the DES components were not affecting absorbance in the ultraviolet range, a blank test of each DES with cyclohexane was completed. None of the DES or individual HBD's had an effect on the cyclohexane baseline at 238.8 nm.

Solvatochromic Dye Measurements of Polarizability
Preparation of a solvatochromic dye was done with N,N, diethyl 4-nitroaniline (98%, Combi-Blocks, San Diego, CA 92126, USA) and acetone (99.5%, Sigma-Aldrich, St. Louis, MO 63103, USA). The dye was first tested by direct addition to several reference solvents including; glycerol (99% Acros Organics, Carlsbad, CA 92008, USA), ethylene glycol (99%, Acros Organics, Carlsbad, CA 92008, USA), deionized (DI) water, cyclohexane (99.9%, Fisher Scientific, Waltham, MA 02451, USA), Dimethyl sulfoxide (99.9%, Fisher Scientific), and acetone. The dye was then used with a solvent delivery method to expedite the experimental process by first preparing a solution of 1.4 × 10 −5 M dye in acetone and then evaporating the solvent to deliver the dye directly to a cuvette. Samples had a concentration of 2 × 10 −5 to 2 × 10 −7 M when measured by the UV-Vis. These measurements were within 0.01 abs of the direct addition measurements. The solvent delivery method was used for all measurements presented in this study.
DES were prepared to be evaluated by the solvatochromic dye. The following chemicals were used for the DES; glycerol (99% Acros Organics), ethylene glycol (99%, Acros Organics), Choline Chloride (99%, Acros Organics), methyltriphenylphosphonium bro-mide (98%, TCI), tetrahexylammonium bromide (98%, Alfa Aesar, Tewksbury, MA 01876, USA), and tetrabutylammonium bromide (98%, Alfa Aesar, Tewksbury, MA 01876, USA). The DES were prepared according to the ratios from Table 1. DES components were first combined at the specified ratios in a sealed flask at room temperature and then heated to 50 • C on a heating plate with gentle mixing. DES were allowed to cool to room temperature (20 • C) under natural cooling and were sealed with Parafilm to prevent any moisture uptake from the environment. Dipolarity was calculated from Equation (6): where V is the wavenumber of the peak of absorbance for 4,4-nitroaniline in either the solvent, cyclohexane, or dimethyl sulfoxide (DMSO) as noted [20]. In this scale the solvent dipolarity is normalized between cyclohexane and DMSO at 0 and 1, respectively.

COSMO Model DES Results: Electron Density
One of the main outputs of the COSMO model is the electron density (resulting from Equation (1)), also known as sigma or σ, that is often plotted for the molecules involved in the model. These plots are referred to as sigma profiles where the area of sigma, ρ(µ), is plotted against surface electron density (µ). The total area under the curve of these profiles is the total area of the molecule. Figure 1A shows the sigma profiles of pure solvents sulfolane, ethylene glycol, and glycerol as well as cyclohexane and benzene.
Benzene has two peaks where the electron rich peak at 0.006 e/Å 2 belongs to the delocalized aromatic π bonds and the electron deficient peak at −0.006 e/Å 2 belongs to the hydrogens. Please note that the sigma profile plot does not sum to zero across all density values. It is not uncommon for there to be more surface area on either the electron rich or deficient side. A solvent can be expected to favorably interact with benzene if it has corresponding peaks that are electron rich or deficient. Sulfolane, for example, has two distinct peaks at −0.006 and 0.013 e/Å 2 . The first peak is the electron deficient that shows the hydrogens that line the cyclic carbons. The electron rich peak is due to the sulfur and oxygen. Notably sulfolane does not have significant surface area that is charge neutral. Benzene occupies the electron rich sulfur-oxygen area of sulfolane when absorbed, while the electron deficient ring portion of sulfolane interacts with adjacent sulfolane [21].
Ethylene glycol and glycerol have very similar sigma profiles as shown in Figure 1A. Despite being the molar minority compound in all of the DES involved in this study, the HBA dominates the sigma profiles of the DES. Figure 1B,C show the sigma profiles of ethylene glycol and glycerol-based DES, respectively. For the quaternary ammonium bromide salts, N4444Br and N6666Br, the non-polar arms of the salts overshadow the polar portions of their corresponding HBD in areas that are electron deficient. N4444Br containing DES and N6666Br containing DES have electron deficient peaks at −0.004 and −0.002 e/Å 2 , respectively, well below the hydrogen bonding donating threshold of −0.0081 e/Å 2 [16,22]. The sharp, electron rich peak in all DES is due to the halide, which is the most electron dense portion of any of the DES. This peak is located at 0.017 e/Å 2 for bromide and 0.019 e/Å 2 for chloride containing DES. The method of separate modeling of the HBA cation and anion increases the total area considered by Equation (1). Most notably the area of available interaction of the halide anion increases as it is no longer at a fixed position relative to the cation.
METPB based DES have an additional peak at 0.002 e/Å 2 when paired with both glycerol and ethylene glycol. This peak is due to π bonds in the phenyl functions groups, which notably are less electron rich than benzene itself. This may be advantageous for promoting HBD-benzene interactions, as the phenyl groups are less electronegative than benzene. Choline chloride containing DES has the most electronegative peak, occurring around −0.009 e/Å 2 , just beyond the hydrogen bonding threshold. Though the nitrogen atom already has four function groups and is positively charged, this shows how uncovered it is compared to the nitrogen atom in either N4444Br or N6666Br. Choline chloride also has the lowest charge neutral characteristic of any of the DES. Given the expectedly large contribution of van der Waals interactions to benzene solubility this may inhibit benzene solubility compared to the other DES [23].
Processes 2021, 9, x FOR PEER REVIEW 7 of 17 ethylene glycol and glycerol-based DES, respectively. For the quaternary ammonium bromide salts, N4444Br and N6666Br, the non-polar arms of the salts overshadow the polar portions of their corresponding HBD in areas that are electron deficient. N4444Br containing DES and N6666Br containing DES have electron deficient peaks at −0.004 and −0.002 e/Å 2 , respectively, well below the hydrogen bonding donating threshold of −0.0081 e/Å 2 [16,22]. The sharp, electron rich peak in all DES is due to the halide, which is the most electron dense portion of any of the DES. This peak is located at 0.017 e/Å 2 for bromide and 0.019 e/Å 2 for chloride containing DES. The method of separate modeling of the HBA cation and anion increases the total area considered by Equation (1). Most notably the area of available interaction of the halide anion increases as it is no longer at a fixed position relative to the cation. Along with the electron density profiles presented, another key output from COSMO is the sigma potential, or µ(σ), profile which shows how favorable an interaction is between a molecule and a surface of a specific charge. A key difference between the sigma profile and sigma potential profile is that while the X-axis in both cases is e/Å 2 in the sigma profile this value is attributed to the molecule itself and in the sigma potential profile this value is attributed to a surface that is interacting with the molecule. These profiles can usually be considered complements of one another. A larger surface area in the electron deficient region on a sigma profile, the more favorable an interaction that is predicted on the electron rich region of a sigma potential profile. This is exactly what is seen with glycerol and ethylene glycol in Figure 2A. The electron deficient regions shown in Figure 1 for these species are due to hydrogens that are available for hydrogen bonding. These corresponds with a negative potential when interacting with electron rich surfaces, which would those that can accept hydrogen bonds. Though these thresholds for charges, above and below ±0.0081 e/Å 2 , are typically framed by considering hydrogen bonding it is important to consider when hydrogen bonding is not available. For example, sulfolane has a large surface area that is electron deficient due to the strong electron withdrawing effects of the SO 2 group, but this area cannot contribute to hydrogen bonding as there are no eligible hydrogens attached to either a halide, oxygen, or nitrogen.

Solvatochromic Dye Measurement Results
While COSMO gives preliminary information on the distribution of positive and negative surface charges it does not give polarizability. This is important as approximately 25% of the driving force for benzene solubility in polar solvents can be attributed to polar  The sigma potential profiles, µ(σ), for pure solvents are shown in Figure 2A. There are three basic types of profiles. Benzene and cyclohexane favor interactions with neutrally charged surfaces, with the minimum potential occurring around 0 e/Å 2 . Sulfolane has a local minimum between the 0 and the hydrogen bonding threshold around 0.081 e/Å 2 , showing a combination of favorable interactions with both neutrally charged surfaces and those with slight electron density. Sulfolane also has an increasingly favorable interaction displaying as a global minimum with electron deficient surfaces that can interact with the electron rich SO 2 group. Glycol and glycerol both have favorable interactions with electron rich and deficient surfaces beyond ±0.081 e/Å 2 and less favorable, but still below a potential of 0 kcal/mol × Å 2 , interactions with neutral surfaces. All of the DES have similar profiles to their respective HBD's, as shown in Figure 2B,C for ethylene glycol and glycerol, respectively. One important distinction between the pure HBD and the DES is that all DES have a peak, unfavorable potentials around the electron deficient hydrogen bonding threshold of −0.081 e/Å 2 and electron rich surfaces with approximately 0.015-0.017 e/Å 2 . The unfavorable interactions around 0.015-0.017 e/Å 2 are due to the HBA anion, which are increased due to the method of modeling cation and anion separately. The key point for these DES is that hydrogen bonding with additional solutes is less favorable than strong electrostatic interactions as the hydrogen bonding is what is allowing the DES to be a stable liquid at room temperature [13].

Solvatochromic Dye Measurement Results
While COSMO gives preliminary information on the distribution of positive and negative surface charges it does not give polarizability. This is important as approximately 25% of the driving force for benzene solubility in polar solvents can be attributed to polar interactions, while the remaining 75% can be attributed to nonpolar forces [23]. To determine dipolarity solvatochromic dyes were used. Table 2 shows the experimentally measured dipolarity of both DES and HBD. The DES all have very similar polarizability with values ranging from approximately 1.10 to 1.20. All glycerol-based DES have polarizability similar glycerol itself. METPB:EG has a slightly higher value than ethylene glycol at 1.22 ± 0.02 and 1.02 ± 0.02, respectively. N4444Br:EG and ChCl:EG are also higher than ethylene glycol at 1.09 ± 0.02 and 1.10 ± 0.04, respectively. N4444Br:Gly has a high polarizability given its HBD. DES made from N4444Br typically have π* values ranging from 0.98-1.08, with the highest values occurring at the eutectic point [24]. Unlike N4444Br DES which have relatively consistent polarizability, ChCl containing DES have been shown to have polarizability that varies significantly with different HBD's from π* values of 0.82 for ammonium thiocyanate to 1.15 with ethylene glycol as measured here [25]. There is a large range of benzene solubility amongst the selected DES, from 0.08 g benzene/g solvent for ChCl:Glycerol to full miscibility for N4444Br:EG. However, the polarizability of these DES are similar, suggesting that the polarizability of the DES is not the dominant factor for aromatic-aliphatic separation. This is consistent with literature finding that nonpolar forces account for nearly 75% of the attractive forces in benzene solubility in polar solvents [23]. Differences in solvent polarizability are still important, as will be discussed in Section 3.3. Glycerol containing DES are shown to have less accurate LLE predictions than EG containing DES. The consistent polarizability of glycerol containing DES may be due to weaker interactions than between the HBA and HBD in ethylene glycolbased DES. This is significant, as polarizability is the main difference between these DES, given that electron density distribution and surface potential (Figures 1 and 2 respectively) are nearly identical for EG and glycerol-based DES.

Liquid-Liquid Equilibrium of DES and Benzene
The equilibrium concentration of benzene in cyclohexane was both experimentally measured and predicted with COSMOtherm using the profiles shown in Section 3.1 with Equation (1). The initial concentration of benzene was between 2 and 10 mass% with a DES:organic phase mass ratio of 1:1.5. A summary of experimental and model predictions is shown in Table 3. Overall, the COSMO predictions for ethylene glycol-based DES were very accurate, with RMSD below 1% for both N4444Br:EG and METPB:EG. The glycerol systems were less accurately modeled by COSMO, with RMSD's of 4% for N4444Br:Glycerol and 6% for METPB:Glycerol. In both cases the glycerol containing DES performed slightly worse than predicted. The lower accuracy for the glycerol containing DES may be due to its consistent polarizability, as discussed in Section 3.2. COSMO may be overpredicting the interactions between glycerol and its HBA which results in higher predicted aromatic extraction than is observed. At the specified solvent:solute ratio of 1:1.5, N4444Br:EG removed the most benzene at 30 ± 1%. METPB:EG and N4444Br:Glycerol performed equally removing 19 ± 2 and 20 ± 2%, respectively. METPB:Glycerol was the worst performing DES only removing 9 ± 1% of the benzene from the organic phase. METPB DES performed worse than their N4444Br counterparts most likely due to the lack of neutrally charged surface, as shown in Figure 1B,C.
All DES solvents performed consistently for initial benzene concentrations below 10%, with small variations within 1% relative benzene removal%. COSMO predicted cyclohexane leaching for N4444Br:EG and N4444Br:Glycerol were 9 mass% of the initial cyclohexane. METPB:EG and METPB:Glycerol only leached 2 and 3 mass% of the initial cyclohexane, respectively.
Under the same conditions pure glycerol and ethylene glycol were able to remove 11 and 16 mass% benzene while leaching 2 and 3 mass% of the initial cyclohexane, respectively. This consistent amount of removal for benzene from cyclohexane regardless of benzene mass% has been observed in other experimental data for DES [18]. Figure 3 shows the predicted ternary diagrams for N4444Br:EG and METPB:Gly with cyclohexane and benzene. Whereas N4444Br:EG has full miscibility with benzene, METPG:Gly has a maximum solubility of 18 mass% benzene. All other ternary diagrams can be found in the Supplementary Information.

Mass Transfer Kinetics
The mass transfer kinetics of benzene moving from the cyclohexane phase into the DES phase were determined by fitting first order equilibrium kinetics described in Equation (5). Figure 4 shows the linear fit of each DES using Equation (5), where the slope is the mass transfer coefficient k. Each coefficient, along with the predicted and measured equilibrium concentration, is shown in Table 4. A sample set of data using an initial mass concentration of 2% benzene is shown in Figure 5 alongside the predicted concentration using the time coefficients of each DES. Overall, first order kinetics fit each examined DES well with each one showing a linear fit in Figure 5. The range of time to equilibrium is consistent with other studies looking at mass transfer into DES. A similar study examining the removal of polyphenol species from solid biomass using glycerol-based DES was shown to have extraction times of 40-60 min to reach equilibrium, but with a second order kinetics model [26].
All other DES had lower amounts of overall benzene removal then N4444Br:EG. METPB:Gly reached equilibrium the fastest with a time constant of 0.180 min −1 .

Mass Transfer Kinetics
The mass transfer kinetics of benzene moving from the cyclohexane phase into the DES phase were determined by fitting first order equilibrium kinetics described in Equation (5). Figure 4 shows the linear fit of each DES using Equation (5), where the slope is the mass transfer coefficient k. Each coefficient, along with the predicted and measured equilibrium concentration, is shown in Table 4. A sample set of data using an initial mass concentration of 2% benzene is shown in Figure 5 alongside the predicted concentration using the time coefficients of each DES. Overall, first order kinetics fit each examined DES well with each one showing a linear fit in Figure 5. The range of time to equilibrium is consistent with other studies looking at mass transfer into DES. A similar study examining the removal of polyphenol species from solid biomass using glycerol-based DES was shown to have extraction times of 40-60 min to reach equilibrium, but with a second order kinetics model [26].

Conclusions
COSMO modeling was used to both determine the electron density distribution o various DES and to predict liquid-liquid equilibrium between benzene, cyclohexane, an DES. Separate modeling of the HBA cation and anion, along with usage of the TZVP bas set resulted in more accurate predictions of LLE. Hydrogen bond donors of ethylene gly col and glycerol were found to dominate the total area considered positive and negativel charged. Meanwhile, the DES containing the hydrogen bond acceptor choline chlorid was far lower in neutrally charged surface area than tetrabutylammonium bromide, te rahexylammonium bromide, and methyltriphenylphosphonium bromide containin DES. This resulted in a far lower solubility of benzene for choline chloride containing DE than all other DES, confirming that nonpolar interactions are more dominant in determin ing benzene solubility than polar forces. Phase equilibrium measurements taken of mix tures between DES and cyclohexane solutions with 2-10 mass% benzene demonstrate that the COSMO model had a very high accuracy in predicting the final LLE if given th initial conditions. Overall, the COSMO predictions for ethylene glycol-based DES wer very accurate, with RMSD below 1% for both N4444Br:EG and METPB:EG. The LLLE o glycerol systems were less accurately modeled by COSMO, with RMSD's of 4% fo N4444Br:Glycerol and 6% for METPB:Glycerol.
Future work in applying DES should utilize the methodologies outlined in this wor to aid in the selection and evaluation of candidate solvents. Given the wide variety of DE available, and new DES that will continue to be discovered, there should be a focus o first developing methodology for evaluating the unique challenges presented by usin DES and then evaluation in specific applications. For aromatic extraction, the selection o DES that have a mixture of polar and non-polar character with high polarizability shoul be a priority. DES such as METPB:EG seem the most promising, with a high single pas removal rate of 20 mass% benzene, and a relatively low uptake of cyclohexane.  All other DES had lower amounts of overall benzene removal then N4444Br:EG. METPB:Gly reached equilibrium the fastest with a time constant of 0.180 min −1 . N4444Br:Gly and METPB:EG had lower time coefficients at 0.096 and 0.084 min −1 respectively. N4444Br:Gly and METPB:EG have similar behavior in most regards from the kinetics shown here and benzene removal (20 and 19% respectively), but had a significant difference in the amount of cyclohexane leached into the DES phase (9 and 2 mass% respectively).
Interestingly, METPB:Gly reaches equilibrium faster than METPB:EG, while N4444Br has the opposite with N4444Br:EG having faster transfer than N4444Br:Gly. METPB:Gly containing DES removed the least amount of benzene, but had the highest time coefficient. It is possible that METPB:Gly having such a low affinity for cyclohexane makes up for the lesser initial and equilibrium concentration difference (shown as the intercept of ln(C n − C e ) in Figure 4). This suggests that mass transfer is not dependent solely on the HBA or the HBD, but rather the DES. While this makes predicting DES performance from just HBA or HBD properties more difficult, it also justifies the continued use of experimental methodologies for evaluating DES even if equilibrium conditions can be accurately estimated by modeling methods such as the COSMO method used in this work.

Conclusions
COSMO modeling was used to both determine the electron density distribution of various DES and to predict liquid-liquid equilibrium between benzene, cyclohexane, and DES. Separate modeling of the HBA cation and anion, along with usage of the TZVP basis set resulted in more accurate predictions of LLE. Hydrogen bond donors of ethylene glycol and glycerol were found to dominate the total area considered positive and negatively charged. Meanwhile, the DES containing the hydrogen bond acceptor choline chloride was far lower in neutrally charged surface area than tetrabutylammonium bromide, tetrahexylammonium bromide, and methyltriphenylphosphonium bromide containing DES. This resulted in a far lower solubility of benzene for choline chloride containing DES than all other DES, confirming that nonpolar interactions are more dominant in determining benzene solubility than polar forces. Phase equilibrium measurements taken of mixtures between DES and cyclohexane solutions with 2-10 mass% benzene demonstrated that the COSMO model had a very high accuracy in predicting the final LLE if given the initial conditions. Overall, the COSMO predictions for ethylene glycol-based DES were very accurate, with RMSD below 1% for both N4444Br:EG and METPB:EG. The LLLE of glycerol systems were less accurately modeled by COSMO, with RMSD's of 4% for N 4444 Br:Glycerol and 6% for METPB:Glycerol.
Future work in applying DES should utilize the methodologies outlined in this work to aid in the selection and evaluation of candidate solvents. Given the wide variety of DES available, and new DES that will continue to be discovered, there should be a focus on first developing methodology for evaluating the unique challenges presented by using DES and then evaluation in specific applications. For aromatic extraction, the selection of DES that have a mixture of polar and non-polar character with high polarizability should be a priority. DES such as METPB:EG seem the most promising, with a high single pass removal rate of 20 mass% benzene, and a relatively low uptake of cyclohexane.