Gautier Viaud
email: gautier.viaud@ecp.fr

Yuting Chen

Benoît Bayol

Paul-Henry Cournède

A Comparison of a Generic MCMC-Based Algorithm for Bayesian Estimation in C++, R and Julia. Application to Plant Growth Modeling

Plant growth is understood through the use of dynamical systems involving many interacting processes and model parameters whose estimation is therefore a crucial issue, all the more so since experimental data obtained from agronomical systems are most of the time characterized by their scarcity and their heterogeneity owing to the complex underlying acquisition. One of the approaches used to solve this kind of problem within a Bayesian paradigm involves Markov Chain Monte Carlo (MCMC) algorithms, one of the drawbacks of the latter being that they lead to some intensive computation because of the statistical framework employed, which is why the efficiency of the implemented computing methods is of particular importance. In this paper, we compare three implementations of a generic MCMCbased algorithm for Bayesian estimation in C++, R and Julia so as to compare the performance and precision of these languages. Here, genericness means that the estimation algorithm can be used for any dynamic model provided that it is implemented in a given modeling template. Such genericness is of crucial importance in a scientific field such as plant growth modeling for which no reference model exists and new models are constantly developed and evaluated. The tests are conducted for the particular cases of Lotka-Volterra model and the Log-Normal Allocation and Senescence model for sugar beet.

.

Generic probability distributions

Parameter estimation algorithms should be designed so as to be easily used with different models, which is why some probability density functions (pdf), necessary for this type of algorithms, need to be expressed in a generic form. In

Introduction

Plant growth modeling has a wide range of potential applications for decision aid in farm management, the most common one being predicting yields in fields at a global scale. Lots of models, whether they be designed at the individual level, see [START_REF] De Reffye | Plant growth models[END_REF] for a review, or the field level (such as the LNAS [START_REF] Chen | Data assimilation to reduce uncertainty of crop model prediction with convolution particle filtering[END_REF] or STICS [START_REF] Brisson | Conceptual basis, formalisations and parameterization of the STICS crop model[END_REF] models), have been developed over the past decades. They can make use of a few or several dozens of parameters, and accurate predictions require correct parametrization of the models -which, because of the complexity and the interactions of the underlying biological processes, can be a difficult task. Two approaches can be chosen for parameter estimation, the frequentist one (with the classical maximum likelihood estimator) and the Bayesian one. In Bayesian estimation, parameters are considered as a random vector whose posterior distribution is deduced from the knowledge of system observations and a prior distribution for the parameters. It offers some interesting perspectives in mechanistic plant growth modeling since some reasonable a priori distributions on the parameters can be deduced from the literature or from the well-known biological processes. The computation of the posterior distribution is generally untractable analytically, and two main families of algorithms can be used for its approximation, Sequential Monte-Carlo (SMC) methods and Markov Chain Monte Carlo (MCMC) methods. In this article, we concentrate on an MCMC-based algorithm for Bayesian estimation and address two main issues often encountered in plant growth models parameter estimation: (i) design both plant growth models and parameter estimation algorithms in a generic form so that each algorithm can be used in concordance with each model and (ii) use the programming language that is the most apropriate for these algorithms that is, ideally, the least time-and memory-consuming and the 1 easiest for writing code.

So far, parameter estimation algorithms were written for a specific model, whence the proliferation of algorithms identical in nature and the rewriting of code. To avoid this issue, we reformulate the theoretical framework of models so as to be able to apply a generic algorithm for different template models. As far as the second issue is concerned, plant growth models describe complex biological processes, generally on a large number of time steps, and a single simulation of the model is therefore by itself already time-consuming. Consequently, it is not surprising within plant growth modeling communities to have, when resorting to analysis or estimation algorithms that require a large number of simulations, computation times up to several hours, sometimes several tens of hours, and code optimization is thus a critical point. This is why we compared the implementation of a generic MCMC algorithm, of primary importance for parameter estimation in plant growth models, in different programming languages: C++ [START_REF]C++ website[END_REF], a common choice for numerical methods when it comes to efficiency, R [3], a long-standing and very used language within the statistics community and with a wide set of libraries, and Julia [START_REF] Cournède | [END_REF], a rather young language developed only a few years ago but whose community is growing rapidly.

We begin by defining in Section 1 the theoretical framework for the genericness of models that will allow parameter estimation algorithms to be applied to different models easily. In Section 2, we recall the basic principles of MCMC algorithms for parameter estimation within a Bayesian perspective and the slight adaptations necessary to deal with the estimation of both parameters and hidden states. A predator-prey model and a plant growth model that serve as test models for MCMC simulations are then presented in Section 3 and the issues related to the programming side and the specificities of the different languages considered are discussed in Section 4. Finally, in Section 5, we present the results of the efficiency for each language in terms of process time and memory.

1 Generic statistical model 1.1 General state space models Plant growth models are often formulated in terms of discrete nonlinear state space models [START_REF] Doucet | Sequential Monte Carlo methods in practice[END_REF]. The plant characteristics are updated at each time step n ∈ 1, N with N the last step of the simulation, where it is made use of a transition function f n and an observation function g n as follows:

X n+1 = f n (X n , E n , Θ, η n ), Y n+1 = g n (X n+1 , Θ, ξ n ),
and where at time n, X n ∈ R d X represents the state variables and Y n ∈ R d Y the observations that are available. The environment influence is taken into account through the variables E n ∈ R d E , and Θ ∈ R d Θ denotes the functional parameters of the model. Two types of noises are considered, the modeling noises η n ∈ R dη and the observation noises ξ n ∈ R d ξ . The former are introduced so as to take into consideration possible model imperfections or environmental stochasticity whereas the latter are because experimental data result from field measurements, which are always a source of errors. The experimental data Y n comprise state variables that are linked to the hidden variables X n via the observation function g n but can differ from them. Owing to practical difficulties to obtain experimental data, these may not be available at all times, especially when dealing with biological systems. The set of the O measurement times is denoted by O = {t i , i ∈ 1, O } and the following notations are introduced: the case of the MCMC algorithm used to run the simulations, these are the pdf of the observations conditional to the parameters and the hidden states p(Y 1→N |Θ, X 1:N ), the transition pdf p(X n+1 |Θ, X n ) and the observation pdf p(Y n+1 |Θ, X n+1 ). For the n-th step of the simulation, the hidden state variables are written as X n = {X i n , i ∈ 1, d X }, the modeling noises as η n = {η i n , i ∈ 1, d η } and the observation noises

ξ n = {ξ i n , i ∈ 1, d ξ }. Let m η : 1, d η → 1, d X and m ξ : 1, d ξ → 1, d X
be applications such that m η ( 1, d η ) and m ξ ( 1, d ξ ) represent the sets of indexes of the state variables on which are set the modeling and observation noises respectively.

X mη(i)+1 n = f (X mη(i) n , η i n ), i ∈ 1, d η , Y i n = f (X m ξ (i) n , ξ i n ), i ∈ 1, d ξ .
In the case of multiplicative noises, this would translate into:

X mη(i)+1 n = X mη(i) n (1 + η i n ), i ∈ 1, d η , Y i n = X m ξ (i) n (1 + ξ i n ), i ∈ 1, d ξ .
It is possible to express the transition pdf by choosing only d η variables from the state variables whose indices are represented by q : 1, d η → 1, d X as:

p(X n+1 |Θ, X n ) = dη i=1 p(X q(i) n+1 |Θ, X 1:d X n , X 1:q(i)-1 n+1 ) × j / ∈q( 1,dη ) δ(X j , m j (X 1:N n , X 1:j-1 n+1 ))
where m j (X 1:N n , X

1:j-1 n+1
) is the value computed for the state variable X j within the model considered. In what follows, since the variables that are deterministic functions of the stochastic variables are computed directly by closure relationships in the simulation program, the Dirac distributions δ(•, •) take values 1, and hence will be omitted in the following. In particular, we can restrain ourselves to the "noised" variables and rewrite the transition pdf as:

p(X n+1 |X n , Θ) = p(η n |Θ) = dη i=1 p(X mη(i)+1 n+1 |Θ, X mη(i) n+1 ). (1)
The observation pdf can be expressed in the same way as:

p(Y n+1 |Θ, X n+1 ) = p(ξ n |Θ) = d ξ i=1 p(Y i n+1 |X m ξ (i) n+1 ). (2)
The generic expression of the pdf of the observations conditional to the parameters and the hidden states p(Y 1→N |Θ, X 1:N ) naturally follows from that of the observation pdf since:

p(Y 1→N |Θ, X 1:N ) = O k=1 p(Y t k |Θ, X t k ) = O k=1 d ξ i=1 p(Y i t k |Θ, X m ξ (i) t k
).

Equation (1) makes it possible to compute the transition pdf as long as are specified the nature of the noises (e.g. additive/multiplicative normal, uniform) and lists of labels corresponding to the parameters necessary to compute the values of these pdfs, such as, for the case of multiplicative normal noises, the modeling noise parameters {σ ηi , i ∈ 1, d η }, the deterministic states {X

mη(i) n+1 , i ∈ 1, d η } and the corresponding stochastic states {X mη(i)+1 n+1 , i ∈ 1, d η }. Likewise, Equation (2
) makes it possible to compute the observation probability density where this time lists of labels must include the observation noise parameters {σ ξi , i ∈ 1, d ξ }, that of the observation states {Y

(i) n , i ∈ 1, d ξ } and that of the corresponding hidden states {X m ξ (i) n , i ∈ 1, d ξ }.
In Section 3, these functions are explicitly written for the case of the LNAS model.

MCMC algorithms for Bayesian estimation

MCMC algorithms have attracted a lot of attention [START_REF] Gilks | Markov chain monte carlo[END_REF] [5] for the analysis of complex statistical models in a Bayesian perspective since they allow for the evaluation of complex and high-dimensional integrals, whence for the estimation of parameters and hidden states. The general principle of MCMC algorithms is to generate a sequence of samples from a probability distribution from which it is not possible to sample directly. This sequence represents an approximation of the target distribution and can be used for the computation of integrals, for instance to estimate the expectation or the standard deviation of parameters. The procedure goes as follows: at each iteration of the algorithm, a sample from a so-called proposal distribution is generated; this sample is either accepted or discarded based on an acceptance ratio which depends on both the actual candidate and the last accepted one.

In plant growth models, we are not only interested in estimating some of the parameters of the model Θ but also the hidden states X 1:N . A joint estimation of (Θ, X 1:N ) is therefore required, according to an update scheme that can be iteratively repeated:

Step

1: Update Θ ∼ p(Θ|X 1:N , Y 1→N ) Step 2: Update X 1:N ∼ p(X 1:N |Θ, Y 1→N )
since there is a high correlation between the parameters and the hidden states. The motivation for the second step is to improve the estimation of the hidden states once the estimation of the parameters have converged. In what follows, the Adaptive Metropolis-Within-Gibbs (AMWG) algorithm [START_REF] Chen | Bayesian Inference in Plant Growth Models for Prediction and Uncertainty Assessment[END_REF] is presented for parameter estimation (Step 1) only.

Step 1: estimating functional parameters

The target distribution is defined as p(Θ, X 1:N |Y 1→N ). The principle of AMWG is to generate an ergodic Markov chain {(Θ (i) , X

:N ), i ∈ 1, N p } whose stationary distribution is precisely the target distribution p(Θ, X 1:N |Y 1→N ) (i) 1
, where N p is the maximum number of iterations or chain length. At the i-th iteration, a candidate (Θ ( ) , X ( )

1:N ) is generated by sampling from a proposal distribution q(Θ ( ) , X ( ) 1:N |Θ (i) , X (i) 1:N ).
As proposed by [START_REF] Fearnhead | MCMC for state-space models[END_REF], we decomposed this proposal distribution into two parts, q(Θ ( ) , X ( )

1:N |Θ (i) , X (i) 1:N ) = q(Θ ( ) |Θ (i) )q(X ( ) 1:N |Θ ( ) ).
This means that a new candidate for the parameters has first to be sampled, and then the state list X 1:N must be generated through the model considered with this set of parameters, i.e. we take q(X ( )

1:N |Θ ( ) ) = p(X ( ) 1:N |Θ ( ) )
where p(X ( ) 1:N |Θ ( ) ) is the pdf corresponding to the model simulation function. The candidate (Θ ( ) , X ( ) 1:N ) is then accepted with probability min(1, α) where

α = p(Θ ( ) , X ( ) 1:N |Y 1→N ) p(Θ (i) , X (i) 1:N |Y 1→N ) × q(Θ (i) , X (i) 1:N |Θ ( ) , X ( ) 1:N ) q(Θ ( ) , X ( ) 1:N |Θ (i) , X (i) 1:N ) .
Accepting the samples with such a ratio ensures that the stationary distribution of the Markov chain will be the target distribution. A common strategy for the sampling of the parameters is random walk: the proposal distribution is therefore taken to be such that q(Θ |Θ) = q RW (Θ -Θ) where q RW is symmetric. A classical choice when it comes to multidimensional parameter estimation is to choose a multivariate normal distribution with zero mean and a covariance matrix Σ (i) for q RW so that

Θ ( ) ∼ N (µ (i) , Σ (i) )
where µ (i) is set to the last accepted candidate Θ (i) . With such a choice and after some rearrangements, the acceptance ratio can therefore be rewritten in a more compact form

α = p(Y 1→N |Θ ( ) , X ( ) 1:N ) p(Y 1→N |Θ (i) , X (i) 1:N ) × p(Θ ( ) ) p(Θ (i) )
where p(Θ) is the prior distribution chosen for the estimated parameters.

The choice of the covariance matrix is of primary importance: if the variance of the proposal distribution is too small, most of the proposed values will be accepted and the difference between two accepted values tend to be rather small, resulting in a slow exploration of the state space. On the other hand, if the variance of the proposal distribution is too large, the candidates will often be rejected and the chain can be stationary for a long time. Adaptive algorithms [START_REF] Haario | An adaptive Metropolis algorithm[END_REF] have been proposed to avoid having to manually change the proposal distribution and employ the following update scheme     

∆ (i) = Θ (i+1) -µ (i) µ (i+1) = µ (i) + γ (i) ∆ (i) Σ (i+1) = Σ (i) + γ (i) [∆ (i) ∆ (i)T -Σ (i) ]
where γ (i) = 1/i is defined so as to make the variations introduced by this adaptive scheme vanish and ensure the ergodicity of the chain. Finally, the covariance matrix is scaled by a factor λ (i) that assures a proper exploration of the target distribution even when the covariance matrix is not well initialized

λ (i+1) = λ (i) exp(γ (i) [α (i) -α ]) Σ (i+1) = λ (i+1) Σ (i+1)
where α is the target acceptance rate. The algorithm is stopped either when the maximum number of iterations is reached or when a specific stopping criterion is verified. Several stopping criteria can be used, they can be based on mean estimates or Monte-Carlo errors, see [START_REF] Chen | Bayesian Inference in Plant Growth Models for Prediction and Uncertainty Assessment[END_REF] for more details. Once the algorithm has stopped, the estimated value of the parameters are computed by averaging. However, it is probable that the stationary distribution of the chain is not reached during the first iterations of the algorithm, which is why a burn-in parameter N bp corresponding to the number of iterations discarded is introduced. Finally, the estimated parameters is computed as

Θ = 1 N p -N bp Np i=N bp +1 Θ (i) .
All this algorithm needs in order to work in a generic way is the computation of the three pdfs presented in Section 1 and the function that generates the states from time step 1 to time step N for a given model M. The different pdfs can be easily computed by specifying the nature of the modeling and observation noises introduced and the names of the variables affected in model M. It is potentially available for the user to specify its own pdf if the standard ones (normal and uniform as of now) do not suffice. Plant growth models are designed outside the algorithm and any can be called to generate a list of states X 1:N with a given set of parameters.

On two models

To illustrate the genericness of our approach, we present two examples of the class of models considered in Section 1.

Model 1: Lotka-Volterra (LV)

The Lotka-Volterra model deals with the evolution of the population of two species, one being a predator and the other a prey. The discrete version we considered reads

X n+1 = N 1 n+1 N 2 n+1 = (1 + a)N 1 n -bN 1 n N 2 n (1 -c)N 2 n + dN 1 n N 2 n
where N 1 n and N 2 n are the number of preys and the number of predators at time n respectively, and Θ = (a, b, c, d) ∈ R 4 + represents the parameters of the model describing the interaction between the two species. The populations N 1 n and N 2 n are assumed to be observed at time n with measurement noises 1 ∼ N (0, (σ 1 ) 2 ) and 2 ∼ N (0, (σ 2 ) 2 ), taken to be additive normal noises,

Y n = Y 1 n Y 2 n = N 1 n + 1 n N 2 n + 2 n Algorithm 1 Adaptive Metropolis-Within-Gibbs Initialization Initialize Θ (0) , X (0)
1:N . Initialize µ (0) , Σ (0) , λ (0) . Parameters estimation Compute the ratio

1: for i = 0 : N p -1 do 2: Sample Θ ( ) ∼ N (Θ (i) , λ (i) Σ (i) ).
α = min 1, p(Y 1→N |Θ ( ) , X ( ) 1:N ) p(Y 1→N |Θ (i) , X (i) 1:N ) p(Θ ( ) ) p(Θ (i) )
.

5:

Sample u ∼ U (0, 1).

6:

if α < u then 7:

Set (Θ (i+1) , X Set (Θ (i+1) , X

1:N ) = (Θ ( ) , X ( ) 1:N ). (i+1)
:N ) = (Θ (i) , X (i+1) 1
1:N ).

10:

Update RW mean vector and covariance matrix

µ (i+1) = µ (i) + γ (i) (Θ (i+1) -µ (i) , Σ (i+1) = Σ (i) + γ (i) [(Θ (i+1) -µ (i) )(Θ (i+1) -µ (i) ) T -Σ (i) ]. 11: Set Θ = 1 N p -N bp Np j=N bp +1 Θ (j) . so that p(Y n+1 |Θ, X n+1 ) = 1 σ 1 √ 2π e (Y 1 n+1 -N 1 n+1 ) 2 2(σ 1 ) 2 × 1 σ 2 √ 2π e (Y 2 n+1 -N 2 n+1 ) 2 2(σ 2 ) 2 .

Model 2: Log-Normal Allocation and Senescence (LNAS)

This model concerns the growth of sugar beet [START_REF] Chen | Data assimilation to reduce uncertainty of crop model prediction with convolution particle filtering[END_REF]. The latter is considered to be made of two main compartments, the root and the foliage. At day n, the root biomass is denoted Q r n and that of the foliage is

Q f n = Q g n + Q s n ,
where Q g n is the biomass of the green leaves and Q s n is the biomass of the senescent leaves. The unit of all biomasses is g.m -2 . The model is discretized according to a daily time step and the environmental variables of day n such as the temperature T n ( • C) and the photosynthetically active radiation PAR n (MJ .m -2 ) are daily averages. To emphasize the variables on which are set modeling noises, deterministic variables are denoted with a superscript det whereas their stochastic equivalent are denoted with a superscript sto. The growth of the plant is led by the thermal time, which represents the accumulated daily temperature above a certain threshold -taken as T b = 0 • C in the case of sugar beet -since germination, corresponding to day i = 0, of the plant:

τ n = n i=0 (T i -T b ).
The thermal time must be higher than a certain value τ init for the plant to reach the emergence stage. When the latter is attained, the plant starts to intercept light and as a consequence to produce biomass through photosynthesis. The biomass produced at day n per unit surface area is denoted by Q det n . It is assumed to follow a Beer-Lambert law:

Q det n = PAR n × µ × (1 -exp(-λ × Q g n )),
where µ (g.MJ -1 ) is an efficiency coefficient, λ (g -1 .m 2 ) a coefficient parameter and (1 -exp(-λ × Q g n ) represents the fraction of radiation intercepted by the foliage. To account for some inaccuracies of the Beer-Lambert law, the variable production of biomass is rendered stochastic by multiplying its deterministic value by a multiplicative normal noise such that

Q sto n = Q det n × (1 + η Q n ),
where η

Q n ∼ N (0, σ Q ).
The biomass produced at day n is distributed between the foliage and root system compartments according to an empirical function γ whose deterministic value is given by

γ det n = γ 0 + (γ f -γ 0 ) × G a (τ n ),
where γ 0 ≥ 0 and γ f ≤ 1 are respectively the initial value and the limit when the thermal time goes to infinity and G a is the cumulative distribution of a log-normal law parameterized by its median µ a and its standard deviation σ a . As for the production of biomass, a modelling noise for the allocation is introduced since this strategy highly depends on environmental conditions. Once again, a multiplicative normal noise is chosen and

γ sto n = γ det n × (1 + η γ n ),
with η γ n ∼ N (0, σ γ ). The biomass of the senescent foliage at day n is denoted by Q s n , it is a proportion of the accumulated foliage biomass making use of the cumulative distribution of a log-normal law parametrized by its median µ s and its standard deviation s s . This process begins once the thermal time has reached a certain threshold τ s . The proportion of the foliage at day n becoming senescent is therefore G s (τ n -τ s )Q f n . The biomass of the whole foliage increases every day: it receives the proportion γ ∈ [0, 1] of the biomass produced on day n,

Q f n = Q f n-1 + γ sto n × Q sto n ,
The biomass of senescent leaves is calculated as a proportion ρ s of the foliage biomass,

Q s n = ρ s n × Q f n ,
from which the biomass of green leaves can be deduced,

Q g n = Q f n -Q s n = (1 -ρ s n ) × Q f n .
Finally, the biomass of the root is increased by what is not allocated to the foliage,

Q r n = Q r n-1 + (1 -γ sto n ) × Q sto n .
For some days, it is assumed that the green leaves biomass Q g and the root biomass Q r are observed with respective measurement noises g and r that are assumed to be multiplicative normal noises. The experimental data potentially available at day n can therefore be written

Y n = Y g n Y r n = Q g n × (1 + g n ) Q r n × (1 + r n ) . Since Q sto n ∼ N (Q det n , (σ Q Q det n ) 2 ) γ sto n ∼ N (γ det n , (σ γ γ det n ) 2
) the transition pdf can be written as

p(X n+1 |Θ, X n ) = 1 σ Q Q det n+1 √ 2π e (Q sto n+1 -Q det n+1 ) 2 2(σ Q Q det n+1 ) 2 × 1 σ γ γ det n+1 √ 2π e (γ sto n+1 -γ det n+1 ) 2 2(σ γ γ det n+1 ) 2
.

In concrete terms, we have to specify into the configuration of the MCMC algorithm before it be run that modeling noises will be multiplicative normal noises as well as the names of the variables affected, i.e. Q det and Q sto with noise standard deviation σ Q on the one hand and γ det and γ sto with noise standard deviation σ γ on the other. A generic function then takes care of computing the value of p(X n+1 |Θ,

X n ). Since Y g n ∼ N (Q g n , σ g Q g n ) and Y r n ∼ N (Q r n , σ r Q r n
), the observation pdf can be written as

p(Y n+1 |Θ, X n+1 ) = 1 σ g Q g n+1 √ 2π e (Y g n+1 -Q g n+1 ) 2 2(σ g Q g n+1 ) 2 × 1 σ r Q r n+1 √ 2π e (Y r n+1 -Q r n+1 ) 2 2(σ r Q r n+1 ) 2
and a similar procedure is used to automatically compute it.

Three different languages

C++ is a procedural, object-oriented, generic programming language. It is general purpose and is used intensively in the industrial and academic fields when it comes to numerical programming since it can deliver high performance in terms of speed and memory consumption. However, the language can be quite complex for writing scientific code; in particular, the standard library does not provide any regular function used in modeling and statistics.

R is a dynamic, lazy, functional, object-oriented language. It is probably the most comprehensive domain specific language for statistical analysis, developed by senior statisticians and researchers and used extensively within the statistics community, which provides an undoubted support. It comes with a platform that allows the loading of many packages, more than 6000 as of October 2014, which can be standard statistical tests, models or data sets. It has immense graphical capabilities, which outperform those of most other statistical languages. For all these reasons, it is very convenient for developing ideas and rapid testing. However, the consumption of time and memory cannot compare to that of C++ or other mainstream languages.

Julia is a fairly new language since it appeared in 2012, [7] [6]. R and Julia are both dynamic languages, which are often considered as highly productive but, as is the case of R, not very efficient. Julia was designed with this idea in mind [START_REF] Cournède | [END_REF] and it was taken advantage of modern techniques such as a rich type information, aggressive code specialization against run-time types and Just-In-Time (JIT) compilation using the LLVM compiler which allows to create functions on the fly. As a result, Julia's performance is claimed to compare to that of statically compiled languages such as C++ while being an interactive dynamic language just like R, providing flexibility and coding productivity.

Possessing this unusual combination of aspects, Julia is a very attractive programming language for plant growth models communities, whence the desire to test its performance compared to those of C++ and R. We implemented the generic MCMC algorithm presented in Section 2 in C++, R and Julia. For researchers unfamiliar with numerical methods, the implementation in R and Julia is much more intuitive than in C++. It must be noted that, since by default R passes all the arguments to functions by copy, we had to embed most of the objects (states, parameters, etc.) inside environment objects to avoid unnecessary copies. On the contrary, variables are always passed by reference in Julia and easily can be in C++, which make them much simpler to design efficient algorithms.

Results

We compared the performance of the different languages on two test cases using the LV and LNAS models.

Test case 1: LV model

We chose a set of parameters (a, b, c, d) = (0.1, 0.02, 0.15, 0.03)

and we generated experimental data Y 1 1:50 and Y 2 1:50 from these values. We then ran the AMWG algorithm with priors

      
a ∼ N (0.12, 0.012), b ∼ N (0.018, 0.0018), c ∼ N (0.17, 0.017) d ∼ N (0.028, 0.0028).

As can be seen from Table 1, the values of the estimated parameters agree with the ones used for simulating the experimental data set.

Test case 2: LNAS model

The environment variables available were the daily average temperature and the photosynthetically active radiation over 160 days. at which were given the values of the biomasses of the green foliage Q g and of the root Q r . We estimated three parameters of the model: the radiation use efficiency µ, the median of the allocation log-normal distribution µ γ and the initial value of the allocation strategy γ 0 . All the other parameters were set to a constant value throughout the simulations. We used as prior distributions for the estimated parameters    µ ∼ N (3.6, 0.1), µ γ ∼ N (0.75, 0.08), γ 0 ∼ N (600.0, 20.0).

We chose this particular test case to ensure that the algorithms were correctly written since the results are known for a non-generic MCMC algorithm specifically designed for the LNAS model within the PyGMAlion platform of the Digiplante team.

General remarks

Because the main concern of these simulations was to compare the efficiency of different programming languages, the same fixed number of iterations was set for each language. We ran N p = 350000 iterations with a burn-in period of N bp = 50000 iterations. The estimated parameters for these test cases are shown in Tables 1 and2 Obviously, there must still be room for improvement in each language as we designed these implementations as researchers who would have common knowledge of these languages. It was envisioned that R had little chance to compete with the other two languages in terms of time, but not with such high ratios: as far as test case 2 is concerned, on Machine 1, R takes roughly 3,700 more time than C++ whereas Julia takes a bit more than twice, which is a decent achievement for a dynamic language compared to a static one. It has to be noted that C++ is advantaged since the compile time is not taken into account. The results obtained on the computing cluster are somewhat identical, except for Julia's performance. A possible explanation for this is that the binary version of Julia was used here, since it was not possible to compile it properly. Rather surprisingly, the ratio of times between test case 2 and test case 1 is much higher for Julia than for C++. This probably arises because of some operations that cannot be time-compressed, such as loading packages, including files, JIT compilation of functions, etc. In R's defence, MCMC-based algorithms are not the best-case scenario for this language since no matrix computation is performed. In terms of memory, once again, C++ is the more efficient, R achieves a decent score and Julia seems a bit more greedy. A last point, rather subjective, has to be mentioned: the time required for the development of the algorithms were much less in R and Julia than in C++. All things considered, whereas R is not a conceivable option because of the process time needed, Julia offers a respectable alternative to C++, admittedly a bit more time-and memory-consuming, but with a higher expressiveness.

Discussion and perspectives

In this paper, we presented an approach to defining general state space models in a generic form, allowing to design algorithms -such as parameter estimation algorithms -making use of this genericness and which would not require to be written for a specific model.

The Digiplante team built a C++ platform called PyGMAlion based on this principle of genericness [START_REF] Cournede | Development and evaluation of plant growth models: methodology and implementation in the pygmalion platform[END_REF]. It allows the user to easily define a dynamic model comprising structures such as state variables, environmental variables, parameters and observations. Some statistics tools for parameter estimation, sensitivity analysis, data assimilation and uncertainty analysis are available and the possibility to use generic complex algorithms such as MCMC or Convolution Particle Filters (CPF), see [START_REF] Campillo | Convolution Particle Filter for Parameter Estimation in General State-Space Models[END_REF], [START_REF] Chen | Data assimilation to reduce uncertainty of crop model prediction with convolution particle filtering[END_REF], algorithms for different plant growth models is currently under work.

One of the main issues for parameter estimation within plant growth models being numerical computing performance, we compared the memory and time consumptions for different programming languages which highlighted a very promising future for this rather new dynamic language that Julia is: with performances comparable to that of C++, being probably easier for statisticians to code with and having rapidly growing community and package libraries, it may just be a matter of time before it gets more attention from the statistics researchers community. It would be interesting to conduct new tests on other parameter estimation algorithms such as CPF or Particle Markov Chain Monte-Carlo (PMCMC), see [START_REF] Andrieu | Particle Markov chain Monte Carlo methods[END_REF], all the more so since these can be parallelized, which would be a good opportunity to test the parallel computing capabilities of Julia.

  through model M.

  The experimental data was made of 14 observations at days

	O = {54, 68, 76, 83, 90, 98, 104, 110, 118, 125, 132, 139, 145, 160}

Table 1 :

 1 . For each language, Final values of the estimated parameters (means and variances) for test case 1.

		a	b	c	d
	Means 9.99 e-02	1.99 e-02	1.50 e-01	3.00 e-02
	Vars	3.41 e-09	2.09 e-09	1.97 e-08	2.54 e-09
			µ	µ γ	γ 0
	Means		3.75 e-00	7.52 e-01	5.73 e+02
	Vars		3.54 e-03	3.10 e-04	3.08 e+02

Table 2 :

 2 Final values of the estimated parameters (means and variances) for test case 2.we measured the time spent with the Unix command time, which gives three different values: real is the wall clock elapsed time, user is the amount of CPU time spent executing the program, sys is the CPU time spent executing operating system services on behalf of the program. The tests were conducted 100 times and the mean total time as well as its standard deviation were computed. It has to be noted that the results of the mean and of the variance indicated for 100 runs in R are based on linear interpolation of the time needed for 1000 iterations, which is still consistent with the time obtained when the full algorithm was run for time issues. We performed these tests both on a personal machine with Intel Core i7-3687U CPU 2.10GHz × 4, 8 GiB of memory on a 64-bit Linux operating system, denoted Machine 1, and on the computing cluster of CentraleSupélec, with CPU Intel Xeon E5-2695 V2 (2.4Ghz), denoted Machine 2. The results for test case 2 (by far the most time-consuming) are presented in Tables3 and 4.

		C++	R	Julia
	Test case 1
	Machine 1
	real	2.87 e+00	3.37 e+04	5.30 e+01
	usr	2.74 e+00	3.76 e+04	5.29 e+01
	sys	1.02 e-02	1.52 e+01	4.36 e-01
	Machine 2
	real	2.38 e+00	3.54 e+04	5.92 e+01
	usr	2.24 e+00	3.95 e+04	5.78 e+01
	sys	8.20 e-03	2.53 e+01	1.50 e+00
	Test case 2
	Machine 1
	real	2.73 e+01	1.02 e+05	6.48 e+01
	usr	2.73 e+01	1.02 e+05	6.45 e+01
	sys	1.28 e-02	4.62 e+01	2.12 e-01
	Machine 2
	real	2.51 e+01	1.10 e+05	8.34 e+01
	usr	2.51 e+01	1.10 e+05	8.33 e+01
	sys	1.01 e-02	7.05 e+01	1.56 e+00

Table 3 :

 3 Times (in seconds) needed for the algorithm in the different languages to end (average based on 100 runs).

		C++	R	Julia
	Test case 2
	Machine 2
	mem	1.99 e+01	5.32 e+01	1.68 e+02

Table 4 :

 4 Memory (in MB) used by the algorithms in the different languages (average based on 100 runs).


Indra Bahadur

Khalid Osman

Christophe Coquelet

Paramespri Naidoo

Deresh Ramjugernath

Solubilities of Carbon Dioxide and Oxygen in the Ionic Liquids Methyl Trioctyl Ammonium Bis (trifluoromethylsulfonyl) imide, 1-Butyl-3-Methyl Imidazolium Bis (trifluoromethylsulfonyl) imide, and 1-Butyl-3-Methyl Imidazolium Methyl Sulphate

Keywords: carbon dioxide, oxygen, gas solubility, gravimetric analysis, ionic liquids

Ionic liquids (ILs) are being considered as solvents for gas absorption processes as they have the potential, in general, for improved efficiency of gas separations, as well as lower capital and operating costs compared to current commercial processes. In this study the solvent properties of ILs are investigated for use in the absorption of carbon dioxide (CO 2 ) and oxygen (O 2 ). The absorption of these gases in ILs was measured in a temperature range of 303.15 to 333.15 K and at pressures up to 1.5 MPa by gravimetric analysis. The ILs used were: methyl trioctyl ammonium bis (trifluoromethylsulfonyl) amide [MOA][Tf 2 N], 1-butyl-3-methyl imidazolium bis (trifluoromethylsulfonyl) amide [BMIM][Tf 2 N] and 1-butyl-3methyl imidazolium methyl sulphate [BMIM][MeSO 4 ]. The measurement technique employed in this study is fast, accurate, and requires small quantities of solvent. The results indicated that absorption of both gases increased with a decrease in operating temperature and an increase in pressure. [MOA][Tf 2 N] had the highest CO 2 and O 2 solubility.

[BMIM][Tf 2 N] was determined to have the highest selectivity for CO 2 absorption.

[BMIM][MeSO 4 ] achieved the lowest CO 2 absorption with a moderate O 2 absorption, revealing this IL to be the least desirable for CO 2 and O 2 absorption. Calculation of Henry's law constants for both gases in all ILs at all isotherms confirmed the deductions made from absorption data analysis. The absorption data was modelled using the generic Redlich-Kwong cubic equation of state (RK-EOS) coupled with a group contribution method.

INTRODUCTION

Recent research has shown that the main contributor to an increased rate and extent of global climate change, are CO 2 emissions from various fossil fuel industries including fossil fuelfired power plants, coal-to-liquids (CTL) and gas-to-liquids (GTL) industries, steel and cement industries. [START_REF]Climate Change 2007: The Physical Science Basis. Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change[END_REF] Thus, numerous studies have been undertaken to investigate different techniques for capturing CO 2 from flue gas emitted by these industries, and thereafter storing the CO 2 underground or by other means, thereby preventing its emission into the atmosphere. Despite significant research however, a technique with low energy penalty to feasibly capture CO 2 still remains to be found.

A common CO 2 capture technique closest to commercialisation is the use of solvents to selectively absorb CO 2 from flue gas emitted by fossil fuel industries. Current and conventional solvents include alkanolamines and carbonate-based solvents. These solvents however possess numerous disadvantages, including high corrosiveness, low absorption capacity and high energy requirements for desorption and recycling. Moreover, these solvents have proven to be volatile under industrial flue gas conditions, which has serious repercussions with regard to human health and the environment. [START_REF] Bretherick | Hazards in the Chemical Laboratory[END_REF]2 In some cases, the adverse effects of these solvents on human health, safety, and the environment, combined with their volatility and flammability, has led to legislation to minimize their use. 4 A relatively new type of solvent known as ionic liquids, has received much attention in recent years. Ionic liquids are molten salts that typically have melting points below 373.15 K. [START_REF] Deenadayalu | Ternary excess molar volumes of {methyltrioctylammonium bis(trifluoromethylsulfonyl)imide + ethanol + methyl acetate, or ethyl acetate} systems at T =[END_REF] The specific advantages of ionic liquids is that many can be tuned to have negligible vapour pressure, [START_REF] Bahadur | Apparent Molar Volume and Isentropic Compressibility for the Binary Systems {Methyltrioctylammonium Bis(trifluoromethylsulfonyl)imide + Methyl Acetate or Methanol} and (Methanol + Methyl Acetate) at T = 298[END_REF] good dissolution characteristics, 10 and are non-flammable 11 and non-toxic 12 , depending on the cation and anion used. It is for these reasons, that many ionic liquids have been regarded as -green‖ solvents, which can possibly replace more hazardous conventional solvents in many chemical processes, including CO 2 capture [START_REF] Tang | Ether-and alcohol-functionalized task-specific ionic liquids: attractive properties and applications[END_REF][START_REF] Davis | Task-specific ionic liquids[END_REF][START_REF] Rodríguez | Temperature and composition dependence of the density and viscosity of binary mixtures of water + ionic liquid[END_REF] . There are literally thousands of ionic liquids which can be synthesized, tested, and optimized for thermal, mechanical and chemical stability, and suitable density for various applications. 14 Gas absorption measurements in ionic liquids can provide information about the interaction such as dipole-dipole interactions that exist between these gases and ionic liquids. [START_REF] Anthony | Anion Effects on Gas Solubility in Ionic Liquids[END_REF] These studies can also provide important information about the original solvent behaviour of ionic liquids. [START_REF] Anthony | Anion Effects on Gas Solubility in Ionic Liquids[END_REF] Absorption data for many gas-ionic liquid systems are scarce compared to systems utilising conventional physical and chemical solvents. Research has however been conducted at an accelerated pace as the potential of ionic liquids becomes known. [START_REF] Shiflett | Solubility of CO 2 in Room Temperature Ionic Liquid[END_REF][START_REF] Almantariotis | Absorption of Carbon Dioxide, Nitrous Oxide, Ethane and Nitrogen by 1-Alkyl-3-methylimidazolium (C n mim, n = 2, 4, 6) Tris (pentafluoroethyl) trifluorophosphate Ionic Liquids (eFAP)[END_REF][START_REF] Gonzalez-Miquel | Anion Effects on Kinetics and Thermodynamics of CO 2 Absorption in Ionic Liquids[END_REF][START_REF] Kumełan | Solubility of Tetrafluoromethane in the Ionic Liquid[END_REF][START_REF] Zhou | Solubilities of CO 2 , H 2 , N 2 and O 2 in ionic liquid 1-n-butyl-3-methylimidazolium Heptafluorobutyrate[END_REF][START_REF] Afzal | Solubilities of some gases in four immidazoliumbased ionic liquids[END_REF][START_REF] Oh | High-pressure phase behavior of carbon dioxide in ionic liquid 1-butyl-3-methylimidazolium bis(trifluoromethylsulfonyl)imide[END_REF][START_REF] Sa´nchez | Solvent properties of functionalized ionic liquids for CO 2 absorption[END_REF][START_REF] Raeissi | Hydrogen Solubilities in the IUPAC Ionic Liquid 1-Hexyl-3-methylimidazolium Bis(Trifluoromethylsulfonyl)Imide[END_REF][START_REF] Ren | What Are Functional Ionic Liquids for the Absorption of Acidic Gases?[END_REF] Numerous sources in literature have investigated ionic liquids for their absorption capacity of CO 2 . CO 2 mole fractions were measured by gravimetric analysis in a variety of fluorinated and non-fluorinated ionic liquids at temperatures of 283.15-333.15 K and pressures up to 90 bar by Anderson et al. [START_REF] Zhou | Solubilities of CO 2 , H 2 , N 2 and O 2 in ionic liquid 1-n-butyl-3-methylimidazolium heptafluorobutyrate[END_REF] and Muldoon et al. [START_REF] Muldoon | Improving Carbon Dioxide Solubility in Ionic Liquids[END_REF] . Significant focus was made on CO 2 mole fraction measurement in non-fluorinated ionic liquids by Palgunadi et al. [START_REF] Palgunadi | Fluorine-free imidazolium-based ionic liquids with a phosphorous-containing anion as potential CO 2 absorbents[END_REF] and Zhang et al. [START_REF] Zhang | Dual aminofunctionalised phosphonium ionic liquids for CO 2 capture‖[END_REF] , at atmospheric pressure and temperatures of 298.15-334.15 K. Shiflett and Yokozeki [START_REF] Shiflett | Solubilities and diffusivities of carbon dioxide in ionic liquids[END_REF] , Arshad [START_REF] Arshad | CO 2 capture using Ionic Liquids[END_REF] , and Hasib-ur-Rahman et al. [START_REF] Hasib-Ur-Rahman | Ionic liquids for CO 2 capture-Development and progress[END_REF] Absorption measurements were performed by gravimetric analysis using an Intelligent Gravimetric Analyser (IGA-001) to determine the absorption isotherm of CO 2 and O 2 in the above three ionic liquids .at temperatures of 303.15, 313.15, 323.15 and 333.15 K and pressures up to 1.5 MPa. The Henry's law constant for both gases in all ionic liquids at all isotherms were calculated. The absorption data was modelled using a generic RK-EOS with the aid of a group contribution method. The 3D structure of the ionic liquids is given in Figure 1. The effect of cation, anion and the operating conditions on absorption of CO 2 and O 2 in fluorinated and non-fluorinated ionic liquids were analysed, together with the CO 2 selectivity of each ionic liquid over O 2 .

EXPERIMENTAL PROCEDURE

Materials

Carbon dioxide (CO 2 ) and oxygen (O 2 ) were obtained from Afrox Ltd (South Africa) with a minimum purity of 99.9 %. The ionic liquids: [MOA][Tf 2 N] was obtained from Fluka,

[BMIM][Tf 2 N] and [BMIM][MeSO 4 ] were obtained from Sigma-Aldrich, with a minimum purity of 98.0 %. Characterization of the chemicals is provided in Table 1.

Experimental Procedure

Density and Refractive Index Measurements

The measurements of density (ρ) and refractive index (n) of the pure ionic liquids were undertaken using a digital vibrating-tube densimeter (Anton Paar DMA 5000) and refractometer (ATAGO, model RX-7000a, Japan), respectively with an accuracy of ± 0.03 K in temperature. The estimated accuracy in density and refractive index was ± 0.00003 g•cm -3 and ± 0.00005, respectively. The experimental density (ρ) and refractive index (n) of pure ionic liquids together with literature [START_REF] Deenadayalu | Ternary excess molar volumes of {methyltrioctylammonium bis(trifluoromethylsulfonyl)imide + ethanol + methyl acetate, or ethyl acetate} systems at T =[END_REF][START_REF] Jacquemin | Density and viscosity of several pure and water-saturated ionic liquids[END_REF][START_REF] Pereiro | Physical Properties of 1-Butyl-3methylimidazolium Methyl Sulfate as a Function of Temperature[END_REF] densities (ρ) and refractive index (n) are given in Table 1. Moreover, the water content was determined using a Metrohm 702 SM Titrino Meter and was found to be 400 ppm in MeSO 4 ]. Samples were degassed before measurement as explained below.

[MOA][Tf 2 N], 600 ppm in [BMIM][Tf 2 N] and 500ppm in [BMIM][

Gas Absorption and Desorption Measurements

Gas absorption and desorption measurements of various gases in ionic liquids were conducted by gravimetric analysis using an Intelligent Gravimetric Analyser (IGA-001) built by Hiden Analytical Ltd. A quantity of approximately 70 to 100 mg of solvent was typically loaded onto a stainless steel sample holder which was then loaded into a stainless steel reactor. The sample holder was attached to tungsten and gold wire leading up to a microbalance which has a resolution of ±0.001 mg. The weight of the sample holder is countered using a counterweight so that only the sample weight may be tracked.

Each sample was degassed in the apparatus to ensure the removal of all volatile impurities using a Vacuubrand GMH-MD1 vacuum pump and an Edwards WRGS-NW35 turbomolecular pump to achieve ultra-high vacuum of 1x10 -4 Pa. The apparatus can operate in the pressure range of (0.05 to 1.50) MPa using two pre-calibrated pressure controllers and in the temperature range from (293.15 to 773.15) K using a Polyscience SD07R-20 refrigerated recirculating water bath and a Severn TF50/3/12/F furnace. A Pt100 probe with a resolution of ±0.01 K is located inside the reactor next to where the sample holder is suspended. Pressure was measured using a pressure transmitter with a resolution of ±1x10 -6 MPa. Once loaded, the sample was degassed. Thereafter, the gas was introduced and controlled at the desired pressure and temperature. The low quantity of sample required and the programmable nature of the apparatus through the use of IGA software provided by Hiden Isochema Ltd., made absorption measurements efficient and less time consuming.

Each data point took 1 to 4 hours to ensure equilibrium. MPa were considered in the case of CO 2 , whilst in the case of O 2 it is 0.050 -0.700 MPa due to limitations of gas regulation.

Due to the sample being of liquid form, there was a significant change in the sample density due to the absorption of the injected gas, thus making the sample more buoyant and altering the actual weight reading of the sample. The weight reading was hence not a true reflection of the actual mass of gas absorbed. To correct this, the buoyancy force was taken into account by first noting the weight reading of the sample at different gas pressures using a nonabsorbent gas such as nitrogen or helium. The method of buoyancy correction is supplied as supporting information. Buoyancy correction is also explained in Macedonia et al. [START_REF] Macedonia | Adsorption Studies of Methane, Ethane, and Argon in the Zeolite Mordenite: Molecular Simulations and Experiments[END_REF] . Absorption measurements were conducted and analysed with the data discussed in this study. Desorption measurements were also conducted however, to confirm the potential recyclability of the ionic liquids as solvents.

RESULTS AND DISCUSSION

The absorption and desorption of CO different temperatures as a function of pressure. Desorption measurements showed no consistent pattern of being higher than absorption measurements. Hysterisis is minimal, indicating that CO 2 can be fully desorbed and recovered using these ionic liquids. From Figures 345678, it can be seen that the absorption of both gases at all partial pressures was the highest at a lower temperature T = 303.15 K, and lowest at a higher temperature T = 333.15 K for each ionic liquid. This behaviour has also been observed in literature for other solvents including ionic liquids, confirming that absorption is an exothermic process. [START_REF] Shiflett | Solubilities and diffusivities of carbon dioxide in ionic liquids[END_REF] Another observation of the absorption results, which is also consistent with literature [START_REF] Tang | Ether-and alcohol-functionalized task-specific ionic liquids: attractive properties and applications[END_REF][START_REF] Rodríguez | Temperature and composition dependence of the density and viscosity of binary mixtures of water + ionic liquid[END_REF][START_REF] Anthony | Anion Effects on Gas Solubility in Ionic Liquids[END_REF] is that an increase with gas partial pressure results in higher gas absorption for each ionic liquid.

Additionally, the effect of cation and anion on absorption of gases in ionic liquids was studied and summarized in Figures 9 and10 respectively. By using the density data, we calculated the molar volumes of these ionic liquids at T = 303.15 K and it was found that the molar volume decreases in the following order:

[MOA][Tf 2 N] (588.64 cm 3 •mol -1 ) > [BMIM][Tf 2 N] (292.87 cm 3 •mol -1 ) > [BMIM][MeSO 4 ] (207.09 cm 3 •mol -1 )
. It is known that the absorption of gases in ionic liquids is mainly affected by two factors: (i) interaction between the gas and the anion of the ionic liquids and (ii) the molar volume of ionic liquids [START_REF] Anthony | Anion Effects on Gas Solubility in Ionic Liquids[END_REF][START_REF] Zhou | Solubilities of CO 2 , H 2 , N 2 and O 2 in ionic liquid 1-n-butyl-3-methylimidazolium heptafluorobutyrate[END_REF] .

Regarding CO 2 absorption, Figures 345 confirming with good certainty that increased fluorination of the cation achieves increased CO 2 absorption and particular cation types, higher cation chain lengths increase CO 2 absorption. This confirms findings from other literature as well [START_REF] Shiflett | Solubility of CO 2 in Room Temperature Ionic Liquid[END_REF][START_REF] Anderson | Solubility of CO 2 , CH 4 , C 2 H 6 , C 2 H 4 , O 2 , and N 2 in 1-hexyl-3-methylpyridinium bis(trifluoromethylsulfonyl)imide: Comparison to other ionic liquids[END_REF] that a larger chain length of the cation in an IL results in higher absorption being achieved.

The effect of the anion on CO 2 absorption can be confirmed by also comparing the measured results between [BMIM][Tf 2 N] and [BMIM][MeSO 4 ] seen in Figures 4 and5, where As far as molar volume was concerned, it is proposed that ionic liquids with higher molar volume would achieve greater absorption of a gas solute. [START_REF] Anthony | Anion Effects on Gas Solubility in Ionic Liquids[END_REF][START_REF] Zhou | Solubilities of CO 2 , H 2 , N 2 and O 2 in ionic liquid 1-n-butyl-3-methylimidazolium heptafluorobutyrate[END_REF][START_REF] Seki | In Situ Attenuated Total Reflection Infrared Spectroscopy of Imidazolium-Based Room-Temperature Ionic Liquids under -Supercritical‖ CO 2[END_REF][START_REF] Kazarian | Combining ionic liquids and supercritical fluids: in situ ATR-IR study of CO 2 dissolved in two ionic liquids at high pressures[END_REF][START_REF] Camper | Gas solubilities in room-temperature ionic liquids[END_REF] This may be because increasing the molar volume of ionic liquids increases the void space (free volume), [START_REF] Finotello | Ideal Gas Solubilities and Solubility Selectivities in a Binary Mixture of Room-Temperature Ionic Liquids[END_REF][START_REF] Shannon | Fractional Free Volume as the Basis of Gas Solubility & Selectivity in Imidazolium-based Ionic Liquids[END_REF] and gases appear to be dissolved by occupying these free spaces available in ionic liquids. [START_REF] Finotello | Ideal Gas Solubilities and Solubility Selectivities in a Binary Mixture of Room-Temperature Ionic Liquids[END_REF][START_REF] Scovazzo | Regular solution theory and CO 2 gas solubility in room temperature ionic liquids[END_REF][START_REF] Huang | Why is the Partial Molar Volume of CO 2 So Small When Dissolved in a Room Temperature Ionic Liquid? Structure and Dynamics of CO 2 Dissolved in [Bmim] + [PF6][END_REF] It is clear from the result that the absorption of CO 2 in the tested ionic liquids decrease in order: [START_REF] Zhou | Solubilities of CO 2 , H 2 , N 2 and O 2 in ionic liquid 1-n-butyl-3-methylimidazolium heptafluorobutyrate[END_REF] Regarding O 2 absorption, Figure 7 shows a favourable comparison of measurements done in this work, with measurements conducted by Anthony et al (2005) [START_REF] Anthony | Anion Effects on Gas Solubility in Ionic Liquids[END_REF] at 323.15K. Data was especially comparable at low pressure while some deviation occurred at higher pressure of 0.7 MPa. This suggests that the technique of gravimetric analysis using N 2 for buoyancy correction may be accurate for systems O 2 especially at pressure below 0.4 MPa. indicating that this ionic liquid may not necessarily be favourable for gas separation processes appreciable amounts of both CO 2 and O 2 , due to its low selectivity. By contrast,

[MOA][Tf 2 N] > [BMIM][Tf 2 N] > [BMIM][MeSO
[BMIM][Tf 2 N] achieved the 2 nd highest CO 2 absorption capacity and the lowest O 2 absorption capacity, indicating that this ionic liquid is more selective towards CO 2 absorption and is hence more applicable in CO 2 capture processes.

An effective way of comparing ionic liquids and solvents for gas absorption is to calculate the Henry's law constant, K H , of the absorbed gas in the binary system at each temperature.

where x gas is the mole fraction of CO 2 and O 2 in IL and P gas is the equilibrium pressure.

CORRELATION OF ABSORPTION DATA

The modelling of gas absorption in ionic liquids has been reported in numerous literature sources. More complex models take into account the asymmetry or cations and anions, as well as dipole and quadrupolar interactions.

Due to the relative simplicity of the ionic liquid species in this study in terms of asymmetry, absorption data of CO 2 and O 2 in all ionic liquids was successfully correlated using a much simpler equation of state (EOS). The generic RK-EOS was used, which was previously applied successfully for ionic liquids containing fluorinated anions and imidazolium cations [START_REF] Shiflett | Solubility of CO 2 in Room Temperature Ionic Liquid[END_REF] :

) ( ) ( b V V T a b V RT P     (6)
where a(T) in eq 6, is a function of temperature and species mole fraction. For the components in this work, the Van der Waals-Berthelot mixing formula was used as successfully applied for similar systems. [START_REF] Shiflett | Solubility of CO 2 in Room Temperature Ionic Liquid[END_REF][START_REF] Yokozeki | Solubility of refrigerants in various lubricants[END_REF]

  j i j 1 , j j i 1 1 ) ( x x k T a a T a i N j i i               where   j ij i ji j i ji j ij x l x l x x l l k i    (7)   T P T R a i ci 2 ci 2 i 427480 . 0   where                3 0 k C C k i k T T T T T   (8)
and l ij , l ji , k ij , and m ij are binary interaction parameters with m ij = m ji and m ii = 0. Coefficients β k are treated as adjustable fitting parameters. As with many other ionic liquids, critical temperature and pressure could not be found in literature for the ionic liquids in this study. This is due to ionic liquids decomposing before reaching their critical point.

          N j i i x x k m b b b 1 , j i ij j j i
Critical temperature and pressure were thus computed using a modified Lydersen-Joback-Reid group contribution method: [START_REF] Alvarez | A modified Lydersen-Joback-Reid method to estimate the critical properties of biomolecules[END_REF][START_REF] Valderrama | Critical properties, normal boiling temperatures, and acentric factors of fifty ionic liquids[END_REF]

         2 M M M M b C T n T n B A T T where     bM b 2 . 198 T n T (10)   2 M M C     P n C M P (11)
where A M = 0.5703, B M = 1.0121, C M = 0.2573, E M = 6.75. n is the number of occurrences of any particular functional group in the molecule. Group contributions of ∆T bM are provided in the work of Alvarez and Valderrama. 58 In this study, we modelled the absorption of CO 2 and O 2 in all three ionic liquids. Model parameters are provided as supporting information in Table S1 while estimated critical properties are shown in Table S2.

Figures 345678, illustrate the model-calculated CO 2 and O 2 partial pressure using parameters obtained by data regression, in comparison to the measured data. It is evident that model estimates compare quite favourably with measured data. The root-mean-square error for each system is provided in Table 9, showing very good accuracy of model estimates, which is consistent with the work of Shiflett and Yokozeki [START_REF] Shiflett | Solubilities and diffusivities of carbon dioxide in ionic liquids[END_REF] for CO 2 absorption in [BMIM][PF 6 ] and

[BMIM][BF 4 ].
It is observed from Figures 678that O 2 absorption decreases significantly with increasing temperature. The temperature effect is far more pronounced than with CO 2 absorption. As a result modelling of the data became more difficult, and while Table 9 contains regressed parameters for CO 2 absorption in all ionic liquids, common for all temperatures, parameters for O 2 absorption differ for each isotherm.

CONCLUSIONS

Absorption measurements of CO 2 and O 2 in the ionic liquids [MOA][Tf 2 N], [BMIM][Tf 2 N] and [BMIM][MeSO 4

] have been reported. The data of each isotherm was confirmed to incur low hysterisis. The ionic liquid [MOA][Tf 2 N] achieved the highest CO 2 and O 2 absorption from all the ionic liquids studied and some other ionic liquids compared in the literature, which indicates good CO 2 solubility but undesirably low CO 2 selectivity The ionic liquid

[BMIM][Tf 2 N], although having not achieved the highest CO 2 absorption of the three ionic liquids studied, was found the most suitable for CO 2 absorption from flue gas due to its higher selectivity of CO 2 over O 2 .. The calculation of Henry's law constants for CO 2 and O 2 in all ionic liquids at all isotherms confirmed the findings made by analysis of the absorption data, and provide an effective numerical comparison of ionic liquids performance as gas absorbents. A further analysis of enthalpy and entropy of absorption revealed solubility in [BMIM][MeSO 4 ] to be highly temperature sensitive. The generic RK-EOS model achieved accurate estimates for all systems studied, with a root mean square error of no higher than 1.98 kPa. The model was found to be more robust for systems containing CO 2 , while parameter regression became more temperature sensitive with systems containing O 2 .

(10) Bahadur 303.15 1.5x10 -1 -2.3x10 -2 0.0 0.0 1.5x10 -1 1.0x10 14 -2.3x10 -2 9.1x10 -1 7.1x10 -6 313.15 1.5x10 -1 -2.3x10 -2 0.0 0.0 1.5x10 -1 1.0x10 14 -2.3x10 -2 9.1x10 -1 7.1x10 -6 323.15 1.5x10 -1 -2.3x10 -2 0.0 0.0 1.5x10 -1 1.0x10 14 -2.3x10 -2 9.1x10 -1 7.1x10 -6 333.15 1.5x10 -1 -2.3x10 -2 0.0 0.0 1.5x10 -1 1.0x10 14 -2.3x10 -2 9.1x10 -1 7.1x10 -6 CO 2 (1) + [BMIM][Tf 2 N] (2) 303.15 6.4x10 -2 8.4x10 -4 0.0 0.0 -2.2x10 -1 1.2x10 1 -8.8x10 2 0.3x10 1 1.4x10 -4 313.15 6.4x10 -2 8.4x10 -4 0.0 0.0 -2.2x10 -1 1.2x10 303.15 1.7x10 -1 -1.7x10 -2 0.0 0.0 9.9x10 -1 9.9x10 -1 7.25x10 5 7.18x10 [START_REF]Climate Change 2007: The Physical Science Basis. Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change[END_REF] 1.5x10 -4 313.15 1.7x10 -1 -1.7x10 -2 0.0 0.0 9.9x10 -1 9.9x10 -1 7.25x10 5 7.18x10 3 1.5x10 -4 323.15 1.7x10 -1 -1.7x10 -2 0.0 0.0 9.9x10 -1 9.9x10 -1 7.25x10 5 7.18x10 3 1.5x10 -4 333.15 1.7x10 -1 -1.7x10 -2 0.0 0.0 9.9x10 -1 9.9x10 -1 7.25x10 5 7.18x10 [START_REF]Climate Change 2007: The Physical Science Basis. Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change[END_REF]

  considered fluorinated imidazolium-based ionic liquids in their studies, at temperatures of 298.15-348.15 K and 0.1-800 bar. Anderson et al. 50 utilised a Rubotherm for high pressure VLE measurement. The aim of this study was to measure CO 2 and O 2 absorption in Methyl trioctyl ammonium bis (trifluoromethylsulfonyl) imide ([MOA][Tf 2 N]), 1-butyl-3-methyl imidazolium bis (trifluoromethylsulfonyl) imide ([BMIM][Tf 2 N]) and 1-butyl-3-methyl imidazolium methyl sulphate ([BMIM][MeSO 4 ]). While CO 2 absorption in [BMIM][Tf 2 N] is well measured 34 , O 2 absorption has not. The measurement of O 2 absorption is important since flue gas from coal power plants may contain up to 5% unreacted O 2 (by volume), and solvents need to be highly CO 2 selective to enable the efficient recovering of pure CO 2 for disposal. Moreover, CO 2 and O 2 absorption in the relatively new ionic liquids had not been measured previously.

  Figure S-1 in supporting information illustrates the sample weight change with time for CO 2 in [MOA][Tf 2 N] at 303.15 K and at 0.05MPa, with that particular point taking approximately 1 hour for equilbirum due to small sample size and efficient temperature and pressure control. The absorption and desorption isotherms at 303.15, 313.15, 323.15 and 333.15 K were measured. Pressures of 0.050 -1.500

  2 and O 2 in [MOA][Tf 2 N], [BMIM][Tf 2 N] and [BMIM][MeSO 4 ] were determined by gravimetric analysis using an Intelligent Gravimetric Analyser (IGA). Systems with a fixed overall equilibrium composition of CO 2 and O 2 in [MOA][Tf 2 N], [BMIM][Tf 2 N] and [BMIM][MeSO 4 ] were measured as a function of temperature and pressure are given in Tables 2-7. Figures 3-8 graphically show CO 2 and O 2 absorption and desorption in [MOA][Tf 2 N], [BMIM][Tf 2 N] and [BMIM][MeSO 4 ] at

  , respectively. Regarding the three ionic liquids studied in this work, two contain the same fluorinated anions [Tf 2 N] but different class of cation such as [MOA], [BMIM] and two ionic liquids have same cation [BMIM] but different anions of [Tf 2 N] and [MeSO 4 ]. The molar masses are 648.85, 419.36 and 250.32 g•mol -1 for [MOA][Tf 2 N], [BMIM][Tf 2 N] and [BMIM][MeSO 4 ],

  show that [MOA][Tf 2 N] achieved the highest CO 2 absorption at each temperature and partial pressures studied, followed by [BMIM][Tf 2 N], with [BMIM][MeSO 4 ] achieving the lowest CO 2 absorption. This is because of the effect of cation, anion and molar volume of ionic liquids.

Figure 9

 9 Figure 9 provides a further comparison of CO 2 absorption in the ionic liquids measured in

Figure 10

 10 Figure 10 illustrates that [BMIM][Tf 2 N] achieves higher CO 2 absorption than

[

  BMIM] [Tf 2 N] achieved higher CO 2 absorption than [BMIM][ MeSO 4 ]. This is because [BMIM][Tf 2 N] has a higher molar volume and greater interaction between CO 2 and [Tf 2 N] than [MeSO 4 ] anion, due to increased fluorination.

  4 ] being possible to achieve maximum mole fractions of absorbed CO 2 corresponding to 0.343, 0.283, and 0.138, respectively at T = 303.15 K and 1.5 MPa. The superior CO 2 absorption achieved by [MOA][Tf 2 N] and [BMIM][Tf 2 N] confirms that ionic liquids with fluorinated anions achieve higher CO 2 absorption than ionic liquids with non-fluorinated anions. The study also found that [MOA][Tf 2 N] achieved higher CO 2 absorption at T = 323.15 K than more conventional ionic liquids such as [BMIM][BF 4 ] and [BMIM][PF 6 ] as measured by (Shiflett and Yokozeki) 49 and [HMIM][Tf 2 N] by Anderson et al.,

Figures 6 - 8 ,

 68 Figures 6-8, show that [MOA][Tf 2 N] achieved the highest O 2 absorption at each temperature

Figure 1 .

 1 Figure 1. 3D structure of ionic liquids: (a) [MOA][Tf 2 N], (b) [BMIM][Tf 2 N] and (c)

Figure 2 .

 2 Figure 2. Experimental Set-up of IGA-001 used for the measurements of Gas

Figure 3 .

 3 Figure 3. Isothermal absorption and desorption of CO 2 in [MOA][Tf 2 N]: (◊) absorption, (□)

Figure 4 .

 4 Figure 4. Isothermal absorption and desorption of CO 2 in [BMIM][Tf 2 N]: (◊) absorption, (□)

Figure 5 .

 5 Figure 5. Isothermal absorption and desorption of CO 2 in [BMIM][MeSO 4 ]: (◊) absorption,

Figure 6 .

 6 Figure 6. Isothermal absorption and desorption of O 2 in [MOA][Tf 2 N]: (◊) absorption, (□)

Figure 7 .

 7 Figure 7. Isothermal absorption and desorption of O 2 in [BMIM][Tf 2 N]: (◊) absorption, (□)

Figure 8 .

 8 Figure 8. Isothermal absorption and desorption of O 2 in [BMIM][MeSO 4 ]: (◊) absorption,

Figure 9 . 47 Figure 10 .Figure 11 .Figure 1 .

 94710111 Figure 9. Figure: Isothermal absorption of CO 2 in [MOA][Tf 2 N] in comparison to other ionic

Figure 3 .

 3 Figure 3.

Figure 4 .

 4 Figure 4.

Figure 5 .

 5 Figure 5.

Figure 6 .

 6 Figure 6.

Figure 7 .

 7 Figure 7.

Figure 8 .

 8 Figure 8.

Figure 9 .

 9 Figure 9.

  Figure 10.

Figure S- 1 :

 1 Figure S-1: Kinetic data of sample weight for CO 2 solubility in [MOA][Tf 2 N] at 303.15 K and at 0.05 MPa

Table 8 contains

 8 Henry's law constant values measured in MPa for CO 2 and O 2 in all three ionic liquids and at each isotherm. A high deviation of 0.4 MPa is noted, due to higher Henry's law constants calculated at pressures of 1.3 and 1.5 MPa. . Nevertheless, a comparison of Henry's law constants confirms the conclusions drawn from the absorption measurements. The high selectivity of [BMIM][Tf 2 N] for CO 2 absorption is ascertained by comparing the Henry's constants of CO 2 and O 2 for all isotherms.

Table 9 contains

 9 Enthalpy and Entropy of absorption for CO 2 and O 2 in all ionic liquids. Results for [BMIM][Tf 2 N] compare favourably with Arshad (2009) 38 . Enthalpy and Entropy of absorption values show that CO 2 absorption in [MOA][Tf 2 N] is the least sensitive to temperature, while absorption in [BMIM][MeSO 4 ] is the most sensitive. The opposite is true in the case of O 2 absorption, although comparatively the high enthalpies and entropies of absorption for O 2 indicate low O 2 solubility in relation to CO 2 in these ionic liquids.

Table 1 : Densities, ρ, and Refractive Index, n, of Pure ILs together with Literatures at

 1 , I.; Deenadayalu, N. Apparent Molar Volume and Apparent Molar Isentropic From ref 8. b From ref 33. c From ref 36.

	Compressibility	for	the	Binary	Systems	{Methyltrioctylammonium
	bis(trifluoromethylsulfonyl)imide + Ethyl Acetate or Ethanol} at Different Temperatures
	under Atmospheric Pressure. Thermochim. Acta 2013, 566, 77-83.
	(11) Deenadayalu, N.; Bahadur, I.; Hofman, T. Ternary Excess Molar Volumes of
	{Methyltrioctylammonium Bis[(trifluoromethyl)sulfonyl]imide + Methanol + Methyl Acetate
	or Ethyl Acetate} Systems at (298.15, 303.15, and 313.15) K. J. Chem. Eng. Data 2010, 55,
	2636-2642.

(12) Wood, N.; Stephens, G. Accelerating the discovery of biocompatible ionic liquids.

Phys

. Chem. Chem. Phys. 2010, 12, 1670-1674. (13) Bahadur, I.; Deenadayalu, N.; Tywabi, Z.; Sen, S.; Hofman, T. Volumetric properties of ternary (IL + 2-propanol or 1-butanol or 2-butanol + ethyl acetate) systems and binary (IL + 2-propanol or 1-butanol or 2-butanol) and (1-butanol or 2-butanol + ethyl acetate) systems. J. Chem. Thermodyn. 2012, 49, 24-38. (14) Dreisbach, F.; Meister, D.; Lösch, H. W.; Weißert, T.; Petermann, M. Ionic liquids for gas separation: Measurement of selective gas solubilities of the gas-mixture CH 4 /CO 2 . http://www.merrowscientific.com/Portals/55675/docs/Ionic_Liquids.pdf 2013. (15) Wasserscheid, P.; Welton, T. Ionic Liquids in Synthesis. Wiley-VCH: Weinheim, Neural Networks. Ind. Eng. Chem. Res. 2008, 47, 4523-4532.

Table 2 : Experimental Results for Absorption and Desorption of CO 2 in [MOA][Tf 2 N]

 2

	P meas /MPa	T/K	absorption	x CO2	P meas /MPa	T/K	x CO2 desorption	P calc /MPa
	0.0499	303.18	0.017	0.0499	303.17	0.017	0.0500
	0.1000	303.21	0.031	0.0998	303.15	0.035	0.1000
	0.3998	303.18	0.122	0.3998	303.21	0.138	0.4003
	0.7000	303.18	0.198	0.7000	303.15	0.197	0.7005
	1.0000	303.18	0.263	0.9998	303.05	0.280	1.0000
	1.2999	303.17	0.319	1.2997	303.17	0.323	1.3000
	1.4999	303.23	0.343	1.4999	303.23	0.344	1.5011
	0.0499	313.25	0.019	0.0500	313.14	0.013	0.0500
	0.1000	313.09	0.027	0.0998	313.13	0.030	0.0999
	0.4000	313.13	0.104	0.3998	313.07	0.110	0.3999
	0.6999	313.15	0.173	0.7000	313.10	0.187	0.6999
	0.9997	313.13	0.235	0.9999	313.12	0.235	0.9993
	1.2999	313.16	0.287	1.2998	313.07	0.293	1.2989
	1.4996	313.17	0.313	1.4996	313.17	0.314	1.4992
	0.0500	323.26	0.0143	0.0499	323.16	0.010	0.0500
	0.1000	323.18	0.0215	0.0998	323.10	0.022	0.0999
	0.4001	323.14	0.0886	0.3998	323.22	0.094	0.3997
	0.7000	323.20	0.1511	0.6999	323.07	0.154	0.6999
	0.9999	323.19	0.1981	1.0000	323.19	0.212	0.9999
	1.3000	323.17	0.2493	1.3000	323.05	0.260	1.2993
	1.5000	323.15	0.2820	1.5000	323.15	0.282	1.4992
	0.0500	333.24	0.0112	0.0499	333.14	0.006	0.0500
	0.0999	333.14	0.0178	0.0999	333.24	0.017	0.0999
	0.3999	333.12	0.0773	0.3998	333.18	0.086	0.3996
	0.6998	333.11	0.1330	0.6999	333.04	0.141	0.6999
	1.0000	333.24	0.1802	0.9999	333.07	0.191	1.0006
	1.2998	333.16	0.2281	1.2998	333.12	0.234	1.3005
	1.4999	333.13	0.2595	1.4999	333.13	0.258	1.5006
	Uncertainty: T = ±0.01 K; P = 1x10 -6 MPa; x = ±0.00005

Table 3 : Experimental Results for Absorption and Desorption of CO 2 in BMIM][Tf 2 N]

 3

	P meas /MPa	T/K	x CO2 absorption	P meas /MPa	T/K desorption x CO2	P calc /MPa
	0.0499	303.22	0.016	0.0499	303.17	0.014	0.0504
	0.1000	303.16	0.027	0.0999	303.13	0.028	0.0996
	0.4000	303.18	0.096	0.3998	303.14	0.104	0.3996
	0.7000	303.18	0.158	0.6994	303.06	0.165	0.7000
	0.9998	303.25	0.202	0.9996	303.06	0.219	1.0009
	1.3001	303.20	0.260	1.2996	303.10	0.263	1.3013
	1.4999	303.22	0.283	1.4999	303.22	0.282	1.5021
	0.0499	313.28	0.015	0.0499	313.17	0.008	0.0503
	0.1000	313.19	0.020	0.0998	313.15	0.024	0.0984
	0.4001	313.13	0.077	0.3998	313.10	0.083	0.3993
	0.7001	313.18	0.131	0.6997	313.05	0.135	0.6996
	1.0001	313.14	0.178	0.9995	313.07	0.182	0.9996
	1.2998	313.14	0.216	1.2999	313.12	0.223	1.299
	1.5000	313.11	0.245	1.5000	313.11	0.245	1.4977
	0.0500	323.27	0.011	0.0499	323.15	0.005	0.0501
	0.1000	323.22	0.016	0.0998	323.15	0.013	0.1019
	0.3999	323.20	0.064	0.3997	323.09	0.051	0.3992
	0.7000	323.14	0.108	0.6998	323.12	0.098	0.6994
	1.0000	323.24	0.147	0.9997	323.17	0.140	1.0000
	1.3000	323.09	0.186	1.2999	323.10	0.179	1.2985
	1.5000	323.17	0.210	1.5000	323.17	0.210	1.4974
	0.0498	333.18	0.004	0.0499	333.07	0.005	0.0501
	0.0999	333.15	0.010	0.1000	333.06	0.010	0.1005
	0.4020	333.21	0.042	0.3999	333.22	0.052	0.3982
	0.7001	333.21	0.094	0.6999	333.14	0.091	0.6995
	0.9998	333.17	0.119	0.9998	333.17	0.124	1.0003
	1.3000	333.19	0.145	1.2999	333.12	0.141	1.3015
	1.5000	333.17	0.155	1.5001	333.14	0.149	1.5028
	Uncertainty: T = ±0.01 K; P = 1x10 -6 MPa; x = ±0.00005

Table 4 : Experimental Results for Absorption and Desorption of CO 2 in [BMIM][MeSO 4 ]

 4

	P meas /MPa	T/K absorption x CO2	P meas /MPa	T/K	x CO2 desorption	P calc /MPa
	0.0499	303.18	0.005	0.0505	303.18	0.008	0.0500
	0.1000	303.20	0.010	0.1001	303.19	0.015	0.1000
	0.3997	303.16	0.045	0.4011	303.16	0.048	0.4000
	0.7001	303.09	0.084	0.7003	303.09	0.085	0.6998
	1.0001	303.20	0.104	0.9996	303.10	0.106	1.0006
	1.3001	303.17	0.130	1.3001	303.09	0.132	1.3003
	1.5000	303.17	0.138	1.5000	303.17	0.140	1.5012
	0.0506	313.25	0.008	0.0509	313.15	0.007	0.0500
	0.0998	313.20	0.010	0.1005	313.14	0.013	0.1000
	0.4005	313.13	0.038	0.4000	313.06	0.040	0.3998
	0.6998	313.11	0.074	0.6999	313.17	0.071	0.6994
	1.0000	313.14	0.096	0.9995	313.05	0.097	0.9994
	1.3000	313.16	0.118	1.2999	313.11	0.117	1.2989
	1.5000	313.16	0.129	1.5000	313.16	0.129	1.4988
	0.0507	323.23	0.007	0.0502	323.14	0.003	0.0500
	0.1003	323.11	0.007	0.1004	323.11	0.007	0.1000
	0.4000	323.15	0.035	0.4004	323.19	0.035	0.3997
	0.7001	323.19	0.058	0.7001	323.08	0.058	0.6998
	1.0001	323.24	0.078	1.0001	323.15	0.077	1.0001
	1.2999	323.23	0.097	1.3000	323.21	0.099	1.3004
	1.4999	323.13	0.112	1.4999	323.13	0.109	1.4991
	0.0504	333.23	0.003	0.0506	333.04	-0.001	0.0500
	0.1011	333.10	0.004	0.1004	333.10	0.003	0.1000
	0.3998	333.22	0.028	0.4000	333.20	0.030	0.3998
	0.6999	333.25	0.047	0.7000	333.17	0.050	0.6999
	0.9998	333.09	0.065	0.9998	333.07	0.069	0.9998
	1.2998	333.13	0.083	1.3000	333.08	0.085	1.3004
	1.5000	333.05	0.092	1.5000	333.05	0.092	1.5009
	Uncertainty: T = ±0.01 K; P = 1x10 -6 MPa; x = ±0.00005

Table 5 : Experimental Results for Absorption and Desorption of O 2 in [MOA][Tf 2 N]

 5

	P meas /MPa	T/K	absorption	x O2	P meas /MPa	T/K desorption	x O2	P calc /MPa
	0.0499	303.11	0.004	0.0498	303.17	0.003	0.0499
	0.0999	303.19	0.008	0.0999	303.11	0.008	0.0999
	0.3997	303.15	0.016	0.3997	303.11	0.016	0.3997
	0.6998	303.15	0.023	0.6998	303.15	0.024	0.6998
	0.0498	313.31	0.002	0.0498	313.14	0.001	0.0511
	0.0999	313.11	0.004	0.0999	313.10	0.002	0.1012
	0.3997	313.21	0.013	0.3998	313.07	0.012	0.3996
	0.7000	313.17	0.020	0.7000	313.17	0.021	0.7000
	0.0498	323.25	0.000	0.0498	323.13	0.000	0.0499
	0.0999	323.18	0.001	0.0999	323.23	0.002	0.0999
	0.3999	323.24	0.008	0.3998	323.16	0.007	0.3999
	0.6999	323.15	0.013	0.6999	323.15	0.013	0.6999
	0.0498	333.28	0.000	0.0498	333.15	0.0002	0.0497
	0.0999	333.18	0.001	0.0999	333.11	0.0003	0.0994
	0.3996	333.23	0.003	0.3997	333.25	0.003	0.3996
	0.7000	333.05	0.004	0.7000	333.05	0.003	0.7000
	Uncertainty: T = ±0.01 K; P = 1x10 -6 MPa; x = ±0.00005

Table 6 : Experimental Results for Absorption and Desorption of O 2 in [BMIM][Tf 2 N]

 6

	P meas /MPa	T/K absorption	x O2	P meas /MPa	T/K desorption	x O2	P calc /MPa
	0.0498	303.19	0.001	0.050	303.162	0.001	0.0500
	0.0999	303.17	0.002	0.100	303.071	0.002	0.0999
	0.3999	303.22	0.006	0.400	303.109	0.005	0.3999
	0.7000	303.25	0.006	0.700	303.248	0.005	0.7000
	0.0498	313.29	0.001	0.050	313.154	0.001	0.0500
	0.0999	313.12	0.002	0.100	313.059	0.001	0.0999
	0.3999	313.12	0.004	0.400	313.054	0.004	0.3999
	0.7001	313.20	0.005	0.700	313.197	0.005	0.7001
	0.0499	323.27	0.000	0.050	323.143	0.001	0.0500
	0.0998	323.09	0.001	0.100	323.220	0.001	0.0999
	0.3998	323.24	0.003	0.400	323.138	0.003	0.3998
	0.6999	323.23	0.005	0.700	323.234	0.004	0.6999
	0.0500	333.29	0.0003	0.050	333.138	0.0003	0.0500
	0.0999	333.19	0.0003	0.100	333.114	0.0002	0.0999
	0.3999	333.15	0.0004	0.400	333.066	0.0003	0.3999
	0.7000	333.19	0.0005	0.700	333.186	0.0005	0.6999
	Uncertainty: T = ±0.01 K; P = 1x10 -6 MPa; x = ±0.00005

Table 7 : Experimental Results for Absorption and Desorption of O 2 in [BMIM][MeSO 4 ]

 7

	P meas /MPa	T/K absorption x O2	P meas /MPa	T/K	desorption	x O2	P calc /MPa
	0.0495	303.18	0.007	0.0495	303.19	0.007	0.0500
	0.0994	303.15	0.008	0.0996	303.06	0.008	0.0999
	0.3994	303.22	0.012	0.3994	303.11	0.012	0.3989
	0.6995	303.06	0.013	0.6995	303.06	0.013	0.6997
	0.0496	313.30	0.007	0.0495	313.14	0.007	0.0486
	0.0994	313.15	0.007	0.0996	313.06	0.007	0.0973
	0.3995	313.08	0.011	0.3995	313.04	0.011	0.3942
	0.6995	313.15	0.012	0.6995	313.15	0.012	0.7000
	0.0495	323.33	0.006	0.0496	323.14	0.005	0.0496
	0.0993	323.11	0.007	0.0996	323.16	0.007	0.0999
	0.3993	323.20	0.009	0.3996	323.17	0.009	0.3992
	0.6995	323.23	0.010	0.6995	323.23	0.010	0.6995
	0.0495	333.35	0.003	0.0497	333.13	0.003	0.0501
	0.0998	333.10	0.003	0.0995	333.13	0.003	0.0999
	0.3993	333.10	0.005	0.3996	333.24	0.005	0.3993
	0.6996	333.15	0.007	0.6996	333.15	0.007	0.6996
	Uncertainty: T = ±0.01 K; P = 1x10 -6 MPa; x = ±0.00005

Table 8 : Henry ś Law C onstants of CO 2 and O 2 (K HCO2 and K HO2 ) in [MOA][Tf 2 N], [BMIM][Tf 2 N] and [BMIM][MeSO 4 ] Estimated from Absorption Data

 8

	T/K	[MOA][Tf 2 N]	[BMIM][Tf 2 N] K HCO2 /MPa	[BMIM][MeSO 4 ]
	303.15	3.7	4.6	9.6
	313.15	4.2	5.6	10.5
	323.15	4.9	6.7	12.9
	333.15	5.5	9.0	15.3
			K HO2 /MPa
	303.15	19.6	65.6	26.6
	313.15	30.9	82.9	28.9
	323.15	59.1	128.6	36.1
	333.15	138.1	223.7	58.9

Table 9 : Enthalpy and Entropy of Absorption for CO 2 and O 2 in the Ionic Liquids Studied in this Work

 9

		[BMIM][Tf 2 N]	[MOA][Tf 2 N] ∆H/kJ•mol -1	[BMIM][MeSO 4 ]
	CO 2 -13.46 ± 1.8	-10.46 ± 1.8	-18.15 ± 4.46
	O 2	-26.80 ± 2.32	-25.53 ± 3.45 ∆S/J•mol -1 •K -1	-9.25 ± 3.29
	CO 2 -45.30 ± 5.74	-34.27 ± 4.96	-47.52 ± 5.47
	O 2	-79.12 ± 19.07 -81.80 ± 11.02 -35.62 ± 2.61

Table S1 : Binary Interactions and Fitting Parameters by using RK-EOS model

 S1

	T/ (K)	β 0	β 1	β 2	β 3	ℓ 12	ℓ 21	τ 12	m 12	σ
					CO 2 (1) + [MOA][Tf 2 N] (2)

Table S2 : Critical Properties of ILs used in this study ionic

 S2 1.5x10 -4 O 2 (1) + [MOA][Tf 2 N] (2) liquid molar mass/(g.mol -1 ) T b /K T C /K P C /bar V C /(cm 3 .mol -1 ) ω [MOA][Tf 2 N]

			648.85		1190.0 1447.35 10.31	2002.30	1.0096
	[BMIM][Tf 2 N]		419.37		784.6 1133.41 25.69	956.02	0.3526
	[BMIM][MeSO 4 ]	250.32		595.0	877.42 35.51	663.94	0.3913
	303.15 -1.7x10 -1 4.3x10 -2	0.0	0.0	-6.4	-8.7x10 1	-1.1x10 -5	2.1	0.0
	313.15 1.2	-2.3x10 -1 0.0	0.0	7.2x10 -8	-3.5x10 -6	6.8x10 1	9.7x10 -1 6.3x10 -4
	323.15 3.3	-4.9x10 -1 0.7	-0.2	0.0	0.0	-3.7x10 5	5.6x10 1	2.5x10 -5
	333.15 1.7	1.4x10 -1	1.2	-0.3	-8.0x10 -3	1.8	-3.0x10 2	9.9x10 -1 1.5x10 -4
				O 2 (1) + [BMIM][Tf 2 N] (2)
	303.15 5.5	-1.4	0.0	0.0	6.3x10 -10	-9.9x10 -8	-5.4x10 2	1.2	5.0x10 -5
	313.15 1.2	-3.0x10 -1 0.0	0.0	-6.2x10 -3	1.2	-1.0x10 2	9.9x10 -1 5.0x10 -5
	323.15 1.2x10 2	-4.9x10 1 1.0	1.0	-4.2x10 -7	-6.1	-2.5x10 4	3.9x10 1	5.0x10 -5
	333.15 1.1x10 2	-4.8x10 1 1.0	1.0	7.8x10 -4	2.8x10 3	8.6x10 6	8.0x10 2	2.5x10 -5
			O 2 (1) + [BMIM][MeSO 4 ] (2)
	303.15 9.9	-3.1	0.0	0.0	-1.2	-3.5	-1.1	1.0x10 -1 1.7x10 -4
	313.15 3.4	-1.1	0.0	0.0	-5.0x10 -2	4.1	-2.0x10 -1	1.0	2.0x10 -3
	323.15 1.0	-1.2x10 1 1.0	1.0	-5.2 x10 -2 -3.1x10 5	-3.2x10 2	1.0	1.5x10 -4
	333.15 2.9	-1.3 x10 1 1.0	1.0	-3.9x10 -9	5.7x10 -7	-6.5x10 1	-2.5	1.8x10 -4

ACKNOWLEDGMENTS

The authors acknowledge University of KwaZulu-Natal for a doctoral scholarship for Mr K. Osman and postdoctoral scholarship for Dr I. Bahadur. This work is based upon research supported by the South African Research Chairs Initiative of the Department of Science and Technology and the National Research Foundation.

AUTHOR INFORMATION

Corresponding Author *E-mail: ramjuger@ukzn.ac.za (D. Ramjugernath). Tel.: +27 (31) 2603128. Fax: +27 (31) 2601118.

Notes

The authors declare no competing financial interest.

Supporting information: Buoyancy Correction

The weight change of the sample as gas absorption/desorption occurs is calculated by taking into account all contributing forces including the mass of the counterweight and its associated attachments (hook and chain), the sample container and its associated attachments (hook and chains), the density of the sample and that of the absorbing gas. The following equation applies:

where W is the weight reading, g is acceleration due to gravity in m.s -2 , m s is the dry sample mass, m a is the mass absorbed which needs to be determined, m c is the mass of counterweight, m I and m II are the mass of the balance components on the sample side and counterweight side respectively, and ρ f is the density of the absorbing gas. All masses and the weight reading W are measured in grams. V I , V II and V C are the volume of balance components on the sample side, the volume of balance components on the counterweight side, and the volume of the counterweight respectively. V as is the volume of the sample and the absorbed gas. All volumes are measured in cm [START_REF]Climate Change 2007: The Physical Science Basis. Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change[END_REF] .

where m a,N2 is the mass of sample and absorbed nitrogen and V N2 is the volume of sample and absorbed nitrogen.

The method of buoyancy correction is explained in detail by Macedonia et al., [START_REF] Macedonia | Adsorption Studies of Methane, Ethane, and Argon in the Zeolite Mordenite: Molecular Simulations and Experiments[END_REF] In this study, the buoyancy of the sample was corrected using N 2 gas, which is negligibly soluble in the samples in question and possesses a molecular mass more comparable to CO 2 than the use of helium gas, thus providing good accuracy is asserting the assumption of V as =V N2 .

The error in weight measurement can be expressed as:

where u(w) [START_REF] Bretherick | Hazards in the Chemical Laboratory[END_REF] and u(w) 2 is the sum of error in mass reading of each component in eqs 1 and 2, respectively.

Due to each mass of each balance component in eq 1 being accurate to ±0.001 mg, and with the need of buoyancy correction eq 2, u(w) may be expressed as: u(w) = u(w)1 + u(w)2 = 0.005 + 0.005 = ±0.01 mg (4)


Fan Zhang
email: fan.zhang@mines-paristech.frchristophe.coquelet@mines-paristech.fr

Elise El-Ahmar

Chien-Bin Soo

Xavier Canet

Christophe Coquelet

Application of GC-PC-SAFT EoS to Organic Sulfur Compounds

Keywords:

Design and optimization of separation processes require accurate knowledge of the thermodynamic properties and phase equilibria of involved pure compounds and mixtures. Thermodynamic models are thus needed to determine these properties. Model development relies on not only appropriate theory but also experimental data. However, for the organic sulfur compounds which are commonly found in diverse industrial sectors, few or even no experimental data exist in the literature. Therefore, models with predictive features may act as an alternative to handle engineering purposes.

. The predictive model (named GC-PC-SAFT) was applied to investigate two series of typical organic sulfur compounds: sulfide (R-S-R') and 1-thiol (R-SH). The group parameters of (S) and (SH) were fitted to vapor pressure and liquid density data (from [3]) of 9 sulfides and 7 1-thiols, respectively. The regression results show that the average deviations on vapor pressure are generally lower than 5%, while those on liquid density are generally lower than 2%.

•

  IntroductionGC-PC-SAFT EoS Reference  Conclusion & PerspectiveApplication of GC-PC-SAFT EoS with a dipolar term to investigate the sulfides and 1-thiols Good correlation and prediction of pure compound properties (AAD generally less than 5%) Satisfactory prediction of mixture VLE and h E data without any binary interaction parameters (k ij =0) Prediction for multi-compound systems Improvement in representing the solvation Extension to other organic sulfur compounds m (S) depends on the position of (S)•{μ, x p m} are applied directly to the entire molecule •1-thiols are considered as non auto


Alice Raeli
email: alice.raeli@math.u-bordeaux1.fr

Mejdi Azaïez

Michel Bergmann

Numerical Modelling for Phase Changing Materials (PCM)

There are several reasons for a numerical study of partial differential equations; many of them describe physical phenomena such as electromagnetism, elasticity, etc... In this work we present an approach for the heat equation with discontinuous conductivity and its relative boundary conditions. We present a finite difference method (FDM) on hierarchical meshes based on octrees (quadtrees for 2D cases) built to solve the mathematical model that describes a Phase Changing Material (PCM). We show an analysis of this method and its building reasons; then we describe the technical tools and necessary observations for a consistent cell centered method. The mesh refinement follows the hybrid material's interface looking for the accuracy on the discontinuity presented by the physical problem.

Problem

ABENGOA Solar is focused on renewable and sustainable energy. In this work we study a hybrid material composed by graphite and salt grains.

Left: Heliostats Area, Solar Tower Plant -ABEN-GOA Solar (Sevilla, Spain). Center: Phase Changing Material. Right: Tank containing the PCM. The heat is conserved in the tank and it causes the salt grains to liquefy, generating the initial energy storage. Later, the solidification transfers this heat to release the stored energy.

Mathematical Hybrid Model

Let G be the graphite and S be the PCM. k = (k S (T ), k G ) is the thermal conductivity. The holes generated by the phase change generate a discontinuity of temperature, depending on a constant contact resistance; however the flux conservation condition ensures the continuity of the normal derivatives.

The enthalpy function is strictly connected to a function that describes the phase change:

H(T ) =      [(ρC) s S (1 -f (T )) + (ρC) l S f (T ))]T +ρ l S f (T )L a in S (ρC) G T in G f (T ) = 1 2 (1 + tanh T -T F δT )

Octrees discretization tools

A quadtree is 2:1 balanced if any two neighboring nodes differ by at most 1 in depth. We suppose here to work on a quadtree-based 2D discretization. The grids we refer to in this work are graded on faces and edges with a maximal jump of two levels on the vertices. The same reasons can be developed in the 3D cases. Definition 3.1 A linear octree is an octree that is represented by a linear array instead of a tree data structure. Its structure allows an efficient memory saving (this advantage is commonly used in adaptive mesh refinement and coarsening (AMR) algorithms). Quadtree type discretizations seem intuitive for spatial distribution topology but not easy to manage. A first necessity for the treated linear octree discretization was the identification of all possible configurations. We used a Z-order method for the grids that allows us to obtain a unique way case by case. We choose a base-5 8-digits numerical key ( resp. 26-digits for 3D case) built looking at neighbours configuration of each octant following the curve orientation (example in figure).

2D key order where F stands for face and C for corner.

Using a proper functional based on the levels [L] = L-nL, with L the level of the concerned octant and nL the level of the neighbour, the value attribution on the key elements follows:

0 neighbours in this side 1 [L] = 0 2 [L] = -1 3 [L] = 1 4 [L] = -2 5 [L] = 2

Elements of an octree based cell centered FDM Definition 4.1 (Consistency)

A finite difference scheme is said to be consistent with the partial differential equation it represents, if for any sufficiently smooth solution u of this equation, the truncation error of the scheme tends uniformly towards zero with respect to x, when h tends to zero, i.e. if:

lim h→0 || h || ∞ = 0
Moreover, if there exists a constant C > 0, such that || h || < C p , with p > 0, then the scheme is said to be accurate at the order p.

First example reference configuration. The concerned point is c 4 .

For a cell centered finite difference method we start by an example (in figure) presented in [START_REF] Bergmann | A supra-convergent finite difference scheme for the variable coefficient Poisson equation on non-graded grids[END_REF] for the interested reader.

Let h be the side lenght of the cell with center c 4 , with a neighbours analysis it is simple to deduce the projected distances on x and y directions from c 4 to the neighbours centers with dependence on h. Let u(P ) be the solution function and a i the weights of the finite difference method of concerned points in a configuration; the a i are obtained such that This result is obtained by solving the sistem of the six relative constraints for the concistency. We can observe that the totality of concerned points with an edges implementation is seven so we can determine infinite solutions of the complete system:

         1 1 1 1 1 1 1 0 -h 0 -h 3h 2 -h 2 3h 2 0 h h 0 h 2 -3h 2 -3h 2 0 h 2 2 0 h 2 2 9h 2 8 h 2 8 9h 2 8 0 -h 2 0 0 3h 2 4 3h 2 4 -9h 2 4 0 h 2 2 h 2 2 0 h 2 8 9h 2 8 9h 2 8                    a 4 a 1 a 2 a 3 a 5 a 6 a 7           =         0 0 0 1 0 1        
For a described mesh the number of constraints is 6 (10 on 3D case) for the consistency, and the number of neighbours in any involved configuration is from 7 to 13; consequently the system requiring the consistency admits a solution for a stencil that involves all the neighbours of a centered point. For a unique solution we need:

• a way to identify uniquely any configuration (the chosen key presented above);

• additional constraints that ensure a solution.

The additional constraints

The idea is to ensure the consistency and minimize the deviation from second order as best as possible. Let M be the constraint matrix, ω the weights, λ the necessary additional unknowns, f the right hand side vector for consistency and F ( ω) a weight function. The problem to minimize has the lagrangian form:

L = F ( ω) -λ(M ω -f ) (1)
⇓ ∂L ∂ω = 0 ⇒ ∂F (ω) ∂ω -M T λ = 0 ∂L ∂λ = 0 ⇒ M ω = f
We write minimization problem (1) in matrix form, let be:

• x the vector that includes the weights and the λ added terms, b the right hand side vector that includes the consistency parameters ( f ) and 0 for the other additional equations.

• A the final problem matrix for searching the configuration's weights, B the matrix obtained by imposing the 6 constraints of consistency, C the submatrix of the 4 constraints with third order derivatives, I a large enough identity matrix.

Ax = b ⇔ ∂F (ω) ∂ω -M T λ = 0 M ω = f
Two cases have to be distinguished for N the number of concerned points in a configuration:

• N ≤ 10 : M = B F (ω) = 1/2((1 -α)C T C -αI)ω T ω Ax = ((1 -α)C T C -αI) -B T B 0 x = b
Ensured 6 constraints for consistency imposition and minimization of second order constraints and weights norm.

• N > 10 :

M = B C F (ω) = 1/2(ω T ω) Ax = I -B C T B C 0 x = b
Ensured 10 constraints for second order imposition and minimization of the weights norm.

Remark 5. [START_REF] Bergmann | A supra-convergent finite difference scheme for the variable coefficient Poisson equation on non-graded grids[END_REF] The consistency constraints for the 3D cases is 10 and the total number of equations for the second order accuracy is 20. Consequently the two cases are N ≤ 20 and N > 20.

Algorithm steps

• to identify the configuration on each octant;

• to build the weights matrix and keep the solution;

• to build the system of PDE solver matrix with the relative weights and data octant by octant line;

• solve the system with an assigned source term function and extract the numerical solution of the problem.

Numerical Results: Application of the Finite Element Method

Test 7.1 (Configuration Test) A first test is built on the configuration test where we looked for consistent weights applying our algorithm.

Order (left) and a result error (right) of the test configuration presented in [START_REF] Bergmann | A supra-convergent finite difference scheme for the variable coefficient Poisson equation on non-graded grids[END_REF] (dark blue stands for -7e -06).

Test 7.2 (Penalization Method) We established consistency and convergence of the method testing random grids (we suppose that we do not know a priori the molten salt disposition in the graphite structure), then we supposed that the salt has a circle composition and we focused on surface and internal part (subject to phase change) ignoring the graphite that has constant properties. So we applied the penalization method to improving a local concentration of variables approaching the heat problem of the material on microscopic way.

Error example (top) and method order (bottom): the order 1 was expected from theory.

Test 7.3 (Discontinuity Simulation)

We simulated the discontinuity of the conductivity k using a mollified heaviside function on more different surfaces for a jump from 1 to 1 million on a small number of cells, we established a control for the jump modifying his width and the related cells.

Error simulation on a bar jump of conductivity (top) and on a circle surface (center) with its relative numerical solution (bottom).

Current Works and Future Plans

Currently the same examples are replicated in 3D cases to confirm the validity of the method and to complete this part. The number of variables concerned on the weights systems is different but the complete algorithm follows the same steps until the Poisson solver. In a second moment we will apply the method to real variables of the material problem (we remark that the material exists and its conductivity properties are known) studying the phase change to optimize the energy storage for the molten salt. In addition with this work we expect to be able to perform a reduction order model for a faster solver dedicated to PCM equations.

  ∂ t H(T ) -∇ • (κ∇T ) = g κ∂ n T = 0 (Γ N ) T = T D (Γ D )[κ∂ n T ] = 0 (κ∂ n T )

7

  i=1 a i u i = u xx (c 4 ) + u yy (c 4 ) where the u i terms are obtained by taylor series expansions centered on the relative point. A first result (Min et al.) is: Theorem 4.1 Consider the discretization of the Poisson equation at the center of the cell c 4 . If only face-adjacent cells are to be used, then there does not exist any locally consistent linear scheme on non-uniform Cartesian grids.