On the modeling and simulation of coupled adsorption and thermosolutal convection in supercritical carbon dioxide

In this paper, we present a numerical study along with an exhaustive adsorption investigation in a binary dilute mixture model nearby the solvent’s critical point in a configuration relevant for soil remediation. By means of this model, mass and heat transfer efficiency were qualitatively and quantitatively discussed through this work. The convergence of the solution was evaluated on the values of the Nusselt and Sherwood numbers. The results reveal intense convection expanding into the cavity close to the critical point, thus enabling homogeneous adsorption of the solute. Moreover, the mass fraction perturbation isolines exhibit the existence, along the adsorbent plate, of a thin boundary layer which becomes thinner when approaching the critical point.


Introduction
Supercritical fluid (SCF) stands for a fluid for which temperature and pressure are set above its gas-liquid critical point. In the supercritical state of matter, the fluid physical properties are intermediate between those of a liquid and a gas, making it a very good solvent. Moreover, SCF are very sensitive to temperature and pressure variations, thus a small change in pressure results in large variations of the fluid density as well as the solubility of solutes. More precisely, in the vicinity of the critical point, these tunable solvents exhibit better mass transfer properties compared to those of liquid solvents due to their higher diffusion coefficient and lower viscosity. SCF potential applications have known an increasing involvement and significant progress in a wide range of industrial technologies, for instance in various pharmaceutics sectors related to natural substances extraction of biological interest [1][2][3], agrofood fields [4][5][6][7], energy [8], and waste treatment sectors [9][10][11][12][13][14]. Supercritical carbon dioxide (scCO 2 ) significance has grown fairly quickly in the past decades. It has been extensively investigated in most SCF deposition (SFD) or crystal growth processes and has been the subject of several researches and publications [15][16][17][18][19]. In fact, SFD is a rising technique for the purpose of depositing metals and metaloxides on surfaces, and is notably relevant for technological processes aiming to prepare supported metallic nanoparticles or metallic films [16,17]. Moreover, SFD peculiarly employing scCO 2 is attractive since carbon dioxide offers many advantages by reason of economic and environmental concerns; it is non-flammable, non-toxic, mainly inexpensive, and remains an attractive solvent since CO 2 critical temperature is near ambient temperature. Additionally, scCO 2 has been employed in most of the depositions because no liquid waste is generated, no solvent residue is left on the substrate and the mass transfer rates are fast compared to those of liquids [17]. In addition, the rate and the velocity of the substrate deposition can be controlled by changing the pressure and temperature of the supercritical solvent.
Furthermore, SCF displays very interesting qualities in several potential processes including SCF extraction (SFE) for soil remediation associated with adsorption on activated carbon [20]. This modus operandi was primarily suitable for high molecular weight hydrocarbons [21]. Actually, numerous organics extraction such as that of hexachlorobenzene, phenol, and naphthalene from soil employing scCO 2 was reported in [21][22][23][24]. However, due to environmental legislation, low contaminant concentration is compulsory. Therefore, an almost contaminant-free process is required when the supercritical solvent is recirculated. As, for volatile pollutants, the solvent and solute disjunction by means of decompression remains inefficacious to attain the requested low concentration criterion, an adsorption supplemental step onto activated carbon remains necessary. In this context, Madras et al. [24] proposed a process coupling supercritical extraction and adsorption for the remediation of soils contaminated by heavy molecular weight organic compounds. In fact, organic pollutants are constantly extracted through scCO 2 before being adsorbed on activated carbon. It was established that if the soil desorption is implemented at low temperature, then the process operates in optimum mode if desorption is followed by adsorption at high temperature [24].
Although the widest application of scCO 2 is in extraction, the adsorption process using scCO 2 exhibits currently more attention and presents a number of potential applications. In fact, several studies investigated the various supercritical adsorption singularities of many systems [12][13][14][25][26][27][28][29][30][31][32][33][34][35], thus allowing a selective transfer of species from fluid phase to a solid adsorbent. Moreover, Jha and Madras [31] showed that the heat of adsorption does not depend on the solute considered. When an adsorbent is exposed to a fluid phase, the molecules of the fluid diffuse on its surface. Isothermal adsorption of a particular system provides a source of information on adsorbent pore properties and surface characteristics. The reversible type isotherm is obtained when the adsorption is limited to a few molecular layers. This condition is usually substantiated for adsorption on microporous materials endowed with small-scale external surfaces such as activated carbons. However, SCF adsorption on surfaces or porous particle inner surfaces is much more complex. In fact, in a very close neighborhood of the critical point, the adsorption equilibrium constant is extremely sensitive to temperature and pressure changes leading, especially, to a decrease of the excess amount adsorbed when pressure is increased [32].
In order to devise an adsorption plant, reliable mass transfer models are still necessary. Notwithstanding, exiguous tryouts have been implemented in order to model the complete supercritical adsorption process, including hydrodynamics [31][32][33][34][35]. That is the reason why this paper aims to investigate adsorption and mass transfer in a binary mixture nearby the solvent's critical point in a configuration relevant for soil remediation. In order to take advantage of the SCF accurate properties and enhance thermodynamics conditions enabling a meliorate adsorption along with a good mass transfer, the influence of critical point vicinity is investigated using numerical simulations in the present work.

Methods
The physical model consists of a dilute mixture of naphthalene (the minority species) and supercritical CO 2 , confined in a rectangular cavity of aspect ratio L/H = 10 with a heated and adsorbent plate on a part of the bottom boundary as shown in Fig. 1. Naphthalene is the species adsorbed on the bottom plate. It was selected as a model solute because of the wide experimental data of its solubility in scCO 2 available in the literature, allowing an easier validation of the thermodynamic model [15]. Moreover, the naphthalene-CO 2 mixture can be considered as a reference system, representative of binary mixtures of interest for SFE and adsorption. The phase diagram of this kind of mixture reveals a discontinuous mixture critical line (Fig. 2). The low-temperature part ceases at the LCEP (lower critical endpoint) that matches with the junction between the solid-liquid-vapor coexistence line and the mixture critical line beginning at the critical point of the solvent (here CO 2 ). The LCEP is usually slightly above the solvent's critical point. In the same way, the high-temperature part starts at the critical point of the solute (naphthalene) and ends at the UCEP (upper critical endpoint). For temperatures ranging between T LCEP and T UCEP , the mixture exhibits equilibrium between a solute solid phase and a mixture fluid phase. In the vicinity of the LCEP, solid solubility and adsorption equilibrium constants are highly sensitive to temperature and pressure variations [15,32], which points out the reason why SFE and adsorption processes are customarily carried out in the aforementioned phase diagram region. Moreover, near the gas-liquid critical point of a pure compound, fast thermal equilibration is achieved by a thermoacoustic effect, named the piston effect, while thermal diffusion becomes extremely slow [36][37][38]. More precisely, when the fluid is locally heated, a very thin thermal boundary layer forms along the heated boundary because of the vanishing thermal diffusivity. Since the thermal expansion coefficient diverges near the critical point, this boundary layer strongly expands and adiabatically compresses the bulk fluid leading to a fast and homogenous increase of pressure and bulk temperature. The same phenomenon occurs in mixtures near the solvent's critical point and it was shown in [15,32] that the strong pressure increase induced by the piston effect has important consequences on the solubility of solids and on adsorption equilibrium constants near the LCEP. The critical properties of naphthalene and CO 2 together with those of the mixture at the LCEP are reported in Table 1. Initially, at t = 0, the fluid is assumed in thermodynamic equilibrium at a uniform temperature T i slightly higher than the temperature  Table 1 Critical parameters, acentric factor ω, and molar mass M for the pure compounds and the mixture of the peculiar mixture critical point LCEP, T cm = 307.65 K. Thus, T i is written as T i = (1 + ε)T cm , where the parameter ε = (T i − T cm )/T cm defines the dimensionless proximity to the mixture critical point. The mixture is saturated, i.e., the initial mass fraction w i corresponds to the solubility of naphthalene in scCO 2 at temperature T i . Additionally, the fluid is stratified along with a mean density equal to the mixture critical density ρ cm = 470 kg m −3 . Since ρ cm is very close to the CO 2 critical density, the very strong variations of thermophysical properties of CO 2 on the critical isochore are still observed. Weak heating of the order of a few mK per second is then progressively applied at the solid bottom plate corresponding to y = 0 and for x extending between 1 and L/H − 1 (see Fig. 1). In order to avoid discontinuities between isothermal and heated parts of the bottom boundary, the following hyperbolic tangent function is used: Thus, for ΔT fixed equal to 10 mK for all the simulations, the bottom plate temperature is gradually raised according to the formula: We consider for the mathematical model the 2D compressible time-dependent Navier-Stokes equations, along with the energy and mass diffusion equations together with the Peng-Robinson equation of state. The Peng-Robinson equation implicitly accounts for the divergence of the thermal expansion coefficient β, of the isothermal compressibility χ, and of the specific heat at constant pressure C P near the liquid-gas critical point of CO 2 . Carbon dioxide and naphthalene are respectively named components 1 and 2. The governing equations are solved in the low Mach number approximation framework [39] allowing computational cost lessening while accounting for density variations due to compressible effects. In this context, the pressure P is split into two parts: a homogeneous time-dependent thermodynamic part P th (t) looming up in the state and energy equations, and a non-homogeneous dynamic part P dyn (x, y, t) appearing in the momentum equations and depending on time and space. We employed the low Mach number approximation amendment proposed by Accary et al [40] near the critical point in order to account for the strong fluid stratification.
Moreover, we designate by T cm , ρ cm , and H respectively the characteristics temperature, density, and length. Additionally, we consider the piston effect time scale [41] t PE = t d /(γ − 1) 2 as being the characteristic time, where t d and γ are respectively the characteristic time of thermal diffusion and the mixture specific heat ratio (computed from the equation of state [15]), while the corresponding velocity V PE = H/t PE denotes the characteristic velocity.
Since we are dealing with dilute mixtures with CO 2 as the majority component, the transport properties such as the dynamic viscosity μ, the specific heat at constant volume C V , and the thermal conductivity λ, are assumed to be equal to those of pure CO 2 . As the wall heating considered in this study is small and the divergence of C V and μ near the critical point is weak, these two properties are supposed to be constant throughout the simulation and equal to their value for the initial condition. The values provided by the NIST are then used. On the other hand, the correlation of Arai et al. [42] is used for the variation of pure CO 2 thermal conductivity λ. The binary mass diffusion coefficient D 21 is modeled with the Wilke-Chang equation [43]. All the transport properties are made dimensionless relatively to their initial value. Thus, the dimensionless governing equations are written: where V denotes the velocity of components u and v respectively in the x and y directions, e y is the unity vector of the Cartesian orthonormal basis, w is the naphthalene mass fraction, γ 0 and C v0 denote the corresponding values for a perfect gas (γ 0 = 1.3, C v0 = 3R/M 1 ), and λ * and D Ã 21 correspond respectively to the dimensionless conductivity and mass diffusion coefficient.
In the above equations, the dimensionless parameters stand for the Reynolds number Re, the Froude number Fr, the Prandtl number Pr, and the Lewis number Le respectively defined by: Two other dimensionless parameters are of interest for the present problem: the Mach and the Rayleigh numbers respectively defined by: denotes the sound speed with R = 8.3145 J mol −1 K −1 the perfect gas constant. In Eq. (5), U Ã k and V Ã k denote respectively the dimensionless partial molar internal energy and volume expressed as follows: where U k and V k represent respectively the partial molar internal energy and volume, computed from the Peng-Robinson equation [15].
The mixture thermodynamic state is described by the Peng-Robinson equation of state in the framework of the one-fluid theory: where where ω is the acentric factor (see Table 1) and the superscript (*) refers to dimensionless quantities. The binary interaction parameters k 12 and l 12 were fitted on experimental solubility data and depend on initial temperature as [15]: k 12 ¼ 0:0395 þ 0:0114ð 308: 15 T i −1Þ and l 12 ¼ −0:1136−0:3103ð 308:15 T i −1Þ On the adsorbent plate, the adsorption reaction is modeled by a first-order reaction. As for temperature, the discontinuity between the non-reactive parts and the adsorbent plate is prevented by the use of the function Φ (Eq. (1)). The boundary condition for w on the bottom wall is then: where K a is the adsorption constant depending on temperature and pressure. In view of small temperature differences ΔT considered, K a is linearized as follows [32]: The derivatives (∂K a /∂T) P and (∂K a /∂P) T depend on those of the adsorption equilibrium constant that diverge near the solvent's critical point [32]. It is notable that this divergent behavior is not specific to naphthalene but it is a universal behavior for dilute mixtures since it is due, in particular, to the divergence of the solute infinite dilution partial molar volume [44]. In dimensionless form, the boundary condition (8) is written as: The mathematical model is completed by the following boundary conditions (Fig. 1): -No-slip conditions are prescribed for velocity on all the boundaries; -A homogeneous Neumann boundary condition is applied for w at the non-reactive vertical and upper walls; -The vertical walls are adiabatic while the horizontal ones are kept at the initial temperature, except for the heated plate where a Dirichlet boundary condition is used for T: The code employed in this study is homemade [32][33][34]. The numerical method is based on a second-order semi-implicit scheme for time discretization and a spectral collocation method for space approximation to ensure a high accuracy (see [12]). The Gauss Lobatto collocation points are used. The mesh is therefore naturally refined near the walls, allowing a better representation of the thin thermal boundary layers. The refined mesh size near the walls is in the order of 3 × 10 −4 mm (about 20% of grid points inside the boundary layers). The convergence of the solution in terms of mesh size was evaluated on the values of the Nusselt and Sherwood numbers. It was found that a mesh 251 × 81 was sufficient.

Results and discussion
In the present work, numerical simulations were carried out for three initial temperatures T i = 308.15 K, 309.15 K, and 318.15 K. These various temperatures correspond to distances to the mixture critical point equal to 0.5 K, 1.5 K, and 10.5 K respectively. The initial mass fraction w i dovetails with the solubility of naphthalene in CO 2 , which means that the mixture is saturated. The transport properties values obtained for these initial conditions are reported in Table 2.
While the critical temperature is nearing, we perceive the strong increase of the thermal expansion coefficient β, the thermal conductivity λ along with the specific heat ratio γ. These behaviors lead to the strong Rayleigh, Reynolds, Froude, and Prandtl numbers increment as T i is decreased from 318.15 to 308.15 K (see Table 3). Moreover, it must be highlighted that the Damköhler number Da was set at 10 −5 for all the simulations, which is a relevant value for adsorption of Naphthalene. Indeed, from available literature data [12][13][14]30], the Damköhler number for adsorption of naphthalene on various solids can be estimated in the range 10 −14 -10 −3 . However, it was shown in [32] that the phenomena observed for the adsorption reaction near the LCEP are qualitatively the same whatever the value of Da. The Damköhler number measures the relative Table 2 Transport coefficients corresponding to various temperature values T i importance of adsorption rate compared to diffusion rate. Therefore, when Da ≪ 1, mass diffusion is much faster than adsorption so that it reaches equilibrium before adsorption.

Stability and convergence to the steady state
The values of different characteristic times for each simulation case are listed in Table  4. We can note the strong increase of the thermal diffusion characteristic time and the strong decrease of the piston effect time scale as the critical temperature is approached.
On the other hand, the mass diffusion characteristic time does not change very much with temperature. The solution state (stationary, periodic, chaotic) is evaluated on the temporal evolution of the temperature perturbation T − T i at x = L/2 and in the hot boundary layer (Fig. 3a) and in the middle of the cavity (Fig. 3b) for the three initial temperatures. Qualitatively, the different curves of Fig. 3 indicate that a stationary solution is obtained. The temperature in the hot boundary layer stabilizes faster for the three initial conditions (around t = 60 s). In the middle of the cavity, the stabilization time increases as the initial temperature decreases to the critical point value (Fig. 3b): around 62 s, 92 s, and 99 s for T i = 318.15 K, 309.15 K, and 308.15 K, respectively. These values are globally consistent with the thermal diffusion characteristic times reported in Table 4. Obviously, convection also influences the temperature homogenization inside the cavity and Table 3 shows that the Rayleigh number strongly increases when T i is closer to the critical temperature. Therefore, convection is much more intense at T i = 308.15 K than at T i = 318.15 K, as discussed in section 5.2. Figure 4 presents the vertical profiles of temperature at x = L/2 for different times and for the three initial conditions. A completely different behavior is observed far and close to the mixture critical point. For T i = 318.15 K, a diffusive profile is almost achieved at t = 6 s and the final increase of the bulk temperature is then induced by convection. For the two smaller initial temperatures, the final temperature profile still exhibits a boundary layer-like variation. Moreover, the temperature profiles reveal that a time larger than 100 s is required to reach convergence. That is the reason why, for the rest of this paper, all the results are represented at a time greater than this reference convergence time.  Table 4 Piston effect, t PE , thermal diffusion, t d , and mass diffusion, t md , characteristic times for the three initial temperatures  Figure 5 shows the temporal evolution of the thermodynamic pressure for the three initial temperatures. In all the cases, we observe the strong pressure increase, typical of the piston effect, from the beginning of the simulation. We note also that the pressure level that reached convergence is higher as the critical temperature is approached, corresponding to the divergence of isothermal compressibility near the critical point. This larger pressure increase obviously has an effect on the adsorption as shown by Eq. (9).

Heat and mass transfer investigation in the vicinity of the lower critical endpoint
We are interested in the following section on the effect of the vicinity to the mixture critical point on temperature and mass fraction fields inside the cavity. Figure 6 displays  the temperature isolines and mass fraction perturbation at two different times corresponding to 2.5t PE and 430t PE for the three initial conditions. The solutions obtained for the various values of T i are similar. The temperature field at this time is governed by the piston effect, as in pure CO 2 [41], with three distinct zones: a hot boundary layer along the heated plate, cold boundary layers along isothermal walls, and the cavity bulk that is homogeneously heated by the piston effect. In addition, as in pure CO 2 , the thermal boundary layers are thinner when approaching the critical point because of the strong decrease of thermal diffusivity. In the same way, since Da ≪ 1, the mass fraction near the adsorbent plate is governed by mass diffusion. Indeed, the isolines of the mass fraction perturbation reveal the existence, along the adsorbent plate, of a thin boundary layer that becomes thinner when approaching the critical point because of the decrease of the mass diffusion coefficient D 21 (Table 2). Figure 6 also presents the steady solution obtained for each initial condition. Concerning the mass fraction, it must be underlined that, although the isoline map does not evolve anymore, the values of w i − w keeps on increasing while the adsorption reaction still goes on. Indeed, since Da ≪ 1, the reaction will reach equilibrium much later than mass diffusion. As expected, the temperature field is very different depending on the proximity to the critical point. Because of the diverging behavior of the thermal expansion coefficient β (Table 2), the Rayleigh number strongly increases as approaching the critical temperature (Table 3): Ra = 3683 for T i = 318.15 K, Ra = 21938 for T i = 309.15 K, and Ra = 34139 for T i = 308.15 K. Consequently, the convective flow is much more intense for T i = 308.15 K (Fig. 6), allowing a better mixing inside the cavity. We have calculated the flow enstrophy defined by: This physical quantity, directly related to the norm of the vorticity ∇ × V, is a kind of measure of the vortices intensity inside the cavity. It was found that E strongly increases as the critical point is approached: E = 3 for T i = 318.15 K, E = 14 for T i = 309.15 K and E = 16 for T i = 308.15 K. Therefore, the enstrophy is more than 5 times larger at 308.15 K than at 318.15 K, which clearly represents a substantial increase. The distribution and the intensity of the convection cells directly affect the mass fraction field as revealed by Fig. 7 showing the vorticity field together with the isolines of w i − w. The jet flow between contrarotative rolls is much weaker for T i = 318.15 K allowing the development of a large mass fraction "plume" rising at the middle of the heated plate (Fig. 6). For T i = 308.15 K, in contrast, the very intense convection flow restricts the diffusion of the solute mainly in the boundary layer along the plate. The mass fraction field obtained for T i = 309.15 K is intermediate between those previously described (Fig. 6).
The effect of the proximity to the mixture critical point on heat and mass transfer was also evaluated on the mean Nusselt, Nu, and Sherwood, Sh, numbers on the heated plate, respectively defined by: with L plate = 8 the dimensionless length of the heated plate. Results are reported in Table 5 for the three initial temperatures. We note a strong increase of both Nu and Sh as the critical temperature is approached (+ 86% for Nu, + 46% for Sh), which means that the neighborhood of the critical point favors heat and mass transfer from the heated plate. In [45] As shown in Fig. 8, the values of the mean Sherwood number that we obtained are consistent with Eq. (12) (with m 3 = 0.182) in the vicinity to the critical point. On the other hand, for the highest initial temperature, the Sherwood number is strongly underestimated by the correlation. Adsorption investigation in the vicinity of the lower critical endpoint The adsorbed solute amount is greatly affected by the convection flow. Actually, Fig. 9a exhibits the relative mass fraction perturbation profile (w i − w)/w i on the lower boundary for the three initial temperatures. Although the maximum amount adsorbed on the plate increases when moving away from the critical point, the profiles disclose that the adsorption of the solute is strongly inhomogeneous for T i = 318.15 K, with a peak of adsorption in the middle of the plate, corresponding to the large mass fraction "plume" noted in Fig. 6. On the other hand, close to the critical point (T i = 308.15 K), the peak in the central zone is less pronounced, leading to a more homogeneous adsorption of the solute and, therefore, a more homogeneous material in the framework of SFD processes for example. The same tendency with respect to the distance to the critical point is observed on the temporal evolution of the relative mass fraction perturbation (w i − w)/w i (Fig. 9b).  As already outlined in section 5.2, since Da ≪ 1, the reaction has not reached equilibrium although the mass fraction field in the cavity shows a stationary solution (Fig. 6), and, consequently, the mass fraction perturbation on the bottom plate keeps on increasing.
In order to quantify the effect of convection on the adsorbed amount, we have performed simulations in the same configuration with no gravity. The mean adsorbed amount of solute, defined by: and obtained with and without gravity is reported in Table 6. First, it can be noted that, even in the absence of convection, the adsorbed amount of solute decreases when approaching the critical point (− 1.35% at T i = 308.15 K compared to the value obtained at T i = 318.15 K). This decrease is due to the increase near the critical point of the pressure rise induced by the piston effect (Fig. 5) and the negative value of the adsorption constant derivative with respect to pressure, leading to a decrease of the adsorption constant K a (Eq. (9)). However, Table 6 also shows that this decrease of q mean with the distance to the critical point is much larger in the presence of convection (− 15% at T i = 308.15 K compared to the value obtained at T i = 318.15 K). This stronger decrease is directly caused by the intense convection flow observed for the smallest initial temperature (Fig. 7). Indeed, the adsorbed amount also depends on the contact time of the solute with the adsorbent plate: the longer the solute is in contact with the adsorbent surface, the greater is the adsorbed amount.
In order to study the contact time of the solute molecules with the adsorbent plate, we represent in Fig. 10 the x profile of the vertical velocity v in the boundary layer along the bottom boundary for the three initial temperatures. These profiles are particularly interesting because they show that the smaller the temperature T i is, the greater the vertical velocity v is and, consequently, the faster the solute molecules are removed from the vicinity of the adsorbent wall. Therefore, Fig. 10 clearly shows that the contact time of the solute decreases with the initial temperature. Moreover, the vertical velocity v for T i = 308.15 K is 6 times larger than for T i = 318.15 K. This large vertical velocity is directly related to the intense convection cells depicted by the vorticity field in Fig. 7. Therefore, we can conclude that, because of the increasing Rayleigh number near the critical point (induced by the divergence of the thermal expansion coefficient), the contact time of the solute with the adsorbent plate is strongly reduced and, consequently, the adsorbed amount of solute is decreased when the critical temperature is approached. Table 6 Mean adsorbed amount of solute, q mean , at t = 128s for the three initial conditions with (gravity) and without (g = 0) convection