Manh Nguyen

Univ Laval

Cuong Le Van, Cermsem

GREQAM & EUREQua Philippe Michel

N M Hung

C Le Van
email: levan@univ-paris1.fr

Non-convex agreggative technology and optimal economic growth

Keywords:

This paper examines a model of optimal growth where the agregation of two separate well behaved and concave production technologies exhibits a basic non-convexity. Multiple equilibria prevail in an intermediate range of interest rate. However, we show that the optimal paths monotonically converge to the one single appropriate equilibrium steady state.

Introduction

Problems in the one-sector optimal economic growth model where the production technology exhibits increasing return at first and decreasing return to scale afterward have received earlier attention. [START_REF] Skiba | Optimal Growth with a Convex-Concave Production Function[END_REF], examined this question in continuous time and provided some results, which were further extended rigorously in Majumdar and Mitra (1982) for a discrete time setting. With a convex-concave production function, it has been shown that the time discount rate plays an important role: when the future utility is heavily discounted, the optimal program converges monotonically to the "low" steady state (possibly the degenerated state characterized by vanishing long run capital stock) while in the opposite case, it tends in the long run to the optimal steady state, usually referred to as the Modified Golden Rule (MGR) state in the literature. [START_REF] Dechert | A complete characterization of optimal growth paths in an aggregated model with a non-concave production function[END_REF] further showed that the corresponding dynamic convergence is also monotonic. They equally pointed out that this convergence now depends upon the initial stock of capital if the rate of interest falls in an intermediate range of future discounting.

In the present paper, we put emphasis on the existence of many technologyblueprint books, where each technology is well behaved and strictly concave, but the aggregation of theses technologies gives rise to some local non-convex range. Consider two Cobb-Douglas technologies depicted in Figure 1 where output per capita is a function of the capital-labor ratio. The intersection of the production graphs is located at point C where k = 1, therefore the α-technology is relatively more efficient when k ≤ 1, but less efficient than the β-technology when k ≥ 1. The two production graphs have a common tangent passing through A and B. Thus, the aggregate production which combines both α-technology and β-technology exhibits a non-convex range depicted by the contour ACB. Beside the degenerated state (0,0) in Figure1, there may exist two MGR long run equilibria b k α and b k β . In this case, we must ask which of these two states will effectively be the equilibrium, and how the latter will be attained over time.

We shall show in that when future discounting is high enough, the equilibrium is the optimal steady state b k α corresponding to the technology relatively more efficient at low capital per head. Conversely, when future discounting is low, the equilibrium is the optimal steady state b k β corresponding to the technology relatively more efficient at high capital per head. For any initial value of the initial capital stock in these cases, the convergence to the optimal steady state equilibrium is monotonic. In contrast , when future discounting is in some intermediate range, there might exist two optimal steady states and the dynamic convergence now depends on the initial stock of capital k 0 . We show that there exists a critical value k c such that every optimal path from k 0 < k c will converge to b k α , and every optimal path from k 0 > k c will converge to b k β .

The paper is organized as follows. In the section 2, we specify our model. In section 3, we provide a complete analysis of the optimal growth paths and in section 4, we summarize our findings and provide some concluding comments. .

Optimal steady states

The Model

The economy in the present paper produces a homogeneous good according two possible Cobb-Douglas technologies, the α-technology f α (k) = Ak α , and the β-technology f β (k) = Ak β where k denotes the capital per head and 0 < α < β < 1. The efficient technology will be y = max © Ak α , Ak β ª = f (k). The convexified economy is defined by cof (k) where co stands for convexhull. It is the smallest concave function minorized by f. Its epigraph, i.e. the set

{(k, λ) ∈ R + × R + : cof (k) ≥ λ} is the convex hull of the epigraph of f, {(k, λ) ∈ R + × R + : f (k) ≥ λ} (see figure 1). One can check that cof = f α for k ∈ [0, x 1 ] , cof = f β for k ∈ [x 2
, +∞[, and affine between x 1 and x 2 . More explicitly, we have

αAx α-1 1 = βAx β-1 2 = Ax α 1 -Ax β 2 x 1 -x 2 which implies x 1 = µ α β ¶ β β-α µ 1 -α 1 -β ¶ 1-β β-α and x 2 = µ α β ¶ α β-α µ 1 -α 1 -β ¶ 1-α β-α .
In our economy , the social utility is represented by P t=+∞ t=0 γ t u(c t ) where γ is the discount factor and c t the consumption. At period t, this consumption is constrained by the net output f (k t )k t+1 , where k t denotes the per head capital stock available at date t.

The optimal growth model can be written as

max +∞ X t=0 γ t u(c t )
under the constraints

∀t ≥ 0, c t ≥ 0, k t ≥ 0, c t ≤ f (k t ) -k t+1 , and k 0 > 0 is given.
We assume that the utility function u is strictly concave, increasing, continuously differentiable, u(0) = 0 and (Inada Condition) u 0 (0) = +∞. The discount factor γ satisfies 0 < γ < 1.

Let V denote the value-function, i.e.

V (k 0 ) = max +∞ X t=0 γ t u(c t )
under the constraints

∀t ≥ 0, c t ≥ 0, k t ≥ 0, c t ≤ f (k t ) -k t+1
, and k 0 ≥ 0 is given.

Remark 1: Before proceeding the analysis, we wish to say that our technology specification used for aggregation purpose in this paper is not restrictive. Indeed, consider the following production function

f (k) = max{Ak α , Bk β }, with A 6 = B. Define e k = k λ , e c = c λ , v(c) = u( c λ ), where λ satisfies Aλ α = Bλ β . Let A 0 = Aλ α = Bλ β .
It is easy to check that the original optimal growth model behind becomes

max +∞ X t=0 γ t v(e c t )
under the constraints ∀t ≥ 0, e c t ≥ 0, e k t ≥ 0, e c t ≤ e f ( e k t ) -e k t+1 , and e k 0 > 0 is given.;

where e f (x) = max{A 0 x α , A 0 x β }.

Analysis of the optimal growth paths

The preliminary results are summarized in the following proposition.

Proposition 1 (i) For any k 0 ≥ 0, there exists an optimal growth path (c * t , k * t ) t=0,...,+∞ which satisfies:

∀t, 0 ≤ k * t ≤ M = max h k 0 , e k i , 0 ≤ c * t ≤ f (M),
where e k = f ( e k).

(ii) If k 0 > 0, then ∀t, c * t > 0, k * t > 0, k * t 6 = 1, and we have Euler equation

u 0 (c * t ) = γu 0 (c * t+1 )f 0 (k * t+1 ).
(iii) Let k 0 0 > k 0 and (k 0 * t ) be an optimal path associated with k 0 0 . Then we have: ∀t, k 0 * t > k * t . (iv) The optimal capital stocks path is monotonic and converges to an optimal steady state. Here, this steady state will be either b k α = (γAα)

1 1-α or b k β = (γAβ) 1 1-β .
Proof. (i) The proof of this statement is standard and may be found in Le [START_REF] Van | Dynamic Programming in Economics[END_REF], chapter 2. (ii) From Askri and Le [START_REF] Askri | Differentiability of the value function of nonclassical optimal growth models[END_REF], the value-function

V is differentiable at any k * t , t ≥ 1. Moreover, V 0 (k * t ) = u 0 (f (k * t ) -k * t+1 )f 0 (k * t )
and this excludes that k * t = 1 since 1 is the only point where f is not differentiable. From Inada Condition, we have c * t > 0, k * t > 0, ∀t. Hence, Euler Equation holds for every t.

(iii) It follows from [START_REF] Amir | Sensitivity analysis in multisector optimal economic models[END_REF]

that k 0 0 > k 0 implies ∀t, k 0 * t > k * t .
From Euler Equation we have

u 0 (f (k 0 ) -k * 1 ) = γV 0 (k * 1 )
and

u 0 (f (k 0 0 ) -k * 1 ) = γV 0 (k 0 * 1 ). If k * 1 = k 0 * 1 then k 0 = k 0 0 : a contradiction. Hence, k * 1 < k 0 * 1 . By induction, ∀t > 1, k 0 * t > k * t . (iv) First assume k * 1 > k 0 . Then the sequence (k * t ) t≥2 is optimal from k * 1 . From (iii), we have k * 2 > k * 1 . By induction, k * t+1 > k * t , ∀t. If k * 1 < k 0 , using the same argument yields k * t+1 < k * t , ∀t. Now if k * 1 = k 0 ,
then the stationary sequence (k 0 , k 0 , ..., k 0 , ...) is optimal.

We have proved that any optimal path (k * t ) is monotonic. Since, from (1), it is bounded, it must converge to an optimal steady state k s . If this one is different from zero, then the associated optimal steady state consumption c s must be strictly positive from Inada Condition. Hence, from Euler Equation, either

k s = b k a or k s = b k b since it could not equal 1.
It remains to prove that (k * t ) cannot converge to zero. On the contrary, for t large enough, say greater than some T, we have u 0 (c * t ) > u 0 (c * t+1 ) since f 0 (0) = +∞. Hence, c * t+1 > c * t for every t ≥ T. In particular, c * t+1 > c * T > 0, ∀t > T. But k * t → 0 implies c * t : a contradiction. We obtain the following corollary:

Corollary 2 If γAα > 1, then any optimal path from k 0 > 0 converges to b k β . If γAβ < 1, then any optimal path from k 0 > 0 converges to b k α .

Proof.

In Proposition 1, we have shown that any optimal path (k * t ) converges either to b k α or to b k β . But when γAα > 1, we have b

k α > 1, f ( b k α ) = A( b k α ) β and f 0 ( b k α ) = βA( b k α ) β-1 6 = 1 γ .
Consequently, b k α could not be an optimal steady state. Therefore, (k * t ) cannot converge to b k α . From the statement (iv) in Proposition 1, it converges to b k β .

Similarly, when γAβ < 1, any optimal path from k 0 > 0 converges to b k α .

In Figure 1, when b k α ≥ 1, α-technology is clearly less efficient than β-technology, thus b k α is not the optimal steady state. Similarly for b k β ≤ 1. In these cases, there will be an unique optimal steady state. But when the discount factor is in an intermediate range defined by γAα ≤ 1 ≤ γAβ, there exists more than one such state. We now give an example where b k α and b k β are both optimal. Since x 1 and x 2 are independent of A and γ, we can choose A and γ such that

αAx α-1 1 = βAx β-1 2 = 1 γ , with 0 < γ < 1.
It is easy to check that x 1 and x 2 are optimal steady states for the convexified technology and hence for our technology. Since x 1 = b k α , x 2 = b k β , we have found two positive optimal steady states.

Let now b k α and b k β , depicted in Figure 1, be two optimal steady states and ask the question which of them will be the long run equilibrium in the optimal growth model. We first get an immediate result in:

Proposition 3 Assume γAα ≤ 1 ≤ γAβ. If γAα is close to 1, then any optimal path (k * t ) from k 0 > 0 converges to b k β . If γAβ is close to 1, then (k * t ) converges to b k α . Proof. First, observe that when γAα ≤ 1 then f ( b k α ) = A ³ b k α ´α and when 1 ≤ γAα, f( b k β ) = A ³ b k β ´β . Now consider the case γAα = 1 < γAβ. We have b k α = 1 and A > 1.
It is well-known that given k 0 > 0, there exists a unique optimal path from k 0 for the β-technology. Moreover, this optimal path converges to b k β .

Observe that the stationary sequence

³ b k α , b k α , ..., b k α , ... ´is feasible from b k α , for the β-technology, since it satisfies 0 ≤ b k α = 1 < A ³ b k α ´β = A. Hence, if ³
e k t ´is an optimal path for β-technology starting from b k α and if (k t ) is an optimal path of our model starting also from b k α , we will have

∞ X t=0 γ t u(f ( b k α ) -b k α ) < ∞ X t=0 γ t u(f ( e k t ) -e k t+1 ) ≤ ∞ X t=0 γ t u(f (k t ) -k t+1 ) = V ( b k α ).
That shows that b k α can not be an optimal steady state. Hence, any optimal path from k 0 > 0 must converge to b k β .

Since b k α is continuous in γ, V continuous and since

P ∞ t=0 γ t u(f ( b k α ) - b k α ) < V ( b k α ) when γAα = 1,
this inequality still holds when γAα is close to 1 and less than 1. In other words, b k α is not an optimal steady state when γAα is close to 1 and less than 1. Consequently, any optimal path with positive initial value will converge to b k β .

Similar argument applies when γAβ is near one but greater than one.

What then happens when the discount factor is within an intermediate range? We now would like to show : Proposition 4 Assume γAα < 1 < γAβ. If both b k α and b k β are optimal steady states then there exists a critical value k c such that every optimal paths from k 0 < k c will converge to b k α , and every optimal paths from k 0 > k c will converge to b k β .

Proof. Consider at first k 0 < b k α . Since b k α is optimal steady state, we have k * t < b k α , ∀t > 0. Since the sequence (k * t ) is increasing, bounded from above by b k α , it will converge to b k α . Similarly, when k 0 > b k β , any optimal path converges to b k β .

Let k = sup

n k 0 : k 0 ≥ b k α o such that any optimal path from k 0 converges to b k α . Obviously, k ≤ b k β , since b k β is optimal steady state. Let k = inf n k 0 : k 0 ≤ b k β o such that any optimal path from k 0 converges to b k β . Obviously, k ≥ b k α , since b k α is optimal steady state. We claim that k = k. It is obvious that k ≤ k. Now, if k < k, then take k 0 , k 0 0 which satisfy k < k 0 < k 0 0 < k.
From the definitions of k and k, there exist an optimal path from k 0 , (k * t ), which converges to b k β and an optimal path from k 0 0 , (k 0 * t ), which converges to b k α . For t large enough, k 0 * t < k * t , which is impossible since k 0 < k 0 0 (see Proposition 1, statement (iii)). Posit k c = k = k and conclude.

Remark 2: The existence of critical value is standard since the paper by [START_REF] Dechert | A complete characterization of optimal growth paths in an aggregated model with a non-concave production function[END_REF]. See also, for the continuous time setting. Askenazy and Le [START_REF] Askenazy | A model of optimal growth strategy[END_REF]. But in these models, the technology is convexconcave. The low steady state is unstable while the high is stable. An optimal path converges either to zero or to the high steady state. In our model, with a technology, say concave-concave, any optimal path converges either to the high steady state or the low steady state. This certainly helps identifying the production function for data aggregation task. If the future discount rate falls in a range defined by γAα < 1 < γAβ, then the computation of the critical capital stock k c is essential in view of the determination of the relevant production technology at stake. Third, when there are several production technologies, it is possible to proceed with pair-wise aggregation in order to determine the relevant technology for long run equilibrium. Assume that we have a third technology, say the ε-technology, to take into account. Pair-wise aggregation of α and β-technology allows us to eliminate the α-technology, say. Therefore, we now have to perform the same analysis with β -technology and ²-technology, and so on so forth when several technologies are at stake. Pair-wise consideration in this way would help determining the relevant technology corresponding to the optimal steady state. Fourth, under the regular conditions of concavity in Ramsey model, the long run equilibrium could be achieved with a decentralized market mechanism. Recall that non-convexity is thought to be the main cause for market failure. But in the case we consider in this paper, the economy attains one (and only one) long run equilibrium corresponding to a well-behaved concave production technology. Therefore, despite the non-convexity arising from technology aggregation, there is no market failure and decentralized allocation will indeed implement the Modified Golden Rule State.

Figure 1 :

 1 Figure 1:

Concluding commentsIt is shown in this paper that when future discounting is high enough, precisely when γAβ < 1, the resulting long run equilibrium is the optimal steady state b k α . For any value of the initial capital stock, the convergence to this equilibrium is monotonic. On the other hand, when future discounting is relatively low, precisely when γAα > 1, the same result will be obtained but with the equilibrium optimal steady state b k β . When future discounting is in some middle rang, i.e. when γAα < 1 < γAβ, there might exist two optimal steady states and the dynamic convergence will depend on the initial stock of capital. We show that there is a critical capital stock k c such that every optimal paths from k 0 < k c will converge to b k α , and every optimal paths from k 0 > k c will converge to b k β .Several useful remarks can be made. First, it is conceivable that the results obtained in this paper are unaffected when either one or both production technologies entails some fixed costs, i.e. positive output is made possible only if the capital per capita exceeds a threshold level, but their aggregation exhibits the kind of non-convexity depicted in Figure1. Second, for the economist-statisticians, this paper hopefully highlights the importance of informations other than those contained in the technology-blueprint book. Under either high or low future discounting, only one technology is relevant in the sense that it is the chosen technology in long run equilibrium.


Nicolás Gañán

Jérôme Morchain

Séverine Camy

Jean-Stéphane Condoret

Rate-based simulation of a high pressure counter-current packed column for supercritical CO 2 extraction of alcohol from dilute aqueous mixtures

In this work, the modeling and simulation of a fractionation packed column for the recovery of isopropanol from dilute aqueous mixtures using supercritical CO 2 is presented. The model is based on the numerical resolution of differential mass balances for each component over the column height. The multicomponent mass transfer between phases is described using a "rate-based" approach and the concept of local mass transfer coefficients. The model was validated by reproduction of experimental steady-state results for the fractionation of 5% isopropanol aqueous solutions obtained in a bench scale counter-current column. The effect of process conditions on the separation performance was satisfactorily described by the model calculations, showing that operation pressure and CO 2 flow rate enhance IPA recovery and extract purity, while operation temperature has a negative effect. Model deviations (AARD) were in all cases lower than 20%.

Introduction

Fractionation of liquid mixtures using supercritical solvents has been proposed and studied as a promising tool for new separation challenges, as well as an alternative to current separation processes [START_REF] Brunner | Counter-current separations (Review)[END_REF][START_REF] Bejarano | Fractionation technologies for liquid mixtures using dense carbon dioxide[END_REF]. In this way, the use of supercritical carbon dioxide (scCO 2 ) as a solvent proposes a more environmentally friendly technology when compared to conventional organic solvents. The main potential appli cations include deterpenation of citrus oils [3 6], purification and removal of solvent residues from edible oils [START_REF] Espinosa | Process analysis and phase equilibria for the removal of chemicals from fatty oils using near-critical solvents[END_REF], fractionation of fish oils [START_REF] Riha | Separation of fish oil ethyl esters with supercritical carbon dioxide[END_REF][START_REF] Catchpole | Fractionation of fish oils using supercritical CO 2 and CO 2 +ethanol mixtures[END_REF], recovery of high value lipid compounds from edible oils or by products [START_REF] Vázquez | Recovery of squalene from vegetable oil sources using countercurrent supercritical carbon dioxide extraction[END_REF][START_REF] Fornari | Countercurrent supercritical fluid extracion of different lipid-type materials: experimental and thermodynamic modelling[END_REF] and extraction of alcohols and other organic sol vents from aqueous solutions [START_REF] Brignole | Supercritical fluid extraction of alcohols from water[END_REF][START_REF] Budich | Supercritical fluid extraction of ethanol from aqueous solutions[END_REF].

In a typical fractionation process, the mixture to be separated and the supercritical solvent are fed to a separation column (generally a packed column) operated in counter current mode. During contact, the components of the mixture distribute between the liquid and super critical phases according to their volatility and affinity for scCO 2 to an extent determined by the prevailing temperature, pressure and com position conditions, as well as the effective mixing and contact area of the phases. The supercritical phase leaving the column at the top ("extract") is then totally or partially expanded in a separator where solutes condense and are recovered. The processed liquid phase leaves the column at the bottom ("raffinate"). The general features of this process are similar to the gas liquid contacting operations, and there fore its mathematical description is usually analogous or directly based on the theories developed for those operations. The main differences are connected to the particular properties of the supercritical phase (which cannot be considered simply as a gas, because of its significantly higher density and viscosity) together with the complexities of high pressure phase equilibria in ternary or multicomponent systems.

Due to the usually large solvent flow rates required by such pro cesses when CO 2 is used, process optimization is fundamental to assess economic feasibility, including not only the fractionation column but also the separators and solvent regeneration cycle. Costly and time consuming experimental work [START_REF] Díaz | Modeling and optimization of supercritical fluid processes[END_REF] can be alleviated through the use of computational tools based on robust and reliable thermodynamic and mass transfer models.

Similarly to gas liquid and liquid liquid operations, two different approaches can be proposed for modeling and simulation of super critical fractionation columns: (i) methods based on the concept of equilibrium stage, and (ii) rate based or "non equilibrium" methods, based upon hydrodynamic and interfacial mass transfer along with a thermodynamic equilibrium assumption at the interface.

In "equilibrium stage" methods the column is modeled as a suc cession of discrete units where the phases mix and reach equilibrium before separating. Concentrations are assumed to be homogeneous in each phase. This approach is applied to packed columns through the concept of "height equivalent to a theoretical stage" (HETS), i.e., the height of packed bed which produces a separation output similar to an equilibrium stage. Many commercial process simulation softwares propose calculation modules for contactors based on this approach. Although easy to implement and very useful to assess and compare performances of a complete process, equilibrium stage methods are somewhat limited due to the poor generalization capability of the concept of HETS or "stage efficiency", which is very dependent on the physical properties of the system and the flow rates, as well as the contacting technology. Therefore their evaluation is complex and highly sensitive to the operating conditions.

On the other hand, rate based methods take into account the mass transfer kinetics between phases and do not assume that the phases are globally in equilibrium after contact but only at the interface. For this reason, they are often called "non equilibrium" methods. In these models, the mass transfer and hydrodynamic particularities of each system are explicitly computed based on physical and geometrical properties. Ruivo et al. [START_REF] Ruivo | Dynamic model of a countercurrent packed column operating at high pressure conditions[END_REF] applied this approach to the dynamic si mulation of the separation of methyl oleate and squalene with scCO 2 in a counter current structured packing column. Their model consists in a set of differential mass balances for each phase and algebraic equations accounting for mass transfer and hydrodynamic aspects. Phase equili brium at the interface was taken into account through an empirical correlation, and the interfacial fluxes were calculated using overall mass transfer coefficients. With this model, steady state but also the initial transient period and the dynamic response to perturbations could be described. In a more recent work, Fernandes et al. [START_REF] Fernandes | Non-isothermal dynamic model of a supercritical fluid extraction packed column[END_REF] have ex tended this model to non isothermal operation by incorporating an energy balance. Martin and Cocero [START_REF] Martín | Mathematical modeling of the fractionation of liquids with supercritical CO2 in a countercurrent packed column[END_REF] have developed a steady state model for the scCO 2 fractionation of liquid mixtures in counter current packed columns with external reflux based on differential mass and energy balances, describing mass transfer by Maxwell Stefan multi component diffusion theory. Their model was validated against ex perimental data on edible oils fractionation and fish oil ethyl esters separation. The model was solved by an iterative method, guessing the outlet gas flow rate and composition and calculating backwards the corresponding inlet values until obtaining the specified or experimental values within a certain tolerance. The only fitted parameter was the effective interfacial area.

In the present work, a dynamic rate based model of a high pressure counter current packed column is presented and applied to the simu lation of the extraction of isopropanol (IPA) from dilute aqueous solu tions using scCO 2 . This kind of separations is of great interest for the recovery of bio alcohols in fermentation processes, and has been ex perimentally studied in bench scale fractionation column in our pre vious work [START_REF] Lalam | Fractionation of aqueous isopropanol mixtures in a countercurrent packed column using supercritical CO 2[END_REF]. Two different types of packing and various operating conditions have been studied. In that work, the fractionation process was described using the equilibrium stage approach. The number of equilibrium stages was found to be between one and two for a 2 m height column, this number being dependent on operating conditions, which reflects the fact that the stage approach is not well suited to that particular case. The same ternary system has been previously studied by Rathkamp et al. [START_REF] Rathkamp | Evaluation of packed columns in supercritical extraction processes[END_REF] who characterized mass transfer performances of a spray column, with and without a packing, and proposed a descrip tion using the conventional NTU HTU approach. In their work they considered droplets of fluid phase flowing through a continuous liquid phase, with transfer of isopropanol only, and described thermo dynamics using a constant partition coefficient, which is likely to be not enough accurate for these complex systems. We propose here to further analyse our previous experimental results making use of a dynamic rate based model. In order to preserve the numerical cost at a reason able expense, the thermodynamic equilibrium data over a large range of operating conditions are calculated first and embedded in the dynamic simulator as correlations. With this approach, the concentration profiles

= Re Q Sa μ L L L eff (6)
The effective interfacial area was assumed to be equal to the wetted packing area, which was estimated with the following correlation from Onda et al. (Eq. ( 7)):

= - ⎡ ⎣ ⎢ - ⎛ ⎝ ⎞ ⎠ ⎤ ⎦ ⎥ a a σ σ a d Ga We 1 exp 1.45 ( ) p c L p p L L eff 0.75 0.35 0.05 0.2 (7)
where σ c is the "critical" surface tension of the packing material and σ L the surface tension of the liquid. Ga L and We L are the Galilei and Weber numbers for the liquid phase, respectively. The calculation of themo dynamical and physico chemical properties is presented in the fol lowing sections. Finally, the liquid hold up (h L ) was estimated with a standard cor relation for randomly packed columns [START_REF] Otake | Liquid hold-up in packed towers, operating and hold-up without gas flow[END_REF], according to Eq. ( 9):

= - h ad G a R e 1.295( ) L pp L L 0.44 0.676 (8)
All symbols and non dimensional numbers are defined in the no menclature. All details regarding the characteristics of the fractionation column can be found in our previous work [START_REF] Lalam | Fractionation of aqueous isopropanol mixtures in a countercurrent packed column using supercritical CO 2[END_REF]. The main information is reported in Table 2.

Thermodynamic model

As mentioned, the model assumes thermodynamic equilibrium at the interface. The IPA water CO 2 ternary phase equilibrium was com puted using the Simulis ® Thermodynamics software (ProSim SA, France). A combined equation of state/activity coefficient model (EoS/ G E ) approach with suitable complex mixing rules was selected. The equation of state is the Soave Redlich Kwong equation [START_REF] Soave | Equilibrium constants from a modified Redlich-Kwong equation of state[END_REF] with the modification proposed by Boston Mathias [START_REF] Boston | Phase equilibria in a third-generation process simulator[END_REF]. The co volume b of the SRK equation is obtained from a standard mixing rule (arithmetic mean of the corresponding parameters of the pure compounds). The term a of the SRK equation is computed using the PSRK mixing rule [START_REF] Holderbaum | PSRK: a group contribution equation of state based on UNIFAC[END_REF] and the UNIQUAC model [START_REF] Abrams | Statistical thermodynamics of liquid mixture: a new expression for the excess Gibbs energy of partly or complete miscible systems[END_REF][START_REF] Anderson | Application of the UNIQUAC equation to calculation of multi-component phase equilibria: 1. Vapor-liquid equilibria[END_REF] is used to calculate the excess energy at zero pressure, G E 0 , involved in the PSRK mixing rule. More detailed in formation about the modeling of the system phase behaviour can be found in our previous work [START_REF] Lalam | Fractionation of aqueous isopropanol mixtures in a countercurrent packed column using supercritical CO 2[END_REF].

Physico chemical properties

The equations involved in the column model require the knowledge or an estimation of several physico chemical properties of the pure compounds. The needed physico chemical properties of pure IPA, water and CO 2 are presented in Table 1, as obtained from NIST.

The density and viscosity of the fluid and liquid phases were as sumed to be constant along the column and equal to those of pure CO 2 and pure water at the specified pressure and temperature conditions, neglecting the effect of the low concentration of solutes and the com position changes along the column. This simplification, limited to dilute liquid solutions, was verified with available experimental data.

For density of the liquid phase, experimental results reported by Hebach et al. [START_REF] Hebach | Density of water + carbon dioxide at elevated pressures: measurements and correlation[END_REF] indicate that the density of CO 2 saturated water within the range of operation conditions covered in this work differs only around 1% with respect to pure water at the same conditions. Regarding the fluid phase, to the best of our knowledge there are no density data for the systems CO 2 + water or CO 2 + isopropanol at su percritical conditions in the open literature. However, as a reference, data of saturated density of the system CO 2 + ethanol at near critical and supercritical conditions (40 and 55 °C, 7 9 MPa) indicate that the difference with pure CO 2 density is below 12% [START_REF] Tsivintzelis | Phase compositions and saturated densities for the binary systems of carbon dioxide with ethanol and dichloromethane[END_REF]. Calculations performed with the Simulis © software using our thermodynamic model at the inlet: flow rate and composition are constant at the outlet: zero gradient for all quantities.

A second order upwind scheme is used for internal nodes (k = 3, …,N -1) as well as for the definition of the Neumann boundary Note that the Reynolds number for the liquid phase has to be cal culated in terms of effective interfacial area (a eff ), as in Eq. ( 6):

for a mixture consisting in 99% CO 2 , 0.5% water and 0.5% IPA (close to the composition of the fluid phase in the column) predict a 9% devia tion with respect to pure CO 2 density, as provided by the NIST data base.

Supercritical fluid phase viscosity was estimated with the method of Lucas [START_REF] Poling | The Properties of Gases and Liquids[END_REF] for dilute mixtures of IPA and water (up to 2%) in scCO 2 and was found similar to the viscosity of pure CO 2 under the same condi tions within a maximum deviation range of 15%. For the viscosity of the liquid phase, data reported by Tanaka et al. [START_REF] Tanaka | Viscosity of (water + alcohol) mixtures under high pressure[END_REF] show that the viscosity of aqueous solutions of IPA under pressure increases with IPA concentration, reaching a maximum at around 60% wt. On the other hand, it is known that the dissolution of supercritical solvents in liquids decreases their viscosity. For example, data of CO 2 + oleic acid and CO 2 + squalene show that the saturated liquid viscosity is reduced with increasing pressure, along with a higher dissolution of CO 2 [START_REF] Brunner | Gas Extraction: An Introduction to Fundamentals of Supercritical Fluids and the Application to Separation Processes[END_REF]. To the best of our knowledge, there is no information regarding the viscosity of water saturated with CO 2 at high pressure conditions, but a similar behaviour is expected to occur, although in a lesser extent considering that scCO 2 solubility in water is lower than in lipidic liquids. Therefore, considering the opposite effects of IPA and CO 2 on the aqueous phase viscosity, the latter was assumed to be equal to that of pure water.

The surface tension of the liquid phase (σ L ) was estimated using the method of Tamura, recommended for aqueous solutions of organic compounds [START_REF] Poling | The Properties of Gases and Liquids[END_REF], corrected to account for the effect of pressure. The details of this method are given in Appendix A. Briefly, the mixture surface tension is calculated by interpolation of the pure compound values as a function of the mixture composition. The model was fitted to experimental data of surface tension for the ternary system IPA + water + CO 2 at high pressure conditions [START_REF] Chun | Interfacial tension in high-pressure carbon dioxide mixtures[END_REF], and allowed es timation of the surface tension variation along the column. This was thought to be necessary due to the fact that this property is sensitive to the presence of IPA and has a strong influence on the value of the wetted or effective interfacial area, a key parameter for mass transfer.

The diffusion coefficients of IPA and water in the supercritical phase (D 13 and D 23 ) and of IPA and CO 2 in the liquid phase (D 12 and D 32 ) were estimated using the conventional method of Wilke and Chang [START_REF] Poling | The Properties of Gases and Liquids[END_REF] for liquids. The calculation details are given in Appendix B. Although this method provides an estimation of binary diffusion coefficients at in finite dilution, the values obtained are expected to hold because of the low concentration conditions that prevail in both phases in the case of dilute mixtures fractionation.

Numerical resolution

The model was implemented in the Matlab © software. The differ ential mass balance equations (Eqs. ( 1) and ( 2)) were solved using a finite difference discretization method based on N grid nodes. The first and last nodes receive the conventional following boundary conditions: condition at node N, according to Eq. ( 9):

∂ ∂ = - + + - - Q x z Q x Q x Q x z O z ( ) 3( ) 4( ) ( ) 2Δ (Δ ) L i k L i k L i k L i k 1 2 2 (9)
A second order centred scheme is used for the second node (k = 2), according to Eq. [START_REF] Vázquez | Recovery of squalene from vegetable oil sources using countercurrent supercritical carbon dioxide extraction[END_REF].

∂ ∂ = - + + - Q x z Q x Q x z O z ( ) ( ) ( ) 2Δ (Δ ) L i k L i k L i k 1 1 2 (10)
In this way, the set of 6 partial differential equations (3 components and 2 phases) is transformed into a set of 6N ordinary differential equations. The resulting system of ordinary differential equations and algebraic equations is integrated over time using a fourth order explicit Runge Kutta (RK4) method. For each set of operating conditions, the time step is automatically adjusted in order to meet the stability cri terion (current number < 1).

In general, concentrations as well as total and partial flow rates are function of time and location along the column. When the solution to be purified is concentrated, the solute extraction by scCO 2 may modify significantly the liquid flow rate (and all associated properties). Also, the dissolution of CO 2 in the liquid phase may be important in some cases, and therefore the fluid flow rate can be reduced in the same way. However, as previously mentioned, we consider in the present work dilute liquid feeds so it is acceptable to consider that the flow rates and the liquid hold up are practically constant along the column. Therefore the resolution of the mass balance equations is simplified: for each time step they are solved considering constant h L , Q L and Q F along the column, and only the initial composition profile of each phase is re quired as initial condition. This assumption holds as long as the net mass flux through the interface is small as compared to the mass flow rate of each phase. This clearly constitutes a limitation of our modeling approach which is indeed only valid for dilute mixtures but actually corresponds to the recovery of dilute alcohols from fermentation broths.

Solving this set of equations requires the calculation of the equili brium concentrations at the interface for each node at each time step. Implementing these thermodynamic equilibrium calculations using Simulis © proved to yield very lengthy computations because of the numerous calls to an external procedure (i.e., the Simulis © software toolboxes). To improve the computational efficiency of this dynamic simulator, ternary equilibrium data were beforehand computed for a wide range of compositions, at the specified operation temperature and pressure, and stored as a set of cubic spline interpolation functions.

Results and discussion

Simulation parameters and preliminary runs

As mentioned, the model was validated by fitting experimental steady state fractionation data of dilute aqueous mixtures of IPA using scCO 2 , obtained by Lalam et al. using a bench scale counter current column, as reported in a previous work [START_REF] Lalam | Fractionation of aqueous isopropanol mixtures in a countercurrent packed column using supercritical CO 2[END_REF]. The adjustable coefficients of Eqs. ( 4) and [START_REF] Budich | Countercurrent deterpenation of citrus oils with supercritical CO 2[END_REF] were obtained in order to reproduce the recovery ratio and mass fraction of IPA in the raffinate stream under different temperature, pressure and solvent to feed ratio conditions.

Table 2 shows the column and packing characteristics used in the experiments, required by the model as input information. Two types of packing have been evaluated: stainless steel springs and a porous me tallic (FeCrAlY) foam.

A first set of preliminary runs were performed in order to determine appropriate discretization and integration parameters. In order to ob tain reproducible and stable numerical results, 100 grid points (each of 2 cm height) and a time integration step of 0.3 s were sufficient in most cases. Calculations were carried out until the outlet concentrations were constant (steady state conditions). However, the accurate prediction of the transient period duration was not an objective of this work. The dynamic simulation approach is used here as a strategy to avoid the initialization problem that usually arises in steady state models, which are solved iteratively and whose performance rely on the quality of the initial guess of the outlet conditions [START_REF] Martín | Mathematical modeling of the fractionation of liquids with supercritical CO2 in a countercurrent packed column[END_REF].

As mentioned, due to the constant flow rate hypothesis, only the initial composition profiles along the column are required as initial conditions. Different initialization strategies were tested. The fluid phase was considered initially as pure scCO 2 in all runs, but different initial profiles for the liquid phase were tested: (i) pure water; (ii) a composition profile equal to that of the liquid feed, and (iii) variable composition profile with a linear gradient of IPA from the liquid feed composition at the top to zero at the bottom. In all cases, the model proved to converge toward the same steady state solution, with dif ferences regarding the duration and behaviour of the transient period, as shown in Fig. 2 for a given set of operating conditions. Notice that the curves corresponding to the initial liquid phase configurations (ii) and (iii) almost completely overlap.

The model also provides information about the calculated con centration profiles along the column, as well as the interfacial mass transfer fluxes and partial flow rates. As an example, Fig. 3 shows the calculated profiles for the fractionation of a 5% IPA solution at 40 °C, 10 MPa and a solvent to feed ratio (Q F /Q L ) equal to 9.5 kg/kg using the springs packing. This kind of information can be very useful for tech nological improvement, as it allows detecting if the mass transfer is low in a certain part of the column. It can be checked that the magnitude of the mass flux between phases remains small compared to the mass flux of each phase (total liquid mass flux is equal to 0.587 kg/(m 2 s)) thus validating the nearly constant mass flow rate assumption.

Influence of pressure and temperature

Table 3 shows the density and viscosity values for the liquid and fluid phase used in the simulations at all the studied temperature and pressure conditions. Table 4 shows the calculated IPA diffusion coeffi cients and mass transfer coefficients in the liquid and fluid phases. Note that, due to the hypothesis of constant fluid properties along the column, single values of k L and k F are calculated for each condition. Fig. 4 shows numerical and experimental results of IPA mass frac tion and recovery ratio of water and IPA in the raffinate phase as a function of operating pressure, at constant temperature (40 °C) and Q F / Q L ratio (6.04 kg/kg), for the fractionation of a 5% IPA solution using springs packing. Experimental results are given in terms of composition and recovery ratio in raffinate, because, on an experimental point of view, this outlet phase is easier to recover and analyse than the extract phase. Recovery ratio in the raffinate is obtained by dividing the flowrate of the considered compound in the raffinate by the flowrate of the compound in the feed.

The experimental reproducibility has been assessed and demon strated in our previous work [START_REF] Lalam | Fractionation of aqueous isopropanol mixtures in a countercurrent packed column using supercritical CO 2[END_REF]. Each run was reproduced in quad ruplicate and results were found reproducible with a mean deviation equal to 0.005% for raffinate IPA mass fraction and 0.35% for raffinate The values of the coefficients of the Onda correlations (α and β) were manually fitted in order to reproduce the experimental results, with a good agreement using α = 1.05 and β = 0.001, valid for all experiments related to the springs packing at 40 °C. With these values, the absolute average relative deviation (AARD) was 16.4% for the IPA mass fraction and 17.6% for the IPA recovery, both in the raffinate. This implies that the mass transfer coefficient values predicted by the ori ginal correlations of Onda are divided approximately by 5 in order to reproduce the experimental data obtained in our column. Onda's cor relations were developed from data obtained in conventional gas liquid systems, at atmospheric pressure, and therefore may not be directly applicable to our specific case, in which the gas phase is replaced by a dense supercritical phase and the packing type and column dimensions may induce specific flow patterns. Therefore, the original coefficients were replaced by adjustable parameters in order to account for these features. Concerning flow patterns, it must be remembered that the correlations of Onda were developed from data obtained in gas ab sorption packed columns with a high ratio between column diameter and packing size. In our fractionation column, this ratio is lower and wall effects are suspected to be much more pronounced, reducing the effective interfacial area and the column performance. Moreover, axial dispersion has not been taken into account in the model. Thus, it has been chosen to reflect these phenomena, including pressure effects, in the values of mass transfer coefficients. In Fig. 4, it can be observed that the model describes correctly the enhancing effect of pressure on the column efficiency within the studied range. Although the fluid phase density and viscosity increase with pressure (see Table 3), with a subsequent decrease of the mass transfer coefficient, this effect is compensated by an IPA partition coefficient (m 1 ) enhancement (almost 2 fold) as CO 2 density increases. This can be seen in Fig. 5, which shows the representative values of m 1 as well as the IPA mass transfer resistance in both phases, defined according to Eqs. ( 13) and ( 14):

= R k m 1 F F 1, 1, 1 (13)
= R k 1 L L 1, 1, (14)
The m 1 values shown in Fig. 5 are mean values calculated for comparison using the column global composition. Although the parti tion coefficient is composition dependent, in the case of dilute mixtures its variation is limited within a narrow range, as can be seen for a ty pical run in Fig. 6.

As can be observed in Table 4, the k L values calculated by the model are practically independent of pressure, because the change in the li quid relevant properties (density and viscosity) with pressure is very low. Therefore, the mass transfer resistance on the liquid side also re mains almost constant, as shown in Fig. 5. On the other hand, the k F It can be seen that the model was able to fairly reproduce the ex perimental data trend (AARD of 21.3% and 20.6% for IPA concentra tion and recovery, respectively), with a good matching for Q F /Q L ≥ 6 but higher deviations at Q F /Q L < 6. In fact, for low Q F /Q L ratios (< 6) the model overpredicts the fractionation performance, yielding to lower concentration values of IPA in the raffinate phase than those observed experimentally. In this low Q F /Q L region, lower values of the para meters α and β should be used in order to reproduce correctly the ex perimental results, as discussed in the next paragraph. This possibly indicates a change in the column hydrodynamic regime, which is not explicitly taken into account in the applied correlations. In particular, as low Q F /Q L ratio correspond to low CO 2 flow rate, a hypothesis could be that the column is operated at conditions below loading point, where liquid and CO 2 phases weakly interact.

Simulations were also performed at higher Q F /Q L ratios (up to 20), predicting that for Q F /Q L > 15 the recovery and concentration of IPA in the raffinate change at a very low rate. Operation in this region becomes more expensive, as the cost of CO 2 recompression increases, with only a minor improvement in the separation efficiency. Besides, it can be observed from raffinate recovery ratio (Fig. 9(b)) that the ex traction of water in the supercritical fluid phase increases, reducing selectivity of extraction and thus the IPA purity in the extract.

Influence of packing type

In order to evaluate the model performance with different packing characteristics, simulations were performed and compared with ex perimental results using a metallic foam packing, as reported by Lalam et al. The main characteristics of this foam are presented in Table 2.

Fig. 10 shows the experimental and numerical results for the frac tionation of a 5% IPA solution at 10 MPa and 40 °C and at different Q F / Q L ratios. The results with springs packing are also presented for comparison. Numerical results for the foam packing were calculated using α = 1.38 and β = 0.0018 as adjusted coefficients. As in the case of springs, a fair agreement was obtained for Q F /Q L > 6 (AARD of 33% for IPA concentration and recovery in the raffinate), with higher de viations for Q F /Q L < 6 (47%).

It is interesting to note that, for both types of packing, there seems to be a change in the experimental data slope around Q F /Q L = 6, which is not correctly reproduced by the model. As mentioned, this may be due to a change in the column hydrodynamic regime and inner flow patterns, which would affect the phase distribution, the effective in terfacial area as well as the mass transfer coefficient values. This sudden change of behaviour can also remind the conventional "loading point" of a packed bed as encountered in gas liquid columns. This analogy where V w , V o are the molar volumes of pure water and the pure organic solute, T is the absolute temperature, q is a parameter which depends on the type and size of the organic solute (for alcohols, it is equal to the number of carbon atoms in the molecule), and Ψ w , Ψ o are the superficial bulk volume fractions of water and organic solute, defined as:

= + ψ x V x V x V w w w w w o o (A.5) = + ψ x V x V x V o o o w w o o (A.6)
This set of equations can be solved provided the values of V w , V o , x w , x o , σ w , σ o , q and T, considering that:

+ = ψ ψ 1 w σ o σ (A.7)
In this work, in order to take into account the effect of scCO 2 , the surface tension of pure water (σ w ) was replaced by the value for water under CO 2 pressure within the experimental range covered here. More details can be found in the work by Tamura et al. [START_REF] Tamura | Practical method for estimating surface tensions of solutions[END_REF]. where M B and μ B are the solvent molar mass (g/mol) and viscosity (cP), V A is the molar volume of the solute at its normal boiling temperature (cm 3 / mol), T is the system temperature (K) and ϕ is the association factor of the solvent (ranging from 2.6 if the solvent is water to 1.0 if it is unassociated).

Fig. 2 .Fig. 3 .

 23 Fig. 2. IPA concentration in liquid (a) and fluid (b) outlet currents vs. time for the fractionation of a 5% IPA solution at T = 40 °C, P = 10 MPa, Q L = 0.48 kg/h and Q F /Q L = 9.5 kg/kg, using different initial liquid phase profiles: (--) pure water; (-) 5% IPA solution; (•••) linear profile.

Fig. 4 .

 4 Fig. 4. (a) IPA wt% in raffinate (in CO 2 -free basis) and (b) raffinate recovery ratio of IPA and water versus operation pressure for the fractionation of dilute IPA solutions (5 wt%, T = 40 °C, Q L = 0.48 kg/h, Q F /Q L = 6.04 kg/kg) using springs packing. Dots: experimental results from Lalam et al. [18]. Lines: simulation results with α = 1.05 and β = 0.001.

Fig. 5 .

 5 Fig. 5. Calculated IPA partition coefficient (dotted line) and mass transfer resistances (solid lines). Column average values as a function of operation pressure for the fractionation of dilute IPA solutions (5 wt%, T = 40 °C, Q L = 0.48 kg/h, Q F /Q L = 6.04 kg/kg) using springs packing.

Fig. 6 .

 6 Fig. 6. IPA partition coefficient profile along the column as a function of operation pressure for the fractionation of a 5% IPA solution at T = 40 °C and Q F /Q L = 9.5 kg/kg using springs packing.

Fig. 9 .

 9 Fig. 9. (a) IPA wt% in raffinate (in CO 2 -free basis) and (b) raffinate recovery ratio of IPA and water versus solvent-to-feed ratio for the fractionation of dilute IPA solutions (5 wt%, T = 40 °C, P = 10 MPa, Q L = 0.48 kg/h) using springs packing. Dots: experimental results from Lalam et al. [18]. Lines: simulation results with α = 1.05 and β = 0.001.

Fig. 10 .

 10 Fig. 10. (a) Raffinate recovery ratio of IPA and water and (b) IPA wt% in raffinate (in CO 2 -free basis) versus solvent-to-feed ratio for the fractionation of dilute IPA solutions (5 wt%, T = 40 °C, P = 10 MPa, Q L = 0.48 kg/h) using stainless steel springs (solid lines and black dots) and metallic foam packing (dotted lines and empty dots).

  Appendix B. Method of Wilke and Chang.This method provides an estimation of the diffusion coefficient of a solute (A) in a solvent (B) at infinite dilution conditions (D AB o )

Table 1

 1 Physico-chemical properties of pure compounds (from NIST a ).

	Property	IPA	water	CO 2
	Molar mass (g/mol)	60	18	44
	Critical temperature (°C)	235.9	373.6	31.1
	Critical pressure (MPa)	4.9	22.1	7.38
	Critical volume (cm 3 /mol)	222.0	55.9	91.9
	Molar volume at normal boiling point (cm 3 /mol)	77	18	-

NIST webbook of chemistry http://webbook.nist.gov/chemistry/.

a National Institute of Standards and Technology.

Table 2

 2 Packed column parameters used in the calculations.

	Parameter	Value
	Column length (m)	2
	Column internal diameter (m)	0.0175
	Packing dry void fraction	0.75 (springs)
		0.92 (foam)
	Packing characteristic length (m)	0.003 (springs)
		0.004 (foam)
	Packing specific surface area (m 2 /m 3 )	942 (springs)
		610 (foam)
	Stainless steel "critical" surface tension (N/m)	0.075

Table 3

 3 Fluid and liquid phase properties used in the calculations.

	40	10	630	5.00	996	6.54
	40	12	718	5.85	997	6.54
	40	15	780	6.75	999	6.54
	40	20	840	7.83	1001	6.55
	50	10	384	2.84	992	5.49
	60	10	290	2.38	987	4.69
	70	10	248	2.25	982	4.06
	80	10	222	2.20	976	3.57

T (°C) P (MPa) ρ F (kg/m 3 ) μ F × 10 5 (Pa s) ρ L (kg/m 3 ) μ L × 10 4 (Pa s)

Table 4

 4 Calculated IPA diffusion coefficient (D) and mass transfer coefficient (k) in fluid and liquid phase. T (°C) P (MPa) D F × 10 8 (m 2 /s) D L × 10 9 (m 2 /s) k F × 10 5 (m/s)

	k L × 10 6
	(m/s)