Supercritical CO₂ Extraction of Extracted Oil from Pistacia lentiscus L.: Mathematical Modeling, Economic Evaluation and Scale-Up

Abstract

In this study, the extracted oil of Pistacia lentiscus L. the Tunis region was extracted using supercritical carbon dioxide (SC-CO₂) extraction containing different major components in the oil such as α-pinene (32%) and terpinene-4-ol (13%). The investigation of the effect of different variables on the extraction yield with 5% level of confidence interval showed that the CO₂ pressure was the main significant variable to influence the oil yield. In order to better understand the phenomena, three parameters were considered to adjust all parameters of broken and intact cell (BIC) model: grinding efficiency (G), the internal mass transfer parameter ($k_S{a}_0$), and the external mass transfer parameter ($k_f{a}_0$), which were estimated by experimental extraction curves to calculate the diffusion coefficient. From an economic point of view, we found out that the high cost of production of the extracted oil was due to the low mass of extracted oil obtained from this type of plant.

1. Introduction

Plant species with medicinal properties are becoming crucially essential in research due to their utilities in many areas such as traditional medicine, food, cosmetics and pharmaceuticals [1]. In Tunisia, a species of plant called Pistacia lentiscus L., commonly known as “Dharw” in the Maghreb and “Elustaka” in the Middle East, has medicinal properties in its oil extract. The extracted oil has numerous medicianal properties, for instance anti-atherogenic [2], anti-inflammatory [3,4,5], anti-oxidant [6,7,8,9,10], antimicrobial [9,11,12,13], hypotensive [14,15], anticancer [5,16], anti-arthritis [17], wound treatment, and anti-asthmatic and anthelmintic activities [18,19,20].

The essential oil can be extracted from specific plants by several extraction techniques: hydro-distillation [21], Soxhlet extraction [22] and supercritical carbon dioxide (SC-CO₂) extraction [23,24,25]. The major objective of the extraction process is to provide a more concentrated form of the desired material. Although the cost never compromises the quality, it can be a decisive factor in choosing an adequate process. However, the extraction effectiveness and the safety process must be priorities. In fact, as the limits of solvent residues are increasingly subjected to review, the supercritical solvent used to extract the oil in supercritical fluid extraction (SFE) can replace toxic solvents. To obtain the desired product which meets the needs of consumers, this can be a key element in defining the extract quality.

Modeling of SC-CO₂ extraction from natural matter is a very important tool and it represents a challenge in the research field. Many models have been developed and are currently used in supercritical fluid extraction [26,27,28,29,30]. The most widely used model was developed by Sovová [28,29]. This model is known as the broken and intact cell (BIC) model. It assumes diffusive and convective transport phenomena during the extraction that occurs in three periods. Also, he developed an analytical solution to estimate the extraction parameters by comparing the results of the extraction curves calculated by the model with the experimental data. Stástova et al. [30] simplified the model developed by Sovová [29] to describe the buckthorn extraction curves to evaluate grinding efficiency, mass transfer coefficients, and flow asymmetry.

To the best of our knowledge, the bulk of the prior research involving the extraction of essential oils of Pistacia lentiscus are focused on hydro-distillation extraction. Few studies have reported the supercritical extraction process for the plant of Pistacia lentiscus (PL). Congiu et al. [31] isolated the essential oil from the leaves and the berries of Pistacia lentiscus collected in the region of Sardinia (Costa Rey and Capoterra) using a supercritical CO₂ extraction technique coupled with the fractional separation technique (SFE). Thus, the authors separated the essential oil from cuticular waxes. The obtained yields of the extracted oil from the leaves and the berries were 0.45% and 0.20%, respectively, with the presence of major compounds such as β-caryophyllene, germacrene, β-myrcene, and α-pinene.

It is well-known that the quantity of the extract yield from supercritical CO₂ extraction is affected by several operating parameters such as the CO₂ pressure, CO₂ mass flow rate, time of extraction, and average particle size. In order to optimize the extract yield of the oil produced from the leaves, an efficient way might be to systematically create models around the key ingredient levels of the product via some type of response surface experimental design [32]. Response surface methodology (RSM) is a collection of mathematical and statistical techniques that make a full description of the effect of independent variables near the optimum conditions [21,33,34]. Several classes of treatment structures can be used as RSM experiments [35].

The main objective of our study is to use the response surface methodology (RSM) to study the effect of three operating conditions (pressure, average particle size, and CO₂ flow rate) and their interaction on the extract yield. We also aim to study the influence of operating parameters on mass transfer by evaluating a process applying the BIC model proposed by Sovová [29] in 1994 which was simplified later by Stástova et al. [30] on the essential oil extraction curves that are acquired from Pistacia lentiscus leaves growing in the northern parts of Tunisia. In the end, we obtained an economic evaluation in the scale-up process for the CO₂ extract of these plant leaves.

2. Materials and Methods

2.1. Materials

The Pistacia lentiscus (PL) plant leaves were obtained from a local market in Tunis (Tunisia). The leaves were air-dried under a controlled temperature of 37 °C for 48 h, then stored in vacuum-sealed plastic bags under refrigeration prior to extraction. Immediately prior to supercritical fluid extraction, only the leaves of the samples were used, and they were rounded in a blender to get particle powder with diameters 220 µm and 650 µm. The chemicals used were absolute ethanol (Carlo Erba, Val de Reuil, France), ultrapure water and carbon dioxide (Messer Group, Nancy, France, 99.95%).

2.2. The Procedure of Supercritical Fluid Extraction

The main purpose of using supercritical fluid extraction (SFE) was to obtain the extracted oils. The SFE apparatus is shown in Figure 1.

The extract was collected as a function of time during the process through valves located at the base of the separators. The samples were weighed after 30 min of collection to avoid measuring CO₂ remaining in the bottle.

The temperature of extraction was maintained constant in all the experiments (40 °C) to prevent the heat degradation of thermolabile components in the extracted oil.

Table 1. Conditions of supercritical fluid extraction.

Experiments	P [bar]	dP [µm]	${\mathit{Q}}_{\mathbf{C}{\mathbf{O}}_2}$ [kg h−1] *	${\mathit{\rho }}_{\mathbf{C}{\mathbf{O}}_2}$ [kg m−3]	${\mathit{\mu }}_{\mathbf{C}{\mathbf{O}}_2}·{10}^5$ [kg m s−1]
1	220	650	0.604	857.20	8.18
2	220	220	0.602	857.20	8.18
3	80	650	1.202	277.90	2.23
4	80	650	1.209	277.90	2.23
5	80	220	0.602	277.90	2.23
6	80	650	0.603	277.90	2.23
7	140	650	1.204	763.27	6.51
8	180	220	0.913	819.51	7.45
9	180	220	0.904	819.51	7.45
10	180	650	0.908	819.51	7.45

* Deviation ± 0.03 kg h⁻¹ of CO₂.

shows the experimental conditions of the SFE unit where P is CO₂ pressure, $Q_{{\mathrm{CO}}_2}$ is CO₂ flow rate, d_P is the average particle size of the leave, ${\rho }_{{\mathrm{CO}}_2}$ is CO₂ density, and ${\mu }_{{\mathrm{CO}}_2}$ is CO₂ velocity. Observing Table 1, we noted that the conditions of extraction used during these experiments show that the CO₂ flow was manually controlled, and the estimated variance of the experiments is between 2.5% to 5% of the average flow.

Referring to a previous publication [29] and in order to ensure the solubility of major compounds [36,37,38,39], the collection of extracted oils and waxes was conducted using separators. The first separator was maintained at a low temperature (−5 °C) with the same extraction pressure as the experiment to precipitate waxes while the second separator was maintained at 30 °C and 40 bar for oil extract collection.

The bulk density of milled Pistacia leaves was about 291 kg m⁻³, and the void fraction of the bed was equal to 0.53. Glass beads were placed on the bottom of the extractor, the powder of Pistacia leaves (23 ± 0.05 g) was placed above them and another layer of glass bead was put on the top. In addition, two filters (frits <15 µm) were used in both the inlet and the outlet of the extracting vessel.

2.3. Analysis

Gas Chromatography-Flame Ionization Detector/Mass Spectrometry (GC-FID/MS)

GC-FID analysis was carried out with a Shimadzu GC2010 Plus (Nancy, France), equipped with an HP-5 capillary column (Shimadzu, Nancy, France, with dimension: 30 m × 0.25 mm, film thickness 0.25 µm). The injector and the detector were set at 250 and 300 °C. The temperature column was programmed at 50 °C for 1 min then gradually increased to 270 °C at 3 °C/min. Next, it was held for 5 min then increased to 300 °C at 20 °C/min and subsequently held for 5 min. The split ratio was 5:1 whereas the split flow was equal to 10 mL/min. Nitrogen was used as a carrier gas with a constant pressure of 100 kPa. The GC was also equipped with an Auto-Injector (Shimadzu, Nancy, France, AOC-20i) and the injected volume was equal to 1 µL.

For GC-FID-MS analysis, a Shimadzu GCMS-GC2010-QP2010 Plus equipped with a DB5-MS capillary column (Shimadzu, Nancy, France, with dimension: 30 m × 0.25 mm, film thickness 0.25 µm) was exploited. The injector and the detector were set at 250 °C. The oven temperature was programmed at 50 °C for 1 min, gradually increased to 250 °C at 5 °C/min. It was held for 10 min then increased later to 270 °C at 5 °C/min then held for 5 min. After that, it was increased to 280 °C at 5 °C/min and held for 10 more minutes. Although the split ratio was 10:1, the split flow was equal to 10 mL/min. In this process, helium was used as a carrier gas with a constant speed of 1.69 mL/min. The GC was also equipped with an Auto-Injector (Shimadzu, Nancy, France, AOC-5000) and the injected volume was equal to 1 µL. Mass units were monitored from 35 to 400 m/z at 70 eV. The mass spectra of the components were identified using data from the NIST Library (NIST08s).

2.4. Response Surface Methodology (RSM)

Response surface methodology was used to study the influence of supercritical operating extraction parameters such as CO₂ pressure (P), CO₂ flow rate ($Q_{{\mathrm{CO}}_2}$), and average particle size of the leaf (d_P), on the extract oil yield. These three response variables were coded as x₁, x₂, and x₃, respectively. The range and levels of independent factors were chosen based upon the results of preliminary tests and are gathered in

Table 2. Conditions of supercritical fluid extraction.

Variable	Symbol	Factor Level
		−1	1
Pressure (bar)	x1	80	220
CO2 flowrate (kg/h)	x2	0.6	1.2
Average particle size (µm)	x3	220	650

. The individual and interactive effects of these parameters on the dependent variable were studied. Equation (1) represents the linear model with interactions for the three operating conditions,

$$Y_D=a_0+a_1{x}_1+a_2{x}_2+a_3{x}_3+a_{12}{x}_1{x}_2+a_{13}{x}_1{x}_3+a_{23}{x}_2{x}_3$$

(1)

where Y_D represents a dependent variable (the yield of extract oil), a₀ is a constant, a₁, a₂, and a₃ are individual linear coefficient, a₁₂, a₁₃, and a₂₃ are the interactive linear coefficient, and x₁, x₂, and x₃ are the coded values of independent factors (pressure, CO₂ flow rate, and average particle size respectively).

Nemrod-w software package was used for the regression analysis of the experimental data obtained [33]. Fit quality of the mathematical model equation was expressed by the determination coefficient R², and its statistical significance was checked by an F-test. The significance of the regression coefficient was tested by a t-test. Significance level was given as *** p < 0.001, ** p < 0.01, * p < 0.05. Differences with p-value superior to 0.05 were not considered significant. For our experimental design validation, optimum conditions were fixed based on the data obtained from the experimental design.

2.5. Modeling of the Supercritical Extraction Process

Stastováet et al. [30] made several simplifications on Sovová’s model [29] by introducing two parameters: the grinding efficiency (G) and the dimensionless time (ψ),

$$\mathsf{\psi }=\frac{tQ{y}_s}{N{x}_0}$$

(2)

where t is the extraction time, Q represents the solvent mass flow rate, $y_s$ is the oil solubility in the solvent, N is the initial mass of the solid, and x₀ is the initial oil concentration in the solid.

The mass of extracted oil (E) can be calculated by the following Equation (3)–(7):

$$\frac{E}{N.x_0}=\{\begin{array}{c}\mathsf{\psi }\left[1-\exp \left(-\mathrm{Z}\right)\right]\mathit{for}\mathsf{\psi }<\frac{G}{Z}\\ \mathsf{\psi }-\frac{\mathrm{G}}{\mathrm{Z}}\exp \left[Z\left(h_k-1\right)\right]\mathit{for}\frac{G}{Z}\le \mathsf{\psi }<{\mathsf{\psi }}_{\mathrm{k}}\\ 1-\frac{1}{\mathrm{Y}}\ln \left\{1+\left[\exp \left(\mathrm{Y}\right)-1\right]\exp \left[\mathrm{Y}\left(\frac{\mathrm{G}}{\mathrm{Z}}-\mathsf{\psi }\right)\right]\left(1-\mathrm{G}\right)\right\}\mathit{for}\mathsf{\psi }\ge {\mathsf{\psi }}_{\mathrm{k}}\end{array}$$

(3)

Inside the interval of dimensionless time, two regions exist between G/Z and ψ_k, where

$${\mathsf{\psi }}_k=\frac{G}{Z}+\frac{1}{Y}\ln \left\{1-G\left[1-\exp \left(Y\right)\right]\right\}$$

(4)

There is a region definition in dimensionless coordinate, namely h_k.

$$h_k=\frac{1}{Y}\ln \left[1+\frac{\left\{\exp \left[Y\left(\mathsf{\psi }-\frac{G}{Z}\right)\right]-1\right\}}{G}\right]$$

(5)

The dimensionless quantities Z and Y are proportional to the mass transfer coefficients according to the first and second extraction period,

$$\mathrm{Z}=\frac{N{k}_f{a}_0{\rho }_f}{Q\left(1-\epsilon \right){\rho }_s}$$

(6)

$$\mathrm{Y}=\frac{N{k}_s{a}_0{x}_0}{Q\left(1-\epsilon \right)y_s}$$

(7)

where $k_f$ and $k_s$ are the external and the internal mass transfer coefficients respectively, a₀ is the particle specific interfacial area, ${\rho }_f$ stands for the solvent density, ${\rho }_s$ represents the solid density, and $\mathsf{\epsilon }$ is the bed void volume.

The model was implemented in MATLAB^TM. Three adjustable parameters were considered: the grinding efficiency (G), the internal mass transfer parameter ($k_s{a}_0$), and the external mass transfer parameter ($k_f{a}_0$). Equilibrium type “A” model was considered according to Sovová [28]. Since the external mass transfer parameter ($k_f{a}_0$) had no sensitivity [40,41], Fiori et al. [39] suggested an approach to determine this parameter referring to the literature correlations—Sherwood (Sh), Reynolds (Re) and Schmidt (Sc) numbers of experimental runs. CO₂ physical properties were evaluated according to NIST database [42] and the oil extract properties were assumed to be related to the major compound found in GC-FID/MS analysis, α-pinene [43,44]. Binary diffusion coefficient (D_AB) between the CO₂ and the major compound was obtained by correlations. In the case of the supercritical fluid extraction, Sherwood number (Equation (8)) is a function of only Reynolds and Schmidt when natural convection is not significant [45,46,47],

$$Sh=c_{0}R{e}^{c_1}S{c}^{c_2}$$

(8)

where, $c_0, c_1$ and $c_2$ are the adjustable parameters. According to the most-proposed correlations in the literature [45,46,47], $c_0$ should be higher than 1, $c_1$ is constrained between 0.5 and 0.8, and $c_2=1/3$.

Catchpole et al. [48] and Lito et al. [49] used the approach to estimate only two adjustable parameters (G, and $k_s{a}_0$) directly using the experimental kinetic curve, and they calculated the other parameter ($k_f{a}_0$) using the Sherwood correlation because they considered that it was not significant. However, in our approach, we estimated not only these two parameters (G and $k_s{a}_0$), but also the external mass transfer parameter ($k_f{a}_0$) using the experimental curve for more accuracy.

The relevance of the model fitting to the experimental data was assessed considering two statistical criteria, namely the coefficient of determination R² detrmined using Equation (9) and the root means square error (RMSE) given by Equation (10),

$$R^2= 1-\frac{{{\displaystyle \sum }}_i^n{\left(y_{\left(i\right)exp}-y_{\left(i\right)model}\right)}^2}{{{\displaystyle \sum }}_i^n{\left(y_{\left(i\right)exp}-\overline{y_{exp}}\right)}^2}$$

(9)

$$RMSE= \sqrt{{\displaystyle \sum }_i^n\frac{{\left(y_{\left(i\right)exp}-y_{\left(i\right)model}\right)}^2}{n}}$$

(10)

where n represents the number of available experimental data, and $y_{exp}$ and $y_{model}$ are the experimental extraction yield and the extraction yield predicted by the model, respectively.

2.6. Cost Estimation of Processes and Scale-Up

The manufacturing cost of the supercritical extract was estimated through methodology proposed elsewhere [50,51]. Concerning material cost, electricity and labor, they were collected from regional information in Tunis (Tunisia 2016). The fixed investment cost was obtained from the literature proposed by Turton et al. [52] to evaluate the cost of manufacturing according to Equation (11) including depreciation (10% of FCI).

$$COM=0.340FCI+2.73C_{OL}+1.23\left(C_{UT}+C_{WT}+C_{RM}\right)$$

(11)

where COM is the cost of manufacturing of supercritical extract of Pistacia, FCI corresponds to the fixed cost of investment, C_OL represents the cost of operational labor, C_UT is the cost of utilities, C_WT the cost of waste treatment, and C_RM is the cost of raw material.

According to Carvalho [53], Equation (12) can be used in order to determine the solvent flow rate required to maintain the same kinetic behavior in different SFE units (scale-up) for a given feed mass and bad geometry. Researchers [54,55,56,57] declared that the extraction time has an influence on an extraction’s COM and the extraction rate increases by increasing the solvent flow rate. They also reported that oil yield in the extraction can be positively influenced by the solvent flow rate increases as the following equation:

$$\frac{Q_{{\mathrm{CO}}_2}^2}{Q_{{\mathrm{CO}}_2}^1}={\left(\frac{F_2}{F_1}\right)}^2\left(\frac{H_{B1}}{H_{B2}}\right)\left(\frac{d_{B1}}{d_{B2}}\right)$$

(12)

The extractor geometry data and the installed supercritical extraction cost were obtained from Núñez and del Valle [58,59,60]. In fact, a plant with two extraction vessels, each one having internal volume varying between 0.2, 0.4, 0.6 and 1.0 m³, was evaluated. For instance, for a plant with vessels of 1.0 m³, the aspect ratio was H/d = 8 (with 0.542 m × 4.334 m of inner diameter and height respectively). The wall thickness withstands 390 bar. The fixed cost of investment (FCI) of each SFE unit was determined in USD based on the values of Rocha-Uribe et al. [59]. All the values are reported in

Table 3. Estimated cost of each supercritical fluid extraction (SFE) unit, including all equipment (Chemical Engineering Plant Cost Index, CPECI, 2014 = 580).

Extractor Vessel (liters)	H (m)	d (m)	H/d	Fixed Cost (FCI) (US $)
100	2.01	0.252	8.0	$853,975
200	2.54	0.317	8.0	$1,378,550
400	3.19	0.399	8.0	$2,225,400
600	3.66	0.457	8.0	$2,944,900
1000	4.34	0.542	8.0	$4,191,250

and are calculated using the Chemical Engineering Plant Cost Index (CPECI) value for 2014 (CPECI 2014 = 580) [61].

According to Experiment (2), the relations of laboratory-scale H/d is 3.6, with the flow rate from 3.36 × 10⁻⁵ kg/s, the bed apparent density ρ = 296 kg/m³, with the operational conditions of 220 bar and 40 °C for the extraction process were taken into consideration.

Figure 2 shows the cycle of solvent (pure CO₂) during the supercritical extraction process in an operating unit. The steps are considered as being primarily solvent collection in reservoir (64 bar and 25 °C), followed by a cooling process (10 °C), pumping and pressurization of the extraction vessel (220 bar and 35 °C), followed by temperature increase (40 °C) until obtaining the desired extraction condition and finally, after this process, reducing the pressure (60 bar and 60 °C) for solute precipitation for reuse.

In this case, it is assumed that during the process of decompression, the solute is separated, and the pure solvent is returned to the system. Based on the findings of this and the mutual values found in other research papers, it is identified that extraction yields (mass extract/Pistacia load) were estimated to reach 0.3%, 0.5%, 0.7%, 1.0% and 1.5%.

The cost of 8000 h per year operational work, with continuous 24 h per day, 8 h daily shifts (2 workers/shift) was considered. The cost of labor was operationally considered to be 6.60 USD/h (1475.60 USD/month tax included). The utility cost was estimated relying on energy consumption involved in the solvent cycle CO₂, cold water, and electricity [51,58]. The specific energies of CO₂ for cooling, heating, and pumping in the solvent cycle were equal to −261.29 kJ/kg, 219.2 kJ/kg and 55.0196 kJ/kg, respectively. These calculations were based upon the work of Rock-Uribe et al. [59].

The electricity cost was equal to 217.10 USD/MWh (price charged in Tunis, Tunisia, tax included). Concerning raw material costs, the considered values were: 1.35 USD/kg of dried and milled Pistacia leaves, 0.15 USD/kg of CO₂ and 0.97 USD/kg of ethanol for cleaning purposes [60]. It was considered that 2% of CO₂ mass was lost during the extraction cycle for all the SFE process scale evaluated in this work. The cost of waste treatment was not considered because CO₂ was fully recycled and Pistacia leaf can be used in soil enrichment or energy generation.

3. Results and Discussion

Table 4. The yield of supercritical extraction, with three different variables (P, dp and $Q_{{\mathrm{CO}}_2}$ ).

Experiment	x1 P [bar]	x2 dP. [µm]	x3 ${\mathit{Q}}_{\mathbf{C}{\mathbf{O}}_2}\left[\mathbf{kg}\mathbf{h}{}^{-1}\right]$	Yield [%]
1	220	650	0.604	0.234
2	220	220	0.602	0.285
3	80	650	1.202	0.093
4	80	650	1.209	0.119
5	80	220	0.602	0.123
6	80	650	0.603	0.117
7	140	650	1.204	0.221
8	180	220	0.913	0.220
9	180	220	0.904	0.229
10	180	650	0.908	0.174

shows the experimental yield results and operating conditions of supercritical extraction for the ten experiments. The experimental oil extract yield was calculated using the following Equation (13):

$$Yield \left(\mathrm{wt}\%\right)=\frac{Mass of leaf extract}{Mass of raw material}\times 100$$

(13)

The yield observed for the tested conditions varied between 0.093% and 0.285%. We found that Experiment (2) with the highest pressure, lowest flow rate, and the lowest average particle diameter gave the greatest tested income extraction conditions. The replications performed at 80 bars (0.10% ± 0.0184%), which were Experiment (3) and (4), demonstrated higher variability than those performed at 180 bar (0.022% ± 0.0057%) (Experiments (8) and (9)); therefore, they have coefficients of variation equal to 17.3% and 2.5% respectively.

By comparing all results, we observe that the present work provided different yields that were in some cases inferior to those reported by other authors [31,60]. In fact, Bampouli et al. [62] obtained an outcome from the leaves of Pistacia lentiscus (PL) var. chia (from Chios, Greece) varying between 1.6% and 5% w/w. The conditions were ranging from 100 to 250 bar and 45 °C with a flow rate of 1.5 to 3.0 kg CO₂/h. Also, Congiu et al. [31] acquired yields between 0.25% and 0.45% for the leaves of Pistacia lentiscus (PL) coming from the regions of Costa Rey and Capoterra (Sardinia, Italy) providing 90 bar and 50 °C with a flow rate of 0.9 kg CO₂/h.

The variation of the obtained yields must be due to the areas of cultivation processing, treatment of raw materials and experimental conditions in the extraction process. Appendix A shows the characterizations of the essential oil obtained from the leaves and carried out using GC-MS-FID. We observed that the major compounds for the extraction of Pistacia lentiscus in the Tunis region are α-pinene (32%), followed by terpinene-4-ol (13%), 1-8-cineole (6%), α-terpineol (4%), β-caryophyllene (4%) and borneol (4%), as summarized in

Table 5. Areas of compounds found in the oils obtained by SFE from leaves.

Compounds	RI	Exp.1	Exp.2	Exp.3	Exp.4	Exp.5	Exp.6	Exp.7	Exp.8	Exp.9	Exp.10
α-pinene	939	33.30	30.00	34.21	30.10	33.21	32.10	34.31	31.10	32.41	32.1
Terpinene-4-ol	1178	13.04	13.24	13.08	13.68	13.02	12.01	12.77	12.06	13.12	12.16
1-8-cineole	1033	5.10	6.10	6.02	6.62	5.85	6.42	5.66	6.11	5.06	7.11
α-terpineol	1189	4.61	4.01	4.58	4.88	4.12	4.21	4.06	4.67	4.68	4.55
β-caryophyllene	1434	4.02	4.92	4.22	4.82	4.88	4.03	4.12	4.43	4.01	4.93
Borneol	1165	3.92	3.12	4.62	4.02	4.22	4.12	4.45	4.16	4.85	4.66
Others		36.01	38.61	31.27	35.88	34.70	37.11	34.63	37.47	35.87	34.49

. Furthermore, as expected, these compositions are not significantly influenced by changing the operating conditions due to the constant extraction temperature.

3.1. Study the Effect of Operating Conditions on the Yield Using RSM

Response surface methodology (RSM) was used to study the individual and the interactive influence of operating extraction parameters on the extract yield to find the optimal operating conditions.

Table 6. Design yield (Y_D) for supercritical extraction of Pistacia lentiscus.

Experiment	x1	x2	x3	YD (%)
1	1.0000	1.0000	−0.9934	0.23
2	1.0000	−1.0000	−1.0000	0.28
3	−1.0000	1.0000	0.9769	0.09
4	−1.0000	1.0000	1.0000	0.12
5	−1.0000	−1.0000	−1.0000	0.12
6	−1.0000	1.0000	−0.9967	0.12
7	−0.1429	1.0000	0.9835	0.22
8	0.4286	−1.0000	0.0247	0.22
9	0.4286	−1.0000	−0.0049	0.23
10	0.4286	1.0000	0.0082	0.17

shows the experimental design yield for the ten experiments.

We used the analysis of variance (ANOVA) to evaluate the statistical significance of the linear model represented in Equation (1). The model can describe the variation of the results because it is significant at <5%. We, also, verified the model efficiency and the adaptability to the experimental data by estimating the coefficient of actual and predicted determination (R² and predicted R² respectively) calculated by the analysis of variance. We found out that the actual determination coefficient indicates that the fitted model explains 91.2% of the variability in the extraction yield. The predicted R² was 0.998 (a good agreement) indicating that our experimental design can be used for modeling the response variables employed, as shown in Figure 3.

Table 7. Coefficient of a linear regression equation (Equation (1)).

Coefficient	Coefficient Value	Test Experiment	p-Value along with Confidence %
a0	0.183	12.43	***
a1	0.084	3.71	*
a2	0.000	0.01	99.1%
a3	−0.002	−0.10	92.2%
a12	−0.013	−0.80	48.5%
a13	0.017	0.71	53.3%
a23	0.019	0.86	45.3%

* p < 0.05; *** p < 0.001.

gathers the statistical results of the constant parameters in Equation (1): the linear intercept constant (a₀), the individual linear effects of the three independent variables (a₁, a₂, and a₃), and their interactive linear effects (a₁₂, a₁₃, and a₂₃). Therefore, the linear regression equation used to evaluate the experimental yield becomes (Y_D),

$${\mathrm{Y}}_{\mathrm{D}}=0.183+0.084x_1$$

(14)

This equation indicates that the main factor that significantly influences the yield was the CO₂ pressure when the confidence of 5% was considered. The equation identifies the best conditions through variation of chosen parameters to maximize the extraction efficiency, which is presented in Experiment (2).

For a better understanding of the statistical results, Figure 4 represents the 2D response surface of the experimental yields by the function of pressure, average particle size and the flow rate of CO₂. As can be observed in Figure 4, the highest yield was obtained around the maximum point of pressure (P = 220 bar) when the flow rate of CO₂ and average particle size were around minimum points.

3.2. Analysis and Validation of Experimental Design

The statistical analysis for the final selected model shows that the effect of the CO₂ pressure is the only variable that has a significant effect on the yield, compared to the other variables that have no effects. For this reason, we analyzed the selection effect of the identification and validation points used for our experimental design. The existence of a correlation between the parameters increases the size of the confidence intervals [63], therefore we need to control the value of the correlation coefficients 2 to 2.

In our investigation, we have used the D-optimality criterion [64] to separate the experiments used for both parametric identification and validation model. This method consists of choosing a set of parametric identification points to obtain the highest determinant of the Fischer matrix information [65]. To investigate the influence of selection on the confidence intervals of each parameter, we studied the correlations between the parameters using the following matrix correlation coefficients summarized in

Table 8. Matrix correlation coefficients.

Parameters	a1	a2	a3	a12	a13	a23
a1	1.0000	−0.3039	−0.3291	0.2608	−0.3623	−0.2981
a2	−0.3039	1.0000	0.3909	−0.0455	−0.2335	−0.2095
a3	−0.3291	0.3909	1.0000	−0.3748	−0.3771	0.4994
a12	0.2608	−0.0455	−0.3748	1.0000	0.2565	−0.2247
a13	−0.3623	−0.2335	−0.3771	0.2565	1.0000	−0.2563
a23	−0.2981	−0.2095	0.4994	−0.2247	−0.2563	1.0000

.

Figure 5, which represents the frequency of the correlation coefficient, mimics 87% of the correlation between the parameters pairs are in the range of 0.2 to 0.4. This indicates the use of the Detmax Fedrov algorithm [66] that has no effect on the correlation between the parameters. Therefore, the used supports of a linear model with interaction, and subsequently the experimental design, are applicable at least in this study.

3.3. Effect of Operating Conditions on the Mass Transfer

All experimental data were used to determine the model parameters (G, and $k_s{a}_0$, and $k_f{a}_0$) (see Appendix B).

Table 9. The adjusted parameters (G and k_sa₀) between two approaches.

Exp.	P [bar]	dP [µm]	${\mathit{Q}}_{\mathbf{C}{\mathbf{O}}_2}$ [kg.h−1]	G			ksa0.105 (s−1)
				1st App.		2nd App.	1st App.		2nd App.
				Lito	Catchpole	AYDI A.	Lito	Catchpole	AYDI A.
1	220	650	0.604	0.36	0.36	0.37	7.48	7.48	7.35
2	220	220	0.602	0.61	0.61	0.61	7.01	7.01	7.01
3	80	650	1.202	0.39	0.39	0.40	0.182	0.182	0.182
4	80	650	1.209	0.36	0.36	0.37	0.107	0.107	0.107
5	80	220	0.602	0.52	0.52	0.52	9.28	9.28	9.30
6	80	650	0.603	0.44	0.44	0.45	7.58	7.58	7.58
7	140	650	1.204	0.35	0.35	0.37	7.83	7.84	7.71
8	180	220	0.913	0.62	0.62	0.66	2.75	2.75	2.31
9	180	220	0.904	0.54	0.54	0.59	2.32	2.32	1.91
10	180	650	0.908	0.23	0.23	0.23	0.108	0.108	0.108

shows that the two adjusted parameters (grinding efficiency (G) and internal mass transfer parameter ($k_s{a}_0$), which is estimated by the experimental kinetic curves, have no significant change between the two approaches. As expected [3,13,31], the value of grinding efficiency (G) increases by decreasing the average particle size. This parameter is not only related to the particle size, but also to the shape of its distribution (normal, bimodal). As a result, it influences the curve shape due to the solvent flow asymmetry effect.

In the first approach, we applied the correlations proposed by Lito, Catchpole, and King using only the binary diffusion coefficient (CO₂-α-pinene) to estimate the external mass transfer parameter ($k_f{a}_0$). On the other hand, the adjusted parameters were not only grinding efficiency (G) and the internal mass transfer parameter ($k_s{a}_0$) but also the external mass transfer parameter, which was estimated by experimental extraction curves in the second approach.

Table 10. Parameters (k_fa₀, and D_AB) between two approaches.

Exp.	P [bar]	dP [µm]	${\mathit{Q}}_{\mathbf{C}{\mathbf{O}}_2}$ [kg.h−1]	kfa0.103 (s−1)			DAB. 109(m2 s−1)
				1st App.		2nd App.	1st App.		2nd App.
				Lito	Catchpole	AYDI A.	Lito	Catchpole	AYDI A.
1	220	650	0.604	1.41	1.46	0.736	8.63	9.06	2.97
2	220	220	0.602	6.56	6.77	1.72	8.63	9.06	1.06
3	80	650	1.202	0.114	0.118	8.70	38.0	40.2	0.73
4	80	650	1.209	0.114	0.119	8.95	38.0	40.2	0.76
5	80	220	0.602	0.355	0.369	3.42	38.0	40.2	1.04
6	80	650	0.603	7.63	7.93	4.03	38.0	40.2	13.3
7	140	650	1.204	2.67	2.77	1.26	10.7	11.3	3.15
8	180	220	0.913	9.19	9.51	0.773	9.4	9.9	6.61
9	180	220	0.904	9.14	9.46	0.929	9.4	9.9	8.79
10	180	650	0.908	1.97	2.04	2.64	9.4	9.9	13.4

shows the parameters ($k_f{a}_0$, and D_AB) for both approaches. The adjusted parameter ($k_f{a}_0$) estimated by the experimental kinetic curves in the second approach is different than the correlated parameter proposed by Lito and Catchpole in the first approach. This parameter affects the values of the diffusion coefficient as shown in Table 10.

In the second approach, the parameter values of the Sherwood number were experimentally obtained from the extraction curves and k_f values. The parameter correlation suggested in the second approach was determined from c₀ and c₁ settings, resulting in the equation as shown in

Table 11. Number and the coefficient of determination r² between two approaches.

Exp.	P [bar]	dP [µm]	${\mathit{Q}}_{\mathbf{C}{\mathbf{O}}_2}$ [kg.h−1]	Sh			r2
				1st App.		2nd App.	1st App.		2nd App.
				Lito	Catchpole	AYDI A.	Lito	Catchpole	AYDI A.
1	220	650	0.604	0.24	0.24	0.96	98.478	98.477	98.472
2	220	220	0.602	0.13	0.13	4.16	99.352	99.352	99.355
3	80	650	1.202	0.45	0.44	0.06	99.716	99.716	99.718
4	80	650	1.209	0.45	0.44	0.05	99.632	99.632	99.626
5	80	220	0.602	0.16	0.16	1.24	99.297	99.297	99.357
6	80	650	0.603	0.30	0.30	0.13	98.302	98.301	98.322
7	140	650	1.204	0.37	0.37	0.40	98.548	98.541	98.663
8	180	220	0.913	0.17	0.16	5.84	98.268	98.268	99.137
9	180	220	0.904	0.17	0.16	5.11	98.421	98.421	99.019
10	180	650	0.908	0.31	0.31	0.29	99.381	99.381	99.381

.

$$Sh=0.0349R{e}^{0.58}S{c}^{1/3}$$

(15)

The adjusted parameters are within the range reported by [19] and their applicability for Reynolds and Schmidt values are in the ranges 2 ≤ Re ≤ 60 and 2 ≤ Sc ≤ 12.

The values of the mass transfer coefficients ($k_f{a}_0$) ranged from 1.4 × 10⁻² to 3.7 × 10⁻¹ for the estimation performed with the correlations of Lito and Catchpole.

Concerning the adjustment made with the three parameters G, $k_f{a}_0$ and $k_s{a}_0$ (the value of $k_f{a}_0$ is 0.020), it ranged from 7.3 × 10⁻³ s⁻¹ to 8.9 × 10⁻² s⁻¹. These values are lower than those obtained by adjustment with the correlations. The parameters $k_s{a}_0$ were in all cases between 1.9 × 10⁻⁵ s⁻¹ and 1.8 × 10⁻⁴ s⁻¹.

3.4. Cost Estimation of Processes and Scale-Up

The greatest impact on the cost of manufacturing extracted oil production of P. lentiscus in Tunisia is represented by the raw material cost (RMC), followed by the fixed cost of investment and the utility cost, as indicated in Figure 6.

Table 12. Manufacturing cost of a supercritical extract of Pistacia lentiscus leaves, in a batch of 60 min (220 bar and 40 °C).

Volume (m3)	COM (US$/kg)
	Expected Oil Yield (%)
	0.3	0.5	0.7	1.0	1.5
0.1	999.63	599.78	428.41	299.89	199.93
0.2	942.63	565.58	403.99	282.79	188.53
0.4	952.04	571.22	408.02	285.61	190.41
0.6	945.12	567.07	405.05	283.53	189.02
1.0	813.95	488.37	348.83	244.18	162.79

shows the manufacturing costs of the supercritical extract of Pistacia lentiscus in the US. We note that for the same yield, the lowest costs were those obtained at the highest production volume, as was expected.

Concerning yield values obtained on a pilot scale (0.30%), the cost of the manufacturing process of the supercritical extract weas between 999.63 USD/kg and 813.95 USD/kg. These values are considered very high when compared to the costs of other vegetable raw materials such as rosemary extract which is worth 49.71 USD/kg [21], ginger oleoresin which costs 99.80 USD/kg [20], Curcuma longa L. extract which is worth 164.4 USD/kg [40] and habanero pepper extract with a cost of 540.19 USD/kg [33].

A yield increase in the extraction process reduces the cost of manufacturing. Thus, the high cost of production is due to the low yields gained from this kind of plant. Prices of traded extracted oils of Pistacia lentiscus sold in 5 mL, 10 mL or 30 mL vials are about 5.83 USD/g (5830.00 USD/kg). These values are obtained from the local market in Tunis, Tunisia.

4. Conclusions

Pistacia lentiscus L. plant from the Tunisian region is appeased of medicinal properties in its extract oil that can be produced using supercritical carbon dioxide (SC-CO₂) extraction. In this study, we observed that the α-pinene (32%) was the major compound of the extracted oil of Pistacia lentiscus in the Tunis region. The experiment (2) having the highest pressure, lowest flow rate, and the lowest average particle diameter gave the greatest tested income extraction conditions.

We investigated the influence of CO₂ pressure, average particle size, and CO₂ flow rate and their interaction on the extract yield using the response surface methodology (RSM). It was observed that the main factor that significantly influences the yield was the CO₂ pressure and, therefore, the best conditions through variation of chosen parameters to maximize the extraction efficiency were presented in Experiment (2).

We studied the influence of operating parameters on mass transfer by evaluating a process applying broken and intact cell (BIC) on the essential oil extraction curves that are acquired from the leaves of Pistacia lentiscus L. The two adjusted parameters (grinding efficiency (G) and internal mass transfer parameter ($k_s{a}_0$)), which were estimated by the experimental kinetic curves, have no significant change between the two approaches (Lito and Catchpole approach, and AYDI A approach). However, the external mass transfer parameter ($k_f{a}_0$) proposed by AYDI A was different from the correlated parameter in the first approach that significantly influences the values of the diffusion coefficient.

The economic evaluation in the scale-up process was obtained for the SC-CO₂ extraction of these plant leaves. We indicated that the lowest costs were obtained at the highest production volume for the same yield. The manufacturing cost of oil production is reduced by a yield increase in the extraction process because of the low yields obtained from this type of plant.