Experimental Determination and Computational Prediction of Dehydroabietic Acid Solubility in (−)-α-Pinene + (−)-β-Caryophyllene + P-Cymene System

The solubility of dehydroabietic acid in (−)-α-pinene, p-cymene, (−)-β-caryophyllene, (−)-α-pinene + p-cymene, (−)-β-caryophyllene + p-cymene and (−)-α-pinene + (−)-β-caryophyllene were determined using the laser monitoring method at atmospheric pressure. The solubility of dehydroabietic acid was positively correlated with temperature from 295.15 to 339.46 K. (−)-α-pinene, p-cymene, and (−)-β-caryophyllene were found to be suitable for the solubilization of dehydroabietic acid. In addition, the non-random two liquid (NRTL), universal quasi-chemical (UNIQUAC), modified Apelblat, modified Wilson, modified Wilson–van’t Hoff, and λh models were applied to correlate the determined solubility data. The modified Apelblat model gave the minor deviation for dehydroabietic acid in monosolvents, while the λh equation showed the best result in the binary solvents. A comparative analysis of compatibility between solutes and solvents was carried out using Hansen solubility parameters. The thermodynamic functions of ΔsolH0, ΔsolS0, ΔsolG0 were calculated according to the van’t Hoff equation, indicating that the dissolution was an entropy-driven heat absorption process. The Conductor-like Screening Model for Real Solvents (COSMO-RS) combined with an experimental value was applied to predict the reasonable solubility data of dehydroabietic acid in the selected solvents systems. The interaction energy of the dehydroabietic acid with the solvent was analyzed by COSMO-RS.


Introduction
Pine resin is an inexpensive and biodegradable natural resource, which is abundant in pine and coniferous trees [1]. It is commonly used as an antimicrobial agent [2], paint, ink toner, and for coatings [3,4]. The most important products of pine resin are rosin and turpentine obtained via distillation [5]. Rosin is non-volatile in normal temperature since the boiling point is 250 • C at 0.66 kPa [6]. Rosin is composed mostly of rosin acids and some neutral matter [7]. The main components of rosin are abietic-type resin acid and primarytype resin acid. Turpentine is a terpene mixture mainly composed of monoterpenes such as (−)-α-pinene, β-pinene, camphene, p-cymene, and sesquiterpenes such as longifolene and (−)-β-caryophyllene. The monoterpenes are volatile components in a boiling point range of 155-175 • C at atmospheric pressure, while the boiling points of sesquiterpenes are in a range of 254-256 • C at atmospheric pressure. About 1.3 million tons of pine resin is produced annually in China [8]. The collection of pine resin is the most significant and fundamental part of the rosin and turpentine industry chain. Resin tapping from live pine trees is normally performed every 10-15 days between April and October by the 'debarking' (bark-peeling) method [9]. Regular cuts are made in the trunk of the pine tree, allowing the pine resin to be secreted and collected in large quantities. The turpentine contained in the pine resin is continuously volatilized and lost in an open environment under the sun and rain, causing the turpentine content to decrease from around 30 to 15%. Due to the loss of turpentine, which can also act as a solvent for rosin acid, the rosin acid solidifies, thereby forming hard lumps that clog the resin tracts of the pine tree. Thus, the solid-liquid equilibrium of rosin acids and turpentine is essential to improve resin collection.
Dehydroabietic acid, one of the abietic-type resin acids contained in rosin, is a monocarboxylic acid with a tricyclic phenanthrene skeleton [10][11][12], as shown in Figure 1. As can be seen, the molecular backbone of dehydroabietic acid is composed of a benzene ring rather than a double bond that is prone to oxidation. As a result, dehydroabietic acid is more stable than the other abietic-type resin acid. The unique structure of dehydroabietic acid offers a wide range of applications in the fields of pharmaceuticals, pesticides, surfactants, preparation of fine chemicals, and in the synthesis of biologically active substances [13]. Currently, the data on the solubility of dehydroabietic acid in solvents are mainly focused on alcohols [14] but have not yet been reported in the literature for turpentine systems. Therefore, this study aims to investigate the solid-liquid equilibrium of the rosin acids and turpentine. Dehydroabietic acid was applied as a model compound for rosin acid, while (−)-α-pinene, p-cymene and (−)-β-caryophyllene were employed as model compounds for turpentine. The solubility of dehydroabietic acid was predicted using the COSMO-RS model. The solubility of dehydroabietic acid in (−)-α-pinene, p-cymene, (−)-β-caryophyllene, (−)-α-pinene + p-cymene, (−)-α-pinene + (−)-β-caryophyllene, and p-cymene + (−)-β-caryophyllene were determined using a laser monitoring method. The corresponding molecular thermodynamic models of solid-liquid equilibrium were established, such as the non-random two liquid (NRTL), universal quasi-chemical (UNIQUAC), modified Apelblat, λh, modified Wilson, modified Wilson-van't Hoff, etc. The compatibility of solute and solvent was described using the Hansen solubility parameter. The van't Hoff equation was set to analyze the thermodynamic properties of the dissolution process.

Solid State Properties of Dehydroabietic Acid
The differential scanning calorimetry (DSC) curve of dehydroabietic acid in its pure form is included in Figure 2. The onset melting temperature of dehydroabietic acid (T m ) is 443.22 K, which is in general agreement with the literature values [15] (442.65 to 443.55 K in the literature). The melting enthalpy (∆ fus H) value for pure dehydroabietic acid is 17.19 kJ/mol as obtained by analyzing the DSC curve with the analysis software. The X-ray diffraction (XRD) patterns of raw dehydroabietic acid and recovered equilibrated dehydroabietic acid from three pure solvents and three binary mixed solvents are graphically shown in Figure 3. It can be seen from Figure 4 that the characteristic peaks of dehydroabietic acid in all of the solvent systems are the same as the raw material. This indicates that no crystal form transition occurred during the solid-liquid phase equilibrium of dehydroabietic acid in all pure and binary mixed solvent systems.

Solubility Correlation
A regression analysis of the solubility of dehydroabietic acid in the three monosolvents was carried out with four thermodynamic models (modified Apelblat, λh, NRTL and UNIQUAC models). Three thermodynamic models (λh, modified Wilson and van't-Hoff-Wilson) were employed to regress the solubility of dehydroabietic acid in the three binary solvents. The regression parameters for six models are presented in Tables 4-10. The reliability and suitability of the regression results were evaluated by using root mean square deviation (RMSD), relative deviation (RD) and average relative deviation (ARD) [18,19]. The mathematical expressions for RD, ARD and RMSD are the Equations (1)-(3), respectively.   The RD, ARD and RMSD values are also presented in Tables 4-10 and S1,S2. The relative deviations of the solubility of dehydroabietic acid in the three monosolvents are well distributed by using the four models. However, when the NRTL model is employed to correlate p-cymene, there is a distinct difference in the relative deviation distribution. The maximum RMSD for the mixture of p-cymene is 12.00745 in the NRTL model, and the maximum value of ARD is 4.89582%. The modified Apelblat model is applied to correlate the solubility of dehydroabietic acid in (−)-α-pinene, p-cymene and (−)-β-caryophyllene with the smallest average relative deviation of the four models. The calculated results, as shown in Table 4, indicate that the values of ARD are 0.29696, 0.185545, and 0.36719%, respectively. Similarly, the 10 3 RMSD has minimal values, which are 0.49179, 0.53532, and 0.61823. The modified Apelblat model equation offers excellent correlation results between the solubility of dehydroabietic acid in pure solvents and binary solvent mixtures in comparison to other model equations.
As can be seen from Tables 8-10, the average relative deviation obtained from the Wilson-van't Hoff model for the solubility of dehydroabietic acid in (−)-β-caryophyllene + p-cymene is lower than that obtained from the Wilson model alone, indicating that the combined model is more effective in correlating the solubility of dehydroabietic acid in (−)-β-caryophyllene + p-cymene. As shown in Table S2, during the three models applied to correlate the determined solubility data of dehydroabietic acid in (−)-β-caryophyllene + p-cymene, the relative deviations of the other models are unevenly distributed, except for the λh model.
For the three models that correlated the solubility of dehydroabietic acid in p-cymene and (−)-α-pinene, the Wilson model alone showed an uneven distribution of relative deviations. However, for the binary solvents of p-cymene + (−)-α-pinene, the deviations of all models are unevenly distributed, with a particularly large gap in the Wilson model. The maximum RMSD for the binary solvents of (−)-α-pinene + (−)-β-caryophyllene is 5.41597 in the modified Wilson model, and the maximum value of ARD is 3.56843%.

Evaluation of Thermodynamic Models
In order to select the best correlation model for dehydroabietic acid in different solvents, the Akaike Information Criterion (AIC) [20] was employed. The AIC is expressed as follows: The L(θ) and k denotes the maximum likelihood value of the model, and the amount of estimable parameters to assess the model, respectively. AIC can also be expressed as: where N is the number of observations; RSS is the residual sum of squares; and x i and x ci are the experimental and calculated values of solute solubility, respectively. The AIC calculation results of each model are displayed in Table 11. The Akaike weights are used to describe the results of the application of the model more clearly and are expressed as follows [21]: where M is the amount of chosen models; AIC min is the lowest AIC value of the chosen models; and AIC i is the AIC value of the ith model. The results of comparing the model correlations are shown in Table 11. In principle, the model with the lowest AIC value is the most appropriate model. It can be seen from Table 4 that, for the monosolvents, the lowest AIC values were obtained when correlating the solubility of dehydroabietic acid in (

Thermodynamic Properties of the Solution
The thermodynamic analysis is contributed in order to understand the dissolution process of dehydroabietic acid. The relationship between solubility and the thermodynamic parameters of dissolution is derived from the van't Hoff equation [22] and expressed as follows.
where T i denotes the experimental temperature; n denotes the number of temperature points; and x 1 means the mole fraction solubility of dehydroabietic acid [23]. The fitted curve lnx 1 (1/T − 1/T hm ) with a slope value of enthalpy of solution (∆ sol H 0 ), entropy of solution (∆ sol S 0 ), and Gibbs free energy of solution (∆ sol G 0 ), can be calculated according to the following equation [24]: The Gibbs free energy was affected by enthalpy and entropy changes. In order to further understand how the thermodynamic parameters affect the Gibbs free energy, it is necessary to introduce the concept of the thermodynamic weighting coefficients δ H and δ TS [25,26], which can be defined as follows: The values of ∆ sol G 0 , ∆ sol H 0 , ∆ sol S 0 , δ H and δ TS of dehydroabietic acid in solvents were calculated according to Equations (8)- (13). The results are shown in Table 12. The values of intercept and slope were obtained via regression of the experimental solubility data. The curves of lnx 1 versus (1/T−1/T hm ) for dehydroabietic acid in three monosolvents and three binary mixed solvents are shown in Figure 6. Table 12. Values of apparent thermodynamic properties of solutions (∆ sol G 0 , ∆ sol H 0 , ∆ sol S 0 ) of dehydroabietic acid in solvents (P = 101.3 KPa) a .

Solvents
Intercept   It can be seen from Table 12 that ∆ sol H 0 is positive in all cases, suggesting that the process of dissolution of dehydroabietic acid is endothermic, which is favorable for further dissolution. Thus, the solubility becomes greater as the temperature increases, which is consistent with the solubility results. It also appears that the compatibility of dehydroabietic acid with monosolvents is basically ranked in the order of (−)-α-pinene < (−)-βcaryophyllene < p-cymene. The order of compatibility of dehydroabietic acid with binary solvents is (−)-α-pinene + (−)-β-caryophyllene < (−)-α-pinene + p-cymene < p-cymene + (−)-β-caryophyllene.
The order of ∆ sol G 0 of dehydroabietic acid in monosolvents is (−)-α-pinene > (−)-βcaryophyllene > p-cymene. The order of ∆ sol G 0 of dehydroabietic acid in binary solvents is (−)-α-pinene + (−)-β-caryophyllene > (−)-α-pinene + p-cymene > p-cymene + (−)-βcaryophyllene. This indicates that the larger the value of ∆ sol G 0 , the smaller the solubility. Conversely, the largest δ H for dehydroabietic acid in monosolvents is p-cymene, and the maximum δ H in binary solvents is p-cymene + (−)-β-caryophyllene. In addition, the value of entropy ∆ sol S 0 is also positive, indicating that the dissolution process of dehydroabietic acid is irreversible. Moreover, in the dehydroabietic acid dissolution processes, δ H is always greater than that of δ TS , which shows that enthalpy is the main contributor to the Gibbs free energy change.

Solubility Prediction
Conductor-like Screening Model for Real Solvents (COSMO-RS) predicts the thermodynamic properties of any mixed solution based on a quantitative calculation without any experimental data or group interaction parameters. The solubility of dehydroabietic acid in three monosolvents ((−)-α-pinene, p-cymene, and (−)-β-caryophyllene) and three binary solvents was predicted by COSMO-RS. The predicted results are shown in Figures 7 and 8 along with the comparison between the predicted and experimental solubility.
During the fitting of solubility for dehydroabietic acid in organic solvents by using the COSMO-RS model, the results obtained were unsatisfactory, with ARD values above 60% (the largest ARD value reaching 69.31549%) in each group when only the experimental ∆G fus (T m = 443.22 K, ∆H fus =19.17 kJ/mol) is used as the conditioning parameter of COSMO-RS. However, the best fit was obtained when using experimental reference solubility data (x i ), with the largest ARD value of 23.75923% in several sets of predictions. This may be attributed to the ∆G fus of the solid solutes, as this was estimated using quantitative structural property correlation methods. In contrast, the ∆G fus parameters obtained during the estimation process were only calculated quantitatively and were not validated by experimental data. As a result, the estimated values deviate significantly from the experimental values.  It can be concluded that the method of predicting solubility with experimental reference solubility data is more reliable. Thus, the experimental data in Table 13 were used to predict the solubility of dehydroabietic acid in binary mixtures with various mass fractions at different temperatures. The predicted results are shown in Figures 9-11. In the binary solvents containing p-cymene, it can be clearly found that the solubility of dehydroabietic acid increased with the increasing temperature and mole fraction of p-cymene. Moreover, in (−)-α-pinene + (−)-β-caryophyllene, the solubility of dehydroabietic acid increased with the increasing temperature and decreasing mole fraction of (−)-α-pinene.

Molecular Interaction Energies Analysis
The σ-profiles can be used to describe the salvation and intermolecular interactions between dehydroabietic acid and the solvents. The shielding charge density distribution (σ-profile) on the molecular surface was calculated by the TURBOMOLE software package under the BP (B88-VWN-P86) and TZVP basis sets. The σ-profiles of the solvent and solute were calculated as shown in Figure 12. As can be seen from the figure, (−)-α-pinene is the least polar compound, which is reflected in the narrow distribution of charge density around zero. In the −0.01 to −0.02 hydrogen bond donor region, two broad peaks are generated by hydroxyl oxygen for dehydroabietic acid. One broad peak is generated by carbonyl in the hydrogen bond acceptor region of 0.01 to 0.02, reflecting the fact that dehydroabietic acid is the most polar compound with both a hydrogen bond donor and a hydrogen bond acceptor. Further, in the non-polar region of −0.01 to 0.01, both p-cymene and dehydroabietic acid have a similar three peaks, as they both share the non-polar group benzene ring.    COSMO-RS calculations show that the hydrogen bond donor moment and the hydrogen bond acceptor moment of the dehydroabietic acid molecule are 2.8643 and 1.1416, respectively, indicating that the acidity of the hydrogen bond in dehydroabietic acid is greater than the basicity. Therefore, dehydroabietic acid is more likely to be solvated with a strong hydrogen basic solvent. The hydrogen bond acceptor moments of (−)-α-pinene and (−)-β-caryophyllene were 0.0008 and 0.0345, respectively. There is a more obvious stretching peak in the hydrogen bond acceptor region of the graph for (−)-β-caryophyllene. In contrast, the hydrogen bond acceptor moment for the highly soluble p-cymene is 0. However, the solubility of dehydroabietic acid in the p-cymene is the strongest. It shows that the acidity and basicity of the non-polar solvent has little influence on the solubility of dehydroabietic acid.
The different interaction energies between the dehydroabietic acid molecules and the solvent molecules have implications for the study of their solubility [27]. The COSMOtherm program can describe the intermolecular interaction forces in terms of molecular surface interactions. When performing solubility calculations with COSMO-RS, the interaction forces between the solvent and the solute were obtained; the results are shown in Figure 13. In this paper, the misfit interaction energy (H MF ), hydrogen-bond interaction energy (H HB ), van der Waals interaction energy (H vdW ), and total mean interaction energy (H tot ) of dehydroabietic acid in solvents were calculated using COSMO-RS. In general, the more positive the misfit interaction energy, the weaker the electrostatic interaction was between molecules [28]. As can be seen from Figure 13b, dehydroabietic acid has the weakest electrostatic interaction with p-cymene during dissolution in monosolvents and the weakest electrostatic interaction with p-cymene + (−)-β-caryophyllene in binary solvents. This trend is the opposite to the experimental solubility. The dissolution process is complicated, and the influencing factors may include van der Waals force, the degree of molecular association caused by the polarity of the molecule, the electrostatic interaction, solvation, the relative molecular mass of the solvent and solute, and the type and number of dissolving active groups, etc. The reason for this phenomenon of anomaly can be attributed to the weak influence of electrostatic effects on the studied system.
In addition, the more negative the hydrogen bonding interaction energy, the stronger the intermolecular hydrogen bonding [29]. As shown in Figure 13c, the strongest intermolecular hydrogen bonding is obtained with p-cymene during the dissolution in monosolvents. After the increase in temperature to 322 K, the values of the hydrogen bonding interaction energies of dehydroabietic acid in (−)-α-pinene and p-cymene are close to each other. Therefore, the strongest intermolecular hydrogen bonding in binary solvents is (−)-α-pinene + p-cymene. This suggests that hydrogen bonding interactions have a certain influence on the dissolution of dehydroabietic acid.
The van der Waals force is an attractive force present between molecules and is much weaker than chemical bonds. The higher the van der Waals force is, the higher the melting and boiling points of substances are. For substances with similar composition and structure, the van der Waals force increases with the increase in relative molecular mass. As can be seen from Figure 13d, the van der Waals interaction energy of (−)-β-caryophyllene is the strongest during the dissolution of dehydroabietic acid in monosolvents. Correspondingly, its boiling point and relative molecular mass are also the largest. In binary solvents, the van der Waals interaction energy of (−)-α-pinene + (−)-β-caryophyllene is the strongest.
As can be seen from Figure 13a, the total average interaction energy for the dissolution of dehydroabietic acid in the studied solvent is negative. It suggests that the dissolution process is thermodynamically favorable.
Further descriptions of the chemicals involved in this work can be found in Table 14.

Purification of Dehydroabietic Acid
The dehydroabietic acid was prepared with disproportionated rosin as raw material using the method reported in Ref. [30]. The disproportionated rosin was completely dissolved using 95% ethanol and then put into a stirred reactor. The disproportionated rosin was reacted with 2-ethanolamine at 35 • C for 50 min by means of a reaction-crystallization coupled with an ultrasonic wave to obtain 2-aminoethanol salt. After the reaction, the 2-aminoethanol salt of dehydroabietic acid crystals were obtained via vacuum filtration. The 2-aminoethanol salt of dehydroabietic acid crystals was dissolved with 50% ethanol and then extracted 3 times by adding isooctane to a water bath at 70 • C. Subsequently, a pure 2-aminoethanol salt of dehydroabietic acid was obtained and recrystallized 3 times from 50% ethanol. Finally, pure dehydroabietic acid with a purity of 99% was obtained via acidification with dilute hydrochloric acid. The specific rotation was +62.0 • .
The mass fraction of the materials mentioned above were analyzed using gas chromatography. Dehydroabietic acid contains carboxyl groups and is a polar substance with a high boiling point; therefore, it needs to be pre-treated by adding a 6% tetramethylammonium hydroxide methanol solution before it can be injected into the chromatograph. First, a small amount of the sample was dissolved in methanol. Then 1-2 drops of phenolphthalein indicator and then 6% tetramethylammonium hydroxide methanol solution were added until the solution was pink, which did not fade after gentle shaking. Finally, the sample was injected for analysis. The samples were analyzed by an Agilent 7820A gas chromatograph system equipped with a DB-5MS capillary column (30 m × 320 µm × 0.25 µm) and a flame ionization detector (FID). N 2 of 99.999% purity was used as the carrier gas, at a constant flow rate of 25 mL·min −1 . The injector temperature was 523.15 K and the detector temperature was 563.15 K. The sample volume was 0.4 µL. The column temperature was 373.15 K at the beginning, ramped to 543.15 K at 5 K·min −1 , and finally held at 543.15 K for 5 min. The mass fraction of 99% for dehydroabietic acid was obtained using gas chromatography (GC) analysis.

Characterization of the Solid Phase
The melting point T m and melting enthalpy ∆ fus H of dehydroabietic acid were determined using differential scanning calorimetry (NETZSCH STA 449F3, Nanining, China). The differential scanning calorimetry (DSC) instrument was first pre-calibrated based on pure indium. Subsequently, 3.31 mg of dehydroabietic acid was tested at a ramp rate of 5 K-min −1 in the range of 30-300 • C under nitrogen protection at a gas flow rate of 40 mL/min.
It is necessary to determine whether the selected solvents changed the crystal morphology of dehydroabietic acid. The crystal structure of dehydroabietic acid in pure form and when recovered from the solvent was determined by X-ray diffraction (XRD). The samples were scanned with Cu-Kα radioactivity. The XRD patterns of the sample were conducted at a scanning speed of 6 • /min from 5 to 40 • .

Reliability Analysis of Devices and Methods
Solubility is the basic thermodynamic data for chemical processes such as crystallization and separation as well as product purification. It requires a high degree of accuracy and, therefore, the reliability of the experimental methods and apparatus used for the determination of solubility needs to be checked. The solubility of potassium chloride and benzoic acid in water has been widely reported. Thus, the organic benzoic acid-water and the inorganic potassium chloride-water systems were selected as standard systems for this paper. The solubility of benzoic acid and potassium chloride in water was determined using a laser monitoring method and compared with literature values.
As can be seen in the Figure 14, the data on the solubility of benzoic acid and potassium chloride in water determined using the laser monitoring method matches the literature data [31,32]. The results show that the experimental equipment and methods used to determine solubility in this paper are reliable.

Solubility Measurements
The solubility of dehydroabietic acid in three monosolvents: ((−)-α-pinene, p-cymene, and (−)-β-caryophyllene), and three binary solvents: (p-cymene + (−)-β-caryophyllene, (−)-α-pinene + (−)-β-caryophyllene, and p-cymene + (−)-α-pinene) were measured using the laser monitoring method at P = 101.3 kPa. The device used in the laser monitoring method consists of a jacketed glass vessel (60 cm 3 ), a temperature controlling system, a laser monitoring system, and a stirring system. The jacketed glass vessel could be maintained at the desired temperature within ± 0.1 K by circulating water in a cryogenic thermostat (HS-205, Beijing, China). The laser monitoring system was a combination of a laser generator, photoelectric transformer, and a light strength display. The jacketed glass vessel and the laser monitoring system were placed in a black box to prevent the laser intensity from the light outside, improving the stability and accuracy of the experiment. In the experiments, a predetermined amount of the dehydroabietic acid and the solvent were precisely weighed with an electronic analytical balance (Mettler Toledo AB135-S, Changshu, China, accuracy: 0.0001 g) and then added to the dissolution kettle. The magnetic stirrer was then turned on to mix the solid and liquid phases. The cryogenic thermostat was also turned on, and the temperature was raised approximately 1 K·h −1 . The laser beam from the laser generator entered the dissolution kettle from one side and was accepted by the photoelectric transformer on the other side. The laser signals were converted into electrical signals and directed by the light strength display. In the early stages of the dissolution, most of the solid solutes of dehydroabietic acid were not dissolved but suspended in the liquid, which made most or even all of the incident laser light to become reflected and obscured. Hence, the value on the light intensity display was meager. With the temperature being slowly increased (less than 0.2 K/h near the equilibrium temperature), the solid particles gradually dissolved into the liquid phase, while the value on the light intensity display gradually increased. The laser intensity increased to the maximum value and tended to stabilize when the last dehydroabietic acid particle of the solute was fully dissolved, and the laser intensity reached the maximum. The temperature that corresponded to the maximum intensity was recorded as the equilibrium temperature.
The mole fraction solubility of dehydroabietic acid (x 1 ) in the selected solvent is calculated by Equations (14) and (15).
x 1 = m 1 /M 1 m 1 /M 1 + m 2 /M 2 + m 3 /M 3 (15) where m 1 indicates the mass of dehydroabietic acid; m 2 and m 3 indicate the mass of the solvents; M 1 represents the molar mass of the dehydroabietic acid; M 2 and M 3 represents the molar mass of the solvent; and x 1 denotes a molar fraction of the dehydroabietic acid in the solvent mixture. The mass fraction of the solvent is expressed by the following equation: where m 1 and m 2 denote the mass of solvent 1 and solvent 2, respectively.

Hansen Solubility Parameter (HSP)
The Hansen total solubility parameter (HSP) consists of a dispersion component (E D ), a dipole component (E P ), a hydrogen bonding component (E H ), and the total solubility parameter of the molecular structure, which can be calculated by using the principle of summation.
Dividing it by the molar volume gives the square of the total solubility parameter (δ 2 ) [17,33].
The "similarity and mutual solubility" theory was used to evaluate the variability of the HSP (∆δ t ) [34], which is expressed as Equation (20).
where δ t1 and δ t2 are the total HSP of solute and solvent, respectively. The HSPs of dehydroabietic acid and its solvent are shown in Tables 15 and S3.

The Modified Apelblat Model
The modified Apelblat equation, which ignores the effect of the solute activity coefficient at atmospheric pressure, has a simpler form with three correlation parameters [35] (A, B, C) (See Table 4). This method is suitable for the interpolation correlation of solubility data and fits the experimental data well. However, for systems with a wide range of solubility values, the correlation results of the modified Apelblat equation are biased. The formula is given as follows: This empirical equation describes the relationship between solubility and temperature, in which x 1 is the molar fraction of dehydroabietic acid dissolved in the solvents at the absolute temperature T. However, for systems with a wide range of solubility values, the correlation results of the Apelblat equation deviate considerably.

λh Model
The λh equation is a semi-empirical equation of the following Equation (22) [36,37]: where λ means the non-ideal nature of the saturated solution, T m can indicate the melting point of the solute, and h is equal to the enthalpy of dissolution per mole of solute divided by the gas constant. In the practical fitting process, the melting point, dissolution temperature, and molar fraction for dehydroabietic acid were used as data. λ, h were regressed as the parameters for binary and multivariate systems (see Table 5). The equation directly relates the solubility to temperature, avoiding the activity coefficient and greatly simplifying the calculation process.

NRTL Model
The NRTL equation can be applied to both miscible and partially miscible systems [38], and the excess Gibbs free energy for binary systems can be described as Equations (23)- (27): The expression for the solute activity coefficient calculated from the NRTL model can be expressed as [39]: The NRTL equation fitted to the experimental data contains three parameters (∆g 12 , ∆g 21 , α 12 ). In the actual calculation of this experiment (α 12 = 0.3), the NRTL contains only two adjustable parameters. ∆g 12 and ∆g 21 represent two different intermolecular interactions for the Gibbs free energy, and the parameter values are given in Table 6.

UNIQUAC Model
The UNIQUAC equation is a general analog chemical model obtained using the concept of partial composition and statistical thermodynamic methods [40]. It can be used for multivariate mixtures of non-polar and polar components to predict equilibrium data. The model is easy to calculate and applicable to partially miscible systems. Although the model does not accurately describe the data for all systems, it has a wide range of applications in the industry that is well established. For a binary mixture system, the activity coefficient of component i can be represented by the Equation (30): The parameters can be obtained from Equations (31)- (34).
z is the number of interacting molecules in close proximity around the central molecule, taken as z = 10. θ i and ϕ i are the average surface area fraction and average volume fraction of component i, respectively. r i and q i are the structural parameters of pure component i, which describe the molecular size and surface area, respectively. r i and q i can be used to calculate the R i and Q i of each group by the group contribution method. ∆u 12 and ∆u 21 are the two adjustable parameters of the UNIQUAC equation (see Table 7). The values of q and r can be calculated by the equations [41]: Q i is the surface area parameter of group i, and R i is the volume parameter of group i. Q i and R i can be obtained from the van der Waals volume A vm and surface area V vm . The results of the calculation of r and q are shown in the Table S3.

The Modified Wilson Model
Comer and Kopecni proposed the equation Wilson for the correction of solute solubility in mixed solvents, which is simple in structure and has only two correlation parameters for binary mixed solvent systems. The equation has the following form for binary mixed solvent systems [42]: where λ 12 and λ 21 are the two correlation parameters of the modified Wilson equation with some predictive capability (see Table 9); x m is the molar fraction of solute; and w 1 and w 2 are the mass fractions of the solvent mixture. This equation is obtained by correlating a portion of the data, which can be used to predict the solubility values of solutes at other compositions and temperatures.

The Modified Wilson-van't Hoff Model
The van't Hoff equation was introduced into the Wilson model, which represents an integrated model for fitting solute solubility data in a co-solvent system. The Wilson-van't Hoff equation can be expressed as: A 1 , B 1 , A 2 , and B 2 are the van't Hoff model's parameters acquired by fitting the solubility of the solute in the pure solvent at different temperatures. The relevant parameters are listed in Table 10.

COSMO-RS Prediction
COSMOtherm software is based on the COSMO-RS theory used to predict saturated vapor pressure, solubility, vapor-liquid phase diagrams, solid-liquid phase diagrams, etc. The software is based on quantitative calculations to forecast the thermodynamic properties of fluids in the absence of experimental data.
The COSMO-RS analysis was performed using COSMOthermX software (version 2020). The theoretical basis of this software is the COSMO-RS model. The molecular conformation of dehydroabietic acid was first optimized by the DMOL3_PBE_20.ctd program package with the PBE DFT function. Then, the COSMO file of dehydroabietic acid was obtained by the TURBOMOLE2021 package, while that of the solvents were obtained from the COSMO-RS database DNP basis set. From the compound information in these COSMO files, the sigma surfaces of the dehydroabietic acid and the solvents could be shown and transformed into corresponding sigma profiles. Finally, the molecular interaction energies and solubility values were calculated by inputting the above COSMO files into the COSMOthermX software.
The calculation of ∆G f us must be considered when predicting the solubility of dehydroabietic acid with COSMOtherm. COSMO-RS predicts the solubility mainly from the following equation [43]: where x i sol is the mole fraction of solute dissolved in the corresponding solvents. µ i pure and µ i s represent the chemical potential of pure solute and the infinite dilution chemical potential of solute in the corresponding solvents, respectively. The value of ∆G f us reflects the Gibbs free energy of fusion, which significantly affects the precision of the estimated solubility values. The ∆G f us was obtained from the DSC data or the reference data.

Conclusions
In this paper, the solubility of dehydroabietic acid was investigated in three monosolvents ((−)-α-pinene, (−)-β-caryophyllene, and p-cymene) and three mixed solvents ((−)-α-pinene + p-cymene, (−)-β-caryophyllene + p-cymene, and (−)-α-pinene + (−)-βcaryophyllene). The experimental results show that with increasing temperature, the solubility of dehydroabietic acid in the selected solvents increases. The solubility of dehydroabietic acid in p-cymene is significantly greater than that in (−)-α-pinene and (−)-β-caryophyllene. In addition, the solubility of dehydroabietic acid in binary solvents containing p-cymene is significantly higher than that not containing p-cymene. An analysis of the solubility of dehydroabietic acid with the solvent using the HSP showed that the best miscibility with p-cymene was achieved. The solubility of dehydroabietic acid in selected solvents systems was fitted using the NRTL, UNIQUAC, modified Apelblat, modified Wilson, modified Wilson-van't Hoff, and λh models; with the modified Apelblat model showing the best correlation with the lowest AIC value for dehydroabietic acid in monosolvents. In addition, the λh equation was the most appropriate to correlate the solubility of dehydroabietic acid in binary solvents.
The solubility of dehydroabietic acid in selected solvents systems was predicted by COSMO-RS. The analysis of the sigma profiles of dehydroabietic acid with solvent molecules, and the intermolecular interactions, showed that the synergetic effect of mul-tiple interaction energies offers a significant contribution to the dissolution process of dehydroabietic acid. ∆ sol H 0 , ∆ sol S 0 , and ∆ sol G 0 were obtained from the thermodynamic analysis according to the van't Hoff equation, demonstrating the dissolution is an irreversible heat-absorbing process.