Qifan Ying

Hui Wang
email: huiwang@dhu.edu.cn

Eric Lichtfouse

Numerical simulation on thermal behavior of partially filled metal foam composite phase change materials

Keywords:

As a kind of material with high thermal conductivity, metal foam embedded in phase change materials (PCM) can significantly improve the thermal energy storage performance of PCM, so as to better serve the latent heat thermal energy storage (LHTES). However, the presence of metal foam will suppress the natural convection of PCM, reduce its LHTES capacity and increase the total cost. In view of the above problems, this study designed LHTES units with metal foam filled at the top or the bottom with different height ratios. The melting heat transfer characteristic was studied by numerical simulations. Results show that the flow rate of liquid paraffin in the bottom filling configuration is one to two orders of magnitude higher than that in the top filling configuration, which plays a positive role in natural convection and is more conducive to heat transfer enhancement. However, the top filling configuration performs better in the heat storage capacity. The metal foam with a filling height ratio of 0.75 can increase the heat storage capacity by 7.47% while saving 25% material. On this basis, we proposed a new criterion, named the comprehensive strengthening goodness factor (CSGF), for evaluating the heat transfer performance of LHTES system. CSGF is different from the general criterion in that it can evaluate the heat storage capacity and the economic performance at the same time. Under the evaluation of CSGF, the optimum metal foam-paraffin configurations for different occasions were selected.

Introduction

Nowadays, saving energy consumption and exploring more efficient ways of energy utilization have become research hotspots [START_REF] Bui | Enhancing building energy efficiency by adaptive façade: a computational optimization approach[END_REF][START_REF] Li | Applications and technological challenges for heat recovery, storage and utilisation with latent thermal energy storage[END_REF]. Renewable energy resources, such as solar energy, is a popular way to store thermal energy. In thermal energy storage systems, the focus is on latent heat thermal energy storage (LHTES). This can be attributed to the fact that phase-change materials (PCM) being employed as heat storage media in the LHTES system, which can achieve a large amount of thermal energy storage within a small temperature difference [START_REF] Gao | Application and research progress of phase change energy storage in new energy utilization[END_REF]. In solar energy storage system, PCM is usually used to improve its intermittency and uneven spatial distribution [START_REF] Javadi | Performance improvement of solar thermal systems integrated with phase change materials (PCM), a review[END_REF][START_REF] Wu | A novel composite PCM for seasonal thermal energy storage of solar water heating system[END_REF]. In the electronic component thermal management system, it is used to absorb a large amount of heat emanating from electronic components and prevent high temperature damage [START_REF] Ren | Thermal management of electronic devices using pin-fin based cascade microencapsulated PCM/expanded graphite composite[END_REF][START_REF] Amin | Experimental investigation of multiwall carbon nanotube/paraffin based heat sink for electronic device thermal management[END_REF].

The common types of PCM are generally organic, inorganic and eutectic. Many researchers have found that organic PCM has stable chemical and thermal properties, and compared with inorganic PCM, it is non-corrosive and has no phase segregation and super cooling phenomenon [START_REF] Heyhat | Battery thermal management with thermal energy storage composites of PCM, metal foam, fin and nanoparticle[END_REF], which is a promising energy storage material in the future. However, due to the low thermal conductivity, the thermal energy storage effect is not up to expectations, so the application is limited. At present, methods to enhance the heat transfer of PCM mainly rely on external additives to strengthen, such as adding fins [START_REF] Sciacovelli | Maximization of performance of a PCM latent heat storage system with innovative fins[END_REF][START_REF] Ji | Non-uniform heat transfer suppression to enhance PCM melting by angled fins[END_REF][START_REF] Yao | Study on solidification performance of PCM by longitudinal triangular fins in a triplex-tube thermal energy storage system[END_REF][START_REF] Arshad | Towards the thermal management of electronic devices: a parametric investigation of finned heat sink filled with PCM[END_REF], Nanoparticles [START_REF] Bashar | Experimental investigation of transient melting and heat transfer behavior of nanoparticle-enriched PCM in a rectangular enclosure[END_REF][START_REF] Al-Jethelah | Charging nanoparticle enhanced bio-based PCM in open cell metallic foams: an experimental investigation[END_REF][START_REF] Kumar | Experimental and theoretical investigations on thermal conductivity of the paraffin wax using CuO nanoparticles[END_REF][START_REF] Wang | Approaches for expedition of discharging of PCM involving nanoparticles and radial fins [J][END_REF], graphene [START_REF] Mehrali | Shape-stabilized phase change materials with high thermal conductivity based on paraffin/graphene oxide composite[END_REF][START_REF] Yang | Hybrid graphene aerogels/phase change material composites: thermal conductivity, shape-stabilization and light-to-thermal energy storage[END_REF][START_REF] Zhu | Graphene/SiO2/n-octadecane nanoencapsulated phase change material with flower like morphology, high thermal conductivity, and suppressed supercooling[END_REF][START_REF] Li | Thermal conductivity enhancement and shape stabilization of phase-change materials using three-dimensional graphene and graphene powder[END_REF], and metal foam [START_REF] Qu | Numerical model of the passive thermal management system for high-power lithium ion battery by using porous metal foam saturated with phase change material[END_REF][START_REF] Yao | A new prediction model for the effective thermal conductivity of high porosity open-cell metal foams[END_REF][START_REF] Wang | Theoretical and experimental advances on heat transfer and flow characteristics of metal foams[END_REF][START_REF] Zhang | Numerical simulation of flow and heat transfer in composite PCM on the basis of two different models of open-cell metal foam skeletons[END_REF]. In these strengthening methods, metal foam is concerned. As summarized in Fig. 1, as a porous medium material, it has many obvious advantages such as small density, high surface area, high porosity and high thermal conductivity. This special skeleton interconnected network structure makes its heat transfer performance excellent. It combines well with PCM for wettability, and this metal foam composite phase change material has also been demonstrated to greatly enhance its heat transfer capability [START_REF] Tian | A numerical investigation of heat transfer in phase change materials (PCMs) embedded in porous metals[END_REF][START_REF] Zhu Z-Q, Huang Y-K, Hu N | Transient performance of a PCM-based heat sink with a partially filled metal foam: effects of the filling height ratio[END_REF][START_REF] Chen | Thermal performance of a pcm-based thermal energy storage with metal foam enhancement[END_REF][START_REF] Tang | Pore-scale numerical simulation of the heat transfer and fluid flow characteristics in metal foam under high Reynolds numbers based on tetrakaidecahedron model[END_REF][START_REF] Yang | Influence of aspect ratios for a tilted cavity on the melting heat transfer of phase change materials embedded in metal foam[END_REF][START_REF] Zhang | Experimental study on heat transfer characteristics of metal foam/paraffin composite PCMs in large cavities: Effects of material types and heating configurations[END_REF].

However, with the further research, it is found that metal foam can effectively enhance the melting heat transfer of PCM, but also bring some negative effects. Tao et al. [START_REF] Tao | Lattice Boltzmann simulation on phase change heat transfer in metal foams/paraffin composite phase change material[END_REF] numerically studied the natural convection during the melting period of metal foam-paraffin composite materials. Results showed that the natural convection of paraffin was suppressed due to the obstruction of porous media, and this effect was enhanced at high pores per inch (PPI) with porosity greater than 0.9, which worsened the heat transfer. Lafdi et al. [START_REF] Lafdi | Experimental study on the influence of foam porosity and pore size on the melting of phase change materials[END_REF] pointed out that natural convection is beneficial to improve the uniformity of the temperature distribution of the composites, and the composites with high porosity and large pore size metal foam can reach the steady-state temperature faster, but the heat conduction decreases at the same time. Therefore, it is necessary to measure the strength of heat conduction and natural convection. Dinesh et al. [START_REF] Dinesh | Effect of foam geometry on heat absorption characteristics of PCM-metal foam composite thermal energy storage systems[END_REF] numerically simulated the total heat absorption and latent heat of composite materials (metal foams with different porosity), and found that the metal foams with low porosity enhanced the thermal conductivity of composite materials due to the presence of more metals, but inhibited the natural convection of paraffin in the pores and reduced the total heat absorption of paraffin. Mahdi et al. [START_REF] Mahdi | Multiple-segment metal foam application in the shell-andtube PCM thermal energy storage system[END_REF] recommended that metal foams with porosity greater than 0.9 in combination with paraffin wax, noting that metal foams with too low porosity will significantly reduce the volume of paraffin wax, resulting in less heat storage. These similar findings have been replicated in other studies [START_REF] Mesalhy | Numerical study for enhancing the thermal conductivity of phase change material (PCM) storage using high thermal conductivity porous matrix[END_REF][START_REF] Caliano | Analysis of a phase change material-based unit and of an aluminum foam/phase change material composite-based unit for cold thermal energy storage by numerical simulation[END_REF]. It is known from previous studies that fully filled with metal foam will cause the natural convection and heat storage of PCM to weaken, and increase the cost. This makes partial filling an ideal way to allow more space for natural convection and heat storage, as well as lower costs. The excellent performance of finite metal foam filling is also confirmed [START_REF] Zhu Z-Q, Huang Y-K, Hu N | Transient performance of a PCM-based heat sink with a partially filled metal foam: effects of the filling height ratio[END_REF][START_REF] Zhang | Role of partial and gradient filling strategies of copper foam on latent thermal energy storage: an experimental study[END_REF][START_REF] Xu | Evaluation and optimization of melting performance for a latent heat thermal energy storage unit partially filled with porous media[END_REF].

Actually, the evaluation of the performance of this composite is based on both heat storage and material cost. However, the current study discusses the two separately. Joshi et al. [START_REF] Joshi | Thermal performance augmentation of metal foam infused phase change material using a partial filling strategy: an evaluation for fill height ratio and porosity[END_REF] analyzed the influence of copper foam with different filling heights on the heat storage

γ thermal expansion coefficient (K -1 ) ε porosity λ thermal conductivity (W m -1 K -1 )
μ dynamic viscosity (kg⋅m [42] carried out numerical simulation of ice storage spheres (metal foam-water) in the ice storage system, and used the above criteria (w, Pc) to determine the porosity of copper foam that makes the optimal cold storage performance and economic performance, respectively. Zuo et al. [START_REF] Zuo | Numerical investigation and optimal design of partially filled sectorial metal foam configuration in horizontal latent heat storage unit[END_REF] also used dimensionless Pc' to investigate the economy of composite materials, and modified it under their research model. However, to the best of our knowledge, there is no evaluation criterion that considers both heat storage and economic performance.

Inspired by the finite metal foam filling strategy, LHTES units with two types of filling forms were designed in this study, top filling and bottom filling. Numerical simulation was carried out for different filling heights of metal foams under the two types of forms. The influence of the configuration on natural convection and heat storage was fully considered, and the performance evaluation of LHTES units was mainly studied. A new criterion was proposed to evaluate the heat storage and economic performance comprehensively. Based on which the optimal filling configuration and filling height ratio was obtained to provide reference for practical applications.

Physical model and numerical approach

Physical model

Fig. 2 shows the physical models of different filling configurations in this study. The size of the model (H × W) is 90 mm × 60 mm, in twodimensional. The paraffin wax is used as the heat storage materials and copper foam is used for the enhanced heat transfer material. The two key parameters of copper foam are porosity and PPI. Porosity is defined as the ratio of the volume of pores in the metal foam to the total volume, which is taken as 90% in this paper. PPI is the number of pores per inch, and 10PPI is used in this paper. The thermophysical properties of paraffin wax and copper foam are listed in Table 1.

The filling ratio ξ is defined as the height of metal foam divided by the total height of the model. The left wall of the model is heated at a constant temperature, 423.15 K, other walls are adiabatic. Fig. 2(a) is the case without metal foam. Fig. 2(b) is the case fully filled with metal foam. Fig. 2(c) -(e) represents the cases partially filled with metal foam from the top and the filling ratio is 0.25, 0.5 and 0.75, respectively. Fig. 2(f) -(h) represents the cases partially filled with metal foam from the bottom. The filling ratio is also set to 0.25, 0.5 and 0.75, respectively. Fig. 3 shows the locations of the four the monitoring points. The distances of O 1 , O 2 , O 3 and O 4 from the bottom are 11.25, 33.75, 56.25 and 78.75 mm.

Equations

The enthalpy porosity approach and the volume averaging approach are employed to construct the numerical model presented in this study based on the following assumptions:

(1) The liquid PCM is considered incompressible and moves laminar, while its buoyancy effect follows the Boussinesq approximation. (2) PCM and metal foam are isotropic.

(3) The thermophysical properties of PCM and metal foam are constant with various time, and they maintain local thermal equilibrium.

Based on the above premise, the governing equations can be written performance of paraffin wax, and chosen the optimal filling strategy of copper foam based on the standard of reaching 90% of the heat storage of pure paraffin unit. Wang et al. [START_REF] Wang | Effect of copper metal foam proportion on heat transfer enhancement in the melting process of phase change materials[END_REF] added metal foam of different mass into paraffin wax, and found that the heat storage and heat storage rate of composite materials showed opposite trends with the increase of the filling rate of metal foam. Taking these two factors into consideration, the filling rate at the intersection of these two curves is chosen as the optimal metal foam content. Xu et al. [START_REF] Xu | Evaluation and optimization of melting performance for a latent heat thermal energy storage unit partially filled with porous media[END_REF] proposed a heat storage performance evaluation index w for metal foam-paraffin composite units with different structures, which was defined as heat storage per unit mass per unit time. In the follow-up study, the economic performance evaluation index P c is further proposed, which is defined as the heat storage of composite material per unit cost per unit time [START_REF] Xu | Melting performance enhancement of phase change material by a limited amount of metal foam: Configurational optimization and economic assessment[END_REF]. Lou et al.

as follows [START_REF] Tian | A numerical investigation of heat transfer in phase change materials (PCMs) embedded in porous metals[END_REF]: Continuity equation:

∂ρ pcm ∂t + ρ pcm ( ∂u ∂x + ∂v ∂y ) = 0 (1)
Momentum equations:

ρ pcm ε ( ∂u ∂t + 1 ε ( ∂(uu) ∂x + ∂(uv) ∂y ) ) = - ∂p ∂x + μ pcm ε ( ∂ 2 u ∂x 2 + ∂ 2 u ∂y 2 ) + S x (2)
ρ pcm ε

( ∂u ∂t + 1 ε ( ∂(vu) ∂x + ∂(vv) ∂y ) ) = - ∂p ∂y + μ pcm ε ( ∂ 2 v ∂x 2 + ∂ 2 v ∂y 2 ) + S y (3)
S x and S y in Eqs. ( 2) and ( 3) are defined as momentum source terms in × and y directions, which can be expressed by the following formula:

S x = - (1 -f ) 2 f 3 + δ A m u - ( μ pcm K + ρ pcm C i ̅̅̅ ̅ K √ |u| ) u (4) S y = - (1 -f ) 2 f 3 + δ A m v - ( μ pcm K + ρ pcm C i ̅̅̅ ̅ K √ |u| ) v + ρ pcm gf (T pcm -T 0 ) (5) f = ⎧ ⎨ ⎩ 0T ≤ T m1 (T -T m1 )/(T m2 -T m1 )T m1 ≤ T ≤ T m2 1T m2 ≤ T (6)
The first term of S x and S y represents the flow resistance of PCM in the mushy zone when it is melted. f is the liquid fraction in the paraffin melting process, which is used to characterize the degree of melting process, and is a constant between 0 and 1. The higher the value, the more paraffin has been liquefied. 0 and 1 represent the completely solid and liquid cases, respectively [START_REF] Esapour | Melting and solidification of PCM embedded in porous metal foam in horizontal multi-tube heat storage system[END_REF]. δ is equal to 0.001 to avoid a zero denominator [START_REF] Shatikian | Numerical investigation of a PCM-based heat sink with internal fins[END_REF]. In addition, the coefficient A m of the mushy zone is generally set to 10 5 [START_REF] Zheng | Thermal performance of copper foam/paraffin composite phase change material[END_REF][START_REF] Chamkha | Phase change heat transfer in an L-shape heatsink occupied with paraffin-copper metal foam[END_REF][START_REF]Numerical study on melting of paraffin wax with Al2O3 in a square enclosure -ScienceDirect[END_REF]. T m1 and T m2 indicate the temperature of PCM in the solid and liquid state, respectively.

The remaining terms of S x and S y represent the changes in viscosity, inertia resistance and buoyancy caused by the metal foam. The parameters related to metal foam are calculated by semi-empirical correlations. Κ and C i are the permeability and inertia coefficient of copper foam porous medium respectively, which can be obtained from the following equations [START_REF] Fourie | Pressure drop modelling in cellular metallic foams[END_REF].

K = (εd k ) 2 36χ(χ -1) (7
)

C i = 0.00212(1 -ε) -0.132 ( d f d p ) -1.63 / ̅̅̅̅ K √ (8)
where χ, d k , d p and d f are the tortuosity coefficient, characteristic length, pore size and fiber diameter of the metal foam respectively [START_REF] Calmidi | Transport phenomena in high porosity fibrous metal foams[END_REF]. These are essential parameters to describe the structural characteristics of metal foams, and they directly determine the porosity and pore density of metal foams, which can be obtained from Eqs. ( 9)- [START_REF] Arshad | Towards the thermal management of electronic devices: a parametric investigation of finned heat sink filled with PCM[END_REF].

χ = 2 + 2cos ( 4π 3 + 1 3 cos -1 (2ε - 1
)

) ( 9
)
d k = χ 3 -χ d p ( 10
)
d p = 22.4 × 10 -3 ω ( 11
)
d f = 1.18 ̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ̅ (1 -ε) 3π √ d p ( 12
)
Energy equation:

[ (1 -ε) ( ρc p ) mf + ε ( ρc p ) pcm ] ∂T ∂t + ( ρc p ) pcm ( u ∂T ∂x + ∂T ∂y ) = λ eff ( ∂ 2 T ∂x 2 + ∂ 2 T ∂y 2 ) -ερL ∂f ∂t (13)
where, ε is the porosity of metal foam, the subscripts pcm and mf represent the phase change material (paraffin wax) and metal foam (copper) respectively, and λ eff is the effective thermal conductivity of metal foam phase-change material composite, calculated by the model proposed by Bhattacharya [START_REF] Bhattacharya | Thermophysical properties of high porosity metal foams[END_REF]:

λ eff = 0.35 ( ελ pcm + (1 -ε)λ mf ) + 0.65 ( ε λpcm + 1-ε λmf ) (14)

Boundary and initial conditions

The computation domain is shown in Fig. 2. During the melting period, the initial temperature of each filling configuration system is set to 293.15 K. All the four walls adopt the principle of none-slip, and the initial velocity in × and y directions is zero. The left wall is the heating surface, and the heating temperature is set to 423.15 K, which is a constant temperature heating method. The remaining three walls are defined as adiabatic to reduce heat loss. The junction between the composite and the pure paraffin regions were set as interior in the CFD program to ensure the heat and mass exchange between the two regions at the junction.

Numerical procedure

The governing equation is discretized by finite volume method (FVM) based on the double-precision solver. The SIMPLE algorithm is employed to handle the pressure-velocity coupling, and the pressure correction equation is discretized by the PRESTO method. Second order upwind algorithm is used in momentum and energy equations, and second order implicit algorithm is used in transient formulation. In addition, the residual convergence criteria for the continuity, the velocity and the energy equations are set at 10 -4 , 10 -4 and 10 -6 , respectively.

Independence tests and model validation

Independence tests

In order to improve the accuracy of the numerical model and method, the number of grids and the time step were independently verified before model validation in this study. Through the analysis of the results, the optimal grid number and time step are selected to eliminate the errors caused by these factors and improve the computational efficiency. The independence test model is the case of a fully filled copper foam (ε = 0.9, 10PPI). The temperature of the left side wall is constant at 423.15 K, which is used as the heating surface, and the remaining wall is treated with adiabatic condition. The initial temperature of the LHTES unit is 293.15 K.

Firstly, the independence test was carried out on the grid numbers. The models with 5400, 10,564, 33,750 and 86,400 elements were selected respectively, and the time step was set as 1 s. The simulation results of the variations of liquid fraction with time under different gird numbers are shown in Fig. 4. With the increase of grid number, the process quickens sharply and then flattens out. In the case of 33,750 and 86,400 elements, although the number of grids is very different, the two curves overlap approximately, and the maximum deviation is only 0.73%, which is negligible. Considering the calculation accuracy and cost, the number of grid is determined as 33750, and the meshing is depicted in Fig. 5.

Secondly, under the condition that the grid number is 33750, the independence test is carried out for the calculation cases with time steps of 2 s, 1 s, 0.2 s and 0.1 s. Fig. 6 shows the simulation results of the variations of temperature with time at different time steps. Convergence can be achieved under each time step. When the time step is set to 0.2 s, the maximum deviation from the result of 0.1 s is only 1.1%. So a time step of 0.2 s is sufficient.

Finally, 33,750 was chosen as the optimum grid number, and the corresponding time step is 0.2 s.

Model validation

In order to verify the accuracy of the current numerical model, the numerical method in this study was used to deal with the PCM melting problem in porous media in the literature of Zheng et al. [START_REF] Zheng Z-J | Effect of metal foam with two-dimensional porosity gradient on melting behavior in a rectangular cavity[END_REF]. And compared with the simulation results. A square cavity of 60 mm × 60 mm was filled with paraffin and metal foam uniformly, and then the four operating conditions are simulated respectively on this basis. In Case-a, surface (constant heat flow: 4000 W/m 2 ), and the other walls are adiabatic. The melting temperature ranges from 319.63 k to 333.54 K, and the initial temperature is 299.15 K. The metal foam is made of copper (ε = 0.9, 10PPI). The model size and given conditions were set according to the literature experiment. In addition, a monitoring point with ×, y and z coordinates of 100 mm, 13 mm and 50 mm was established. The results are shown in the Fig. 8. It can be seen that the simulation results are in good agreement with the experimental data, and the maximum deviation is only 1.5%, which proves the reliability of the numerical method in this paper.

Results and discussion

Melting characteristics 4.1.1. Melting interface and temperature field

In order to better study the melt heat transfer mechanism under different filling configurations, the computational domain was divided into pure paraffin region and metal foam-paraffin region, and the liquid fraction, melting interface, temperature field and velocity field in each region was discussed.

The liquid fraction changes and liquid-solid interface evolution of the two filling configurations are shown in Figs. 9 and 10, respectively. Fig. 9 shows the similar degree of liquefaction of the two configurations at the same metal foam filling height ratio (ξ). However, with the increase of ξ, the heat conduction was also enhanced, which accelerated the melting of paraffin in the composite region, so that more paraffin was melted at the same time and showed a higher liquid fraction. At the beginning of melting, the liquid fraction in each case develops linearly and has a high slope. This is because the paraffin in the composite region melts rapidly and thermal conductivity dominates heat transfer. For the two filling configurations, since ξ is the same, the liquid fraction curves almost coincide at the initial stage, which is proved in Fig. 10. Before 200 s, the liquid-solid interface of the two types of configurations in the pure paraffin and composite regions are similar. However, when the paraffin in the composite region melts, the development of liquid fraction takes a turning point. It can be seen that the liquid fraction of the bottom filling configuration has a steeper trend than that of the top filling configuration. This is because when the bottom is filled with metal foam, it acts as a heat source, and the melted hot paraffin moves from the left and bottom of the pure paraffin region to the top at the same time. As shown in Fig. 10, the melting interface shows a C-shaped movement to the right, which is similar to the phenomenon observed by Wang et al. [START_REF] Wang | Effect of copper metal foam proportion on heat transfer enhancement in the melting process of phase change materials[END_REF] in the experiment. When the metal foam is filled at the top, it provides a top heat source for the pure paraffin region. However, due to natural convection, the high temperature liquid generated from the pure paraffin region will flow to the composite region. This, in turn, causes the pure paraffin region to lose heat and thus melt more slowly. As can be seen from Fig. 10, the melting interface of the top filling configuration slowly moves to the lower right corner.

Figs. 11 and 12 respectively show the influence of different configurations on temperature change and the temperature distribution at the time of melting completion. At the early melting stage, only trace paraffin close to the hot wall in the pure paraffin region melts, while in the composite region, due to the high thermal conductivity of the metal foam, the paraffin in this region rapidly heats up and liquefies. So in each case, temperatures rose quickly when the melting started. And with the increase of the metal foam, the more hot paraffin obtained after melting, so the unit reaches a higher average temperature. It can be seen from Fig. 11 that the average temperature of the top-filled configuration is higher than that of the bottom-filled configuration. This is because when the metal foam is located at the bottom, the heat of the hightemperature paraffin is absorbed by the pure paraffin region, making the average temperature lower. When the metal foam is at the top, the liquid paraffin accumulates in the composite region and gradually heats up as it melts. The top filling and bottom filling configurations exhibit different temperature variations as the composite region melts. The pure paraffin region in the top filling configuration melts mainly by heat conduction, so the temperature rise rate is slow. However, natural convection was enhanced in this region with the bottom filling configuration, and the temperature rose at a higher rate in the later stage. However, due to the greatly shortened melting time, the temperature eventually reached was lower, which can be verified in Fig. 12.

It is found in Fig. 12 that the bottom temperature of each configuration is low, which is because the natural convection of paraffin transfers heat to the top portion, making it difficult for the bottom portion to melt. For bottom filling configuration, the temperature distribution becomes more uniform with the increase of metal foam. However, this is not obvious in the top filling configuration. Due to the longer melting time, there is little difference in the overall temperature. High temperature paraffin accumulates in the upper portion, and the temperature decreases from top to bottom with the development of melting. It can be seen from the liquid fraction and temperature change diagram that the melting rate and temperature rise rate of the bottom filling configuration are slightly higher than that of the top filling configuration. However, it is also obvious from Fig. 12 that the high temperature region of the top filling configuration is much greater, which can store more sensible heat.

Velocity field

The dimensionless time is defined as follows:

t * = t t m ( 15
)
Where t m is the complete melting time of the system. t is an arbitrary time during the melting process. The value of t* is between 0 and 1.

Figs. 13-15 present the velocity field at different times for pure paraffin and the two types of filling configurations. In the pure paraffin region, each case had a vortex in the upper left corner at the initial stage of melting (before t* = 0.5). That's because, the paraffin near the heating surface melts first. Under the action of buoyancy, the liquid paraffin flows upward, forming a vortex. However, the difference in the number of vortices in Figs. 14 and 15 is due to the location of the metal foam embedment. For the bottom filling configuration, because the metal foam accelerates the phase transformation process, there is a difference in temperature between the pure paraffin and the composite regions. This temperature imbalance results in two small vortices at the junction of the two regions, which located in the lower left and lower right corners of the pure paraffin region. The existence of the vortex zone is explained by the stable circulating flow of liquid, and the more vortices represent the stronger heat transfer. As the paraffin on the top portion is completely melted (after t* = 0.5), the vortex originally located in the upper left corner moves down and joins the vortex in the lower left corner of the pure paraffin region, accelerating the melting process. However, for the top-filled configuration, despite the temperature difference at the junction of the two regions, liquid paraffin moved toward the composite region due to buoyancy, so there were no small vortices in the pure paraffin region as there were in the bottom filling configuration.

Table 2 shows the order of magnitude of velocity of liquid paraffin in different regions, which is similar to the order of magnitude of 10 -4 m/s observed by Tian et al. [START_REF] Tian | Thermal and exergetic analysis of metal foam-enhanced cascaded thermal energy storage (MF-CTES) [J[END_REF] in the case of full metal foam filling. However, due to the different metal foam filling positions in this paper, some interesting phenomena were found. It is observed that the bottom filling configuration has a higher velocity in any region. The density of paraffin in the solid and liquid phases is different due to the different temperature, which flows upward under the action of buoyancy. With no metal foam on top, the liquid paraffin produced in the composite region flows upward easily. Therefore, the flow velocity of the paraffin in the composite region is increased, which also strengthens the natural convection in the pure paraffin region and increases the flow velocity in this region. So it is reasonable to observe an order of magnitude of 10 - 3 m/s, greater than 10 -4 m/s, in the composite region. The velocity in the pure paraffin region is on the order of 10 -4 m/s, slightly larger than the 10 -3 m/s observed by zhang et al. [START_REF] Zhang | Melting heat transfer characteristics of a composite phase change material fabricated by paraffin and metal foam[END_REF] in pure paraffin, which can also be considered reasonable. Similarly, when the top is filled with metal foam, the liquid paraffin remains at the top. Due to the obstruction of the porous structure of the metal foam, the upward flow of the newly generated liquid paraffin in the pure paraffin region is limited. This does not help natural convection, so the overall velocity is slower. In the composite region, the order of magnitude is 10 -5 m/s, less than 10 -4 m/s, and the cycle is weak or even disappears (Fig. 15(a), t* = 0.5, 0.75; Fig. 15(b), t* = 0.75). In addition, it can be found that the velocity of liquid paraffin in the composite region is lower under any configuration, meaning that the metal foam impeded the flow of paraffin. The bottom filling configuration is advantageous for the enhancement of natural convection.

In addition, some interesting phenomena have been observed in this study of velocity field. As shown in Fig. 16, the distance between the center of the vortex in the composite region and the heating surface was defined as θ, and Δθ was defined as the deviation rate.

Δθ = θ ξ 2 -θ ξ 1 θ ξ 1 (16)
The subscript ξ 1 , ξ 2 indicates the different height ratio of metal foam. When the metal foam filling ratio is small (ξ = 0.25), the vortex center is close to the heating surface. This is because the paraffin near the heating surface melts first and the heat transfer is strong. However, the vortex center gradually moves away from the heating surface with the increase of the filling ratio. This indicates that the position of intense heat transfer moves from the wall to the material center, which also verifies that the metal foam can effectively enhance the heat transfer process of PCM. For the top-filled configuration, due to the accumulation of high temperature liquid paraffin in the composite region, the position of intense heat transfer is more centered, that is, the θ is higher than that of the bottom-filled configuration. However, it can be seen from the value of Δθ that when the same metal foam height ratio is increased, the deviation rate of the bottom-filled configuration is larger, which is more conducive to heat transfer enhancement.

Table 2

Order of magnitude of velocity of liquid paraffin in different regions.

Thermal energy storage performance

In this paper, thermal energy storage of metal foam-paraffin composite LHTES units is defined as Eq. ( 17)-( 21):

Q = Q pcm + Q mf ( 17
)
Q pcm = Q la + Q se ( 18
)
Q la = m pcm L (19)
Q se = c p,ps m pcm (T m -T 0 ) + c p,pl m pcm (T -T m ) (20)
Q mf = c p,mf m mf (T -T 0 ) (21)
where Q la is latent heat; Q se , is sensible heat; c p,ps , is the heat capacity of solid PCM; c p,pl is the heat capacity of liquid PCM; L is the latent heat of PCM; T 0 is the initial temperature; T m is the phase transition temperature.

The heat storage rate is defined as Eq. ( 22):

ν = Q/t m (22)
Figs. 17 and 18 presents the complete melting time and thermal energy storage rate for different cases. It can be seen from Fig. 17(a) that for the bottom-filled configuration, the melting time with ξ = 0.25, 0.5 and 0.75 are 1500 s, 999 s and 467 s, and the thermal energy storage rate are 1.26 kJ/s, 1.93 kJ/s and 4.08 kJ/s in Fig. 18(a), respectively. As ξ increases, the melting time decreases and thermal energy storage rate increases. The same pattern of variation can be observed in the top-filled configuration (Figs. 17(b) and 18(b)). However, it is worth noting that the melting time of ξ = 0.25, 0.5 and 0.75 are 2437 s, 1704 s and 913 s, and the thermal energy storage rate are 0.83 kJ/s, 1.20 kJ/s and 2.26 kJ/s, respectively. For partial filling configuration, when ξ is the same, the bottom-filled configuration melts faster and the thermal energy storage rate is higher. For partially filled configurations, when ξ = 0.75, both bottom-filled and top-filled configurations reach the shortest melting time and the highest thermal energy storage rate.

The influence of different filling configurations on thermal energy storage and average temperature are portrayed in Fig. 19. When ξ is the same, the mass of paraffin in each case is equal. So when the paraffin is completely melted, the latent heat is the same. The difference of heat storage capacity between the bottom-filled and top-filled configurations comes from the sensible heat of paraffin and metal foam. It can be seen from Eqs. ( 20) and ( 21) that sensible heat is related to temperature difference, so the change of sensible heat should be consistent with the change of the trend of the average temperature (Fig. 19). As depicted in Fig. 19(a), for the bottom-filled configuration, the thermal energy storage of ξ = 0.25, 0.5 and 0.75 are 1961.4594 kJ, 1925.26366 kJ and 1905.0215 kJ, respectively. As ξ increases from 0.25 to 0.75, heat storage decreased by 2.29%, 4.08% and 5.09%, respectively. The main reason is that with the increase of ξ, the mass of paraffin decreases, which leads to the decrease of thermal energy storage. On the other hand, with the shortening of melting time, the final temperature reached after the complete melting of composite PCM decreases, which also leads to the decrease of thermal energy storage. For ξ = 1, the thermal energy storage loss is up to 5.53% of that of the case when ξ = 0.

For the top-filled configuration, as depicted in Fig. 19(b), when ξ = 0.25, 0.5 and 0.75, the TES are 2028.1684 kJ, 2044.9843 kJ and 2064.4558 kJ, respectively. When ξ increases from 0.25 to 0.75, the heat storage increases by 0.97%, 1.72% and 2.78%. It is found that the thermal energy storage increases slightly with ξ. This is caused by its melting characteristics are different from the bottom-filled configuration. The enhanced heat transfer effect of metal foam only affects the composite region, resulting in a much longer melting time than that of the bottom-filled configuration. When ξ is the same, the final temperature reached at the completion of melting is higher in the top-filled configuration. Consequently, its thermal energy storage is higher than that of the bottom-filled configuration, and with the increase of ξ, the final temperature will further rise, leading to the increase of thermal energy storage. When ξ = 1, the melting time and thermal energy storage decrease sharply at the same time because there is no pure paraffin region.

To sum up, from Figs. 17-19, it can be concluded that the shorter melting time is at the expense of thermal energy storage, and those with more advantages in thermal energy storage have longer melting time. This is a choice that needs to be evaluated.

Economic performance

In practical engineering applications, in addition to the thermal energy storage performance of materials, economic performance is also a vital index, which determines the project cost.

The economic index P c is defined in Eq. ( 23), represents the heat storage performance of materials per unit time per unit cost.

P c = Q t m C ≈ m pcm L t m (im mf + jm pcm ) (23)
C is the total cost of PCM and metal foam, i and j are the unit price of metal foam and PCM respectively.

In dimensionless P c ' is defined as Eq. ( 24), where P c,N represents the case when ξ = 0, M is the ratio of the respective unit prices of metal foam and PCM. The complete melting time, thermal energy storage and total material cost determine the performance of the LHTES units. Consequently, a comprehensive strengthening goodness factor (CSGF) is proposed in this study to compare the performance of different filling configurations. The definition of CSGF is shown in Eq. ( 25) -( 28):

η = aη Q + bη t + cη Pc , a + b + c = 1; 0 ≤ a, b, c ≤ 1 (25) η Q = Q Q N (26) η t = t N t m ( 27
)
η pc = P c ' (28)
where η Q is the dimensionless goodness factor of thermal energy storage, and is the ratio of Q to Q N (ξ = 0); η t is the dimensionless goodness factor of the complete melting time, and is the ratio of t N (ξ = 0) to t m ; η pc is the dimensionless goodness factor of P c , which is defined in Eq. ( 23); a, b and c are the weighting coefficients. η is a constant greater than zero. The larger η is, the higher the comprehensive performance is.

When a = b = 0 and c = 1, it means that only economic factors are considered, η will be equal to P c ', which is included in CSGF as a particular condition shown in the previous section. According to the working conditions in this study, the CSGF when only considering the complete melting time (η t ), thermal energy storage (η Q ), and the same weight (a = b = c = 1/3) are given below. The η t and η Q in different cases are illustrated in Fig. 21(a) and 21(b), respectively. If the complete melting time is the only criterion in practical application, the bottomfilled configuration (ξ = 0.75) with the largest η t can be selected as the optimum configuration. If only TES is considered, the top-filled configuration (ξ = 0.75) with the largest η Q can be chosen as the most suitable configuration. When the three are equally focused in practice and there is no emphasis, as presented in Fig. 21(c), it is clear that the bottom-filled configuration (ξ = 0.75) has an absolute advantage. Consequently, CSGF can be used to determine the optimum filling configuration for practical applications more quickly.

Conclusion

In the current study, for the heat storage unit of LHTES system, two types of filling modes, top-filled and bottom-filled configurations, are designed. Meanwhile, the different filling height ratio of the metal foam is considered. The melting characteristics of different LHTES units were analyzed and a new criterion CSGF was proposed to evaluate their comprehensive performance. Conclusions can be summarized as follows:

(1) When the metal foam is filled at the bottom, the liquid paraffin has a higher flow rate, which helps to strengthen the natural convection in the pure paraffin region, improve the overall melting rate and temperature rise rate, and accelerate the melting process. However, the internal flow rate of the top filling configuration is slower, and the heat transfer enhancement is not as good as that of the bottom filling configuration, but it can achieve a higher final temperature. (2) The top-filled configuration will increase the thermal energy storage, while the bottom-filled will significantly weaken it. The metal foam filled at the top (ξ = 0.75) has the highest thermal energy storage, compared with the fully filled one, is 7.47% higher with about 25% of the metal foam saved, which is of great practical significance for improving economic performance.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

(3) CSGF is a new dimensionless evaluation criterion of the complete melting time, thermal energy storage and cost. Unlike previous criterions, it takes both heat storage and economy into account. Through its analysis, it is concluded that the bottomfilled configuration has the advantage in complete melting time and economy, while the top-filled configuration has the advantage in thermal energy storage. (4) The strategy of partially filling metal foam has been proved to be effective in enhancing heat transfer, increasing thermal energy storage or reducing cost. Both the bottom and top filling config-urations have their own advantages. In future studies, the struc-ture can be further optimized and the metal foam layer with porosity gradient, material gradient and thickness gradient can be set to improve heat transfer and economic performance.

  Mushy coefficient b the weight coefficient of η t c the weight coefficient of η pc c p Specific heat capacity (J⋅kg -1 ⋅K -1 ) C i form coefficient (m -1 ) d f fiber diameter (m) d p pore size (m) d k characteristic length (m) f liquid fraction g acceleration of gravity (m⋅s -2 ) H LHTES unit height (m) H mf metal foam height (m) i unit price of metal foam ($⋅kg -1 ) j unit price of PCM ($⋅kg -1 ) K permeability (m 2 ) L latent heat of fusion (kJ⋅kg -1 ) m mass (kg) M unit price ratio of the metal foam to the PCM Q TES capacity (kJ) S source term t time (s) t m melting time (s) T temperature (K) T m melting point (K) u, v x,y direction velocity (m⋅s -1 ) W LHTES unit width (m) x,y Cartesian coordinates (m) Greek symbols α thermal diffusivity (m 2 ⋅s -1 )

Fig. 1 .

 1 Fig. 1. Structural and performance advantages of metal foam.

Fig. 2 .

 2 Fig. 2. Physical model of different configurations of metal foam in metal foam-paraffin systems.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Position of monitoring points.

Fig. 5 .

 5 Fig. 5. Meshing system (ξ = 1).

Fig. 6 .

 6 Fig. 6. Variations of temperature with time for different time steps.

Fig. 7 .

 7 Fig. 7. Model validation results under the configuration of fully metal foam filled.

Fig. 8 .

 8 Fig. 8. Comparison of current numerical results with the previous experimental results.

Fig. 9 .

 9 Fig. 9. The influence of different fillng configurations on liquid fraction.

Fig. 10 .

 10 Fig. 10. Liquid-solid interface evolution of bottom-filled and top-filled configurations at different.

Fig. 11 .Fig. 12 .

 1112 Fig. 11. Temperature distribution with melting time for different filling configurations.

Fig. 13 .

 13 Fig. 13. Velocity field of pure paraffin(Case-N) at different times. The dimensionless time is 0.25, 0.5 and 0.75, respectively.

Fig. 14 .

 14 Fig. 14. Velocity field of the bottom filling configuration at different times. The dimensionless time is 0.25, 0.5 and 0.75, respectively.

Fig. 15 .

 15 Fig. 15. Velocity field of the top filling configuration at different times. The dimensionless time is 0.25, 0.5 and 0.75, respectively.

  The order of magnitude of velocity (m/s)

Fig. 16 .Fig. 17 .

 1617 Fig. 16. The deviation distance between the vortex center and the heating surface varies with the height ratio of metal foam.

Fig. 18 .

 18 Fig. 18. Influence of the various filling heights on the thermal energy storage rate with bottom filling (a) and top filling (b) configurations.

Fig. 19 .

 19 Fig. 19. Influence of the various filling heights on the TES and average temperature with bottom filling (a) and top filling (b) configurations.

Fig. 20 .

 20 Fig. 20. Variations of the economic index P c ' with M for different filling configurations.

Fig. 21 .

 21 Fig. 21. Influence of the different filling configurations on the (a) η t , (b) η Q and (c) CSGF with M.

(

  bottom-filled, ξ = 0.75) is superior.4.4. Comprehensive evaluation indexCRediT authorship contribution statementQifan Ying: Writingoriginal draft, Conceptualization, Methodology, Software, Visualization. Hui Wang: Writingreview & editing, Conceptualization, Project administration, Supervision. Eric Lichtfouse: Writingreview & editing.

  -1 ⋅s -1 )

	ρ	density (kg m -3 )
	ξ	height ratio of metal foam
	η	value of the CSGF
	δ	small constant number
	χ	tortuosity coefficient
	ω ν	pores per inch thermal energy storage rate (kJ⋅s -1 )
	θ	distance of deviation (m)
	Subscripts
	0	initial time
	eff	effective value
	pcm	phase change material
	mf	Metal foam
	l	liquid state
	la	latent heat
	s	solid state
	se	sensible heat
	w	wall
	Abbreviations
	BTC	bottom filling configuration
	CSGF	comprehensive strengthening goodness factor
	LHTES	latent heat thermal energy storage
	PCM	phase change material
	PPI	pores per inch
	TES	thermal energy storage
	TFC	top filling configuration

Table 1

 1 Thermo-physical properties of paraffin and copper foam.

	Paraffin	Copper

Acknowledgments

The study was funded by State Key Laboratory for Modification of Chemical Fibers and Polymer Materials, Donghua University (KF2123), Shanghai Sailing Program (19YF14011700).


Nicolas Papaiconomou
email: nicolas.papaiconomou@univ-cotedazur.fr

Jean-Pierre Magnin

Isabelle Billard

Dr Billard

Dr Magnin

Treatment of Pulp Black Liquor with an Ionic Liquid: Optimisation by Using Response Surface Methodology

In this study, selective extraction towards ionic liquid trihexyltetradecylphosphonium chloride ([P66614]Cl) of lignin residues, polysaccharides and organic acids present in black liquor (BL), the main principal wastewater of pulping industry was studied. With the objective of finding optimized conditions allowing to extract lignin resideus while polysaccharides and organic acids remain in aqueous solution, a design of experiments approach based on a response surface methodology was used. Three continuous factors, namely initial pH varying from 9 to 13.5, dilution of BL varying from 5 to 20 and volumetric ratio of black liquor vs. ionic liquid RV., varying from 1 to 19, were investigated. Concentration of ligni residues, polysaccharides and organic acids were measured using Folin-Ciocalteu method, the anthrone method and HPLC, respectively. Results showed that a multi-response optimization led to the extraction of 84.8 % of lignin residues, 66.0 % of polysaccharides, and no extraction of OA under optimised conditions.

Introduction

Ionic liquids (IL) are salts with melting temperatures below 100°C and their versatile physico-chemical properties render them interesting in different industrial and commercial applications such as biotechnology. [START_REF] Shiflett | Commercial Applications of Ionic Liquids[END_REF][START_REF] Itoh | Application of Ionic Liquids in Biotechnology[END_REF][START_REF] Kuzmina | Application, Purification, and Recovery of Ionic Liquids[END_REF] Recently, ionic liquids have been used as new solvents for pretreatment of lignocellulosic biomasses (switch-grass, wheat straw, corn stover, bagasse, ….) [START_REF] Tu | [END_REF][5][6][7][8][9] but few studies have already appeared on the use of ILs for the selective extraction of compounds present in such complex systems as pulping factory black liquor, the most prominent aqueous waste issued from pulping factories. [10,11] Several processes, such as well-known Kraft process are currently available for depolymerizing lignin from wood in order to obtain a paper pulp. [START_REF] Biermann | Essentials of Pulping and Papermaking[END_REF][START_REF] Gellerstedt | [END_REF] In all cases, a black liquor is formed along the process, with significant hazard and safety issues related to the volume and composition of waste produced. Black liquor represents up to 1200 Mt produced each year in the world. [14] It generally consists in an alkaline solution with pH close to 12 or 13, composed with polysaccharides (PS), organic acids (OA), and residues and fragments of lignin, a polyphenolic polymer that will be quantified here in polyphenols equivalents (PP). The presence of such polyphenolic polymer fragments givs its deep black color to so-called black liquor. [15,16] Polysaccharides and organic acids represent two classes of compounds with a high biodegradability, while lignin residues are not. [START_REF] Biermann | Essentials of Pulping and Papermaking[END_REF]17] Consequently, the extraction of lignin residues from the black liquor solution is a necessary step in order to increase the accessibility to polysaccharides and organic acids for further biological treatments, without losing the reuse potential of lignin, for instance under the form of polyphenols after depolymeriztion process, like adhesive, [18] composite materials, [19] or as active ingredients in cosmetics or food applications. [20,21] In a previous study, the efficiency of two phosphonium ionic liquids containing Cl -or Br -anions, namely [P66614]Cl and [P66614]Br, in order to extract polyphenolic lignin residues, polysaccharides, and organic acids from pulping factory black liquor was examined. [10] Extraction was found to be influenced by the dilution of black liquor (BLD), volume ratio of aqueous phase vs. ionic liquid (Rv) and initial pH. Response Surface Methodology (RSM), introduced in the early 1950s, is an efficient and practical tool for designing experiments, fitting response surface models and determining optimum operating conditions. [22][23][START_REF] Myers | Response Surface Methodology: Process and Product Optimization Using Designed Experiments[END_REF] Simple, interaction and quadratic effects of selected factors on the studied response can be studied simultaneously starting from an experimental matrix with limited number of runs. Design of experiments (DOE) presents, in comparison to the classical OFAT methodology (One Factor At a Time), many advantages such as time and money saving, limited number of experiences, thereby making it possible to investigate factor interactions. It is widely used in different industrial fields such as food industry, [START_REF] Granato | Mathematical and Statistical Methods in Food Science and Technology[END_REF][START_REF] Yolmeh | [END_REF] analytical chemistry, [27] food flavor and sensor analysis, [28] in energy application, [29] in biotechnological processes [30] and in extraction of plant materials. [START_REF] Aydar | Statistical Approaches With Emphasis on Design of Experiments Applied to Chemical Processes[END_REF][START_REF] Morais | [END_REF] More specifically, RSM has been successfully used in several applications using ILs, such as chemical analysis, [33] in food industry, [34] catalyst applications, [35] biomass pre-treatment [36] or organic synthesis. [37] In a previous article [10] , phosphonium ionic liquids were proven to be promising solvents for extracting lignin residues from black liquor. OFAT (one factor at a time) methodology was previously used for studying influence of black liquor dilution and volume ratio on extraction. None combined effect between parameters were however studied. In addition, initial pH effect was examined for one unique combination of BL dilution and Rv, 0.1% and 10, respectively. Therefore and in order to go beyond the proof of concept reported previously, [10] the aim of this study is first to model the individual extraction percentages and distribution coefficients for lignin residues, hereafter studied as polyphenol equivalent (PP), PS and OA towards [P66614]Cl, using Response Surface Methodology.

That way, all three parameters, BL dilution, initial pH and Rv are take extensively into account over expanded ranges, so that all possible interactions between these three parameters are described. In a second step, with the objective of further treating such waste, multiresponse optimization is established in order to maximize extraction of lignin residue while minimizing PS and OA extraction percentages. Finally, possibilities to extract selectively PP towards [P66614]Cl, starting from different initial conditions is discussed.

Results and Discussion

DOE observations

According to the DOE methodology, 35 experiments were carried out with varying pH, BLD an RV. Because black liquor exhibits a pH of 13.5 and because lignin residues are less soluble in acidic water, [15] pH was used as an extraction parameter. Black liquor dilution and aqueous vs. IL phase volumes are common parameters in extraction experiments, the former being related to the initial concentration of extractable species while the latter correspond to the amount of extracting phase used compared to that of black liquor. Reasons for the value ranges selected for each parameter, namely pH, BLD and Rv varying between 9 and 13.5, 0.2 and 0.05 or 1 and 19, accordingly, are detailed in the experimental section.

For clarity reasons, all results are collected in the support information file. Extensive details and results for the experimental CCD design and extraction percentage and distribution yields for lignin parts, discussed here as polyphenolic equivalents, polysaccharides and organic acids are collected in the supplementary information file (Table S1). In agreement with previous results [6] , black liquor solution studied here contains 30 g.L -1 of PP equivalents, 23 g.L -1 of PS and 19 g.L -1 of OA. In all extraction experiments, upon contacting the ionic liquid with a black liquor sample, a marked change in color is observed within 30 s. The initially black aqueous solution becomes yellow or translucent while the ionic liquid turns brown or black depending on the experimental conditions. For clarity reasons, only two snapshots of a typical extraction experiment are shown in Fig. 1. Snapshots for all extraction experiments are collected in Fig. S1 of the SI file. Dark color of the black liquor is due to the presence of lignin residues. The transfer of this black liquor towards the ionic liquid implies that these lignin residues are significantly extracted out from the black liquor. Results reveal that generally, lignin residues are well extracted towards the ionic liquid, exhibiting distribution coefficients range from 5 to 500. PS are somehow less extracted, with distribution coefficient (DPs) raging from 0. Because the experiments presented in table S1 were done following recommendations of the experimental design software, discussing the data and find specific trends in the influence of parameters on the extraction of PP, PS and OA is not straightforward. Nevertheless, it first appears that worse extraction results for both PP and PS occur when initial pH is left at 13.5, while best results are obtained at pH 9. Second, for the same pH and BLD, an increase in RV yields an increase in distribution coefficients. The influence of BLD on distribution coefficients is less clear and not linearly correlating the extraction yields or distribution coefficients. Besides, the trend in extraction for PP and PS appear to be similar. That is, low and high values for DPS or DPP are obtained under the same experimental conditions, namely with an initial pH value of 11.25 or 9 and a high RV value.

Depending on the experimental conditions, a small dark precipitate appears both all along a few tube walls and in the bottom of the tubes after centrifugation, as showed in figure S1 of the supporting information file. It is known that lignin solubility depends on pH, the lower the pH, the lower the solubility of lignin. [15] However, only a few tubes, namely tubes 10 and 23 exhibit a low pH of approximately 6.6 and 6.3, respectively, and a precipitate. One possible explanation might be in some cases related to the solubility of lignin in acidic solutions. But it is also related to an interaction between acid groups of lignin and ionic liquid. Such an interaction can between, say, an acid group of lignin residues and an ionic liquid cation, might lead to a deprotonation of this acidic group, yielding formation of ion pair with an increased molecular weight. The solvation energy of these ions is thus expected to increase, as well as hinder the interaction between such acid group and water. This interaction, or ion association phenomenon, might yield to an extraction towards the ionic liquid phase as mentioned above, but also to some precipitation when the ionic liquid is saturated or the resulting pair is too large. Such result was previously obtained with a variety of ionic liquids and extracted compounds. [38][39][40] A more thorough investigation of the pH is presented in Figure 2, plotting the pH measured after extraction of black liquor as a function of initial pH for various RV factors 2, 10 and 19. Final pH measured in black liquor after extraction ranged from 1 to 13, depending on initial pH. It appears that with initial pH values of 9, final pH is significantly lower than the initial value of 9. Variation between final and initial pH (|𝛥𝑝𝐻|) thus range from 0.7 to 7.75 with seven experiments out of ten yielding |𝛥𝑝𝐻| above 2. Such a decrease is less pronounced starting from a pH value of 11.25: values range from 0.34 to 7.73, with four experiments out of fifteen yielding |𝛥𝑝𝐻| values above 2. Starting from a value of 13.5 for the initial pH value, hardly any variation is observed. This reveals a certain formation of H + ions upon extraction, a phenomenon that can be related to the above-mentioned explanation, no linear nor evident correlation was found between |𝛥𝑝𝐻| and RV or BLD. Therefore, influence of both factors on both responses was studied using the RSM methodology.

Modeling of pH, extraction percentage and distribution coefficient

Experimental results of CCD were analyzed with the help of Design Expert software according to the polynomial quadratic response function given in eq. (1): 𝑌(𝐵𝐿𝐷, 𝑝𝐻, 𝑅 𝑣 ) = 𝑎 -𝑏(𝐵𝐿𝐷) + 𝑐(𝑝𝐻) + 𝑑(𝑅 𝑣 ) + 𝑒(𝐵𝐿𝐷). (𝑝𝐻) + 𝑓(𝐵𝐿𝐷). (𝑅 𝑣 ) + 𝑔(𝑝𝐻). (𝑅 𝑣 ) + ℎ(𝑝𝐻) [START_REF] Itoh | Application of Ionic Liquids in Biotechnology[END_REF] + 𝑖(𝑅 𝑣 ) [START_REF] Itoh | Application of Ionic Liquids in Biotechnology[END_REF] + 𝑗(𝐵𝐿𝐷) [START_REF] Itoh | Application of Ionic Liquids in Biotechnology[END_REF] (1)

with Y a linear, logarithmic or square root function of each experimental result obtained in this work, that is either pH, the extraction percentages or distribution coefficients of PP, PS or OA, accordingly. BLD, RV and pH are the black liquor dilution, volume ratio and the final pH, respectively. a to j values are the adjustable parameters of the model. Their values are collected in Table S2 of the supplementary information file. Each individual response model, namely, a response to be compared to experimental data, is analyzed in a first step by Analysis of Variance (ANOVA). During this step, response function Y from eq. ( 3) is fitted to experimental results as presented in Table S2 and the software selects various models for describing experimental results. Response function was allowed to describe experimental responses, namely |pH|, DPP, DPS, DOA, PPP, PPS and POA respectively, under the form of either Y=X, 𝑌 = 𝑙𝑛(𝑋 + 𝜀) or 𝑌 = √𝑋 + 𝜀 , X being an experimental response as detailed above and  a numerical value. The latter value was necessary in case some experimental data exhibited null values. From that analysis, detailed in the supplementary information file, the best suited models for the extraction percentage and distribution coefficients of PP, PS and OA, as well as the pH variation are collected in Table S2.

It is worth noticing that run 14 shows unusual data points that deviate significantly from other observations. Because this is most probably due to incorrect experimental run, according to ref. [32], the run is withdrawn from the analysis of externally studentized residuals after model refining and eventual response transformation. Therefore, run 14 is an outlier point located outside standard deviation limits (data not shown) and was consequently eliminated for future analysis. Fitting of response functions was carried out individually by analyzing distributions of residual plots according to the residual normal probability plots for each response, residual values versus predicted plot, and predicted values versus actual values (not shown). In all cases, a linear distribution of residuals along a straight line was observed, and a disordered organization of residuals and predicted values against predicted response values and values of runs is obtained. Along with the results obtained by the ANOVA analysis detailed in the supplementary information file, all diagnostic tests are thus satisfactory, implying that the models can be used in a wide spectrum of design.

Modeling of pH variation

According to the fitting procedure described in the preceding subsection, |pH| response was successfully described as a function ln(|pH|+0.00775) and was found to depend on all three factors, RV, BLD, and pH (Table 1). In this work, because three factors are investigated, response, such as |pH|, extraction percentages or distribution coefficients, as detailed later in the manuscript, will be discussed following two representations. One is a classical, socalled 3D representation, corresponding to the evolution of response as a function of two mains parameters out of the three parameters studied here, namely pH, RV and BLD. The two main parameters are selected according to the perturbation graph available in the complementary material. The third parameter is fixed at its central value, for example, value of 11.25 for pH. In such a representation of the surface response, the bottom of the 3D-plot (shown in yellow in the plots) is representing the corresponding 2D iso-response curves. The second representation is a so-called cubic representation, with each axis of the cube corresponding to a parameter ranging from -1 to 1.

The main objective of the 3D and cubic graphs is to evaluate the behavior of the response related to two or three factors and thus to better know about the space of study.

In the case of pH variation response, perturbation graph reveals that BLD is the factor influencing most ln(pH+0.00775), while influence of RV and initial pH on ln(pH+0.00775) are similar. 3D representation of pH response (shown in Highest values for polyphenol or polysaccharide extraction percentages as well as distribution coefficients are obtained when black liquor is highly diluted, that is when a BLD value close to 0.05 -0.066 is used. An opposite result is found for organic acid extraction. In the latter case, a higher BLD value close to 0.15-0.17 yielded highest distribution coefficients.

3D and cubic model representation of extraction

The previous validation for models collected in Table 1 allowed to describe the evolution of each response, namely extraction percentage and distribution coefficient for all three families of compounds over the whole experimental range of values for RV, BLD and pH. Representation of extraction model for organic acids Black liquor principally contains lactic acid as an organic acid.

Representation of lignin residue extraction model

Results indicate that extraction percentage of organic acid varies from 0 to 63 %. According to perturbation graphs, extraction of lactate depends mainly on RV and BLD (Fig. 6-I). Highest values for the extraction of organic acids are located at the borders of the studied experimental field, namely at low and high values for RV and BLD, respectively. Value of 63 % for the extraction yield of organic acids (POA), the highest value obtained here, is found with values for pH, BLD and RV of 13.5, 0.13 and 1, respectively. Furthermore, at a given RV value, extraction yield of organic acids decreases monotonously with BLD. The higher the RV value, the lower the influence of BLD on POA is, the latter becoming nearly negligible when RV is equal to 19. No extraction of organic acid is observed with a low BLD value, either at pH 9 whatever RV, and at high RV whatever pH. This absence of organic acid extraction is represented by the blue area in Fig. 6.I. Plotting the distribution coefficient of organic acids (DOA) as a function of BLD and RV, as shown in Figure 6.III yields a steep and well defined peak centered around the optimal value calculated here, namely 20.7 obtained at values of 9, 0.17

Optimization of multi responses

Because lignin residues might be broken down to smaller polyphenolic units, valuable compounds with known biological activity, currently the most widely used antioxidant compounds in cosmetics for instance, [41] and because a black liquor only containing polysaccharides and organic acids would exhibit a much higher biodegradability and be much less of a hazardous solution, extracting lignin residues from black liquor towards an ionic liquid, while polysaccharides and organic acids remain in the black liquor, was investigated. Furthermore, in order to use black liquor in further biological processes, for example to degrade biologically such liquor, lignin fragments need be removed from BL, while organic acids, and polysaccharides remain in the liquor, as it provides a necessary feed material for living organisms. The objective of this study is thus to obtain an extraction of polyphenolic part of the black liquor as high as possible as opposed to extraction of polysaccharides and organic acids. To that end, optimization of all responses was carried out simultaneously in a so-called multi-response optimization procedure. In such a procedure, values for factors RV, pH and BLD were adjusted in order to yield optimized values for the two sets of responses, namely extraction percentage and distribution coefficients for all extracted compounds.

Under our experimental conditions, two constraints on pH and BLD were initially set. First, for convenience reasons, treating this aqueous phase without acidifying it was preferred. Furthermore, because black liquor exhibits very high concentrations of lignin residues, it is usually diluted several times with water. However, in order to avoid the production of too large volumes of effluent, a low dilution, that is, diluting black liquor 5 times, corresponding to a BLD value of 0.2 was chosen for the optimization process.

Under those conditions and keeping in mind that the objective of the procedure is to obtain an extraction of PP as high as possible, unlike for PS and OA, multi-response optimization procedure was carried out in order to adjust RV factor with set values of 13.5 and 0.2 for pH and BLD, respectively. Simultaneous optimization of all three extraction percentages, namely PPP, PPS and POA was first carried out, followed by a maximization of DPP and minimization of DPS and DOA under the same conditions. Table 2 summarizes the results obtained. Optimum values for extraction percentages and distribution coefficients were obtained at values of 11.6 and 15.6 for RV, respectively. These values for RV indicate that one volume of IL was able to treat efficiently more than 10 times its volume in black liquor. On the one hand, in agreement with our objective, PP was found to be extracted significantly, exhibiting a value of 84.8 % for the extraction percentage, while only a small fraction of organic acids was extracted towards the ionic liquid. On the other hand, lowest extraction value for polysaccharides obtained in this multiresponse optimization procedure, was 66 % for PPS and 31 for DPS, implying a significant extraction of these compounds towards the ionic liquid. Please note that a value of 31 for DPS implies that for an equal volume of black liquor and ionic liquid, the concentration of polysaccharides in ionic liquid will be 31 times higher than in black liquor, corresponding to an extraction of approximately 97 % of PS extracted. This important co-extraction between PP and PS could be explained by the fact that a fraction of PS is linked to PP as hemicellulose fraction, existing in original black liquor. Furthermore, because PP and PS are compounds with high molecular weight, with large hydrophobic domain in the case of PP, they are more easily transferred towards a hydrophobic ionic liquid. Satisfyingly, extraction experiments carried out under the conditions mentioned above, namely, with BLD and pH of 0.2, 13.5, respectively, and an RV value of 11.6 or 15.6, accordingly, gave results very close to the ones calculated here, and in all case within the optimized value range calculated by the software.

In order to optimize further the extraction conditions, for instance in order to obtain a lower valued for PPS, multiple-response optimization procedure was carried out by setting pH and BLD to values of 13.5 and 0.05, respectively. The value of 0.05 for BLD was selected because in agreement with the results detailed in section 3.6.2, values for PPS significantly below 66% are obtained with low BLD values. Results shown in Table 2 confirm that selecting a lower BLD value yields a much lower extraction of polysaccharides, namely %, a result in better agreement with our objectives. Finally, when letting the optimization procedure adjust all three factors, the multiple-response optimization procedure yields values of 94, 47 and 4 % for the extraction yields of PP, PS and OA, respectively, with values of 11.5, 0.05 and 1 for pH, BLD and RV, respectively. This implies that despite a higher extraction of PP, nearly half of polysaccharides are also extracted towards the ionic liquid phase. In addition, RV was adjusted to a value of 1, implying that an equal value of ionic liquid and BL is needed in such a case. This cannot be considered as a favorable condition compared to those presented in Table 2, where RV values are always above 10. These results show that lignin residues can be efficiently removed from a black liquor using an ionic liquid such as [P66614]Cl.

Polysaccharides are also extracted, though to a lesser extent. The amount of polysaccharides extracted can be significantly reduced by a factor of 2, by diluting four times the aqueous phase, that is, by changing BLD from 0.2 to 0.05. Organic acids, unlike lignin residues and polysaccharides, remain mostly quantitatively in the aqueous phase.

Conclusion

In order to propose new ways of recycling wastewater from pulping factory, so-called black liquor, extraction of lignin residues, polysaccharides and organic acids from black liquor, towards an phosphonium-based ionic liquid, [P66614]Cl, was studied. The influence of the initial pH of the black liquor, the dilution in water of black liquor and the volume ratio of black liquor vs. ionic liquid on the extraction of three family of compounds, namely polyphenolic polymers, polysaccharides and organic acids, was considered. It appears that polyphenols are generally extracted along with polysaccharides, unlike organic acids that remain preferably in the aqueous phase.

Results were then modeled according to a design of experiment approach. Response Surface Methodology was used to describe extraction percentage and distribution coefficient of all three compounds. Single responses were successfully optimized, describing accurately experimental data gathered in this work.

In order to provide a black liquor mostly depleted from lignin residues, hence further treatable by simple recycling processes, multi-response optimization procedure was carried out with varying initial conditions. When setting pH, RV and BLD to 13.5, 15.2 and 0.2, respectively, values simplifying the pretreatment of black liquor prior to extraction, 85, 66 and 0 % of lignin residues, polysaccharides and organic acids, respectively, were found to be extracted from the aqueous phase. When black liquor was diluted 20 times, (BLD 0.05), 88, 31 and 1.6 % of PP, PS and OA were extracted, respectively. Finally, highest extraction of lignin residues was possible by reducing the initial pH of black liquor down to 11.25. Good agreement was found between the expected values optimized using a design of experiment and experimental ones.

The work presented here therefore proposes three different ways of pretreatment of black liquor, by removing efficiently polyphenols. Further recycling of resulting black liquor, depending in these different optimized conditions proposed here will be carried out shortly. So far, ionic liquid was successfully reused as an extracting phase four times by contacting it with fresh black liquor, without any treatment for removing PP or PS, without any observable los of extraction efficiency. Finally, ways to recover the lignin fragments within the ionic liquid, for example using contacting the ionic liquid with an acidic aqueous solution, as well as ways to break down lignin part to exploitable polyphenolic units is currently under investigation and should be reported shortly.

Experimental Section

Chemicals

Ionic liquid trihexyl(tetradecyl)phosphonium chloride ([P66614]Cl) was purchased from Solvionic, and used as received. Stock solution of black liquor resulted from a heat-pressure-basic treatment of pine wood (30% NaOH g/g dry wood, 1% anthraquinone, 175°C, 60 min). It was stored in Regarding extraction experiments, diluted samples were prepared from this stock solution using distilled water, with three dilution factors: 0.050, 0.125, and 0.200. Dilution factor will hereafter be named as black liquor dilution (BLD) and is defined as the volume of initial black liquor used and completed to 1 ml of distilled water. A dilution factor of 0.05 of black liquor, therefore, corresponds to 0.05 mL of BL diluted in 1 mL of distilled water. pH was adjusted with diluted sulfuric acid before final dilution.

Analytical methods

Analyses are exclusively carried out on the aqueous phase. pH is estimated with an Electrode Inlab Semi Micro, connected to a pH-meter (Crison GLP21, sensitivity of 0.01 pH). Measurement is performed in a small working volume (200µl). Variation of pH is calculated as the difference between initial pH and final pH measured after extraction by ionic liquid.

Because of an impossibility to quantify separately all residues resulting from the lignin degradation within the Kraft process, concentration of these water soluble lignin residues is measured in accordance to the Lowry-Ciocalteu method, [42] classically used in order to quantify polyphenolic compounds. Therefore, in the following manuscript, these residues will be discussed in terms of polyphenolic equivalent (PP).

Concentration of polysaccharides (PS) is measured using the anthrone method. [43] Optical density is monitored using a double-beam molecular adsorption spectrophotometer (Shimadzu UV-2501PC, Jasco V-530 software).

with Ci, the initial compound concentration (PP, PS or OA) in initial BL phase, Cf the final compound concentration in BL phase after extraction run, and RV, the volume ratio of aqueous phase vs. ionic liquid. In the following, DPP, DPS and DOA correspond to the distribution coefficients for PP, PS and OA respectively. PPP, PPS, and POA refer to the global extraction percentages for PP, PS and OA respectively.

General extraction protocol

All extraction experiments are carried out in sterile polypropylene tubes (2 ml) at room temperature (26±2°C). Black liquor is manipulated using normal micropipettes and the ionic liquid is manipulated with a positivedisplacement micropipette (Microman R). Typically, a sample of 100-1000 µL of aqueous solution is added to a microtube, followed by the ionic liquid.

The microtube is mechanically shaken during 30 s with a Vortex mixer. Emulsion is formed after stirring. It is broken by centrifugation as previously described in the case of L-lactic extraction by phosphonium ionic liquid [44 ] (i.e. centrifugation at 13,000 rpm during 5 minutes). In all extraction experiments, the bottom phase is the aqueous phase because the density of diluted black liquor is 1.096 g.mL-1, a value higher than that of 0.882 g.mL-1 obtained at 25 °C for undiluted [P66614]Cl. [45] An aliquot of the aqueous phase at the bottom of the tube is taken by making a hole in the microtube with a sterile needle and transferred in a new sterile microtube for further analysis.

Experimental Design

Table 3 summarizes the experimental ranges for the three studied continuous factors: dilution factor of black liquor, hereafter named BLD, pH, and volume ratio of liquor vs. ionic liquid, so-called RV. Ranges for the parameter values studied here were selected according to the following general objectives: i) to modify as little as possible the initial black liquor prior to its extraction using an ionic liquids and ii) to reduce the amount of chemical and liquid phases needed for extracting lignin residues, polysaccharides and organic acids. The amount of ionic liquid was thus always considered smaller than that of BL. RV values ranging between 1 and 19 were thus studied. The initial pH value for black liquor is known to be around 13. In order to study the influence of pH without introducing too large quantities of acids, values for the initial pH ranging between 9 to 13.5 were thus studied. Finally, because black liquor is extremely concentrated and is always diluted, BLD range was fixed between 0.2 and 0.05, corresponding to BL diluted 5 and 20 times, respectively.

A Central composite design (CCD) of experiments is applied using the statistical software Design Expert 8.0 (Stat Ease, Minneapolis, USA). A 35-experiment containing matrix is obtained with 7 repeated experiments to determine experimental error (runs 16/22/24/25/27/30/32 of Table S1 of the SI file). Repeated experiments are situated at the center of the studied area. This central location on the design of experiment corresponds to experimental values of 0.125, 11.25 and 10 for BLD, pH and RV, respectively and coded values of 0 within the design of experiment.

For each compound under study, two different final objectives were pursued. First, one may look for a high extraction percentage, no matter the volume of ionic liquid which would be required for this or one may look for a high distribution ratio, a parameter which depends on RV. In order to match these two different goals, six different responses were individually examined: three global extraction percentages (PPP, PPS, and POA) and three distribution coefficients (DPP, DPS, and DOA). Additionally, variation of pH of the aqueous solution after IL extraction was analyzed.

Entry for the Table of Contents

Insert graphic for Table of Contents here. ((Please ensure your graphic is in one of following formats))

Selective removal of polyphenolic polymers, polysaccharides and organic acids present in pulping factory wastewater, so-called black liquor, was studied using an ioniliquid tetrahexylphosphonium chloride. Optimised conditions were obtained using a design of experiment methodology with three variables, namely black liquor pH, dilution factor and water to ionic liquid ratio.

Figure 1 .

 1 Figure 1. Selected snapshots for the extraction of diluted BL using ionic liquid [P66614]Cl. Top phase is [P66614]Cl. Left: vials prior extraction. Right: Vials after vortex and centrifugation, as detailed in the experimental section. BLD: 0.12, pH 9, RV: 10.

Figure 1 .

 1 Figure 1. A: Evolution of final pH as a function of initial pH and RV. •: RV=1. ■ : RV = 10. ○: RV=19.

Figure 3 .

 3 Figure 3. I: 3D representation for |ΔpH| as a function of RV and BLD at the central point value of pH, corresponding to a pH value of 11.25. Blue, green, and red areas represent lower, intermediary and high values of responses, with 25, 50, 75 % of the maximal value, respectively. Delta pH stands for |ΔpH|. Dilution factor of black liquor stands for BLD. II: Cubic representation for |ΔpH| as a function pH, RV and BLD. Values for |ΔpH| given at the cube corners.

Figure 4 summarizes

 4 3D and cubic representations for the extraction percentage (PPP) and distribution coefficient (DPP) of lignin residues, calculated as polyphenolic equivalents. Extraction of lignin residues (PP) was plotted as a function of pH and RV (Fig.4.I and 4-III) because perturbation graph revealed that PPP and DPS depended mostly on these two parameters. Both maximal responses, PPP and DPP, are situated at the limit of the experimental studied field as indicated in Fig.4-I and 4-III. As shown in Fig. 4-I and 4-II, whatever the factors, lignin residues are significantly extracted, as shown by PPP values. Extraction efficiency ranges from 82.1 to 99.6 %, respectively (Fig. 4.II).10.1002/cplu.202200348 Accepted Manuscript ChemPlusChem This article is protected by copyright. All rights reserved. 21926506, ja, Downloaded from https://chemistry-europe.onlinelibrary.wiley.com/doi/10.1002/cplu.202200348 by Universite Cote D'azur, Wiley Online Library on [28/02/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons LicenseHighest value for PPP (99.6 %) is obtained at 9, 0.07 and 6 for pH, BLD and RV values, respectively. Furthermore, PPP generally decreases with pH, BLD, or RV, accordingly. The highest value of 99.7 is very close to that of 99.4 obtained at the lowest pH, BLD and RV values. Distribution coefficient for lignin residues described as polyphenol-equivalent (DPP) varies from 4.8 to 658 and greatly depends on all three parameter value, as observed in the 3D graph and cubic representation. The maximal value of 658 is obtained at values of 9, 0.05 and 14 for pH, BLD and RV, respectively. According to Fig.4-II, DPS generally decreases with pH and BLD.

Figure 4 .

 4 Figure 4. 3D (left side) and cubic (right side) representation for extraction percentage of polyphenol (I-II) and corresponding distribution coefficient (III-IV). Parameter values are indicated in Table 1. Blue, green, and red areas represent lower, intermediary and high values of responses, with 25-50-75 % of the maximal value, respectively.

Figure 5 .

 5 Figure 5. 3D (left side) and cubic (right side) representation for extraction of polysaccharides (I-II) and its distribution coefficient (III-IV). Parameter values are indicated in Table 1. Blue, green, and red areas represent lower, intermediary and high values of responses, with 25, 50 and 75 % of the maximal value, respectively.

  and 11 for pH, BLD, pH and RV, respectively. The peak covers over half of the response surface of Fig.6.III. DOA is close to 0 as long as BLD is lower than approximately 0.11. Cubic representation shown in Fig.6-IV is characterized by presence of null or very low response values in four out of eight corners: at low pH and BLD whatever RV, and at high pH and RV whatever the BLD value. In agreement with Fig.6-III, higher values for the distribution coefficient are located within the cube. It appears that at high pH values, no extraction is found to occur. This can be explained by the fact organic acids are under their basic form at high pH, the latter being well solvated in water, hence poorly extracted towards an ionic liquid phase. 10.1002/cplu.202200348 Accepted Manuscript ChemPlusChem This article is protected by copyright. All rights reserved.21926506, ja, Downloaded from https://chemistry-europe.onlinelibrary.wiley.com/doi/10.1002/cplu.202200348 by Universite Cote D'azur, Wiley Online Library on [28/02/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Figure 6 .

 6 Figure 6. 3D (left side) and cubic (right side) representation for extraction of lactate (I-II) and its distribution coefficient (III-IV). Parameter values are indicated in Table 1. Blue, green, and red areas represents lower, intermediary and high values of responses, with 25-50-25 % of the maximal value, respectively.

  protected by copyright. All rights reserved. 21926506, ja, Downloaded from https://chemistry-europe.onlinelibrary.wiley.com/doi/10.1002/cplu.202200348 by Universite Cote D'azur, Wiley Online Library on [28/02/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

  protected by copyright. All rights reserved. 21926506, ja, Downloaded from https://chemistry-europe.onlinelibrary.wiley.com/doi/10.1002/cplu.202200348 by Universite Cote D'azur, Wiley Online Library on [28/02/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License a fridge at low temperature (4-6°C).

  protected by copyright. All rights reserved. 21926506, ja, Downloaded from https://chemistry-europe.onlinelibrary.wiley.com/doi/10.1002/cplu.202200348 by Universite Cote D'azur, Wiley Online Library on [28/02/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Table 1 .

 1 Optimized experimental factors for each individual response, leading to individual maximized responses in the studied experimental area.

	and

3.I), at the central point value

(pH 11.25)

, indicates that variation of pH is only significant when RV and BLD exhibit low values. Maximal value of 5.5 for |pH| is obtained when RV and BLD exhibit values of 1 individually or by coupling responses in the case of multiresponse analysis. In order to obtain the optimum value for a given response, such as extraction percentage or distribution coefficient, each factor, namely RV, BLD and pH, was individually optimized in its value range. Different extraction goals can be established in the case of a solution containing multiple compounds to extract or remove, such as those found in a black liquor solution. Because lignin residues within black liquor present a certain toxicity and are hardly biodegradable at the concentration they are in, and because they can yield valuable compounds, the main goal of this work was thus to remove polyphenols while keeping polysaccharides and organic in solution in the diluted black liquor. In this study, both goals were studied either individually or in multi-response analysis. Optimized experimental conditions are summarized in Table

1

.

[a] Low and high values of the 95% confidence interval in which one could expect to find the average response of a group of runs.

Validation of each optimized model, as detailed in Table

1

, was carried out by comparing maximal responses to experimental results carried out under optimized conditions, as shown in Table

1

. It appears that experimental values are close to mathematical values, that straightforwardly validates each model. Extraction percentage or distribution coefficient can thus be predicted by each model in the experimental field studied by the DOE.

Table 2 .

 2 Optimization of multi-responses for extraction percentage and distribution factor with fixed values for both BLD and initial pH

	Responses		Optimized	Optimized	Exp.
			value	value range	value
	pH = 13.5; BLD = 0.2	RV = 11.6
	PPP	Maximize	84.8	82.1 -99.5	88.2
	PPS	Minimize	66.0	22.0 -79.1	72.3
	POA	Minimize	16.7	0 -68.0	0
	pH = 13.5; BLD = 0.2	RV = 15.6
	DPP	Maximize	104	34 -317	125
	DPS	Minimize	31.5	13.3 -74.7	22.5
	DOA	Minimize	0.14	0 -7.34	0
	pH = 13.5; BLD = 0.05	RV = 15.6

Table 3 .

 3 Chemical and physical factor ranges studied in the Central Composite Design (CCD).

	Factor	Name	Min.	Max.	Coded	Coded	Coded
					values	values	values
	A	BLD	0.050	0.200	-1 =0.050	0 =0.125	1 =0.200
	B	pH	9.00	13.50	-1 =9.00	0 =11.25	1 =13.50
	C	RV	1	19	-1 =1	0 =10	1 =19

Acknowledgements 10.1002/cplu.202200348 We thank Pr. G. Mortha (Laboratory of Paper Processes, LGP2, UMR CNRS 5518, Grenoble, France) for providing the wood issued black liquor. Keywords: response surface methodology • ionic liquid • black liquor • polyphenols • polysaccharides

Accepted Manuscript

ChemPlusChem This article is protected by copyright. All rights reserved.