The amount of Thymol and carvacrol (mg/g dried sample) of the essential oil of *Z.multiflora*, extracted by SWE, hydrodistillation and Soxhlet extraction [7].

## 1. Introduction

Extraction always involves a chemical mass transfer from one phase to another. The principles of extraction are used to advantage in everyday life, for example in making juices, coffee and others. To reduce the use of organic solvent and improve the extraction methods of constituents of plant materials, new methods such as microwave assisted extraction (MAE), supercritical fluid extraction (SFE), accelerated solvent extraction (ASE) or pressurized liquid extraction (PLE) and subcritical water extraction (SWE), also called superheated water extraction or pressurized hot water extraction (PHWE), have been introduced [1-3].

SWE is a new and powerful technique at temperatures between 100 and 374^{o}C and pressure high enough to maintain the liquid state (Fig.1) [4]. Unique properties of water are namely its disproportionately high boiling point for its mass, a high dielectric constant and high polarity [4]. As the temperature rises, there is a marked and systematic decrease in permittivity, an increase in the diffusion rate and a decrease in the viscosity and surface tension. In consequence, more polar target materials with high solubility in water at ambient conditions are extracted most efficiently at lower temperatures, whereas moderately polar and non-polar targets require a less polar medium induced by elevated temperature [5].

Based on the research works published in the recent years, it has been shown that the SWE is cleaner, faster and cheaper than the conventional extraction methods. The essential oil of *Z. multiflora* was extracted by SWE and compared with two conventional methods, including hydrodistillation and Soxhlet extraction [7]. The total extraction yields found for the total essential oil of *Z. multiflora* were 2.58, 1.51 and 2.21% (w/w) based on the dry weight for SWE, hydrodistillation and Soxhlet extraction, respectively.

The comparison among the amount of thymol and carvacrol (milligram per gram dried sample) by SWE, hydrodistillation and Soxhlet extraction is shown in Table 1 [7]. The amount of valuable oxygenated components in the SWE method is significantly higher than hydrodistillation and Soxhlet extraction. As hexane is a nonpolar solvent, non-oxygenated components are enhanced compared to subcritical water. On the other hand, in general, non-oxygenated components present lower vapor pressures compared to oxygenated components, and in this sense, its content in hydrodistillated extracts are increased. Because of the significant presence of the oxygenated components, the final extract using the SWE method was relatively better and more valuable.

## 2. Equipment of subcritical water extraction

No commercial SWE equipment is available, but the apparatus is easy to construct in the laboratory. SWE is performed in batch or continue systems but continue system is current. In this system, extraction bed is fixed and flow direction is usually up to down for easily cleaning of analytes. Already, we worked SWE of *Z. multiflora* with first generated system (Iranian Research Organization for Science and Technology [IROST], Tehran, Iran) [7] (Fig. 2(a)). After it, laboratory-built apparatus of SWE equipment was designed by studying of different systems and have advantages versus first generation (Semnan University, Semnan, Iran) (Fig. 2(b)).

The second generated system is presented in Fig. 2(c). The main parts of a dynamic SWE unit are the following: three tanks, two pumps, extraction vessel, oven for the heating of the extraction vessel, heat Exchanger for cooling of exteract, pressure restrictor and sample collection system. One of the pumps is employed for pumping the water (and extract) and another pump is employed for flushing the tubings. Also, one of tank is employed for organic solvent that is as a main solvent or co-solvent. A pressure restrictor is needed to maintain the appropriate pressure in the equipment. It was constructed of stainless steel.

## 3. Effective parameters in subcritical water extraction

### 3.1. Effect of temperature

One of the most important parameters affecting SWE efficiencies is the extraction temperature. As the temperature rises, there is a marked and systematic decrease in permittivity, an increase in the diffusion rate and a decrease in the viscosity and surface tension. SWE must be carried out at the highest permitted temperature. It should be mentioned that increasing the extraction temperature above a certain value gives rise to the degradation of the essential oil components. The maximum permitted extraction temperature must be obtained experimentally for different plant materials. Regarding the extraction of essential oils, it has been shown that temperatures between 125 and 175^{o}C will be the best condition. The extraction temperature for *Z. multiflora* was optimized in order to maximize the efficiency of thymol and carvacrol (structural isomers) as key components (more than 72%) [7]. Its influence was studied between 100 and 175^{o}C, and the mean particle size, flow rate, extraction time and pressure were selected to be 0.5 mm, 2 ml/min, 60 min and 20 bar pressure, respectively (Fig. 3).

It was seen that the efficiency of thymol and carvacrol increased generally with increase in temperature up to 150^{o}C. At 175^{o}C, it decreased, and an extract with a burning smell was produced. It may be the result of degradation of some of the constituents at higher temperatures. Because of the highest efficiency of thymol and carvacrol essential oil at 150^{o}C and the disagreeable odor of the extract at higher temperatures, further experiments were carried out at this temperature.

### 3.2. Effect of particle size

The effect of the mean particle size on the efficiency of thymol and carvacrol at 150^{o}C temperature, 2ml/min flow rate, 20 bar pressure and 150 min extraction time is shown in Fig. 4 [7]. The mean ground leaf particles were selected to be 0.25, 0.5 and 1.0 mm. The final amount of thymol and carvacrol extracted from 0.5-mm-size particles was near to 0.25-mm particles. It shows that, at least in the selected range of mean particle sizes (0.25-0.5 mm), the extraction process may not be controlled by the mass transfer of thymol and carvacrol. It was expected that the rate of the 0.25-mm-size particles was more than the 0.5-mm-size particles, but it did not happen.

A possible explanation for this observation could be that the particles were close fitting at initial times, and the extraction was done slowly. After the expired time, the close fitting particles opened from each other, and at final extraction value of the 0.25-mm-size particles was more than the 0.5-mm-size particles. Regarding the larger 1.0-mm-size particles, the efficiency is substantially lower. It shows that the process may be controlled by the mass transfer of thymol and carvacrol for larger particle sizes.

To prevent the probable vaporization of the essential oils during the grinding of the leaves and also to make the work of the filters easier, for further experiments, the best value of mean particle size was selected as 0.50 mm.

### 3.3. Effect of flow rate

The effect of water flow rate on the efficiency of thymol and carvacrol at 150^{o}C temperature, 0.5-mm-particle size, 20 bar pressure and 150 min extraction time is shown in Fig. 5 [7]. The water flow rate was studied at 1, 2 and 4 ml/min. As can be seen, the rate of the essential oil extraction was faster at the higher flow rate. The rate is slower at 2 ml/min and even slower at 1 ml/min. It is in accordance with previous works [8, 9]. It means that the mass transfer of the thymol and carvacrol components from the surface of the solid phase into the water phase regulated most of the extraction process. Increase in flow rate resulted in increase in superficial velocity, and thus, quicker mass transfer [10]. The main disadvantage of applying higher water flow rates is increasing the extract volume and consequently, lower concentration of the final extracts. In practice, the best flow rate must be selected considering two important factors, including the extraction time and the final extract concentration. It is clear that a shorter extraction time and more concentrated extracts are desirable. To prevent a slower extraction rate and longer extraction times, despite the larger amount of the final extracts, a flow rate of 2 ml/min was selected as the optimum value.

While temperature, particle size and flow rate extraction are the main parameters affecting SWE, type of analyte, extraction vessel characteristics and use of modifiers and additives are also important. Although matrix and other effects play a role, many of these are less critical in SWE than in SFE because of the harsh extraction conditions (high temperature) typical in SWE, particularly for non-polar analytes.

## 4. Extraction mechanism

The SWE process can be proposed to have six sequential steps: (1) rapid ﬂuid entry; (2) desorption of solutes from matrix active sites; (3) diffusion of solutes through organic materials; (4) diffusion of solutes through static ﬂuid in porous materials; (5) diffusion of solutes through layer of stagnant ﬂuid outside particles; and (6) elution of solutes by the ﬂowing bulk of ﬂuid (Fig. 6).

As we know, the extraction rate is limited by the slowest of these three steps. The effect of step (1) is typically small and often neglected. Although the diffusion of the dissolved solute within the solid is usually the rate limiting step for most botanicals, partitioning of solute between the solid matrix and solvent have been reported as the rate-limiting mechanism for SWE of essential oil from savory [10].

The plots the amount of compound extracted versus solvent flow rates and versus solvent volume can determine the relative importance of these steps. For example, if the rate of extraction is controlled by intra-particle diffusion or kinetic desorption, the increase in bulk fluid flow rate would have little effect on extraction rate. On the other hand, if the extraction is controlled by external film transfer diffusion, extraction rates increase with solvent flow rate. In the case where the extraction rate is controlled by thermodynamic partitioning, doubling the bulk fluid flow rate would double the extraction rate, while the curves of extraction efficiency versus the volume of water passed for all flow rates would overlap. In one of our previous work, four proposed models have been applied to describe the extraction mechanisms obtained with SWE of *Z. multiflora* essential oil. These were included (1) partitioning coefficient model, (2) one-site (3) two-site desorption models and (4) thermodynamic partition with external mass transfer model [11]. In other studying unsteady state mass balance of the solute in solid and subcritical water phases (two-phase model) was investigated [12]. Also Computational Fluid Dynamics (CFD) modeling of extraction was considered [13].

## 5. Modeling of SWE

### 5.1. Thermodynamic model (Partitioning coefficient (K_{D}) model)

Partitioning coefficient model, adopted from Kubatova et al. [10], describes the extraction process that is controlled by partitioning of solute between matrix and solvent similar to elution of solute from a partition chromatography column. For extraction, this type of behavior occurs when the initial solute concentration in the plant matrix is small. This model assumes that the initial desorption step and the subsequent fluid-matrix partitioning is rapid. Here the thermodynamics partitioning coefficient, *K*_{D}, is defined as:

Hence the Extraction with subcritical water can be fitted using this simple thermodynamic model. The mass of analyte in each unit mass of extraction fluid and the mass of analyte remaining in the matrix at that period in the entire extraction time is based on the *K*_{D} value determined for each compound. The thermodynamic elution of analytes from matrix was the prevailing mechanism in SWE as evidenced by the fact that extraction rate increased proportionally with the subcritical water flow rate. Therefore, if the K_{D} model applies to a certain extraction, the shape of an extraction curve would be defined by:

M_{a}: cumulative mass of the analyte extracted after certain amount of volume V_{a} (mg/g dry sample)

M_{b}: cumulative mass of the analyte extracted after certain amount of volume V_{b} (mg/g dry sample)

M_{i}: total initial mass of analyte in the matrix (mg/g dry sample)

M_{b}/M_{i} and M_{a}/M_{i}: cumulative fraction of the analyte extracted by the fluid of the volume V_{b} and V_{a} (ml)

K_{D}: distribution coefficient; concentration in matrix/concentration in fluid

ρ: density of extraction fluid at given condition (mg/ml)

e: exponential function

m: mass of the extracted sample (mg dry sample).

The model eq. (1) and the experimental data for *Z. multiflora* from all volumetric flow rate, were used to determine the K_{D} value by minimizing the errors between the measured data and the K_{D} model using Matlab curve fitting solver. The values of K_{D} are shown in Table 2 for different flow rates [11]. It was demonstrated that individual essential oil compounds have a range of K_{D} values from ~4 to ~250 [10].

### 5.2. Mass transfer models

#### 5.2.1. Diffusion model

Mass transfer can be defined as the migration of a substance through a mixture under the influence of a concentration gradient in order to reach chemical equilibrium. The diffusion coefficient (D_{e}) is the main parameter in Fick’s law, and application of this mathematical model to solid foods during solid-liquid extraction is a common way to calculate the effective diffusion coefficient (Crank, 1975 [14]). However, Gekas (1992) noted, values of D_{e} can vary by several orders of magnitude for the same material which may be due to structural changes in the food material during different stages of the process [15]. Therefore, it is important to keep a constant particle size as breakage of cell wall or grinding can reduce the particle size and hence decrease the distance for solute to travel from inside to surface of particle.

Fick derived a general conservation equation for one-dimensional non-steady state diffusion when the concentration within the diffusion volume changes with respect to time, known as Fick’s second law (Cussler, 1984; Mantell et al., 2002) [16, 17]:

Where C is the solute concentration (mg/ml) at any location in the particle at time t (s); C_{o} is the initial solute concentration (mg/ml); D_{e} is the effective diffusion coefficient (m^{2}/s) assuming that D_{e} is constant with the concentration; t is extraction time (s); r is the radial distance from the centre of a spherical particle (m); R is radius of spherical particle (m).

Various solutions of Fick’s second law have been presented for the diffusion of a compound during solid-liquid extraction depending on the shape of the particle. An approximate numerical solution to Fick’s second law (eq. 2) for a spherical particle was given by Crank (1975) and Cussler (1984):

Where M_{t}: total amount of solute (mg/g) removed from particle after time t, M_{∞}: maximum amount (mg/g) of solute extracted after infinite time. M_{t}/M_{∞}: ratio of total migration to the maximum migration concentration, R: average radius of an extractable particle. When time becomes large, the limiting form of Eq. 3 becomes:

To determine the effective diffusion coefficient values two methods were used. The first method was a linear (graphical) solution in which D_{e} was determined from the slope of the ln (1-M_{t}/M_{∞}) vs. time plot (Dibert et al., 1989) [18]. Thus, eq. 4 can be solved by taking the natural logarithm of both sides. It shows that the time to reach a given solute content will be directly proportional to the square of the particle radius and inversely proportional to D_{e}:

Where slope =

The second method of solution used involved nonlinear regression with effective diffusivity (D_{e}) as a fitting parameter. In this method, the effective diffusivity D_{e} was estimated from eq. 4 using a Microsoft Excel Solver program. The program minimizes the mean square of deviations between the experimental and predicted ln(1-M_{t}/M_{∞}) values (Tutuncu and Labuza, 1996) [19]. The first 10 terms of the series solution are taken into consideration by the program as the solution to the series becomes stable after 10 terms (n=10).

Most researchers in this area have adopted diffusion models based on solutions to Fick’s second law for various defined geometrical shapes of the solid given in Crank method, usually based on infinite or semi-infinite geometries of a slab (plane sheet), a cylinder or a sphere. For example, a slab can be used to describe an apple slice or a sheet of herring muscle; a sphere for the description of coffee beans or particles of cheese curd, and a cylinder for the description of cucumber pickles. There is considerable variation in the solutions adopted by various researchers. The starting point for modeling a particular diffusion process is to consider the shape of the solid and the nature of the process itself: uptake of solute into the food, leaching of solute from the food or diffusion of solute through the food, and the experimental conditions in terms of initial and equilibrium solute concentrations.

The solution models usually consider a uniform initial solute concentration throughout the food, no resistance to mass transfer in the diffusion medium and no chemical reaction; but vary for a particular geometry depending on the solute concentrations at the surface of the solid, the volume of the solution (and therefore the relative change in solute concentration in the external solvent), and the time period of the experiment. A representative selection of solutions is given in Table 3. For example, Bressan et al. investigated solute diffusional loss from coffee beans and cottage cheese curd, respectively [20]; they chose solution models to Fick’s second law from Crank which assumed the geometry of the solids approximate to that of infinite spheres. However, the former researchers considered a system in which the solute concentration in the external solvent remains effectively constant and zero throughout the extraction. This can be the case for large infinite solvent volumes and small solute solid concentrations. The latter group of researchers considered a process where the concentration of the solute in the surrounding medium changes significantly. This can be the case for small solvent volumes and high solid phase solute concentrations in the case of leaching, and high external solvent concentrations in the case of solute uptake processes; and consequently adopted a different model.

#### 5.2.2. One-site kinetic desorption model

One-site kinetic desorption model describes the extractions that are controlled by intra-particle diffusion. This occurs when the flow of fluid is fast enough for the concentration of a particular solute to be well below its thermodynamically controlled limit. The one-site kinetic model was derived based on the mass transfer model that is analogous to the hot ball heat transfer model [29, 30]. The assumptions are that the compound is initially uniformly distributed within the matrix and that, as soon as extraction begins, the concentration of compound at the matrix surfaces is zero (corresponding to no solubility limitation). For a spherical matrix of uniform size, the solution for the ratio of the mass, M_{r}, of the compound that remains in the matrix sphere after extraction time, t, to that of the initial mass of extractable compound, M_{i} is given as:

In which n is an integer and D_{e} is the effective diffusion coefficient of the compound in the material of the sphere (m^{2}/s).

Diffusion process | Diffusion equation | Experimental measurements | Calculation Diffusivity |

proteins Diffusion though a potato disk [21] | Protein concentration on the source side initially; at intervals on receiving side | Fitting the experimental values for Mt to the equation by non-linear regression | |

Analytes diffusion though an apple disk [22] | A simple lumped parameter equation model | soluble analytes content in the limited volume of solvent | A single effective diffusivity was calculated for each set of data directly from the formula |

Analytes diffusion though cheese curd [19] | Solution of Fick’s 2nd law for solute loss or uptake for sphere geometry but for finite volume of solvent with attainment of equilibrium | Samples of solvent surrounding curd withdrawn at intervals | Experimental data fitted to the equation model |

Analytes diffusion though carrot cylinders [23] | Solution of Fick’s 2nd law for solute loss or uptake for an infinite cylinder for short time periods. Authors then modify values to apply to finite cylinder | Carrot samples withdrawn for analysis at intervals | Dimensionless time values calculated for each data point using equation, linear regression of Fourier relationship yields De |

Analytes diffusion though potato tissue [24] | Solution of Fick’s 2nd law similar to flow through a membrane | Samples of solvent surrounding slices withdrawn and assayed at intervals, material balance used to calculate Mi and Mt | Equation model modified by omitting terms n"/>0 and natural logarithms, by non-linear regression of ln Mi/Mt against t |

Analytes and pectic substances diffusion though apple tissue [25] | Solution of Fick’s 2nd law for diffusion from an infinite slab developed by Newman [59] | Apple slices withdrawn for analysis at intervals | A series of dimensionless time values found for data using Newman’s tables. De calculated by linear regression of Fourier relationship |

Analytes diffusion though apple tissue [26] | Equation model adopted considers diffusion in a slab in three planes. | Apple slices withdrawn for analysis at intervals | Modification of equation model by natural logarithms enables regression of lnE against t give slope with De term |

Salt and acetic acid diffusion though herring [27] | Solution of Fick’s 2nd law for diffusion from/or uptake by infinite solution. | Fish withdrawn for analysis at intervals | Experimental data fitted to equation model by successive approximations |

Analytes diffusion though potato tissue [28] | Solution of Fick’s 2nd law for diffusion from/or uptake by an infinite slab for uptake from a solution of finite volume. Solution of Fick’s 2nd law for diffusion from/or uptake by infinite slab. | Initial solute content of potato strip and solute contents of solvent surrounding strip at intervals | Minimizing the residuals between experimental data and theoretical values for Mt/M∞ obtained from the appropriate equation models |

The curve for the above solution tends to become linear at longer times (generally after t > 0.5 t_{c}), and ln (M_{ r}/M_{i}) is given approximately by:

Where t_{c} (min) is a characteristic time quantity, defined as:

An alternative form of eq. 7, or so called a one-site kinetic desorption model, can be written for the ratio of mass of analyte removed after time *t* to the initial mass, M_{i}, as given by:

In which M_{t} is the mass of the analyte removed by the extraction fluid after time t (mg/g dry sample), M_{i} is the total initial mass of analyte in the matrix (mg/g dry sample), M_{t}/M_{i} is the fraction of the solute extracted after time t, and k is a first order rate constant describing the extraction (min^{-1}).

Matlab curve fitting solver was used to determine the desorption rate constant, k, from the data for all flow rates. The values for *Z. multiflora* SWE are show in Table 4 [11]. As mentioned, the kinetic desorption model does not include a factor describing extraction flow rate, k should be the same value for all flow rates if the model is said to fit the experimental data. However, this was not the case (Table 4, the average error 3%-17%). The kinetic desorption rate increased for the volumetric flow rate of 1 to 4 ml/min. This indicated that the kinetic desorption model may not be suitable for describing the data at different flow rates of *Z. multiflora*.

#### 5.2.3. Two-site kinetic desorption model

Two-site kinetic model is a simple modification of the one-site kinetic desorption model that describes extraction which occurs from the “fast” and “slow” part [10]. In such case, a certain fraction (F) of the analyte desorbs at a fast rate defined by k_{1}, and the remaining fraction (1-F) desorbs at a slower rate defined by k_{2}. The model has the following form:

The two site kinetic model does not include solvent volume, but relies solely on extraction time. Therefore, doubling the extractant flow rate should have little effect on the extraction efficiency when plotted as a function of time. On the contrary, the thermodynamic model is only dependent on the volume of extractant used. Therefore, the extraction rate can be varied by changing the flow rate. Hence, the mechanism of thermodynamic elution and diffusion kinetics can be compared simply by changing the flow rate in SWE. If the concentration of bioactive compounds in the extract increases proportionally with an increase in flow rate at given extraction time when the solute concentration is plotted versus extraction time, the extraction mechanism can be explained by the thermodynamic model. However, if an increase in flow rate has no significant effect on the extraction of the bioactive compounds, with the other extraction parameters being kept constant, the extraction mechanism can be modeled by the two site kinetic model [10, 31]. The mechanism of control and hence the model valid for SWE may be different depending on the raw material, the target analyte and extraction conditions.

For the two-site kinetic desorption model, the values of k_{1} and k_{2} were determined by fitting the experimental data with the two-site kinetic desorption models by minimizing the errors between the data and the model results. In the two-site model, the extraction rate should not be dependent on the flow rate. The k_{1} and k_{2} values for *Z. multiflora* SWE shown in Tables 5 and 6 demonstrated that the extraction rates were not completely independent of flow rate (the average error 11%-20%).

#### 5.2.4. Thermodynamic partition with external mass transfer resistance model

This model describes extraction which is controlled by external mass transfer whose rate is described by resistance type model of the following form:

in which C is the fluid phase concentration (mol/m^{3}), C_{s} is the solid phase concentration (mol/m^{3}), k_{e} is the external mass transfer coefficient (m/min) and a_{p} is specific surface area of particles (m^{2}/m^{3}) [32]. If the concentration of the solute in the bulk fluid is assumed small and the solute concentration in the liquid at the surface of solid matrix is described by partitioning equilibrium, K_{D}, the solution of eq. 11 for the solute concentration in the solid matrix, C_{s}, becomes:

eq. 11 can be rewritten as the ratio of the mass of diffusing solute leaving the sample to the initial mass of solute in the sample, M_{t}/M_{i}, as given by the following equation.

Because a_{p} is difficult to be measured accurately, a_{p} and k_{e} are usually determined together as k_{e}a_{p}, which is called overall volumetric mass transfer coefficient. The factors that influence the value of k_{e}a_{p} include the water flow rate through the extractor and the size and shape of plant sample.

The values for the model parameters, K_{D} and k_{e}a_{p} in eq. (9) determined by Matlab curve fitting solver from the experimental data obtained at 150°C are summarized for *Z. multiflora* SWE in Tables 7 and 8 for different mass flow rates (Q, mg min^{-1} ) [11]. Linear regresion of the plot between ln(k_{e}a_{p}) and lnQ gives the following correlation for k_{e}a_{p} and Q:

To quantitatively compare the extraction models, the mean percentage errors between the experimental data and the models were considered. Based on the result in fitting from experimental data, the K_{D} model was generally suitable for the description of extraction over all the volumetric flow rates tested. On the other hand, one-site and two-site kinetic desorption models describe the extraction data reasonably at lower volumetric flow rates. Of all the models considered, however, the thermodynamic partition with external mass transfer model could best describe the experimental data.

#### 5.2.5. Two-phase model

A mathematical model can be developed to predict optimal operating parameters for SWE in a packed-bed extractor. Three important steps consist of diffusion of solutes through particles, diffusion and convection of solutes through layer of stagnant fluid outside particles and elution of solutes by the flowing bulk of fluid are assumed. Unsteady state mass balance of the solute in solid and subcritical water phases led to two partial differential equations. The model can be solved numerically using a linear equilibrium relationship. The model parameters were mass transfer coefficient, axial dispersion coefficient, and intraparticle diffusivity. The last parameter was selected to be the model tuning parameter. The two other parameters were predicted applying existing experimental correlations.

#### 5.2.5.1. Model description

The more precise method is based on differential mass balances along the extraction bed. A two-phase model comprising solid and subcritical water phases can be used. Extraction vessel is considered to be a cylinder filled by mono-sized spherical solid particles. The overall scheme of system was like Fig. 7.

The major assumptions used to describe the SWE process for deriving the essential oils extraction model were:

Packed bed extractor was isothermal and isobaric,

The physical properties of subcritical water were constant,

The hydrodynamics of a fluid bed was described by the dispersed plug-flow model,

The radial concentration gradient in the bulk fluid phase was assumed to be negligible,

The volume fraction of bed was not influenced by the weight loss of plant during the extraction,

The essential oil was assumed as a single component.

Under these assumptions, the differential mass balance equation for any component in the particle and bulk liquid phase and associated initial and boundary conditions can be written as following dimensionless forms:

Solid phase:

Subcritical water phase

Eqs. (17) and (21) may be solved using a linear equilibrium relationship between concentrations in the solid phase and SW phase at the interface, as follows [33]:

Where C_{fp} is the solute concentration in the fluid phase at the particle surface, _{p} is the volumetric partition coefficient of the solute between the solid and the fluid phase. Therefore, there are three eqs. (17), (21) and (25) which can be solved, simultaneously, for three unknowns C_{p}, C_{f} and C_{fp}.

The finite difference equations are a set of simultaneous linear algebraic equations must be solved for implicit method to obtain the concentration distribution at any time. In both cases, tridiagonal systems arise which are conveniently solved at each time step by the Thomas algorithm [34]. A computer code was written using MATLAB simulation software.

#### 5.2.5.2. Parameter identification and correlations

The possible control of mass transfer was assayed by estimating the diffusion coefficient in the liquid. We can use the correlation proposed by Wilke and Chang (1955) to estimate this coefficient [35]:

Where φ is 2.26 for water and 1.5 for ethanol and V_{A} was estimated by the Tyn Calus equation:

Where V_{c} is the estimated by the method of Joback and Reid (1987), with the aid of Molecular Modeling Plus software (Norgwyn Montgomery Software, USA) [36]. The axial dispersion coefficient in the supercritical phase was approximated as follows [37] and it may be used for the subcritical phase:

Where the average void volume fraction of the fixed bed was ε = 0.4 and:

Intraparticle diffusivity of the essential oils in the solid phase, D_{m}, was selected to be the tuning parameter of the model. Mass transfer between liquids and beds of spheres, k_{f}, can be represented by Wilson and Geankoplis in two cases [38]:

Density of water at high pressures and temperatures from 273 to 473 K was assumed to be the density of saturated water (kg). It was calculated as follows [39]:

Where ρ is density (kg/m^{3}) and T is temperature (K). The water viscosity at temperatures from 300 to 450 K was calculated by the following equation:

Where μ is the viscosity (Pa.s) and T is temperature (K). Viscosity of water was supposed to be independent of pressure. The model was verified successfully using the SWE data for *Z.multiflora* leaves at 20 bar and 150^{o}C. The optimum value of 2 × 10^{-12} m^{2}/s was obtained for the intraparticle diffusivity [12].

#### 5.2.6. CFD model

There is no scientific literature about the application of CFD modeling approach in the SWE processes. In our previous work, we have tried to do CFD modeling of essential oils of *Z. multiflora* leaves [13].

In CFD modeling in packed bed reactor, there are two cases. In first one, when the reactor to particle diameter ratio is less than10, it is needed to build the exact geometry and the location of the particles and solving the governing equations in meshes and can see the field of flow between particles inside reactor. But in another one, when the reactor to particle diameter ratio is higher than10, this system is defined as porous media. For modeling of these reactors, the Navier-stocks with one additional term, which is contained viscous force and inertia loss, in porous media is solved.

In this work, because of the reactor to particle diameter ratio is higher than10; the system can define as porous media. The momentum equation is defined as:

Where

In this work the porosity of bed is 0.4, so we have following value for C and β coefficients.

The energy equation is defined as:

Where the Effective thermal diffusion coefficient is defined as:

The spices transfer equation is defined as:

S is source term. Comparing the default of spices transfer equation with spices transfer equation of this system can define the

C_{fs} was the concentration of Thymol & Carvacrol in the fluid face at the particle surface and was defined as:

and

In such works, the Reynolds number is defined as:

Re_{d}<1 Darcy regime or creeping flow

1<Re_{d}<150 inertia regime

150<Re_{d}<300 unsteady laminar flow regime

Re_{d}>300 unsteady Turbulent flow regime

The superficial velocity in this packed bed is defined as dividing flow rate by cross section.

Dividing superficial velocity by porosity, real velocity can calculate.

and for calculating D_{p}

In this work the mean particle diameter was 0.6 mm.

In this work the Reynolds number is 0.2, so there is laminar flow in this system.

The material of system was mixture of Thymol & Carvacrol and water. The inlet flow was contained water with mole fraction 1 and the outer flow was consisted of water and Thymol & Carvacrol [13]. The packed bed had 103 mm height and the diameter of bed was 16 mm. The uniform mesh was used for this domain.

The governing equations are solved by a finite volume method. At main grid points placed in the center of the control volume, volume fraction, density and spices fraction are stored. The conservation equations are integrated in space and time. This integration is performed using first order upwind differencing in space and fully implicit in time. For a first-order upwind solution, the value at the center of a cell is assumed to be an average throughout the cell. The SIMPLE algorithm is used to relate the velocity and pressure equations.

The simulation geometry is shown in Fig.8. The 2D calculation domain is divided into 22*146 grid nodes, in the radial and axial directions, respectively. The grid and mesh size are chosen to be uniform in the two directions. The inlet of system was pure water and the outlet of system was extracted Thymol & Carvacrol.

The extraction values versus time for 150^{o}C were shown in Fig. 9. It can be seen, the extraction values were increased as exponentially. After 60 min, the change of extraction values with time was very little.

In order to investigate the applicability of the CFD model, the theoretical results are compared with experimental measurements obtained at optimum conditions (20 bars, 150˚C, and 2 ml/min). Fig. 9 shows that the modeling extraction values profile is increasing rapidly in the period of 0-60 min and thereafter in the second region (60-120 min) the slope reduces until reaches a constant trend in the third period of 120-150 min. In the first region, because of high Thymol & Carvacrol concentrations in the *Z. multiflora* leaves and therefore, high mass transfer driving force, high desorption rate of Thymol & Carvacrol from solid matrix occurs.

## 6. Conclusion

It was tried to give overall view about subcritical water extraction. Effective parameters, mechanism and modeling of extraction were surveyed. Overall by considering mean average errors of models, a mathematical model base on the combination of partition coefficient (*K*D) and external mass transfer gave a good description of subcritical water extraction of *Z*. *multiflora*, while the kinetic model reasonably described the extraction behavior at lower flow rates [11].

On the other side, model was developed by introducing differential mass balances using two phase model, and applying a linear equilibrium relationship. Because of considering the effect of variation of the concentration profile in the SW phase, it seems that the proposed model is more significant from the physical point of view.

The Last model was CFD modeling of extraction from *Z. multiﬂora* leaves using subcritical water. It was concluded that CFD is poised to remain at the forefront of cutting edge research in the sciences of fluid dynamics and mass transfer. Also, the emergence of CFD as a practical tool in modern engineering practice is steadily attracting much interest and appeal. The results of CFD model have been agreed well with experimental data. As shown, along of extractor, Thymol was extracted and was in outflow.

## Notations

specific surface, m^{2}/m^{3} ( = | a or a_{p} | |||

Biot number ( = | Bi | |||

solute concentration in the solid phase, kmol/m^{3} | C | |||

solute concentration in the SW phase, kmol/m^{3} | C_{f} | |||

solute concentration in the fluid phase at the particle surface, kmol/m^{3} | C_{fp} | |||

solute concentration in the solid phase, kmol/m^{3} ( | C_{p} | |||

initial solute concentration in the solid phase, kmol/m^{3} | C_{po} | |||

diffusivity of solute (A) in liquid (B), m^{2}/s | D_{AB} | |||

effective diffusion coefficient (m^{2}/s) | D_{e} | |||

axial dispersion coefficient, m^{2}/s | D_{L} | |||

diffusivity in the solid, m^{2}/s | D_{m} | |||

particle diameter, m | d_{p} | |||

exponential function | e | |||

a certain fraction of the analyte desorbs at a fast rate by k_{1} | F | |||

remaining fraction desorbs at a slower rate by k_{2} | (1-F_{1}) | |||

thermodynamics partitioning coefficient | K _{D} | |||

external mass transfer coefficient (m/min) | k_{e} | |||

mass transfer between liquids and beds of spheres | k_{f} | |||

volumetric partition coefficient of the solute between the solid and the fluid phase | k_{p} | |||

cumulative mass of analyte extracted after certain amount of volume V_{a} (mg/g dry sample) | M_{a} | |||

cumulative mass of the analyte extracted after certain amount of volume V_{b} (mg/g dry sample) | M_{b} | |||

total initial mass of analyte in the matrix (mg/g dry sample) | M_{i} | |||

cumulative fraction of the analyte extracted by the fluid of the volume V_{b} and V_{a} (ml) | M_{b}/M_{i} and M_{a}/M_{i} | |||

total amount of solute (mg/g) removed from particle after time t | M_{t} | |||

maximum amount (mg/g) of solute extracted after infinite time | M_{∞} | |||

ratio of total migration to the maximum migration concentration | M_{t}/M_{∞} | |||

mass of the extracted sample (mg dry sample) | m | |||

Peclet number of the bed ( = | Pe_{b} | |||

Peclet number of the solid ( = | Pe_{p} | |||

average radius of an extractable particle | R | |||

Reynolds number ( = | Re | |||

Schmidt number ( = | Sc | |||

Sherwood number ( = | Sh | |||

Stanton number ( | St | |||

Temperature, K | T | |||

superficial SW fluid velocity, m/s | u | |||

molar volume of the solute at its normal boiling point, cm^{3}/mol | V_{A} | |||

critical volume, cm^{3}/mol | V_{c} | |||

dimensionless axial coordinate along the bed, z/L | X | |||

dimensionless radius ( = | Y |

## Greek symbols

ε | void volume fraction |

μ | viscosity, Pa.s |

ρ | density, kg/m^{3} |

φ | association factor for the solvent |

τ | dimensionless time ( = |