Michel Raous
email: raous@lma.cnrs-mrs.fr

Gianpietro Del Piero
email: gianpietro.del.piero@unife.it

Adhesive Contact : a Survey and a Unified Formulation

This unified formulation is based on : -general laws, typically, energy conservation and dissipation principles, that is, mechanical versions of the first two laws of thermodynamics, -a set of state variables, that is, an array of independent variables which fully determine the response to all possible deformation processes, -a set of elastic potential and dissipation potentials, which are functions of state in terms of which the general laws take specific forms, -a set of constitutive assumptions. The behavior of the interface is first characterized by two given loading curves f n (u n ) and f t (u t ) , which are supposed to be star-shaped with respect to the origin, where u n and u t are the normal and tangential components of the relative displacement on the interface and R n and R t are respectively the components of the contact force. A state variable α is introduced to measure the current intensity of damage. The variables {u n , u t , α} have to satisfy a set of inequalities which defines the state space. To exemplify this, for the special case where there is no viscosity, the loading curve R n =f n (u n ) and the state space are given for the normal components in Fig. 1 (a) and(b). The response due to a deformation process of loadingunloading starting from the origin is represented by the dashed line in the force-displacement plane (a) and in the state space (b). As described in [START_REF] Del Piero | A unified model for adhesive interfaces with damage, viscosity, and friction[END_REF], the strain energy is then given by :

Ψ(u n , u t , α) = 1 2 g n (α)u 2 n + 1 2 g t (α)u 2 t (1)
where g n (u n ) = fn(un) un and g t (u t ) = ft(ut) ut , are the current stiffness of the interface. The dissipation potential is then deduced from the energy balance evaluated on the forcedisplacement response (see Fig. 1) :

Φ d (α, . α) = - 1 2 (g n (α) + 1 2 g t (α))α 2 . α. (2)

When the viscosity is included in the evolution of the damage, a dissipation potential is added :

where h(v) is given and defines the viscous dependence (h(α) > 0).

In order to introduce the coupling of the adhesion to the friction, the following friction dissipation power is added :

where µ is the friction coefficient which is generally a function of α .

Along with the definition of the state of admissible states for (R, u n , u t , α), one writes the energy conservation under the form of the power equation : Using the power equation (and its first or second derivative when the power equation is identically satisfied) and the previous inequalities, one can obtain the behavior law of the interface including the Signorini conditions and the equations of the evolution equation of α.

The laws presented in the survey could be obtained by choosing appropriated forms of f n , f t , h(α), µ(α). This is illustrated on the RCCM model (Raous et al [2] [3]).


Nunziante Valoroso
email: nunziante.valoroso@uniparthenope.it

Michel Raous
email: raous@lma.cnrs-mrs.fr

GRADIENT OF DAMAGE ENHANCEMENT FOR A COHESIVE MODEL

Keywords: Damage mechanics, gradient enhancement, regularization

Gradient enhancements have become increasingly popular in the last decades for dealing with problems in mechanics suffering from spurious mesh sensitivity induced by strain softening. Many proposals exist in this sense and various regularization techniques have been presented and successfully applied to study localization and fracture.

In short, the idea underlying almost all such techniques is that of using some extended constitutive equations in which information about the material microstructure is represented through a length scale-related parameter. The physical interpretation of this quantity on a micromechanical basis is still the object of an open debate, whereby its interpretation as a mere numerical regularization parameter is certainly more appropriate.

From a computational standpoint, once spatial gradients and/or length scales are introduced in the constitutive equations the latter are no longer defined at the local (quadrature point) level but they are established at a larger scale, i.e. the scale of the structural model, in a form that could be rephrased in an integral format. Basically, for usual local models stresses, strains and internal variables are defined in a point-wise fashion whereby, as outlined in [START_REF] Alfano | A displacement-like finite element model for J 2 elastoplasticity: Variational formulation and finite-step solution[END_REF], their values can be regarded as the parameters of a piece-wise constant interpolation. Hence, variables computed at the Gauss point level in classical displacement-based finite element methods can be understood as fields that are in general discontinuous across elements boundaries and inside elements as well. This discontinuous pattern is indeed one of the most striking consequences of the strictly local character of the constitutive law.

Contrariwise, for nonlocal and gradient-enhanced models the presence of gradient or averaging operators in the constitutive equations enforces a greater regularity of strains, stresses and internal variables and the resulting solution will be globally smoothed through elements. Moreover, increased regularity in general characterizes the entire solution and is not limited to the gradient-enhanced variables.

In this work reference is made to an interface cohesive model originally proposed in [START_REF] Valoroso | A damage-mechanics-based approach for modelling decohesion in adhesively bonded assemblies[END_REF]. When considering the simplest two-parameter version for the mode I case, damage evolution is governed by the following relationships:

f (D) ≤ 0; Ḋ ≥ 0; f (D) Ḋ = 0 (1)
for the damage function:

f (D) = Y + G c log(1 -D) (2)
being D ∈ [0, 1] the scalar damage variable, Y the (local) damage-driving force and G c the mode-I fracture energy. The gradient enhancement can be introduced using different arguments, see e.g. [START_REF] Frémond | Damage, gradient of damage and principle of virtual power[END_REF][START_REF] Nguyen | The non-local generalized standard approach: a consistent gradient theory[END_REF][START_REF] Pham | Approche variationnelle de l'endommagement : II. les modèles à gradient[END_REF] among others, and the evaluation of the damage state now amounts to solving the complementarity problem (1) with

g(D) = c∇ 2 D + Y + G c log(1 -D) (3)
in place of the local damage function (2). In the above relationship c is a model parameter having the dimensions of a force that can be related to the characteristic length scale. Moreover, owing to the presence of spatial derivatives in the damage function (3), the computation of the damage state requires the solution of a boundary value problem.

Numerical results that are being presented in this communications will show the superior performances of the damage enhanced model with reference to typical fracture propagation tests, where the gradient model provides meaningful answers quite independently from the mesh size.


Martha Hajiw

Antonin Chapoy
email: antonin.chapoy@pet.hw.ac.uk

Christophe Coquelet
email: christophe.coquelet@mines-paristech.fr

Hydrocarbons -Water Phase Equilibria using the CPA Equation of State with a Group Contribution Method

Keywords: Equation of state, Vapourliquid equilibrium, Group contribution method, Binary interaction parameters

It is proposed in this paper to extend the original group contribution method PPR78 to systems containing water, by combining it to the Cubic-Plus-Association (CPA) equation of state (EoS).

Applying simple geometric combination rules, the binary interaction parameter k ij (T) can be calculated from interaction parameters between hydrocarbon groups and water. This model, called the GC-PR-CPA is applied to predict hydrocarbonswater mutual solubilities over a wide temperature and pressure range, depending on available literature data.

Group interaction parameters, here CH 4 , C 2 H 6 , CH 3 , CH 2 , CH, C, CH aro , CH 2,cyclic , CH cyclic /C cyclic , C 2 H 4 , CH 2,alkene /CH alkene with H 2 O have been defined with solubility data. Predictions of the developed model have been validated against independent solubility data as well as water content in hydrocarbon rich phase. Predictions of the new model are in good agreement for light and medium hydrocarbons, however some deviations are observed for heavier hydrocarbons.

Introduction

During production, transportation and processing, changes in temperature and pressure can lead to water condensation, ice and/or gas hydrates formation. Therefore accurate knowledge of mixtures phase equilibria are important for safe and economical design and operation of pipelines and production/processing facilities. Phase equilibria are predicted with thermodynamic models. Cubic equations of state (EoS), such as the Soave-Redlich-Kwong (SRK EoS) [1] and the Peng-Robinson (PR EoS) [2] equations of state are widely used in the industry. Moreover, industry requires simple predictive models for phase diagram understanding and process design. However classical cubic EoS cannot handle correctly mixtures containing associating compounds, i.e. forming hydrogen bonding, such as water or alcohols. The Cubic-Plus-Association equation of state (CPA-EoS) [3] gives better predictions of such systems, according to Kontogeorgis et al. [3,4] . It combines a cubic EoS, usually the SRK EoS, and the Wertheim's term which describes the physical interactions between associating compounds. It is proposed in this work, to replace the SRK EoS with the PR EoS, and to use the PPR78 model developed by Jaubert and co-workers to quantify binary interactions between two hydrocarbon groups [5][6][7][8][9] . Although the PPR78 model was recently extended to aqueous systems [10] , a new temperature-dependent function is proposed in this study to predict the interactions between water and an hydrocarbon group. We indeed noticed that the addition of a Wertheim's term to the PR EoS was not compatible with the use of the existing group interaction parameters. The advantage of the PPR78 model is its simplicity and its computation time efficiency. Indeed, only component properties are needed: the critical pressure P C , the critical temperature T C and the acentric factor ω. By combining the CPA EoS and the PPR78 model it is expected to have, in a simple way, better and physical reliable predictions for mixtures with water and alcohols. Thus it has been evaluated by predicting mutual solubilities of hydrocarbonswater and compared to the SRK-CPA EoS with adjusted binary interaction parameters or set to zero and two other SAFT type models: the GCA EoS [11] and GC-PPC-SAFT EoS [12] .

The GC-PR-CPA model

Species forming hydrogen bonds (water, alcohols and glycols), also called associating compounds, present unusual thermodynamic behaviours, due to the strong attractive interactions between molecules. Classical equations of state combined with classical mixing rules are well known not to handle correctly the phase behaviour of mixtures containing water and/or alcohols. In this work the Cubic-Plus-Association (CPA) EoS, originally developed by Kontogeorgis and co-workers [3] , is used.

In the original version this model combines the SRK cubic equation of state and the Wertheim's theory. In this work, the SRK equation of state is replaced by the Peng Robinson (PR) EoS [13] , justified by its better predictions of liquid densities and phase diagram in the critical region [14] . The binary interaction parameter calculated with the PPR78 EoS can be used with any other cubic EoS by using the mathematical relationship developed by Jaubert et al. which makes it possible to switch the binary interaction parameters from a cubic EoS to another one [15] .

The PR-CPA pure compounds parameters

The PR-CPA EoS is expressed here, in term of pressure for a pure compound (Eq. ( 1

where T C,i and P C,i are respectively the critical temperature and the critical pressure of the component i, ω is the acentric factor, a i is the energy parameter and b i the co-volume.

In the Wertheim's equation g(d) simpl (Eq. ( 3)) is the simplified expression of the radial distribution function as suggested by Kontogeorgis et al. [3] , η the reduced fluid density, ρ the molar density, x is the mole fraction of the component i and X Ai the mole fraction of the molecule i not bonded at a site A (Eqs. ( 4) and ( 5)).

𝑔 𝑑 𝑠𝑖𝑚𝑝 = (3)

𝑋 𝐴 𝑖 = 1 1+𝜌 𝑛 𝑗 𝑋 𝐵 𝑖 ∆ 𝐴 𝑖 𝐵 𝑖 𝐵 𝑖 𝑗 (4) ∆ 𝐴 𝑖 𝐵 𝑖 = 𝑔 𝑑 𝑠𝑖𝑚𝑝 𝑒𝑥𝑝 𝜀 𝐴 𝑖 𝐵 𝑖 𝑅𝑇 -1 𝛽 𝐴 𝑖 𝐵 𝑖 𝑏 (5)
Therefore, parameters a 0,i , b i , C 1,i , the association energy ε A i B i and the association volume β A i B i are the five parameters of the CPA EoS. For non associating compounds, the three PR parameters are calculated with the critical temperature, critical pressure and acentric factor and ε A i B i and β A i B i are equal to zero. However, to improve predictions for associating compounds, these parameters have been fitted to vapour pressure and saturated liquid density data, using the objective function defined by Eq. ( 6). X i is either the vapour pressure or the saturated liquid density and NEXP the number of experimental points:

𝐹 = 1 𝑁𝐸𝑋𝑃 𝑋 𝑖,𝑒𝑥𝑝 -𝑋 𝑖,𝑐𝑎𝑙𝑐 0.01× 1+𝑋 𝑖,𝑒𝑥𝑝 2 𝑁𝐸𝑋𝑃 𝑖=1 (6)
PR-CPA parameters for water are given in Table 1.

According to the terminology of Huang and Radosz [16] , the four-site (4C) association scheme has been applied for water.

Mixing rules

For mixtures, the equation parameters are defined with mixing rules. In this work, the classical mixing rule (van der Waals one fluid theory) has been applied for a(T) and b (Eqs. (7) and ( 8)).

𝑎 𝑇 =

𝑥 𝑖 𝑥 𝑗 𝑎 𝑖 𝑇 𝑎 𝑗 𝑇 1 -𝑘 𝑖𝑗 𝑁 𝑗 =1 𝑁 𝑖=1 (7)

𝑏 =

𝑧 𝑖 𝑏 𝑖 𝑁 𝑖=1 (8) where N is the number of components in the mixture.

Folas et al. [17] mentioned that the combining rule CR1 ((Eqs. ( 9) and ( 10)) is one of the best mixing rules for the determination of ε and β. CR1 is used in our GC-PR-CPA EoS.

𝜀 𝐴 𝑖 𝐵 𝑗 = 𝜀 𝐴 𝑖 𝐵 𝑖 +𝜀 𝐴 𝑗 𝐵 𝑗 2 (9)
𝛽 𝐴 𝑖 𝐵 𝑗 = 𝛽 𝐴 𝑖 𝐵 𝑖 𝛽 𝐴 𝑗 𝐵 𝑗 (10)

The group contribution method

As seen in Eq. ( 7), the PR-CPA EoS applied for mixtures involves a binary interaction parameter k ij .

It is usually set to zero or defined by correlations specific to binary systems. In this work, k ij is calculated with the PPR78 model. Between 2004 and 2013, Jaubert and co-workers [5][6][7][8][9][10][18][19][20][21][22][23][24][25] proposed a group contribution method and defined parameters of 21 groups to estimate temperature dependent Binary Interaction Parameters (Eqs. ( 11) and ( 12)) for mixtures involving hydrocarbons, water, mercaptans and permanent gases (N 2 , CO 2 , H 2 S, H 2 ).

𝑘 𝑖𝑗 𝑇 = - (12) where N g is the number of groups, k and l are groups, A kl and B kl the interaction parameters between groups (A kl =A lk and B kl =B lk ) and α ik is defined by Eq. ( 13):

𝛼 𝑖𝑘 = 𝑛𝑢𝑚𝑏𝑒𝑟 𝑜𝑓 𝑔𝑟𝑜𝑢𝑝 𝑘 𝑖𝑛 𝑚𝑜𝑙𝑒𝑐𝑢𝑙𝑒 𝑖 𝑡𝑜𝑡𝑎𝑙 𝑛𝑢𝑚𝑏𝑒𝑟 𝑜𝑓 𝑔𝑟𝑜𝑢𝑝𝑠 𝑖𝑛 𝑚𝑜𝑙𝑒𝑐𝑢𝑙𝑒 𝑖 (13)
The group H 2 O has been added to the model. As shown in Figure 1, the binary interaction parameters, adjusted on literature data temperature by temperature, presents a polynomial tendency for hydrocarbons -water systems, similar to a Henry's law for such systems [26] . The trend of this parameter cannot be reproduced using Eq. (12). Therefore the term Sum PPR78 has been modified for the group H 2 O with the introduction of three new interaction parameters C kH2O , D kH2O and E kH2O (Eq. 14). Table 2 presents k ij values for different temperatures between 283.15 K and 507.15 K.

𝑆𝑢𝑚 𝑎𝑠𝑠𝑜 = 𝛼 𝑖𝑘 -𝛼 𝑗𝑘 𝛼 𝑖𝐻 2 𝑂 -𝛼 𝑗 𝐻 2 𝑂 𝐶 𝑘𝐻 2 𝑂 𝑇 2 + 𝐷 𝑘𝐻 2 𝑂 𝑇 + 𝐸 𝑘𝐻 2 𝑂 + α 𝑖𝐻 2 𝑂 - 𝐻 2 𝑂 𝑙=1 𝐻 2 𝑂 𝑘=1 𝛼 𝑗 𝐻 2 𝑂 𝛼 𝑖𝑙 -𝛼 𝑗𝑙 𝐶 𝐻 2 𝑂𝑙 𝑇 2 + 𝐷 𝐻 2 𝑂𝑙 𝑇 + 𝐸 𝐻 2 𝑂𝑙 (14) With C k,H2O =C H2O,l , D k,H2O =D H2O,l , E k,H2O =E H2O,l .
The model proposed has therefore five parameters: two parameters (A kl and B kl ) for non-associating groups and three parameters (C kl , D kl and E kl ) for associating groups such as water. Therefore, if a binary system without associating compounds is studied, e.g. ethanepropane, the model is reduced to the PPR78 EoS: the associating part of the PR-CPA EoS is not used and the binary interaction parameter is the one defined by Eqs. ( 11) and (12). However, for a system with water, e.g. methanewater, the PR-CPA EoS is applied and the binary interaction parameter is calculated with Eq. (11) together with the new term Sum asso (Eq. ( 14)).

In this work, the interaction parameters C kH2O , D kH2O and E kH2O are defined for CH 4 , C 2 H 6 , CH 3 , CH 2 , CH, C, CH aro , CH 2,cyclic , CH cyclic /C cyclic , C 2 H 4 , CH 2,alkene /CH alkene groups in order to define k ij (T) for hydrocarbons-water mixtures. They have been calculated using a flash calculation, by minimizing the objective function defined in Eq. ( 15):

𝐹 = 1 𝑁𝐸𝑋𝑃 𝑥 𝑖,𝑒𝑥𝑝 -𝑥 𝑖,𝑐𝑎𝑙𝑐 0.01× 1+𝑥 𝑖,𝑒𝑥𝑝 2 𝑁𝐸𝑋𝑃 𝑖=1 (15)
where N EXP is the number of experimental data, x i,exp the experimental liquid mole fraction of hydrocarbon and x i,calc the calculated one. Table 3 presents the group interaction parameters obtained with Eq. ( 15). The objective function was limited to hydrocarbon solubility data, because they are more widely available in the literature. Indeed, in the case of VLE, measuring gas solubility is easier than measuring the water content [27] . Furthermore, the binary interaction parameter for systems with water has limited impact on water content predictions (see Figure 3). However, as shown on Figure 2, for the system methanewater, the SRK-CPA EoS with a k ij set to zero over predicts methane solubility, while the SRK-CPA EoS with adjusted binary parameters and the GC-PR-CPA EoS are in good agreement with experimental data. Therefore, adjusting the binary interaction parameter, using correlations or a group contribution method, improves solubility predictions.

Validation of the model

To evaluate the accuracy of the GC-PR-CPA model, solubilities of different hydrocarbons, such as normal alkanes, branched alkanes, aromatics naphthenics, and alkenes in water are compared to experimental data. The model has also been compared to the SRK-CPA EoS with adjusted binary interaction parameters or set to zero. Adjusted binary interaction parameters (BIP) for the SRK-CPA EoS are temperature dependent polynomial functions. They have been adjusted on VLE or LLE data and used directly from our software (Hydraflash™1 ) without any retuning for this work. When results were available, the model has been compared to two SAFT type models: the GCA EoS [11] and the GC-PPC-SAFT EoS [12] . As mentioned before, the GC-PR-CPA EoS has five parameters for pure associating compounds to be adjusted. These parameters are the same as those given by Kontogeorgis et al. in their paper [3] for the SRK-CPA EoS. It is also the case for the SAFT EoS [28] : the energy parameter, the segment diameter, the chain length parameter, the association energy and the association volume parameter are the five parameters to be adjusted for pure compounds as well as the binary interaction parameter for mixtures. Furthermore, it is interesting to point out that heavier hydrocarbons, such as n-pentane, n-hexane, benzene, toluene and ethylbenzene present also a minimum of solubility in water. Tsonopoulos [29] defined correlations for alkanes and aromatics solubility in water as well as the temperature corresponding to the minimum of solubility. In this paper, the ability of the GC-PR-CPA model to reproduce the minimum of solubility is evaluated for aromatics. The model is compared to the original correlation defined by Tsonopoulos (Eqs. ( 16) and ( 17)) and to the modified correlation with readjusted parameters 𝒜, ℬ and 𝒞.

𝑙𝑛 𝑥 𝑖 = 𝒜 𝑖 + ℬ 𝑖 𝑇 + 𝒞 𝑖 𝑙𝑛 𝑇 (16) 𝑇 𝑚𝑖𝑛 = ℬ 𝑖 𝒞 𝑖 (17) Finally, in process simulation, EoS are not only used to determine phase diagram but must be able to correctly estimate thermodynamic properties like densities, enthalpies or estimate Henry's law constant [26] at infinite dilution. As an example, the capability of the GC-PR-CPA EoS to predict the mixing enthalpy h M is evaluated here for the benzene-water system.

Database

Several solubility data have been found for lighter alkanes and aromatics with water. Comparison between data has been made, and data which did not follow the major trend of other experimental data have been deleted. Smoothed data for methane [30] and ethane [31] , and data covering the larger range of temperature and pressure have been used for calculations. The model has been validated with other VLE and/or LLE data, if available. There are fewer data for branched alkanes and naphthenics, but the same process has been applied to screen the data. Finally for heavier components (n-decane and heavier), very few data are existing, but their accuracy is not proved, since their solubility in water is very small (lower than 1.10 -9 ). In literature, the uncertainty on experimental data is estimated to be between 3 % and 5 % for lighter hydrocarbons, and up to 16 % for alkanes greater than C 9 .

Results and discussion

Methanewater

Group parameters C CH4,H2O , D CH4,H2O and E CH4,H2O for the binary system methanewater have been calculated with smoothed solubility data [30] for temperatures from 283.15 to 423.15 K and up to 200 MPa. Figure 2 shows the methane solubility in water for three isotherms. The AAD (Average Absolute Deviation) for each isotherm (473.15 K, 423.15 K, and 310.93) are 4 %, 2 % and 2 %, respectively.

For smoothed data used in calculations, the AAD is 3 % and the Average Absolute Error AAE is 8.2.10 -5 . Thus, with error on experimental data between 3 % and 5 %, it can be concluded that the model GC-PR-CPA reproduced with accuracy the methane solubility in water. And, as shown on Figure 3, the water solubility in the methane rich vapour phase is also well reproduced.

Ethanewater

In their paper, Mao and co-workers [31] presented ethane solubility data in water for temperatures from 273.15 to 473.15 K and pressures up to 100 MPa, with an experimental uncertainty of 7 %.

Parameters C C2H6,H2O , D C2H6,H2O and E C2H6,H2O have been defined with these smoothed data As seen on Figure 4, the model reproduces the ethane solubility in water, with reasonable accuracy (the AAD is 3 % and the AAE 3.9.10 -5 ). Moreover compared to the two other predictive models (Figure 5), the GC-PR-CPA EoS follows the experimental trend with temperature as well as the GCA EoS [11] (7.7 % and 6 % deviations respectively). The GC-PPC-SAFT [12] presents the highest deviation with 25%, even if dipole-dipole interactions are taken into account in this equation. Despite the deviations at high pressures, inherent to the CPA EoS, the model is in good agreement with the experimental data. The deviation at higher pressures could be attributed to the neglect of "crossassociation" between ethane and water.

Other alkanes with water

 Phase equilibrium predictions

Parameters have been obtained using experimental data for the following binary systems: VLE data for propanewater [32] and n-butanewater [33] , LLE data for n-hexanewater [34] and n-nonanewater [35] . For lighter compounds, such as propane, predictions are in good agreement with experimental data, as shown on Figure 6 and Figure 7. The AAD on all data for propane is 12 % and the average of the AAE is 1.7.10 -5 (the maximum being 1.4.10 -4 ). The GC-PR-CPA EoS is better than the GC-PPC-SAFT model (22% deviation) but presents deviations with experimental data (4 % deviation against 7 % for the GCA EoS). The SRK-CPA EoS with adjusted k ij is the closest to experimental data with 2% of deviation. As a group contribution method, the GC-PR-CPA EoS looks to predict with higher accuracy than the SAFT type EoS. Propane solubility is better predicted with the GC-PR-CPA EoS at temperatures higher than 320 K.

N-hexane solubility is given on Figure 8. The AAD for all data is 15% and the AAE is 6.7.10 -7 (the maximum being 5.0.10 -6 ). As seen on Figure 9, the GC-PR-CPA model is better than the two other models, with 14 % deviation against 20 %. The SRK-CPA EoS presents 15 % for these data. For other alkanes up to n-heptane, the AAD is lower than 20 %. For these binary systems, the model is accurate and represents the minimum of solubility, but it shows some discrepancies at high temperatures.

Due to a lack of data and because of the uncertainty of experimental data, predictions are poor for all alkanes from n-decane to eicosane. It should be pointed out the solubility of heavy alkanes is extremely small, i.e. lower than 1.10 -9 . Therefore, it is very difficult to represent their behaviour.

 Discussions

As seen on Figure 5, Figure 7 and Figure 9, other SAFT group contribution models do not reproduce alkanes behaviour in water as well as the CPA EoS (SRK-CPA or GC-PR-CPA). Referring to two studies of Privat et al. [36] and Polishuk et al. [37] on SAFT EoS, inconsistent predictions for pure compounds were highlighted: they have shown that a SAFT type equation of state, in their paper the PC-SAFT EoS, predicts two fluid-fluid saturation curves, two critical points and a LLV triple point [36] .

Moreover, some combinations of fitted SAFT parameters seem to lead to unrealistic phase behaviour at temperatures and pressures of industrial processes. It is thus mandatory to implement a specific algorithm to avoid these problems. These phenomena are not observed with the GC-PR-CPA EoS and the model is accurate at ambient conditions and even at low temperatures. The binary interaction parameters being adjusted for the binary system propane-water and being temperature dependent, it is evident that the SRK-CPA EoS gives better predictions. But it is interesting to see that the GC-PR-CPA EoS is more accurate than what is expected for a group contribution method.

Results for branched alkanes

Group interaction parameters C CH,H2O , D CH,H2O , E CH,H2O have been determined by combining solubility data of i-butane, i-propane and 2,3-dimethylbutane and C C,H2O , D C,H2O , E C,H2O with 2,2dimethylpropane, 2,2-dimethylbutane, 2,2-dimethylpentane and 3,3-dimethylpentane. Limited data are available for the solubility of branched alkanes in water and generally can only be found at 0.1 MPa. So, considering the scarcity of data, higher deviations should be expected. The model represents well i-butane solubility in water at 0.1 MPa (5 % deviation). It gives the right behaviour but over predicts i-pentane solubility: the AAD is 22 % and the AAE is 3.1.10 -6 (the maximum being 4.6.10 -6 ).

Results for 2,2-dimethylpropane, 2,2-dimethylbutane and 2,2-dimethylpentane are given in Table 4.

Predictions for 2,3-dimethylbutane [34] and 3,3-dimethylpentane [38] are shown on Figure 10, the AAD being respectively 10 % (AAE=9.0.10 -7 ) and 11 % (AAE=3.3.10 -7 ). Binary interaction parameters have been adjusted for the SRK-CPA EoS for 2,3-dimethylbutane but not for 3,3-dimethylpentane (Figure 10), explaining why the k ij is equal to zero and why the SRK-CPA gives 8% deviation for 2,3dimethylbutane but 79% for 3,3-dimethylpentane. Regarding to predictions, the low solubility values and the quantity of data available, the model can be considered accurate for these binary systems.

Results for naphthenics

Limited data are available for cyclic components with water. Only cyclohexane-water and methylcyclohexane-water binary systems present reliable data. As seen on Figure 11, the behaviour is well represented for cyclohexane and methylcyclohexane at 0.1 MPa. The AAD are respectively 8 % and 19.5 % and the AAE 2.1.10 -6 and 1.1.10 -6 . In comparison, the SRK-CPA EoS with adjusted binary interaction parameters presents 9.5% deviation for cyclohexane and 7% for methylcyclohexane. Predictions for cyclohexane are similar for both models, but the group contribution method degrades slightly the results. But considering the small values of solubilities, the model can be considered accurate for these binary systems.

Results for alkenes

Ethylene is a group in the PPR78 model. Parameters for C 2 H 4 group with water have been adjusted with VLE data. Four isotherms (308.15 K, 310.93 K, 327.59 K, and 348.15 K) are presented on Figure 12 and Figure 13. The AAD are respectively 1 %, 4 %, 8 % and 2 %. For all data, the AAD is less than 7 % and the AAE is about 6.8.10 -5 . The SRK-CPA EoS presents 8% deviations for ethylene solubility. Considering both models, they are similar, with better predictions at higher pressures for the GC-PR-CPA EoS. Therefore, it can be considered that the model GC-PR-CPA represents with good accuracy the ethylene solubility in water.

Parameters for the group CH 2,alkene / CH alkene have been defined using VLE data for propylene. Figure 14 shows propylene solubility for two isotherms, 310.93 K and 377.59 K. The AAD are respectively 4 % and 6 %, and 8 % for all data. The AAE for propylene solubility is 3.3.10 -5 . Both GC-PR-CPA and SRK-CPA EoS represent well the propylene solubility in water, with a small advantage for the SRK-CPA EoS. Water content data in propylene have also been found in the literature. Results are given on Figure 15 for two isotherms as well. The AAD for all data is 27 % and the AAE 5.7.10 -3 . Water content values vary from 3.65.10 -3 to 2.59.10 -1 . As seen on Figure 15, the model does not represent water content as well as the propylene solubility, but it gives the right shape, contrary to the SRK-CPA EoS which has only 7.5% deviation with experimental data. Higher deviations for the GC-PR-CPA EoS could be explained by the interactions between all groups present. Here, for propylene, only parameters for CH 2,alk / CH alk with H 2 O and CH 3 with H 2 O are determined. But there are also interactions between CH 2,alk / CH alk and CH 3 and it can affect predictions. Indeed, referring to the paper on alkenes by Jaubert et al. [8] , both liquid and vapour phases could not be accurately predicted some for alkene-alkane systems and it could therefore explain deviations in water content predictions.

Results for aromatics  Phase equilibrium predictions

To define the group interaction parameter for aromatics, the interaction between the group CH aro and water has been considered. Thus, parameters C CHaro,H2O , D CHaro,H2O and E CHaro,H2O have been calculated using liquid-liquid equilibrium data for benzenewater. Data available in the literature are mostly given at 0.1 MPa. So parameters obtained give good benzene solubility in water at this pressure (Figure 16). The AAD for these data is 2 % and the AAE is 7.5.10 -6 . Data have also been measured for pressures from 0.5 to 10 MPa for a Gas Processors Association project [START_REF] Valtz | Book Amine BTEX Solubility[END_REF] and the AAD is 5 % and the AAE 3.7.10 -5 . These results show the good accuracy of the model for this system. However, the model under predicts water solubility in the hydrocarbon phase (Figure 17), with 43% deviation at atmospheric pressure.

Groups C Haro and C aro are present in toluene, ethylbenzene and xylenes. After defining C CHaro,H2O , D CHaro,H2O and E CHaro,H2O , it is possible to adjust these three parameters for C aro -H 2 O with liquid-liquid equilibrium data. Deviations for toluene (Figure 18) and ethylbenzene (Figure 19) are respectively 9 % and 10 % on all data, and the AAE 1.0.10 -5 and 3.4.10 -6 . For data presented on figures 18 and 19, the AAD are lower with 0.9% and 7% deviations respectively. There are more discrepancies for xylenes (11 % for m-xylene and p-xylene and 40 % for o-xylene of deviations and the AAE are between 3.6.10 -6 and 1.2.10 -5 ) but the minimum of solubility is always represented. Higher deviations for o-xylene are explained by the scatter in the experimental data: no agreement between references is observed.

Therefore, the GC-CPA-PR model can represent aromatics solubilities in water at 0.1 MPa, but it under predicts water solubility in hydrocarbon liquid phase.

 Minimum of solubility

Using the experimental data presented on Figure 16, Figure 18 and Figure 19, parameters 𝒜, ℬ and 𝒞 from Tsonopoulos correlation (Eq. ( 16) have been readjusted for benzene, toluene and ethylbenzene.

The temperature of minimum of solubility has been calculated for the GC-PR-CPA model and using Eq. ( 17) with readjusted parameters ℬ and 𝒞. Predictions have been compared to experimental data (Figure 16, Figure 18 and Figure 19) and to the original correlation published by Tsonopoulos. All parameters are given in Table 5. The AAD between experimental data and the readjusted correlation are 1.5% for benzene, 3% for toluene and 2.5% for ethylbenzene. In their paper, Bohon and Claussen [START_REF] Bohon | [END_REF] define the minimum of solubility near 291.15 K for all liquid hydrocarbons studied: benzene, toluene, ethylbenzene, m-xylene, p-xylene, biphenyl and naphtalene. As for Tsonopoulos [29] , the temperature of the minimum of solubility is 291 K, but according to Tsonopoulos and Wilson [41] the minimum of solubility for benzene is 283.74 K. Referring to Table 5, temperatures of minimum of solubility, calculated with the readjusted correlation present some discrepancies with their value. The temperature of the minimum of solubility has also been calculated for the GC-PR-CPA EoS, being 285.69 K for benzene, 288.91 K for toluene and 294 K for ethylbenzene. The calculated minimum of solubility of benzene is closer to the one presented by Tsonopoulos and Wilson [41] than the one by Bohon and Claussen [START_REF] Bohon | [END_REF] or Tsonopoulos [29] . The difference in the values could be explained by the experimental data used for parameters adjustment. To show the importance of temperature dependency, a constant BIP has been adjusted on solubility data for the SRK-CPA EoS for the system benzene-water. As seen on Figure 16, the SRK-CPA EoS with a constant BIP does not predict a minimum of solubility and predictions are far from experimental data.

 Mixing enthalpy

Mixing enthalpies, as defined by Privat and Jaubert [42] , were calculated for the binary system benzene water and compared to experimental data at 592 K and 16.4 MPa [43] , measured by Wormald and Slater and mistakenly called excess enthalpies. As seen on Figure 20, the mixing enthalpy behaviour is well reproduced by both the SRK-CPA and GC-PR-CPA EoS and it has to be pointed out that the BIPs have not been adjusted using enthalpy data. The SRK-CPA EoS slightly under predicts the mixing enthalpies, while the GC-PR-CPA EoS under predicts only for a water mole fraction between 0.6 and 0.86. At 592 K and 16.4 MPa, the maximum experimental value is 13.7 kJ.mol -1 . The SRK-CPA EoS predicts a slightly lower value of 12.5 kJ.mol -1 and the GC-PR-CPA EoS 13.1 kJ.mol -1 . The maximum absolute deviation is 3.0 kJ.mol -1 for the SRK-CPA EoS and 1.2 kJ.mol -1 for the GC-PR-CPA EoS. To improve the predictions, interaction parameters can also be adjusted on mixing enthalpies data.

Conclusion

A model combining the PR-CPA EoS and a group contribution method has been developed. To be extended to systems with water, the original equation defining k ij (T) has been modified. Three new group interaction parameters had to be defined between groups already existing and the last group H 2 O. In general, the model gives accurate predictions for lighter alkanes (lower than n-decane) and gives the right behaviour for heavier compounds. Considering the low solubility values for heavy alkanes, it is difficult to have correct predictions. It can also predict other hydrocarbons solubilities such as aromatics, naphthenics and alkenes with a good representation of the minimum of solubility, and properties of mixing, such as the enthalpy. The GC-PR-CPA EoS does not degrade predictions as much as it would be expected from a group contribution method, but it is in general slightly less accurate than the SRK-CPA EoS with adjusted binary interaction parameters.

List of symbols

List of abbreviations

1 𝜂

 1 of the PR EoS are expressed below (Eq. (2)). R = 8.314472 J.mol -1 .K -= 𝑎 𝑐,𝑖 1 + 𝐶

  a(T) temperature dependent energy term in the PR term (MPa.L².mol -²) a 0 parameter in the energy term (MPa.L².mol -²) A i site A in the molecule i A kl , B kl constant parameters in the PPR78 expression of binary interaction parameters b co-volume parameter C kl , D kl , E kl constant parameters in the GC-PR-CPA expression of binary interaction parameters

http://www.hydrafact.com/

Acknowledgments

This work was funded by Chevron, GDF-Suez, Petrobras, Statoil and Total, whose support is gratefully acknowledged.

The authors would like to acknowledge also Dr Eric Auger for the utilisation of the GC-PPC-SAFT.

Martha Hajiw would like to thank Heriot Watt University for the use of Hydraflash™ and the funding.


Alfonso Gonzalez

Luis Pereira

Patrice Paricaud

Christophe Coquelet

Antonin Chapoy

Modeling of Transport Properties Using the SAFT-VR Mie Equation of State

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Nowadays Carbon capture and storage (CCS) is being studied as one of the most promising process to reduce the CO 2 emissions from burning fossil fuels. Several research projects have been focused in improving the safety and reducing both the costs and the environmental risks associated with CCS operations [START_REF] Li | Viscosities, Thermal Conductivities and Diffusion Coefficients of CO2 Mixtures: Review of Experimental Data and Theoretical Models[END_REF]. From a thermodynamic point of view, an accurate description of thermophysical properties is crucial to the design of capture, transport and storage of CO 2 . Therefore, the research and the assessment of new equation of state (EoS) models are important for the deployment of more accurate models.

In the oil and gas industry, with regard to the upstream, there is no interest in looking for more accurate EoSs, because in reservoir simulations the cubic EoSs provide enough accuracy due to uncertainties in the geology and the flow in porous media. Likewise, considering the computational time, thanks to their simplicity, cubic EoSs are scarcely replaced by more complex PVT models in reservoir simulations. However in the downstream industry and CCS, more complex thermodynamic models can be applied in order to increase the accuracy of PVT modeling: for example non-cubic EoSs can provide a better density description at high pressures and high temperatures (HPHT) [START_REF] Yan | -PVT Modeling of Reservoir Fluids Using PC-SAFT EoS and Soave-BWR EoS[END_REF].

This paper focuses on one of the latest versions of the family of SAFT EoS, the SAFT-VR Mie. The adequacy of this EoS to describe transport properties like viscosity and interfacial tension (IFT), when coupled with appropriated models, has been investigated with results showing a very good agreement predictions and experimental data. The investigated systems include five binaries (CO 2 /O 2 , CO 2 /N 2 , CO 2 /Ar, CO 2 /n-C 4 and CO 2 /n-C 10 ) and two multicomponent mixtures (90%CO 2 / 5%O 2 / 2%Ar / 3%N 2 and 90%CO 2 / 6%n-C 4 / 4%n-C 10 ) in a wide range of pressure and temperature conditions.

Model Description

Equation of state

The SAFT-VR Mie EoS proposed by Lafitte et al. 2013 is one of the latest versions of the SAFT equations of state. The SAFT-VR Mie is featured by the use of the Mie potential to describe the attraction-repulsion interactions between the segments which build the molecules. The Mie potential is defined as:

                                a r a r r r r r u a r a r r Mie              ) ( (Eq. 1)
where r is radial distance, ζ the temperature-independent segment diameter, ε the potential depth; λ r and λ a are the repulsive and attractive ranges, respectively. The Mie parameters for the substances studied in this work are reported in the Table 1. The SAFT-VR Mie EoS can be expressed in terms of the reduced Helmholtz energy as the summation of several contributions:

ASSOC CHAIN MONO IDEAL a a a a NkT A a      (Eq. 2)
where a IDEAL is the ideal gas contribution, a MONO the monomer contribution, a CHAIN the contribution due to chain formation and a ASSOC the association contribution. The expressions of each contribution can be found in the studies of Lafitte et al. 2013. Regarding to the phase equilibria, the thermodynamic criterion for phase equilibrium is the equality of chemical potentials of each component in all co-existing phases. Using the SAFT EoS, the chemical potential can be derived as µ=(∂A/∂N) V,T and then used to calculate the fugacity of each component [START_REF] Chapoy | Effect of Impurities on Thermophysical Properties and Phase Behaviour of a CO2-Rich System in CCS[END_REF]. In addition to the pure molecular parameters, temperature independent binary interaction parameters k ij have been adjusted against solubility data (database [START_REF] Knapp | Vapor-Liquid Equilibria for Mixtures of Low Boiling Substances[END_REF] to improve the description of vapor-liquid equilibria. The optimum k ij parameters used in Eq. 3 are reported in Table 2. The VLE predictions using SAFT-VR Mie with the adjusted k ij are presented in the Fig. 1. The percentage absolute average deviation (%AAD) between the predicted saturation pressure and vapor composition and experimental data of Nagarajan are 2.8%Psat and 0.9%y.

 

j i ij j i ij ij k       3 3 3 1  (Eq. 3)

Viscosity model TRAPP (TRAnsport Properties Prediction

) is a predictive model based on the extended corresponding states theory (ECS) (K. C. M. and K. E. [START_REF] Gubbins | Conformal Solution Theory for Viscosity and Thermal Conductivity of Mixtures[END_REF][START_REF] Hanley | Analysis of the Transport Coefficients for Simple Dense Fluids: The Diffusion and Bulk Viscosity Coefficients[END_REF], used to estimate the viscosity and thermal conductivity of pure fluids and their mixtures over the entire phase range. In the original TRAPP, propane was employed as reference fluid [START_REF] Ely | Prediction of Transport Properties. 2. Thermal Conductivity of Pure Fluids and Mixtures[END_REF], although several reference fluids can be selected [START_REF] Huber | Prediction of the Thermal Conductivity of Refrigerants and Refrigerant Mixtures[END_REF]. Accordingly, the residual viscosity of the mixture at a corresponding state point (T 0 and ρ 0 ) is given by:

ENSKOG R R m m m F T                ] [ ) , ( 0 0
(Eq. 4)

where η m 0 is the viscosity of the mixture at low pressure evaluated by [START_REF] Herning | Calculations of the Viscosity of Technical Gas Mixtures from the Viscosity of Individual Gases[END_REF] approximation. The term η R -η R 0 is the residual viscosity of the reference fluid which is calculated as (Younglove and Ely 1987):

1 3 , 5 . 0 0 2 1 . 0 0 1 0 ] ) 1 ( exp[ G G G G R r R R           (Eq. 5)
where η R is the real viscosity at T 0 and ρ 0 of the reference fluid, η R 0

is the viscosity at low pressure and T 0 , and ρ r,R is the reduced density calculated as ρ r,R = ρ 0 / ρ c,R . The corresponding state (T 0 and ρ 0 ) is calculated as T 0 =T/f m and ρ 0 =ρh m . G 1 , G 2 and G 3 are parameters that can be calculated as:

  T E E G 2 1 1 exp   (Eq. 6) 5 . 1 4 3 2 T E E G   (Eq. 7) 2 7 6 5 3 T E T E E G    (Eq. 8)
where the E i parameters have been correlated for propane as reference fluid and are shown in Table 3. In order to calculate the terms F ηm and ∆η ENSKOG , the following mixing rules were applied:

  i j ij j i m h y y h (Eq. 9) ij i j ij j i m m f h y y h f   Eq. 10)   8 ) ( ) ( 3 / 1 3 / 1 j i ij h h h   (Eq. 11) 2 / 1 ) ( j i ij f f f  (Eq. 12)
where y i is the mole fraction of component i; f i and h i are functions of the critical parameters and acentric factor defined as:

)] ln 7498 . 0 05203 . 0 )( ( 1 [ r R R C C i T T T f       (Eq. 13) )] ln 2822 . 0 1436 . 0 )( ( 1 [ r R C R C C R C i T Z Z h         (Eq. 14)
The viscosity dimensional scaling factor can be calculated as:

3 / 4 2 / 1 2 2 / 1 ) ( ) ( ) ( ) ( ij ij i j ij j i m R m h M f y y h M F      (Eq. 15) with j i j i ij M M M M M   2 (Eq. 16)
where M i denotes the molecular weight.

The expression ∆η ENSKOG is the residual viscosity which takes into account the differences in the molecules size based on the hard sphere assumption [START_REF] Poling | The Properties of Gases and Liquids[END_REF]:

ENSKOG x ENSKOG m ENSKOG      
(Eq. 17)

with i i i hs ij ij i j ij j i ENSKOG m Y g y y          0 6 2 (Eq. 18)
where the α parameter is α=9.725x10 -7 and ζ i is the hard-sphere diameter calculated as

ζ i =4.771 h i 1/3
. The radial distribution function of the hard-sphere fluid is given by the following equation [START_REF] Gubbins | Perturbation Theory for the Angular Pair Correlation Function in Molecular Fluids[END_REF]:

      2 3 2 2 1 2 1 3 1 1 ij ij hs ij g              (Eq. 19) where 3 2 2 k k k k k k ij j i ij y y          (Eq. 20) 3 6 k k k y       (Eq. 21)
To calculate the expression Σβ i Y i it is necessary to solve the linear system of equations of the form:

i i i ij Y B    (Eq. 22) with             hs ij ij j i j j j i i g M M M y y Y 3 15 18 1    (Eq. 23)                            jk k i ij k i k k i k ik ik k i ij M M M M M M M g y y B    3 2 3 5 1 2 2 0 (Eq. 24)
where δ ij is the Kronecker function which is a delta function (i.e., δ ij =1 if i=j and δ ij =0 in other cases).

Interfacial tension model

The density gradient theory (DGT) has been used to compute the interfacial tension. This model is based on the square gradient of van der Waals [START_REF] Rowlinson | Translation of J. D. van Der Waals' ‗The Thermodynamik Theory of Capillarity under the Hypothesis of a Continuous Variation of Density[END_REF]) and on the reformulation of Cahn and Hilliard 1958 to compute interfacial tension values from bulk phase properties such as density and composition. The DGT, when coupled with appropriate thermodynamic models, has been successfully applied in the prediction of interfacial properties of a wide class of systems and interfaces. The reader is referred to the study of [START_REF] Pereira | Study of the Impact of High Temperatures and Pressures on the Equilibrium Densities and Interfacial Tension of the Carbon Dioxide/water System[END_REF] and references within for more details. In this work, the SAFT-VR Mie EoS was used to estimate the bulk equilibrium of the investigated systems and DGT used to predict interfacial properties. The main equations within the DGT framework are given below.

In summary, by applying the minimization criterion of the Helmholtz energy to planar interfaces, the interfacial tension values with respects to the density of a reference component is given by [START_REF] Miqueu | Modelling of the Surface Tension of Binary and Ternary Mixtures with the Gradient Theory of Fluid Interfaces[END_REF][START_REF] Miqueu | Modeling of the Surface Tension of Multicomponent Mixtures with the Gradient Theory of Fluid Interfaces[END_REF]:

      L ref V ref i j ref ref j ref i ij d d d d d c IFT         2 (Eq. 25)
where ρ L ref and ρ V ref are the bulk phase densities and the ref subscript denotes the reference component of the mixture. ∆Ω is the variation of the grand thermodynamic potential which is related to the Helmholtz free energy by the following equation:

  p f i i i         0 (Eq. 26)
where f 0 is the Helmholtz free energy density of the homogeneous fluid at local density, µ i are the chemical potential of each component and p is the pressure at equilibrium. The methodology described by [START_REF] Miqueu | Modelling of the Surface Tension of Binary and Ternary Mixtures with the Gradient Theory of Fluid Interfaces[END_REF] was here followed for determining the density distribution of each component across the interface.

In Eq. 25, c ij is the cross influence parameter and the mixing rule used based on the geometric mean of the pure component influence parameters [START_REF] Carey | The Gradient Theory of Fluid Interfaces[END_REF]) and it is given by

  j i ij ij c c c    1 (Eq. 27)
where β ij is the binary interaction parameter and c i and c j are the pure component influence parameters. In this work, the binary parameter β ij has been fixed to 0, making the calculation of the interfacial tension of the systems investigated fully predictive. The influence parameters c i and c j can be derived from theoretical expressions [START_REF] Freeman | Molecular Theory of Surface Tension[END_REF]. Nevertheless, the parameters are generally correlated using Surface Tension (ST) data from pure substances, by rewriting Eq. 25 for c i as follows:

2 0 2 1                   d p f ST c L V i (Eq. 28)
There are different approaches for calculating the influence parameters of each pure component [START_REF] Miqueu | Modelling of the Surface Tension of Binary and Ternary Mixtures with the Gradient Theory of Fluid Interfaces[END_REF]. In this work, c i has been taken as a constant value calculated from ST data far from the critical point [START_REF] Lafitte | Interfacial Properties of water/CO2: A Comprehensive Description through a Gradient Theory-SAFT-VR Mie Approach[END_REF]. Hence, the influence parameters used were adjusted against ST data (NIST webbook 2013) of the pure components at a reduced temperature t r =0.7. The parameters used are reported in Table 4.

Results

Viscosity

The viscosity of CO 2 systems was calculated using the TRAPP model and the densities computed from the SAFT-VR Mie EoS. A literature review of the available experimental data on the viscosity of some mixtures of relevance for CCS operations was carried and the results listed in Table 5. As can be seen in Table 5, experimental data on this property are available in a broad range of temperatures but, in most studies, pressure is limited to either atmospheric pressure or to a maximum pressure of 2.6MPa. Only in the study of [START_REF] Chapoy | Effect of Impurities on Thermophysical Properties and Phase Behaviour of a CO2-Rich System in CCS[END_REF] this property has been investigated for temperatures in the range (273-423) K and pressures up to 150MPa. The viscosity studies were conducted in the single phase region, i.e. Gas (G) phase at low pressures, Liquid (L) above saturation pressure or in the supercritical (SC) region, as described in Table 5. The modeling results for the viscosity of the CO 2 -rich systems are shown in Table 6, along with the percentage absolute average deviation (%AAD) between predicted and experimental viscosity values for each studied system. A total of 264 experimental points were considered and an overall %AAD of 2.68 % was calculated for the viscosity values predicted with the TRAPP+SAFT-VR Mie model. This value represents half of the deviation reported by [START_REF] Huber | The Corresponding-States Principle: Dense Fluids[END_REF] and [START_REF] Poling | The Properties of Gases and Liquids[END_REF] and obtained with the TRAPP model combined with a modified BWR-EoS. Nonetheless, the largest deviations with the TRAPP+SAFT-VR Mie model were obtained at high pressures and low temperatures, as depicted in Fig. 2, with a maximum deviation of 9.3% for the multicomponent system (90%CO 2 , 5%O 2 , 2%Ar and 3%N 2 ).

Interfacial tension

The saturated densities and IFT of CO 2 -rich systems were modeled by coupling the DGT with the SAFT-VR Mie EoS. The results for the binaries CO 2 -n-C 4 and CO 2 -n-C 10 and the ternary CO 2 -n-C 4 -n-C 10 are plotted in Fig. 3 and Fig. 4, respectively. As can be seen in Fig. 3 and Fig. 4, the DGT+SAFT-VR Mie predictions are in excellent agreement with the measured saturated density and IFT data of the systems considered. Largest deviations were observed close to the critical point, where the SAFT-VR Mie fails to capture the impact of pressure on the density of the bulk phases. Nonetheless, as can be seen in Figure 5, the maximum absolute deviation between predicted and measured IFT data was of 0.26mNm -1 . In general, the IFT values for the CO 2 -n-C 4 system and the ternary system are slightly overestimated, whereas an inverse behavior was observed for the CO 2 -n-C 10 system. In addition, the density profiles through the interface were computed with the DGT approach for the ternary system (90%CO 2 , 6%n-butane and 4%n-decane) and the results plotted in Fig. 5. As depicted in Fig. 5(a), the density profiles show a local enrichment of the interface with CO 2 molecules as a peak was noticeable in the density profile of carbon dioxide. The location of peak closer to the vapor side of the interface suggests that the absorption of CO 2 occurred at the surface of the liquid hydrocarbon phase. Furthermore, as can be seen in Fig. 5(b), the pressure increase lead to an increase of the interface thickness in this mixture and to a reduction of the local accumulation on light components, such as CO 2 .

Comparable results were obtained for the binary mixtures and the results are in agreement with those obtained by [START_REF] Miqueu | Modeling of the Surface Tension of Multicomponent Mixtures with the Gradient Theory of Fluid Interfaces[END_REF].

Conclusions

The TRAPP model and the density gradient theory have been combined with the SAFT-VR Mie EoS to model the viscosity and IFT of CO 2 -rich systems. The low deviations observed to measured density data showed the SAFT-VR Mie EoS good capabilities to describe the equilibrium of the studied systems with k ij presented in this work. Furthermore, by using the computed bulk phase properties, the TRAPP and DGT models were able to predict the viscosity and IFT of CO 2 -rich mixture with relative low deviations to experimental data, endorsing the capability of the methods for estimating key transport properties of fluids in CCS streams.

  Fig. 1: Predicted VLE by SAFT-VR Mie for CO2-nC4 (grey) (Nagarajan 1985) and CO2-nC10 (black) (Nagarajan 1986). Symbols: () T=344.3 and (▲) T=377.6K.

Fig. 3 :Figure 4 :Fig. 5 :

 345 Fig. 3: Predicted density (a) and IFT (b) of CO2-nC4 (black) (Nagarajan 1985) [ref] and CO2-nC10 (grey) (Nagarajan 1986). Symbols: () T=344.3 and (▲) T=377.6K.

  Fig. 5: (a) Density profiles calculated for each component in the CO2-nC4-nC10 mixture at 344.3K and 9.31MPa: CO2-(blue), n-butane (green) and n-decane (red). Total interface length of 7.36 nm. (b) Density profiles of carbon dioxide calculated for the CO2-nC10 mixture at 344.3K at different pressures against the results of molecular simulation at same conditions (Müller and Mejía 2009).

Table 1 : Molecular parameters of the compounds used in this work (Lafitte et al. 2013)

 1

	Substance	ms	σ/Å	(ε/k)/K	λr	λa
	Carbon dioxide	1.5000	3.1916	231.88	27.557	5.1646
	n-butane	1.8514	4.0887	273.64	13.650	6
	n-decane	2.9976	4.5890	400.79	18.885	6
	Nitrogen	1.4214	3.1760	72.438	9.8749	6
	Oxygen	1.4283	2.9671	81.476	8.9218	6
	Argon	1.0000	3.4038	117.84	12.085	6

Table 2 : Binary interaction parameters between O 2 , Ar, N 2 and CO 2 , and between CO 2 , n-butane and n- decane for the SAFT-VR Mie used in this work.

 2

		N2	O2	Ar	CO2	n-C4	n-C10
	N2	0	-0.009	-0.017	-0.127	-	-
	O2	-0.009	0	-0.023	-0.033	-	-
	Ar	-0.017	-0.023	0	-0.014	-	-
	CO2	-0.127	-0.033	-0.014	0	0.053	0.057
	n-C4	-	-	-	0.053	0	-0.004
	n-C10	-	-	-	0.057	-0.004	0

Table 3 : E i parameter TRAPP Model (Poling et al. 2001).

 3

	E1	-14.113294896
	E2	968.22940153
	E3	13.686545032
	E4	-12511.628378
	E5	0.01168910864
	E6	43.527109444
	E7	7659.4543472

Table 4 : Density gradient theory influence parameters used in Eq. 26.

 4

	Substance	ci / J•m 5 •mol -2

Table 5 : Literature experimental data for the viscosity of CO 2 -mixtures with N 2 , O 2 , and Ar.

 5

	Source	System	Phase	T/K	P/MPa	N	Uncertainty
	Kestin and Leidenfrost (1959)	CO2/N2, CO2/N2	G	293	0.1-2.17	28	±0.05%
	Kestin et al. (1966)	CO2/N2, CO2/Ar	G	293-304	0.1-2.6	83	±0.1%
	Gururaja et al. (1967)	CO2/N2, CO2/O2	G	298	0.1	17
	Kestin and Ro (1974)	CO2/N2, CO2/Ar, CO2/N2/Ar	G	297-773	0.1	64	±0.15%
	Kestin et al. (1977)	CO2/O2	G	298-674	0.1	10	±0.3%
	Hobley et al. (1989)	CO2/Ar	G	301-521	0.1	24	<0.7%
	Chapoy et al. (2013)	CO2/N2/O2/Ar	G/L/ SC	273-423	1.5-150	38
					Total	264

Table 6 : Absolute average deviation of the viscosity calculations.

 6

	System	xCO2	N	T/K	P/MPa	%AAD
	CO2/N2	0.19-0.90	95	293-773	0.1-2.6	3,12
	CO2/Ar	0.22-0.92	96	293-773	0.1-2.6	1,28
	CO2/O2	0.19-0.92	19	298-674	0.1	3,25

Experimental (Chapoy et al. 2013) and predicted viscosity of a CO2-rich system (90%CO2, 5%O2, 2%Ar and 3%N2).


Marc Lavielle
email: marc.lavielle@inria.fr

Leon Aarons

What do we mean by identifiability in mixed effects models?

Keywords: Model identifiability, Practical identifiability, Structural identifiability, Parameter estimation, Mixed effects model, Pharmacokinetics

We discuss the question of model identifiability within the context of nonlinear mixed effects models. Although there has been extensive research in the area of fixed effects models, much less attention has been paid to random effects models. In this context we distinguish between theoretical identifiability, in which different parameter values lead to non-identical probability distributions, structural identifiability which concerns the algebraic properties of the structural model, and practical identifiability, whereby the model may be theoretically identifiable but the design of the experiment may make parameter estimation difficult and imprecise. We explore a number of pharmacokinetic models which are known to be non-identifiable at an individual level but can become identifiable at the population level if a number of specific assumptions on the probabilistic model hold. Essentially if the probabilistic models are different, even though the structural models are non-identifiable, then they will lead to different likelihoods. The findings are supported through simulations.

Introduction

A statistical model is said to be identifiable when, given an infinite amount of data, it is possible to uniquely estimate the true values of the model parameters. The uniqueness property implies that different values of the model parameters generate different probability distributions of the observable variables. Conversely if two or more sets of parameters generate identical distributions of the observed values the model is not identifiable. However it still may be possible to uniquely identify a subset of parameters: that is, not all the parameters of an unidentifiable model are unidentifiable. In this case we say that the model is partially identifiable.

It should also be noted that even when the parameters of the model are unidentifiable, the model itself may perfectly describe the observed variables, in other words perfectly fit the data.

A distinction is made in the literature between structural identifiability and practical identifiability. Structural identifiability is related to the structure of the underlying mathematical model. It is concerned with whether the parameters of a model can be exactly identified from a given experiment with perfect input-output data [START_REF] Bellman | On structural identifiability[END_REF][START_REF] Cobelli | Parameter and structural identifiability concepts and ambiguities: a critical review and analysis[END_REF][START_REF] Shivva | An approach for identifiability of population pharmacokineticpharmacodynamic models[END_REF][START_REF] Walter | On the identifiability and distinguishability of nonlinear parametric models[END_REF][START_REF] Xia | Identifiability of nonlinear systems with application to hiv/aids models[END_REF]. Chappell et al. [START_REF] Chappell | Global identifiability of the parameters of nonlinear systems with specified inputs: a comparison of methods[END_REF] compare methods for analyzing the global structural identifiability of the parameters of a nonlinear system with a specified input function. Gargash and Mital [START_REF] Gargash | A necessary and sufficient condition of global structural identifiability of compartmental models[END_REF] study the problem of global and local structural identifiability for fixed structure compartmental deterministic model of the biological system. Frohlich et al. [START_REF] Fröhlich | Uncertainty analysis for non-identifiable dynamical systems: Profile likelihoods, bootstrapping and more[END_REF] investigate the effect of structural non-identifiability on the performance of frequentist methods for standard uncertainty analysis. They observe that the profile likelihood approach is the only one that properly identifies structural non-identifiability of parameters. On the other hand practical identifiability is related to the study design, i.e. the limited amount of information that can be obtained from a given experiment [START_REF] Guedj | Practical identifiability of hiv dynamics models[END_REF][START_REF] Petersen | Practical identifiability of model parameters by combined respirometric-titrimetric measurements[END_REF][START_REF] Shivva | An approach for identifiability of population pharmacokineticpharmacodynamic models[END_REF]. A link exists between practical identifiability and sensitivity analysis. Brun et al. [START_REF] Brun | Practical identifiability analysis of large environmental simulation models[END_REF] propose a systematic approach for tackling the parameter identifiability problem of large models based on local sensitivity analysis.

Model identifiability is closely related to model indistinguishability. The objective of indistinguishability analysis is to determine if different models are capable to fit the available input-output data [START_REF] Godfrey | Identifiability and indistinguishability of nonlinear pharmacokinetic models[END_REF][START_REF] Yates | Structural identifiability and indistinguishability of compartmental models[END_REF]. Identifiability and distinguishability of parametric models are important properties when the parameters to be estimated have a biological meaning or when the model is to be used to reconstruct physiologically meaningful variables that cannot be measured directly [START_REF] Walter | On the identifiability and distinguishability of nonlinear parametric models[END_REF].

The question of identifiability of pharmacokinetic (PK) and pharmacokinetic-pharmacodynamics (PKPD) models has been previously studied. Evans et al. [START_REF] Evans | An identifiability analysis of a parent-metabolite pharmacokinetic model for ivabradine[END_REF] consider the identifiability of a parentmetabolite pharmacokinetic model for ivabradine and one of its metabolites. Shivva et al. [START_REF] Shivva | Parameterisation affects identifiability of population models[END_REF] use a simple one compartment population pharmacokinetic model to show that identifiability of the variances of the random effects are affected by the parameterisation of the fixed effects. An approach for assessment of identifiability for fixed and mixed effects PK or PKPD models is proposed in [START_REF] Shivva | An approach for identifiability of population pharmacokineticpharmacodynamic models[END_REF]. Guedj et al. [START_REF] Guedj | Practical identifiability of hiv dynamics models[END_REF] study the identifiability of parameters in a model of HIV dynamics based on a system of non-linear Ordinary Differential Equations. Garcia et al. [START_REF] Garcia | Identifiability of PBPK models with applications to dimethylarsinic acid exposure[END_REF] discuss different types of identifiability that occur in physiologically-based pharmacokinetic (PBPK) models and give reasons why they occur.

Identifiability in mixed effects models has received much less attention. Identifiability of linear mixed effects models is considered in [START_REF] Wang | Identifiability of linear mixed effects models[END_REF]. In the context of mixed effects models, identifiability is a fundamental prerequisite for model identification. It concerns uniqueness of the population parameters estimated from a given set of observations obtained from several individuals of a same population [START_REF] Bonate | Pharmacokinetic-pharmacodynamic modeling and simulation[END_REF][START_REF] Lavielle | Mixed Effects Models for the Population Approach: Models, Tasks, Methods and Tools[END_REF][START_REF] Wu | Mixed effects models for complex data[END_REF]. It was the purpose of the current investigation to examine some issues in identifiability of mixed effects models.

We define the concept of identifiability of mixed effects models in Section 2. Identifiability is usually defined for continuous data models but this concept can be easily extended to categorical, count or time to event data models. We show that a model can be identifiable even if it is not structurally identifiable. Indeed the probability distribution of the individual parameters plays an important role for characterizing the identifiability property of the model. Furthermore, a distinction is made between identifiability of the population parameters (the population parameters can be estimated successfully and unequivocally from the observed data) and identifiability of the individual parameters (the individual parameters can be estimated successfully and unequivocally from the observed data).

Sections 3 and 4 present some specific results taken from the field of pharmacokinetics. We consider first a situation where different parameterizations of a PK model are algebraically undistinguishable. The model is shown to be identifiable as soon as some hypothesis are made on the correlation structure of the PK parameters. It is usually claimed that the bioavailability F cannot be estimated using only PK measurements obtained from an oral administration: only the ratio F/V can be properly estimated. We show that the model is identifiable under some assumptions on the probabilistic model. Then, both F and V can be simultaneously estimated.

Identifiability of mixed effects models

Preliminary remarks

The definition of identifiability for mixed effects models is not always very precise. As an example, the definition given in [START_REF] Wu | Mixed effects models for complex data[END_REF] reduces to: "parameters or models are called non-identifiable if two sets of different parameters lead to the same probability distribution". This definitions remain quite ambiguous and needs to be clarified. In particular, we will need to distinguish the identifiability of the population parameters and the identifiability of the individual parameters.

We need also to make a clear distinction between identifiability, structural identifiability and practical identifiability, defining what these properties mean in the context of mixed effects models.

The concept of structural identifiability was introduced first in the area of systems and control (see [START_REF] Bellman | On structural identifiability[END_REF]), where the systems are deterministic and depend on nonrandom parameters. The situation is quite different in a population approach context where the individual parameters are random variables. In order to analyse the properties of a statistical model it is necessary to take into account the algebraic properties of the structural model as well as the properties of the probabilistic model. In particular, the choice of the parameterization may have a strong impact on the probabilistic properties of the model. For instance, parameterizations of a PK model using volume V and clearance Cl, or using volume V and elimination rate constant k, where k = Cl/V , are absolutely equivalent from a purely algebraic point of view: we can use both interchangeably. That's not the case if we put a distribution on these parameters.

Assume, for instance, that V and Cl are independent and log-normally distributed:

log(V ) ∼ N (log(V pop ), ω 2 V ) log(Cl) ∼ N (log(Cl pop ), ω 2
Cl ) corr (log(V ), log(Cl)) = 0 where corr (log(V ), log(Cl)) is the linear correlation between log(V ) and log(Cl).

Then, we are implicitly assuming that k is also log-normally distributed with a variance ω 2 k larger than the variance of V and that log(k) and log(V ) are negatively correlated. Indeed,

log(k) = log(Cl) -log(V )
Then, since log(Cl) and log(V ) are both normally distributed, log(k) is also normally distributed and log(k) ∼ N (log(Cl pop /V pop ), ω 2 V + ω 2 Cl ) Furthermore, the covariance between log(V ) and log(k) is given by

cov (log(V ), log(k)) = -var (log(V )) = -ω 2 V Then, corr (log(V ), log(k)) = -ω V /ω k
On the other hand, if we use parameters (V, k) in the model, assuming independent distributions, then, we implicitly assume that V and Cl are dependent.

Structural identifiability

Structural identifiability is related to the structure of the underlying mathematical model, for example as discussed above, the PK model.

For sake of simplicity, we will only consider models for univariate data, i.e. for a single outcome.

We need some notation. Let N be the number of individuals. Then, for i = 1, 2, . . . , N , let y i = (y ij , 1 ≤ j ≤ n i ) be the n i observations for individual i collected at times (t ij , 1 ≤ j ≤ n i ).

Let us start with a basic model for continuous data:

y ij = f (t ij ; ψ i ) + ε ij
Here, f is the structural model, which is a parametric function of time, ψ i is a vector of individual parameters, and

(ε ij , 1 ≤ j ≤ n i
) is a sequence of residual errors. We will assume that ε ij is a sequence of random variables with mean 0 and finite variance σ 2 .

Structural identifiability of the model is directly related to the properties of the structural model f . We don't take into account possible differences between individuals.

Without any loss of generality, we will assume that f is defined for t ≥ 0. Let f ( • ; ψ) be the function f defined for a given vector of parameters ψ (we could equivalently use the notation f ψ ).

Then, we will say that the model is structurally identifiable if there exists a one-to-one mapping between the parameter ψ and the function f ( • ; ψ), i.e.

ψ = ψ ⇔ f (t; ψ) = f (t; ψ ) for any t ≥ 0
This definition can be easily extended to any other parametric mixed effects models:

-Continuous data model with non constant residual error model: Assume now that there exists a parametric function g such that

y ij = f (t ij ; ψ i ) + g(t ij ; ψ i )ε ij
Here, the structural model is the pair (f, g).

-Time-to-event data model: The structural model is the hazard function h(t; ψ).

-Count data model: Consider for instance a Poisson model,

y ij ∼ Poisson(λ(t ij , ψ i ))
The structural model is the Poisson intensity λ(t; ψ)

-Categorical data model: Assume a Bernoulli model for binary data as an example,

P(y ij = 1) = 1 -P(y ij = 0) = π(t ij ; ψ i )
Here, the structural model is the function π(t; ψ).

Practical identifiability

Practical non-identifiability is less clearly defined in the literature compared to structural nonidentifiability.

Practical identifiability not only depends on the model structure, but is also related to the experimental conditions together with the quality and quantity of the measurements [START_REF] Petersen | Practical identifiability of model parameters by combined respirometric-titrimetric measurements[END_REF]. Then, a parameter that is structurally identifiable may be practically unidentifiable with a limited amount and quality of experimental data.

Deriving some relationship between identifiability and confidence interval of parameter estimates seems natural. Nevertheless, since practical identifiability is a non-asymptotic property, it is not possible to propose a rigorous definition based on asymptotic confidence intervals.

Raue et al. [START_REF] Raue | Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood[END_REF] propose an appealing definition of practical identifiability based on a likelihoodbased confidence region instead of asymptotic confidence intervals. They suggest an approach that exploits the profile likelihood and enables the detection of both structural and practical nonidentifiabilities.

Unfortunately, such an approach is cumbersome to adopt for (nonlinear) mixed effects models since it requires the estimation of the population parameters and computation of the likelihood many times.

Methods that can be used in practice for detecting some identifiability issues remain quite empirical:

-We can for instance run the estimation algorithm with different initial values. Convergence to different solutions may be suspicious.

-Even if it is an asymptotic criteria, the observed Fisher Information Matrix can also be used. Indeed, the inverse of this matrix provides an approximation of the variance-covariance matrix of the estimated parameters. A large condition number of this variance-covariance matrix (i.e. the ratio of its largest and smallest eigenvalues) reflects a strong correlation structure between estimates and may indicate some identifiability issue.

Identifiability of the population parameters

Identifiability of the population parameters θ is related to the properties of the statistical model of the observations y = (y i , 1 ≤ i ≤ N ):

p(y; θ) = N i=1 p(y i ; θ) (1)
= N i=1 p(y i , ψ i ; θ) dψ i (2)
We are here in a "classical" situation where the statistical model is identifiable if the mapping between θ and the probability distribution p(y; θ) is one-to-one (we use indiscriminately p(y; θ) for the probability distribution function (pdf) and for the distribution of y).

By definition, the Maximum Likelihood (ML) estimate of θ maximizes p(y; θ). If the model is identifiable, then, the ML estimate converges to the true value of θ when N increases, under very general regularity conditions [START_REF] Lavielle | Mixed Effects Models for the Population Approach: Models, Tasks, Methods and Tools[END_REF]. That means that it is possible to estimate θ as precisely as required as soon as the size N is large enough. On the other hand, if the model is not identifiable, then, maximizing p(y; θ) with respect to θ will not lead to a unique solution. We can find, for instance, θ(1) and θ(2) such that, for any vector of population parameters θ, p(y; θ) ≤ p(y; θ(1) ) = p(y; θ(2) )

A number of particular cases have been reported. For instance, Wang et al. [START_REF] Wang | Identifiability of linear mixed effects models[END_REF] study the identifiability of the covariance parameters in a linear mixed effects model. They focus on those models that are not over-parameterized and derive conditions of identifiability and study commonly used covariance structures. In an unpublished work, Nuñez and Concordet consider the identifiability problem in a nonlinear mixed effects model for continuous data, assuming a nonparametric distribution for the individual parameters. They provide several explicit conditions on f which ensure the identifiability of the model.

Since p(y; θ) = p(y, ψ; θ), identifiability of the complete model p(y, ψ; θ) is a necessary condition to ensure the identifiability of the observed model p(y; θ). Unfortunately, it is not a sufficient condition.

For example consider the linear model

y ij = (a i + b i )t ij + ε ij (3)
where a i and b i are normally distributed with unknown means m 1 and m 2 . Here, the vector of population parameters θ includes m 1 and m 2 . The model p(y; θ) is unidentifiable since the population parameters m 1 and m 2 cannot be estimated from a sequence of observations (y ij ): only the sum m 1 + m 2 can be estimated since a i + b i is normally distributed with mean m 1 + m 2 . On the other hand, the joint model p(y i , a i , b i ; θ) is identifiable since m 1 and m 2 can be estimated using sequences

(a i , 1 ≤ i ≤ N ) and (b i , 1 ≤ i ≤ N ).
Analyzing the structural identifiability of the model is important and useful, but it is not sufficient for concluding if the model is identifiable or not. For instance, a statistical model may be identifiable even if it is not structurally identifiable. Indeed, when data coming from various individuals is available, and under some hypothesis, we can take advantage of the probability distribution of the individual parameters to estimate the population parameters. It is then the combination of algebraic relationships and probabilistic relationships that make the model identifiable.

Consider again the linear model [START_REF] Brun | Practical identifiability analysis of large environmental simulation models[END_REF]. The structural model f (t; a i , b i ) = (a i + b i )t is not identifiable since only a i +b i can be estimated. We have seen in the previous example that the statistical model of the observations is not identifiable when a i and b i are both normally distributed. Surprisingly, the model becomes identifiable in most cases when a i and b i are not both normally distributed.

Assume for instance that a i and b i have exponential distributions with parameters λ 1 and λ 2 . Then, assuming that λ 1 = λ 2 , the sum a i + b i is a random variable with pdf

p(z; λ 1 , λ 2 ) = λ 1 λ 2 λ 1 -λ 2 e -λ2z -e -λ1z
There is a one-to-one mapping between (λ 1 , λ 2 ) and this pdf: the model is identifiable. Identifiabililty of the model remains a theoretical property: in practice, accurate estimation of λ 1 and λ 2 will require a huge amount of data.

Identifiability of the individual parameters

We can see the problem of estimating the individual parameters (ψ i ) as an inverse problem: we aim to recover these unobserved vectors of parameters using the observations (y i ).

Structural unidentifiability means that the problem is ill posed: we can find different vectors ψ i and ψ i that produce the same structural predictions:

ψ i = ψ i and f (t; ψ i ) = f (t; ψ i ) for any t ≥ 0
And then,

p(y i |ψ i ) = p(y i |ψ i )
Each individual conditional model p(y i |ψ i ) is therefore unidentifiable, but in a mixed effects model, we usually don't estimate each individual parameter ψ i by maximizing this conditional distribution p(y i |ψ i ). Indeed, each individual parameter ψ i is a random vector with distribution p(ψ i ; θ). Then, for a given vector of population parameter θ (given or previously estimated), we rather consider the posterior distribution

p(ψ i |y i ; θ) = p(y i |ψ i )p(ψ i ; θ) p(y i ; θ)
Thus, ψ i is identifiable as soon as p(ψ i ; θ) = p(ψ i ; θ).

In other words, the population distribution of the individual parameters plays the role of a prior distribution that now makes the problem well posed.

Example 1: Consider a bi-exponential model:

f (t; A i , B i , α i , β i ) = A i e -αi t + B i e -βi t
This structural model is locally identifiable but globally unidentifiable since (A i , α i ) and (B i , β i ) are interchangeable. The model becomes identifiable if we introduce some information about α i and β i , assuming for instance that P(α i > β i ) > 0.5. The use of this prior leads to select the solution where α i > β i and discard the other solution where α i < β i . This structural model f (t; a i , b i ) = (a i + b i )t is clearly locally unidentifiable since only the sum a i + b i can be estimated maximizing the conditional pdf of y i . Let c i = a i + b i and define ĉi as

ĉi = Argmax ci p(y i |c i ) = Argmin ci j (y ij -c i t ij ) 2
Thus, âi + bi = ĉi = j y ij t ij / j t 2 ij , but it is impossible to decompose this sum and compute âi and bi without any additional information.

Assume for instance that

a i b i ∼ N m 1 m 2 , Ω
Let T i be the n i × 2 matrix

T i = t i,1 t i,2 . . . t i,ni t i,1 t i,2 . . . t i,ni
The posterior distribution of a i and b i is now a well defined normal distribution:

a i b i |y i ∼ N µ i,1 µ i,2 , Γ i
where

Γ i = T i T i σ 2 + Ω -1 -1 and µ i,1 µ i,2 = Γ i T i y i σ 2 + Ω -1 m 1 m 2
The maximum a posteriori (MAP) estimates of a i and b i are, respectively, µ i,1 and µ i,2 . They are well defined and unique: the individual parameters of the model are now identifiable.

Remark: introduction of a Gaussian prior information for estimating the parameters of a linear ill-posed problem is equivalent to introduce a Tikhonov regularization term [START_REF] Tikhonov | Numerical methods for the solution of ill-posed problems[END_REF].

3 Illustration: the flip-flop phenomenon

In this example, the structural model is not identifiable without some physiological constraint (we can find two different sets of PK parameters that produce identical PK profiles). Nevertheless, identifiability is ensured under some assumptions on the probabilistic model.

For a sake of simplicity, we will consider a single individual and omit the subscript i in the notation.

The structural PK model

Consider a basic PK model for a single oral administration at time 0,

f (t; k a , V, k) = Dk a V (k a -k) e -k t -e -ka t . (4)
It is easy to see that, for any t ≥ 0,

f (t; k a , V, k) = f (t; k a , V , k )
where

k a = k , k = k a and V = (k/k a )V .
Without any assumptions on the parameter values, the two solutions are indistinguishable. Some assumptions on the parameter space may make the model identifiable (e.g. k a > k).

The probabilistic model

Without any assumption on the parameter space the structural model is not identifiable but the probabilistic model may be identifiable under some assumptions.

If we assume log-normal distributions for (k a , V, k), then (k a , V , k ) are also log-normally distributed since log(V ) = log(V ) + log(k) -log(k a ).

Furthermore, if we assume that log(k a ), log(V ) and log(k) are uncorrelated, then log(V ) and log(k a ) are implicitly correlated, as well as log(V ) and log(k ).

Thus, the two following models are distinguishable:

M 1 : corr (k a , V ) = corr (k a , k) = corr (k, V ) = 0 M 2 : corr (k a , V ) = corr (k a , k ) = corr (k , V ) = 0
In other words, if we assume that the variance-covariance matrix Ω of (log(k a ), log(V ), log(k)) is diagonal, then the model is identifiable since we cannot have simultaneously both Ω and Ω diagonal, where Ω is the variance-covariance matrix of (log(k a ), log(V ), log(k )).

On the other hand, if we don't make any assumption on Ω, then the model is not identifiable: we can find two sets of fixed parameters (ka pop , V pop , k pop ) and (ka pop , V pop , k pop ) and two covariance matrices Ω and Ω such that p(y, ψ; ka pop , V pop , k pop , Ω) = p(y, ψ; ka pop , V pop , k pop , Ω )

where y = (y j , 1 ≤ j ≤ n) is a vector of observed concentrations and ψ = (k a , V, k) is the vector of individual PK parameters.

Statistical implications

Assume now that we have some PK data and we want to use the PK model defined in (4) to fit this data. Then, it is expected to obtain different values of the likelihood under M 1 and M 2 since the two probabilistic models are different (we cannot have simultaneously corr (k, V ) = 0 and corr (k , V ) = 0).

Simulation:

Data with N = 100 subjects and n = 12 measurements per subject were simulated with the following parameter values:

k a,pop = 1, V pop = 10, k pop = 0.1, ω ka = 0.25, ω V = 0.3, ω k = 0.15.
We can then try to estimate the population parameters of the model, using any software tool such as Monolix or NONMEM. According to the initial value, the SAEM algorithm [START_REF] Lavielle | Mixed Effects Models for the Population Approach: Models, Tasks, Methods and Tools[END_REF] in Monolix converges to two solutions: ( k(1) a , V (1) , k(1) ) = (1, 10.4, 0.099) and ( k(2) a , V (2) , k(2) ) = (0.097, 0.99, 1.04). The estimated log-likelihood for these two solutions are, respectively, -1126 and -1145. We see that, in this example, the first solution, which is the "right" solution (i.e. close to the "true" values used for the simulation), correspond to the global maximum of the likelihood, while the other one correspond to a local maximum.

In other words, the model is identifiable if we assume a diagonal matrix Ω, and maximizing the likelihood allows to select the "right" solution.

Here, it is not the values of the population PK parameters that allows one to select a model, but the (strong) hypothesis that we make concerning the covariance structure of the random effects. Indeed, if we simulate uncorrelated PK parameters using the other parameterization, then, the likelihood criteria will select this solution.

On the other hand, if we don't make any assumption about the covariance matrix, i.e. if we estimate a full variance-covariance matrix, then the model is not identifiable and the likelihood criteria cannot select a model. Indeed, the log-likelihood still exhibits two maxima, (1, 10.4, 0.1) and (0.098, 1, 1.03), but with very close values: -1125.7 and -1126.2, respectively.

Remark 1: As expected, even if independent random effects were simulated using the first parameterization, estimated variance-covariance matrix associated to the second solution is not diagonal: estimated correlations between log(k a ) and log(V ) and between log(k ) and log(V ) are, respectively, 0.41 and -0.62.

Remark 2: Similar results are obtained with NONMEM and FOCE:

k(1) a = 1.01, V (1) = 10.4, k(1) = 0.099 using the first set of initial estimates and k(2) a = 0.098, V (2) = 0.995, k(2) = 1.04 using the second one. We have simulated rich data in this example. Then, up to the exchangeability issue between k a and k, the individual PK parameters can be estimated accurately from the individual PK data. FOCE works very well here because this algorithm is precisely based on the estimation of the individual parameters.

4 Illustration: Identifiability of the bioavailability We show with this example that, even if the structural model only allows the identification of the ratio V /F , the probabilistic model, under some assumptions, makes the model identifiable and allows the estimation of both V and F .

We will consider again a single individual and omit the subscript i in the notation.

The model and its properties

-The PK model is a one compartment model for oral administration with parameters (F, k a , V, k)

where F is the bioavailability, i.e. the fraction of administered dose which is absorbed. Let f (t; F, k a , V, k) be the predicted concentration given by this model at time t:

f (t; F, k a , V, k) = D F k a V (k a -k) e -k t -e -ka t (5)
-The residual error model is an exponential model, i.e. observed concentrations are log-normally distributed: log(y j ) = log(f (t j ; F, k a , V, k)) + ε j k a , V and k are log-normally distributed while F has a logit-normal distribution:

logit(F ) ∼ N (logit(F pop ), ω 2 F )
where logit(x) = log(x/(1 -x)) for 0 < x < 1.

First of all, it is easy to see that the structural model is not identifiable. Indeed, let (F, k a , V, k) and (F , k a , V , k ) be two set of individual parameters such that

k a = k a , k = k, V /F = V /F , then f (t; F, k a , V, k) = f (t; F , k a , V , k ) for any t > 0.
The structural model is therefore partially identifiable since only k a , k and R = V /F are identifiable.

Even if the structural model is not identifiable, the model itself is identifiable. We will use a "two step procedure" to demonstrate that we can derive a consistent estimator for all the population parameters of the model. Here, consistency means that this estimator converges to the true values of the population parameters when both the number of individuals and the number of observations per individual tend to infinity.

1. The structural model is partially identifiable. Then, for each individual i = 1, 2, . . . , N , the set of identifiable individual parameters k ai , k i and R i = V i /F i can be perfectly recovered when the number of measurements n i for individual i tends to infinity, by maximizing the conditional distribution p(y i |k ai , k i , R i ).

2.

The ratio R = V /F is a random variable defined as the ratio of a logit-normal and a lognormal variable. This distribution depends on parameters F pop , V pop , ω F and ω V and there exists a one-to-one mapping between these four parameters and the distribution of R. Thus, the maximum likelihood estimator of these four parameters, derived from a N -sample R 1 , R 2 , . . . R N of R, is consistent: it converges to the true values of these parameters when the number of individual N tends to infinity. ML estimators of k a,pop and ω ka (resp. k pop and ω k ) obtained from k a1 , . . . , k a N (resp. k 1 , . . . , k N ) are also consistent.

Remark 1: The model is not identifiable if both V and F are log-normally distributed. Indeed, the log-ratio log(R) is normally distributed:

log(R) ∼ N (log(V pop ) -log(F pop ), ω 2 V + ω 2 F )
and there exists an infinity of possible decompositions leading to the same probability distribution.

Remark 2: This interesting result remains an asymptotical result. In practice, it means that we can expect to estimate all the population parameters of the model with a desired precision, if we have enough data for that. When the number of individuals and measurements is finite, the properties of the ML estimator cannot be derived analytically. A Monte-Carlo study can be used to evaluate these properties for a given design.

Simulation study

We simulate PK data from this model for N = 5000 individuals. A single dose D = 100 is administrated at time 0 and n = 23 measurements are collected at times 0.5, 1, 3, 5, . . . , 21, 23.

The standard deviation of the residual errors

(ε ij ) is σ = 0.1.
Values of the population PK parameters are k a,pop = 1, V pop = 10, k pop = 0.1, ω ka = 0.25, ω V = 0.3 and ω k = 0.15. We will consider two logit distributions for F .

Model A: logit(F ) ∼ N (logit(0.9), 1)

The pdf of F is displayed in Figure 1. We see that this distribution is very different from a lognormal distribution. We can then expect to be able to estimate the population parameters.

We used the SAEM algorithm implemented in Monolix for computing the ML estimate of the population parameters and their standard errors. Table 1 shows that population parameters are indeed very well estimated in this example.

Even if the population parameters are "almost" perfectly estimated, the individual parameters cannot be estimated very precisely. Figure 2 compares the simulated individual parameters, considered here as the "true" values, with the Maximum a Posteriori (MAP) estimates, i.e. the modes of the conditional distributions p(V i |y i , θ) and p(F i |y i , θ) for i = 1, 2, . . . , N . On the other hand, the ratio F i /V i is estimated very accurately.

In this example, the likelihood has a maximum which is very well defined and SAEM converges easily even with poor initial guesses. Figure 3 displays the convergence of 5 runs of SAEM obtained with different initial values. Thus, thanks to the design (i.e. a very large number of subjects and a large number of measurement per subject) and thanks to the probability distribution of the individual parameters, the model can be considered as "practically identifiable".

Estimations obtained with NONMEM FOCE are F = 0.764, ka = 1.00, V = 8.88, k = 0.100. We can see that FOCE introduces some bias in the estimation of F and V . Indeed, it is not possible to estimate correctly the individual parameters since the model is not structurally identifiable. Then, any method based on the estimation of the individual parameters cannot work as well as maximum likelihood estimation.

Model B: logit(F ) ∼ N (logit(0.4), 0.2 2 )

Things will change with this distribution for F . Indeed, we can see Figure 4 that this distribution is now very close to a log-normal one. Even if the model remains identifiable in theory, we cannot expect anymore a good estimation of the population parameters.

Table 2 displays the results obtained with a single run of SAEM. We see that population parameters F pop and V pop are poorly estimated with this run. The ratio F pop /V pop remains very well estimated (0.401 instead of 0.4) as well as the total variance ω 2 F + ω 2 V (0.139 instead of 0.13). We also notice a clear degradation of the results obtained with FOCE ( Fpop = 0.830, Vpop = 20.8, ωF = 0.56, ωV = 0.3).

This poor estimation of the population parameters leads to a misspecified population distribution and a bias in the estimation of the individual parameters. Figure 5 shows that the V i and the F i are underestimated. Nevertheless, ratios R i = V i /F i are correctly estimated.

In this example and because of the lack of identifiability of some parameters, the likelihood does not exhibit a unique isolated maximum. Figure 6 shows that the convergence of SAEM strongly depends on the initial value. The very high correlation (0.9989) between the estimates of F pop and V pop confirms that we should not rely on the estimated values of these parameters.

When such lack of practical identifiability is revealed, a solution with this example consists in fixing either F pop or V pop . A less radical solution may consist in introducing a prior information on F pop or V pop . If we introduce, for instance, a logit-normal distribution for F pop , with mean logit(0.4) and standard deviation 0.2, then the population parameters are correctly estimated ( Fpop = 0.41, Vpop = 10.1) as well as the individual parameters.

We should notice that removing the inter-individual variability of F i is not a solution since that makes the model non identifiable. Indeed, if F i = F pop , then V i /F i follows a log-normal distribution with mean log(V pop /F pop ) (in the log domain) and standard deviation ω F . Then, only the ratio V pop /F pop is identifiable in this model.

Conclusions

We have shown that models that are non-identifiable at an individual level may become identifiable at the population level under conditions in which the probabilistic models differ between alternate models. This requires strong assumptions about the probabilistic models which may be difficult to validate in practice. Similarly, even if the models are identifiable at the individual level it may prove difficult to estimate the parameters of the model unless supported by good experimental design. From a pharmacokinetic point of view it means that the differences between individuals can break the non-identifiability seen at the population level and this may allow better mechanistic understanding of the interindividual differences in pharmacokinetics.

We have mainly considered here the most theoretical aspects of the identifiability of a model. For the numerical examples, we have been using an EM-like algorithm, assuming that the maximum likelihood estimate of the population parameters could be computed. Our first partial results suggest that linearization methods (FO, FOCE) are more sensitive to a lack of identifiability than maximum likelihood methods with no approximation on the model. A detailed discussion around the impact of of the estimation method -and its implementation in a software tool -on the results is beyond the scope of this paper. Such discussion as well as practical suggestions on how to proceed when the model shows signs of un-identifiability clearly deserve to be the subject of further work.

Example 2 :

 2 Consider again model (3) where ε ij ∼ i.i.d. (0, σ 2 ).

Fig. 1

 1 Fig.1Model A: pdf of the logit-normal distribution with parameters (logit(0.9), 1).

Fig. 2

 2 Fig. 2 Model A: estimated values versus true values of log(V i ), logit(F i ) and R i = V i /F i .

Fig. 3

 3 Fig. 3 Model A: convergence of 5 runs of SAEM obtained with different initial values.

Fig. 4

 4 Fig. 4 Model B: pdf of the logit-normal distribution with parameters (logit(0.4), 0.2 2 ).

Fig. 5 Fig. 6

 56 Fig. 5 Model B: estimated values versus true values of log(V i ), logit(F i ) and R i = V i /F i .

Table 1

 1 Model A: ML estimates of the population parameters.

	parameter	true value	initial value	estimation relative s.e. (%)
	Fpop	0.9	0.4	0.909	0.8
	ka,pop	1	0.5	1.006	0.4
	Vpop	10	5	10.188	0.9
	kpop	0.1	0.3	0.100	0.2
	ω F	1	3	1.007	7.0
	ω ka	0.25	1	0.249	1.3
	ω V	0.30	1	0.306	1.4
	ω k	0.15	1	0.148	1.1
	σ	0.10	1	0.100	0.3
	parameter	true value	initial value	estimation relative s.e. (%)
	Fpop	0.4	0.2	0.245	9.9
	ka,pop	1	2	0.999	0.40
	Vpop	10	5	6.11	9.8
	kpop	0.1	0.05	0.10	0.2
	ω F	0.2	1	0.289	24.0
	ω ka	0.25	1	0.249	1.3
	ω V	0.30	1	0.236	20.5
	ω k	0.15	1	0.148	1.1
	σ	0.10	1	0.01	0.3

Table 2

 2 Model B: ML estimates of the population parameters.

Acknowledgments

The research leading to these results received support from the Innovative Medicines Initiative Joint Undertaking under grant agreement 115156, resources of which are composed of financial contributions from the European Union's Seventh Framework Programme (FP7/2007-2013) and EFPIA companies' in kind contribution. The DDMoRe project is also supported by financial contribution from Academic and SME partners.