Adeline Lach
email: a.lach@brgm.fr

Faïza Boulahya
email: f.boulahya@brgm.fr

Laurent André
email: l.andre@brgm.fr

Arnault Lassin
email: a.lassin@brgm.fr

Mohamed Azaroual
email: m.azaroual@brgm.fr

Jean-Paul Serin
email: jean-paul.serin@univ-pau.fr

Pierre Cézac
email: pierre.cezac@univ-pau.fr

Thermal and volumetric properties of complex aqueous electrolyte solutions using the Pitzer formalism -The PhreeSCALE code

Keywords: PHREEQC, PhreeSCALE, Pitzer equation, heat capacity, density, Na 2 SO 4, MgSO 4

est destinée au dépôt

Thermal and volumetric properties of complex aqueous electrolyte solutions using the Pitzer formalism -The

PhreeSCALE code

Introduction

Knowledge of the heat capacity and the density of electrolyte solutions is of interest in many fields. These parameters are necessary for geothermal process studies [START_REF] Schröder | Design and test of a new flow calorimeter for online detection of geothermal water heat capacity[END_REF] and in order to evaluate desalination process performance [START_REF] Lin | Direct contact membrane distillation with heat recovery: Thermodynamic insights from module scale modeling[END_REF]. Density data are also very important in CO 2 sequestration studies [START_REF] Bachu | Sequestration of CO 2 in geological media in response to climate change: capacity of deep saline aquifers to sequester CO 2 in solution[END_REF]. They can all be obtained experimentally either by calorimetric measurements [START_REF] Perron | Apparent Molal Volumes and Heat Capacities of some Sulfates and Carbonates in Water at 25 °C[END_REF] or with a densimeter [START_REF] Al Ghafri | Densities of Aqueous MgCl 2 (aq), CaCl 2 (aq), KI(aq), NaCl(aq), KCl(aq), AlCl 3 (aq), and (0.964 NaCl + 0.136 KCl)(aq) at Temperatures Between (283 and 472) K, Pressures up to 68.5 MPa, and Molalities up to 6 mol•kg -1[END_REF], or calculated from the temperature or the pressure derivative of the excess Gibbs energy. Some of the currently available software products that provide these values are:

-HSC Chemistry, developed by Outotec (2014), contains a module named AQUA that uses the Pitzer model to calculate the osmotic coefficient, the enthalpy, and the heat capacity of homogenous (single phase) aqueous electrolyte solutions.

-AQSOL001 (Aqueous Solutions ApS) can calculate the heat capacity of aqueous solutions. It is based on the eUNIQUAC model described by Thomsen (1997). The free downloadable version from their website includes Na + , H + , OH -, Cl -, SO 4

. A version containing more species can be purchased.

- Wang et al. (2002) described a model that is the sum of three interaction terms: a longrange electrostatic term, a short-range term and a middle-range term. These are represented by a Pitzer-Debye-Hückel equation, the UNIQUAC model, and a symmetrical second virial coefficient-type expression, respectively. This model is included in OLI systems software that calculates thermodynamic properties (Gibbs free energy, enthalpy, entropy, heat capacity and volume) and derived thermodynamic properties (e.g. density, osmotic pressure).

-Thermo-Calc [START_REF] Andersson | Thermo-Calc & DICTRA, computational tools for materials science[END_REF] calculates thermodynamic properties such as enthalpy and heat capacity. Two of this software's databases, AQS2 and TCAQ2, are dedicated to aqueous solutions. AQS2 contains 82 elements and the non-ideality of the aqueous solution phase is described using the complete revised HKF model (Helgeson-Kirkham-Flowers, Tanger and Helgeson, 1988), i.e. taking into account the Debye-Hückel Limiting Law term, ionic solvation, ionic association, as well as the binary, ternary and higher-order interaction terms. TCAQ2 contains 76 elements and the non-ideality of the aqueous solution phase is described using the extended SIT (Specific Interaction Theory) model, i.e. taking into account the Debye-Hückel

Limiting Law term, as well as the binary, ternary and higher-order interaction terms.

-VOPO [START_REF] Monnin | Density calculation and concentration scale conversions for natural waters[END_REF], for VOlumetric PrOperties, calculates the density of natural waters. The major species (Na, K, Ca, Mg, Cl, SO 4 , HCO 3 and CO 3 ) are taken into account in this program, which is based on the Pitzer model and uses the specific interaction parameters determined by [START_REF] Monnin | An ion interaction model for the volumetric properties of natural waters: Density of the solution and partial molal volumes of electrolytes to high concentrations at 25°C[END_REF] at 25 °C and 1 bar.

Most of the software listed above includes either thermal or volumetric properties. Only the software developed by OLI systems calculates both. However, this software contains its own database and it is not possible to use it to define new parameters for unknown systems. We therefore developed software that is able to calculate thermal and volumetric properties of aqueous species and can be used to optimize interaction parameters. PHREEQC V.3 [START_REF] Parkhurst | Description of Input and Examples for PHREEQC Version 3-A Computer Program for Speciation, Batch-Reaction, One-Dimensional Transport, and Inverse Geochemical Calculations[END_REF], a geochemical calculation code widely used by geochemists, offers possibilities for new applications. In particular, studies that were limited to surface and shallow subsurface environments can now extend to deeper geological and geothermal environments and to industrial processes (by including pressure effects).

PHREEQC code can calculate various water properties as activity coefficients of aqueous species and activity of water according to different formalisms including the Pitzer formalism.

However, additional thermodynamic properties that are not yet included in PHREEQC, namely apparent relative molar enthalpy and heat capacity of solutions, are required for these new fields of application. On the other hand, the calculation of density was already implemented in the last version of the PHREEQC code. The authors used the HKFmoRR (Helgeson-Kirkham-Flowers-modified-Redlich-Rosenfeld) equation to compute the molar volume [START_REF] Appelo | Equations for calculating hydrogeochemical reactions of minerals and gases such as CO2 at high pressures and temperatures[END_REF]:

V 𝑚,𝑖 = 𝑉 𝑖 °+ 𝐴 𝑉 0.5𝑧 𝑖 2 𝐼 0.5 1 + 𝑎i𝐵 𝛾 𝐼 0.5 + (𝑏 1,𝑖 + 𝑏 2,𝑖 𝑇 -228 + 𝑏 3,𝑖 (𝑇 -228)) 𝐼 𝑏 4,𝑖 (1)
Where 𝑎̇𝑖 is the ion-size parameter, 𝐵 𝛾 is the Debye-length parameter and 𝑏 𝑖 are for fitting the ionic strength dependence. And the apparent molar volume is calculated by:

V 𝜙 = ∑ 𝜈 𝑖 V 𝑚,𝑖 𝑖 (2)
In this study, the apparent molar volume is calculated by:

V 𝜙 = 𝑉 2 °̅̅̅ + 𝑉 𝑒𝑥 (3)
Where 𝑉 𝑒𝑥 is the excess molar volume calculated using the speciation and the volumetric interaction parameters following the Pitzer formalism. We chose to implement this approach in addition to that proposed by PHREEQC. Indeed, the thermal properties are calculated with the Pitzer's equations and the implementation of the Eq. ( 3) allows maintaining consistency in the formalism of properties calculation.

The aim of our work was to add to this software (PHREEQC) the calculation of apparent relative molar enthalpy, heat capacity, and density of complex aqueous solutions using the Pitzer formalism. The new numerical tool resulting from the present work is called PhreeSCALE.

The principal equations for calculating apparent relative molar enthalpy, heat capacity, and density of aqueous solutions are described below, followed by the modifications to PHREEQC. Then, the PhreeSCALE code is tested in two ways. First, the code is used to calculate the excess properties of several more or less complex chemical systems: all the apparent thermal and volumetric properties of the Na 2 SO 4 -H 2 O and MgSO 4 -H 2 O binary systems, the heat capacity of the Na-K-Cl-SO 4 -H 2 O system and the density of the Na-K-Ca-Mg-Cl-SO 4 -HCO 3 -H 2 O system. During this first series of tests, the code is applied using Pitzer interaction parameters (𝑌, 𝑌 𝐿 , 𝑌 𝐽 and 𝑌 𝑉 ) and their standard properties (𝑌 0 ) taken from literature in order to show that the new equations are correctly implemented in the code.

In a second part of this work, we demonstrate that the Pitzer equations can be used in combination with the HKF model which calculates the standard partial molar thermodynamic properties of aqueous species. But, the sets of interaction parameters coming from literature are generally not consistent with the calculated standard molar properties. Consequently they must be revised when necessary. PhreeSCALE is then coupled with the PEST optimization software [START_REF] Doherty | PEST: Model-Independent Parameter Estimation[END_REF]

Theory

The thermodynamic properties of aqueous solutions, such as the osmotic coefficient, heat capacity or density, can be obtained by deriving the excess Gibbs free energy with respect to temperature or pressure. We use the Pitzer model [START_REF] Pitzer | Activity coefficients in electrolyte solutions[END_REF] to describe the excess Gibbs energy of aqueous solutions, G ex (Eq. 4) and all the aforementioned thermodynamic properties are calculated using the same theoretical framework.

G 𝑒𝑥 (𝑤 𝑤 𝑅𝑇) ⁄ = { - 4𝐼𝐴 𝜙 𝑏 ln(1 + 𝑏𝐼 1 2 ⁄ ) + 2 ∑ ∑ 𝑚 𝑐 𝑚 𝑎 (𝐵 𝑐𝑎 + (∑ 𝑚 𝑐 𝑧 𝑐 𝑐 )𝐶 𝑐𝑎 ) 𝑎 𝑐 + ∑ ∑ 𝑚 𝑐 𝑚 𝑐 ′ (2Φ 𝑐𝑐 ′ + ∑ 𝑚 𝑎 𝜓 𝑐𝑐′𝑎 𝑎 ) 𝑐 ′ 𝑐 + ∑ ∑ 𝑚 𝑎 𝑚 𝑎 ′ (2Φ 𝑎𝑎 ′ + ∑ 𝑚 𝑐 𝜓 𝑎𝑎′𝑐 𝑐 ) 𝑎 ′ 𝑎 + ∑ 𝑚 𝑛 2 𝜆 𝑛𝑛 𝑛 + ∑ 𝑚 𝑛 3 𝜇 𝑛𝑛𝑛 𝑛 +2 ∑ ∑ 𝑚 𝑛 𝑚 𝑛 ′ 𝜆 𝑛𝑛 ′ 𝑛 ′ 𝑛 + 3 ∑ ∑ 𝑚 𝑛 2 𝑚 𝑛 ′ 𝜇 𝑛𝑛𝑛′ 𝑛 ′ 𝑛 +6 ∑ ∑ ∑ 𝑚 𝑛 𝑚 𝑛 ′ 𝑚 𝑛 " 𝜇 𝑛𝑛′𝑛" 𝑛" 𝑛 ′ 𝑛 +2 ∑ ∑ 𝑚 𝑛 𝑚 𝑐 𝜆 𝑛𝑐 𝑐 𝑛 + 3 ∑ ∑ 𝑚 𝑛 2 𝑚 𝑐 𝜇 𝑛𝑛𝑐 𝑐 𝑛 +2 ∑ ∑ 𝑚 𝑛 𝑚 𝑎 𝜆 𝑛𝑎 𝑎 𝑛 + 3 ∑ ∑ 𝑚 𝑛 2 𝑚 𝑎 𝜇 𝑛𝑛𝑎 𝑎 𝑛 + ∑ ∑ ∑ 𝑚 𝑛 𝑚 𝑐 𝑚 𝑎 𝜁 𝑛𝑐𝑎 𝑎 𝑐 𝑛 + ∑ ∑ ∑ 𝑚 𝑛 𝑚 𝑐 𝑚 𝑐 ′ 𝜂 𝑛𝑐𝑐′ 𝑐 ′ 𝑐 𝑛 + ∑ ∑ ∑ 𝑚 𝑛 𝑚 𝑐 𝑚 𝑎 ′ 𝜂 𝑛𝑎𝑎′ 𝑎 ′ 𝑎 𝑛 +6 ∑ ∑ ∑ 𝑚 𝑐 𝑚 𝑛 𝑚 𝑛 ′ 𝜇 𝑐𝑛𝑛′ 𝑛 ′ 𝑛 𝑐 +6 ∑ ∑ ∑ 𝑚 𝑎 𝑚 𝑛 𝑚 𝑛 ′ 𝜇 𝑎𝑛𝑛′ 𝑛 ′ 𝑛 𝑎 } (4)
The different parameters (𝐴 𝜙 , 𝐵 𝑐𝑎 …) are not detailed here but they are summarized in the appendix of [START_REF] Lach | A new Pitzer parameterization for the binary NaOH-H2O and ternary NaOH-NaCl-H2O and NaOH-LiOH-H2O systems up to NaOH solid salt saturation[END_REF].

The temperature dependence of every interaction parameter is assumed to have the following expression (which is a combination of equations from [START_REF] Møller | The prediction of mineral solubilities in natural waters: A chemical equilibrium model for the Na-Ca-Cl-SO4-H2O system, to high temperature and concentration[END_REF] and [START_REF] Christov | Chemical equilibrium model of solution behavior and solubility in the H-Na-K-OH-Cl-HSO 4 -SO 4 -H 2 O system to high concentration and temperature[END_REF]

where Y designates any of the adjustable specific interaction parameters.

The heat capacity (𝑐 𝑝 ) and the density of the solution (ρ) are measurable properties and must therefore be described formally. They can be calculated from the apparent molar heat capacity (𝐶 𝑝,𝜙 ) and the apparent molar volume (𝑉 𝜙 ), respectively, which are obtained from the derivatives of excess Gibbs energy.

Apparent properties

The apparent property (𝑍 𝜙 ) is defined by the following equation [START_REF] Hubert | Aqueous solution of sodium sulfate. Determination of the dissolution enthalpy at 25, 27.5 and 45°C[END_REF]:

𝑍 𝜙 = (𝑍 (𝑛 𝑤 , ∑ 𝑛 𝑖 𝑖≠𝑤 ) -𝑍(𝑛 𝑤 , 0)) ∑ 𝑛 𝑖 𝑖≠𝑤 ⁄ (6)
where 𝑍(𝑛 𝑤 , ∑ 𝑛 𝑖

𝑖≠𝑤

) is the property value of a solution containing 𝑛 𝑤 moles of solvent and ∑ 𝑛 𝑖 𝑖≠𝑤 moles of solutes and 𝑍(𝑛 𝑤 , 0) is the value of the property in the pure solvent.

The apparent relative molar enthalpy

The relative enthalpy (𝐿) of an electrolyte solution is related to 𝐺 𝑒𝑥 by Eq. ( 7) [START_REF] Pitzer | Activity coefficients in electrolyte solutions[END_REF].

𝐿 = -𝑇 2 [𝜕 ( 𝐺 𝑒𝑥 𝑇 ) 𝜕𝑇 ⁄ ] 𝑃,𝑚 (7)
Using Eq. ( 6) and since the relative enthalpy is null for the pure solvent solution:

𝐿 𝜙 = 𝐿 ∑ 𝑛 𝑖 𝑖≠𝑤 ⁄ (8)
Therefore, applying Eq. (7-8) to Eq. ( 4), the apparent relative molar enthalpy is calculated by the following equation:

𝐿 𝜙 = ( 𝑅𝑇 2 ∑ 𝑚 𝑖 𝑖 ) { 𝐼 𝐴 𝐿 𝑅𝑇 2 ln(1+𝑏𝐼 1 2 ⁄ ) 𝑏 -2 ∑ ∑ 𝑚 𝑐 𝑚 𝑎 (𝐵 𝑐𝑎 𝐿 + (∑ 𝑚 𝑐 𝑧 𝑐 𝑐 )𝐶 𝑐𝑎 𝐿 ) 𝑎 𝑐 -∑ ∑ 𝑚 𝑐 𝑚 𝑐 ′ (2Φ 𝑐𝑐 ′ 𝐿 + ∑ 𝑚 𝑎 𝜓 𝑐𝑐′𝑎 𝐿 𝑎 ) 𝑐 ′ 𝑐 -∑ ∑ 𝑚 𝑎 𝑚 𝑎 ′ (2Φ 𝑎𝑎 ′ 𝐿 + ∑ 𝑚 𝑐 𝜓 𝑎𝑎′𝑐 𝐿 𝑐 ) 𝑎 ′ 𝑎 -∑ 𝑚 𝑛 2 𝜆 𝑛𝑛 𝐿 𝑛 -∑ 𝑚 𝑛 3 𝜇 𝑛𝑛𝑛 𝐿 𝑛 -2 ∑ ∑ 𝑚 𝑛 𝑚 𝑛 ′ 𝜆 𝑛𝑛 ′ 𝐿 𝑛 ′ 𝑛 -3 ∑ ∑ 𝑚 𝑛 2 𝑚 𝑛 ′ 𝜇 𝑛𝑛𝑛′ 𝐿 𝑛 ′ 𝑛 -6 ∑ ∑ ∑ 𝑚 𝑛 𝑚 𝑛 ′ 𝑚 𝑛 " 𝜇 𝑛𝑛′𝑛" 𝐿 𝑛" 𝑛 ′ 𝑛 -2 ∑ ∑ 𝑚 𝑛 𝑚 𝑐 𝜆 𝑛𝑐 𝐿 𝑐 𝑛 -3 ∑ ∑ 𝑚 𝑛 2 𝑚 𝑐 𝜇 𝑛𝑛𝑐 𝐿 𝑐 𝑛 -2 ∑ ∑ 𝑚 𝑛 𝑚 𝑎 𝜆 𝑛𝑎 𝐿 𝑎 𝑛 -3 ∑ ∑ 𝑚 𝑛 2 𝑚 𝑎 𝜇 𝑛𝑛𝑎 𝐿 𝑎 𝑛 -∑ ∑ ∑ 𝑚 𝑛 𝑚 𝑐 𝑚 𝑎 𝜁 𝑛𝑐𝑎 𝐿 𝑎 𝑐 𝑛 -∑ ∑ ∑ 𝑚 𝑛 𝑚 𝑐 𝑚 𝑐 ′ 𝜂 𝑛𝑐𝑐′ 𝐿 𝑐 ′ 𝑐 𝑛 -∑ ∑ ∑ 𝑚 𝑛 𝑚 𝑐 𝑚 𝑎 ′ 𝜂 𝑛𝑎𝑎′ 𝐿 𝑎 ′ 𝑎 𝑛 -6 ∑ ∑ ∑ 𝑚 𝑐 𝑚 𝑛 𝑚 𝑛 ′ 𝜇 𝑐𝑛𝑛′ 𝐿 𝑛 ′ 𝑛 𝑐 -6 ∑ ∑ ∑ 𝑚 𝑎 𝑚 𝑛 𝑚 𝑛 ′ 𝜇 𝑎𝑛𝑛′ 𝐿 𝑛 ′ 𝑛 𝑎 } (9)
where:

𝐴 𝐿 = 4𝑅𝑇²(𝜕𝐴 𝜙 𝜕𝑇 ⁄ ) (10)
and

𝑌 𝐿 (𝑇) = ∂Y 𝜕𝑇 (11)

The apparent molar heat capacity

The heat capacity is related to the relative enthalpy by the following equation [START_REF] Pitzer | Activity coefficients in electrolyte solutions[END_REF]:

𝐶 𝑝 ′ = 𝐶 𝑝 °+ (𝜕𝐿 𝜕𝑇 ⁄ ) 𝑃,𝑚 (12)
where 𝐶 𝑝 ° is a linear combination of the standard molar heat capacities of the solvent and the solutes. Therefore:

𝐶 𝑝 ′ = 𝑛 𝑤 𝐶 𝑝,𝑤 °+ ∑ 𝑛 𝑖 𝐶 𝑝,𝑖 0 𝑖 + (𝜕𝐿 𝜕𝑇 ⁄ ) 𝑃,𝑚 (13)
Following Eq.( 6), the apparent heat capacity of a multicomponent solution, 𝐶 𝑝,𝜙 , is defined by the following equation:

𝐶 𝑝,𝜙 = (𝐶 𝑝 ′ -𝑛 𝑤 𝐶 𝑝,𝑤 °) ∑ 𝑛 𝑖 𝑖 ⁄ (14)
Combining Eq. ( 13) and Eq. ( 14), we obtain:

𝐶 𝑝,𝜙 = (∑ 𝑛 𝑖 𝐶 𝑝,𝑖 0 𝑖 ) ∑ 𝑛 𝑖 𝑖 ⁄ + (𝜕 ( 𝐿 ∑ 𝑛 𝑖 𝑖 ) 𝜕𝑇 ⁄ ) 𝑃,𝑚 (15)
Using Eq. ( 8), the heat capacity is related to the apparent relative molar enthalpy by the following equation:

𝐶 𝑝,𝜙 = 𝐶 𝑝2 °̅̅̅̅̅ + (𝜕𝐿 𝜙 𝜕𝑇 ⁄ ) 𝑃,𝑚 (16)
Here, 𝐶 𝑝2 °̅̅̅̅̅ is the standard partial heat capacity of all of the solutes at infinite dilution:

𝐶 𝑝2 °̅̅̅̅̅ = (∑ 𝑚 𝑖 𝐶 𝑝,𝑖 0 𝑖 ) ∑ 𝑚 𝑖 𝑖 ⁄ (17)
where 𝐶 𝑝,𝑖 0 is calculated with the equation of [START_REF] Helgeson | Theoretical prediction of the thermodynamic behavior of aqueous electrolytes by high pressures and temperatures; IV, Calculation of activity coefficients, osmotic coefficients, and apparent molal and standard and relative partial molal properties to 600°C[END_REF] summarized in [START_REF] Johnson | SUPCRT92: A software package for calculating the standard molal thermodynamic properties of minerals, gases, aqueous species, and reactions from 1 to 5000 bar[END_REF].

Therefore, applying Eq. ( 16) to Eq. ( 9) gives the following equation for the apparent heat capacity:

𝐶 𝑝,𝜙 = 𝐶 𝑝2 °̅̅̅̅̅ + ( 𝑅𝑇² ∑ 𝑚 𝑖 𝑖 ) { 𝐼 𝐴 𝐽 𝑅𝑇 2 ln(1+𝑏𝐼 1 2 ⁄ ) 𝑏 -2 ∑ ∑ 𝑚 𝑐 𝑚 𝑎 (𝐵 𝑐𝑎 𝐽 + (∑ 𝑚 𝑐 𝑧 𝑐 𝑐 )𝐶 𝑐𝑎 𝐽 ) 𝑎 𝑐 -∑ ∑ 𝑚 𝑐 𝑚 𝑐 ′ (2Φ 𝑐𝑐 ′ 𝐽 + ∑ 𝑚 𝑎 𝜓 𝑐𝑐′𝑎 𝐽 𝑎 ) 𝑐 ′ 𝑐 -∑ ∑ 𝑚 𝑎 𝑚 𝑎 ′ (2Φ 𝑎𝑎 ′ 𝐽 + ∑ 𝑚 𝑐 𝜓 𝑎𝑎′𝑐 𝐽 𝑐 ) 𝑎 ′ 𝑎 -∑ 𝑚 𝑛 2 𝜆 𝑛𝑛 𝐽 𝑛 -∑ 𝑚 𝑛 3 𝜇 𝑛𝑛𝑛 𝐽 𝑛 -2 ∑ ∑ 𝑚 𝑛 𝑚 𝑛 ′ 𝜆 𝑛𝑛 ′ 𝐽 𝑛 ′ 𝑛 -3 ∑ ∑ 𝑚 𝑛 2 𝑚 𝑛 ′ 𝜇 𝑛𝑛𝑛′ 𝐽 𝑛 ′ 𝑛 -6 ∑ ∑ ∑ 𝑚 𝑛 𝑚 𝑛 ′ 𝑚 𝑛 " 𝜇 𝑛𝑛′𝑛" 𝐽 𝑛" 𝑛 ′ 𝑛 -2 ∑ ∑ 𝑚 𝑛 𝑚 𝑐 𝜆 𝑛𝑐 𝐽 𝑐 𝑛 -3 ∑ ∑ 𝑚 𝑛 2 𝑚 𝑐 𝜇 𝑛𝑛𝑐 𝐽 𝑐 𝑛 -2 ∑ ∑ 𝑚 𝑛 𝑚 𝑎 𝜆 𝑛𝑎 𝐽 𝑎 𝑛 -3 ∑ ∑ 𝑚 𝑛 2 𝑚 𝑎 𝜇 𝑛𝑛𝑎 𝐽 𝑎 𝑛 -∑ ∑ ∑ 𝑚 𝑛 𝑚 𝑐 𝑚 𝑎 𝜁 𝑛𝑐𝑎 𝐽 𝑎 𝑐 𝑛 -∑ ∑ ∑ 𝑚 𝑛 𝑚 𝑐 𝑚 𝑐 ′ 𝜂 𝑛𝑐𝑐′ 𝐽 𝑐 ′ 𝑐 𝑛 -∑ ∑ ∑ 𝑚 𝑛 𝑚 𝑐 𝑚 𝑎 ′ 𝜂 𝑛𝑎𝑎′ 𝐽 𝑎 ′ 𝑎 𝑛 -6 ∑ ∑ ∑ 𝑚 𝑐 𝑚 𝑛 𝑚 𝑛 ′ 𝜇 𝑐𝑛𝑛′ 𝐽 𝑛 ′ 𝑛 𝑐 -6 ∑ ∑ ∑ 𝑚 𝑎 𝑚 𝑛 𝑚 𝑛 ′ 𝜇 𝑎𝑛𝑛′ 𝐽 𝑛 ′ 𝑛 𝑎 } (18)
where:

𝐴 𝐽 = (2 𝑇 ⁄ )𝐴 𝐿 + 4𝑅𝑇²(𝜕²𝐴 𝜙 𝜕𝑇² ⁄ ) (19)
And:

𝑌 𝐽 = 𝜕𝑌 𝐿 𝜕𝑇 + 2 𝑇 𝑌 𝐿 (20)
The mass heat capacity of the solution (𝑐 𝑝 ) can then be calculated from the apparent heat capacity using the following equation:

𝑐 𝑝 = 𝐶 𝑝,𝜙 ∑ 𝑚 𝑖 𝑖 + 1000𝑐 𝑝,𝑤 1000 + ∑ 𝑚 𝑖 𝑀 𝑖 𝑖 (21)

The apparent molar volume

The volume of an electrolyte solution (𝑉) is related to 𝐺 𝑒𝑥 by Eq. ( 22):

𝑉 = 𝑉°+ [𝜕𝐺 𝑒𝑥 𝜕𝑃 ⁄ ] 𝑇,𝑚 (22)
where V° is a linear combination of the standard volume of the solvent and the solutes.

Therefore:

𝑉 = 𝑛 𝑤 𝑉 𝑤 °+ ∑ 𝑛 𝑖 𝑉 𝑖 °𝑖 + [𝜕𝐺 𝑒𝑥 𝜕𝑃 ⁄ ] 𝑇,𝑚 (23)
Following Eq. ( 6), the apparent molar volume is expressed by the equation:

𝑉 𝜙 = (V -𝑛 𝑤 𝑉 𝑤 °) ∑ 𝑛 𝑖 𝑖 ⁄ (24)
Therefore:

𝑉 𝜙 = 𝑉 2 °̅̅̅ + [𝜕𝐺 𝑒𝑥 𝜕𝑃 ⁄ ] 𝑇,𝑚 ∑ 𝑛 𝑖 𝑖 ⁄ (25)
where 𝑉 2 °̅̅̅ is the standard partial volume of all solutes at infinite dilution:

𝑉 2 °̅̅̅ = (∑ 𝑚 𝑖 𝑉 𝑖 °𝑖 ) ∑ 𝑚 𝑖 𝑖 ⁄ (26)
Applying Eq. ( 25) to Eq. ( 4) gives the following equation for the apparent molar volume of any aqueous solution:

𝑉 𝜙 = 𝑉 2 °̅̅̅ + ( 𝑅 ′ 𝑇 ∑ 𝑚 𝑖 𝑖 ) { 𝐼𝐴 𝑉 𝑅 ′ 𝑇 ln(1+𝑏𝐼 1 2 ⁄ ) 𝑏 + 2 ∑ ∑ 𝑚 𝑐 𝑚 𝑎 (𝐵 𝑐𝑎 𝑉 + (∑ 𝑚 𝑐 𝑧 𝑐 𝑐 )𝐶 𝑐𝑎 𝑉 ) 𝑎 𝑐 + ∑ ∑ 𝑚 𝑐 𝑚 𝑐 ′ (2Φ 𝑐𝑐 ′ 𝑉 + ∑ 𝑚 𝑎 𝜓 𝑐𝑐′𝑎 𝑉 𝑎 ) 𝑐 ′ 𝑐 + ∑ ∑ 𝑚 𝑎 𝑚 𝑎 ′ (2Φ 𝑎𝑎 ′ 𝑉 + ∑ 𝑚 𝑐 𝜓 𝑎𝑎′𝑐 𝑉 𝑐 ) 𝑎 ′ 𝑎 + ∑ 𝑚 𝑛 2 𝜆 𝑛𝑛 𝑉 𝑛 + ∑ 𝑚 𝑛 3 𝜇 𝑛𝑛𝑛 𝑉 𝑛 +2 ∑ ∑ 𝑚 𝑛 𝑚 𝑛 ′ 𝜆 𝑛𝑛 ′ 𝑉 𝑛 ′ 𝑛 + 3 ∑ ∑ 𝑚 𝑛 2 𝑚 𝑛 ′ 𝜇 𝑛𝑛𝑛′ 𝑉 𝑛 ′ 𝑛 +6 ∑ ∑ ∑ 𝑚 𝑛 𝑚 𝑛 ′ 𝑚 𝑛 " 𝜇 𝑛𝑛′𝑛" 𝑉 𝑛" 𝑛 ′ 𝑛 +2 ∑ ∑ 𝑚 𝑛 𝑚 𝑐 𝜆 𝑛𝑐 𝑉 𝑐 𝑛 + 3 ∑ ∑ 𝑚 𝑛 2 𝑚 𝑐 𝜇 𝑛𝑛𝑐 𝑉 𝑐 𝑛 +2 ∑ ∑ 𝑚 𝑛 𝑚 𝑎 𝜆 𝑛𝑎 𝑉 𝑎 𝑛 + 3 ∑ ∑ 𝑚 𝑛 2 𝑚 𝑎 𝜇 𝑛𝑛𝑎 𝑉 𝑎 𝑛 + ∑ ∑ ∑ 𝑚 𝑛 𝑚 𝑐 𝑚 𝑎 𝜁 𝑛𝑐𝑎 𝑉 𝑎 𝑐 𝑛 + ∑ ∑ ∑ 𝑚 𝑛 𝑚 𝑐 𝑚 𝑐 ′ 𝜂 𝑛𝑐𝑐′ 𝑉 𝑐 ′ 𝑐 𝑛 + ∑ ∑ ∑ 𝑚 𝑛 𝑚 𝑐 𝑚 𝑎 ′ 𝜂 𝑛𝑎𝑎′ 𝑉 𝑎 ′ 𝑎 𝑛 +6 ∑ ∑ ∑ 𝑚 𝑐 𝑚 𝑛 𝑚 𝑛 ′ 𝜇 𝑐𝑛𝑛′ 𝑉 𝑛 ′ 𝑛 𝑐 +6 ∑ ∑ ∑ 𝑚 𝑎 𝑚 𝑛 𝑚 𝑛 ′ 𝜇 𝑎𝑛𝑛′ 𝑉 𝑛 ′ 𝑛 𝑎 } (27)
where:

𝐴 𝑉 = -4𝑅 ′ 𝑇(𝜕𝐴 𝜙 𝜕𝑃 ⁄ ) (28)
Technically, the volumetric interaction parameters are the pressure derivatives of the interaction parameters:

𝑌 𝑉 = (𝜕𝑌 𝜕𝑃 ⁄ ) (29)
However, so far, the interaction parameters as given by Eq. ( 5) are considered to be independent of pressure. Therefore, despite the fact that it is not consistent with Eq. ( 5), an additional independent set of parameters is included in the database to account for temperature-variable pressure effects. Many authors have proposed different temperature dependences for these parameters [START_REF] Krumgalz | Volumetric Ion Interaction Parameters for Single-Solute Aqueous Electrolyte Solutions at Various Temperatures[END_REF][START_REF] Saluja | Apparent Molar Heat Capacities and Volumes of Mixed Electrolytes: [NaCl(aq) + CaCl2(aq)], [NaCl(aq) + MgCl2(aq)], and [CaCl2(aq) + MgCl2(aq)][END_REF]. We chose the one proposed by Phutela and Pitzer (1986b) which can be formally made compatible with Eq. ( 5):

𝑌 𝑉 = 𝑏 0 /𝑇 + 𝑏 1 + 𝑏 2 𝑇 + 𝑏 3 𝑇 2 + 𝑏 4 𝑇 3 + 𝑏 5 𝑃_𝑏𝑎𝑟 (30)
Integrating Eq. ( 30), makes it possible to link the coefficients 𝑏 𝑖 with the coefficient 𝑎 𝑖 (Eq.

(5)) by the following expression:

𝑎 1 = 𝑏 1 𝑃 + 𝑏 5 2
𝑃 2 + 𝑐 1 ; 𝑎 2 = 𝑏 2 𝑃 + 𝑐 2 ; 𝑎 3 = 𝑏 3 𝑃 + 𝑐 3 ; 𝑎 4 = 𝑏 4 𝑃 + 𝑐 4 ; 𝑎 5 = 𝑏 0 𝑃 + 𝑐 5

(31) However, our interaction parameters are determined at 1 bar below 100 °C and at saturation pressure above 100 °C, as is very often done. Therefore, there is an implicit pressure dependence which is not compatible with Eq. ( 31). Indeed, knowing 𝑎 𝑖 and 𝑏 𝑖 coefficients, it is not possible to determine one 𝑐 𝑖 coefficient that is valid at both 1 bar and saturation pressure. The entire set of interaction parameters must be re-optimized, which is beyond the scope of this work. Therefore, we continued to use pressure-independent interaction parameters.

From the apparent molar volume, the density of the solution can be calculated with the equation:

ρ = (1000 + ∑ 𝑚 𝑖 𝑀 𝑖 𝑖 ) (𝑉 𝜙 ∑ 𝑚 𝑖 𝑖 + 10 6 (𝜌 𝑤 ) -1 ) ⁄ (32)
All of the above equations, and specifically Eq. ( 9), (18), and ( 27), can theoretically be applied to any complex aqueous solution. Indeed, all possible aqueous speciation, including neutral species, is considered. As such, these equations represent an extension of the theoretical work of [START_REF] Felmy | The prediction of borate mineral equilibria in natural waters: Application to Searles Lake, California[END_REF], [START_REF] Pabalan | Apparent molar heat capacity and other thermodynamic properties of aqueous potassium chloride solutions to high temperatures and pressures[END_REF], [START_REF] Monnin | An ion interaction model for the volumetric properties of natural waters: Density of the solution and partial molal volumes of electrolytes to high concentrations at 25°C[END_REF][START_REF] Pitzer | Activity coefficients in electrolyte solutions[END_REF].

Computer program: PHREEQC modification

PHREEQC V3 [START_REF] Parkhurst | Description of Input and Examples for PHREEQC Version 3-A Computer Program for Speciation, Batch-Reaction, One-Dimensional Transport, and Inverse Geochemical Calculations[END_REF] was modified to include the equations described above, and the new version is called PhreeSCALE.

The sources

The temperature dependence of interaction parameters was modified in accordance with Eq.

(5) and the temperature dependence of the derived interaction parameters were added in accordance with Eq. ( 16) and Eq. ( 25) together with the volumetric interaction parameters (Eq. ( 30)). The apparent properties and the equations necessary for computing were also added (the temperature and pressure derivatives of the dielectric constant and the density of pure water, the heat capacity of water, the standard partial molar heat capacity and the volume of aqueous species at infinite dilution, etc.). The equation for calculating the osmotic coefficient was also modified to take into account partial dissociation [START_REF] Lach | A new Pitzer parameterization for the binary NaOH-H2O and ternary NaOH-NaCl-H2O and NaOH-LiOH-H2O systems up to NaOH solid salt saturation[END_REF] according to the following equation:

𝜙 𝑃ℎ𝑟𝑒𝑒𝑆𝐶𝐴𝐿𝐸 = 𝜙 𝑃ℎ𝑟𝑒𝑒𝑞𝑐 ∑ 𝑚 𝑖 𝑖 ∑ 𝜈 𝑖 𝑚 𝑖 𝑖 ⁄ (33)
where 𝜈 𝑖 is the number of moles of ions formed by the complete dissociation of one molecule of solute.

The database

In order to run PhreeSCALE, parameters used in the new equations must be readable from the associated thermodynamic database. One line is therefore added to each species with the keyword HKF in the SOLUTION_SPECIES block. It gathers 8 parameters: the 7 HKF parameters (a 1 , a 2 , a 3 , a 4 , c 1 , c 2 , ω Tr,Pr ) used to calculate the standard partial molar heat capacity and the volume of any aqueous species i at infinite dilution [START_REF] Helgeson | Theoretical prediction of the thermodynamic behavior of aqueous electrolytes by high pressures and temperatures; IV, Calculation of activity coefficients, osmotic coefficients, and apparent molal and standard and relative partial molal properties to 600°C[END_REF] and one additional parameter corresponding to 𝜈 𝑖 in Eq. ( 33). The first seven parameters can be found in the Thermoddem database [START_REF] Blanc | Thermoddem a database devoted to waste minerals[END_REF]. The energy unit is the Joule. The volumetric interaction parameters in the PITZER block were also added, which required the creation of the corresponding keywords B0V, B1V, B2V, C0V, THETAV, LAMDAV, ZETAV, PSIV, MUV and ETAV.

Selected output

Several new output key words were created and correspond to each property. They can be called in the SELECTED_OUTPUT block (declared in the input file) in order to write the corresponding properties in a selected output file. These keywords are described in Table 1.

Results

PhreeSCALE validation using interaction parameters from literature

Before any practical application of the code, the first task to achieve is checking the program.

To do so, we used PhreeSCALE with sets of interaction parameters available in the literature and compared our results to those obtained by the same authors. The test covers four chemical systems with respect to various interaction parameters (Y, Y 𝐿 , Y 𝐽 and Y 𝑉 ) and standard properties (C 𝑝 0 and V 0 ). The first chemical system is the Na 2 SO 4 -H 2 O binary system for which [START_REF] Rogers | High-temperature thermodynamic properties of aqueous sodium sulfate solutions[END_REF] gave interaction parameters to calculate the osmotic coefficient, the apparent relative molar enthalpy and the apparent molar heat capacity. In addition, Phutela and Pitzer (1986b) measured the volumetric properties and determined the volumetric interaction parameters for the same system. The second system is the 2-2 type MgSO 4 -H 2 O binary system, for which Phutela andPitzer (1986a, 1986b) measured the volumetric properties and determined the volumetric interaction parameters. The third system considered for checking the program implemented in PhreeSCALE is the Na-K-Cl-SO 4 -H 2 O system for which [START_REF] Conti | Mixed interaction terms in the use of pitzer equations: application to the heat capacities of the reciprocal system Na+/K+/Cl-/SO42-/H2O at high temperatures[END_REF][START_REF] Conti | Heat capacities of aqueous mixed electrolyte solutions at high temperatures. Application of the Pitzer equations to the mixed system K-Na-Cl-SO4[END_REF] determined the interaction parameters for calculating the heat capacity. Finally, the system Na-K-Ca-Mg-Cl-SO 4 -HCO 3 -H 2 O is considered in order to check the calculation of density at 25 °C using the parameters of [START_REF] Monnin | An ion interaction model for the volumetric properties of natural waters: Density of the solution and partial molal volumes of electrolytes to high concentrations at 25°C[END_REF].

To compare the values of a property x for a data point i calculated by the different authors (superscript "cal") and by PhreeSCALE (superscript "PhreeSCALE"), the absolute value of average relative deviation is calculated with the following equation:

∆ 𝐴𝐴𝐷 = (1 𝑁 ⁄ ) ∑|(𝑥 𝑖 𝑐𝑎𝑙 -𝑥 𝑖 PhreeSCALE ) 𝑥 𝑖 𝑐𝑎𝑙 ⁄ | 𝑁 𝑖 ( 34
)
where N is the total number of data points.

If the comparison concerns apparent properties, it can be necessary to correct the values computed by PhreeSCALE in order to make them compatible. Indeed, in this study the apparent properties are computed considering the ionic species. Authors like Phutela andPitzer (1986a, 1986b) or [START_REF] Rogers | High-temperature thermodynamic properties of aqueous sodium sulfate solutions[END_REF] consider the electrolyte instead the ionic species. The correction to make results from the unicity of a solution property (𝑐 𝑝 or 𝜌):

𝑐 𝑝 𝑙𝑖𝑡𝑡 = 𝑐 𝑝 𝑃ℎ𝑟𝑒𝑒𝑆𝐶𝐴𝐿𝐸 ( 35
)
Where 𝑐 𝑝 𝑙𝑖𝑡𝑡 correspond at the heat capacity provides by the literature.

Using Eq. ( 21) gives:

𝐶 𝑝,
where the subscript 𝑘 refers to an electrolyte and the subscript 𝑖 to an ionic species. As the denominator is identical in the two members of Eq. ( 36), it comes:

𝐶 𝑝,𝜙 𝑙𝑖𝑡𝑡 = 𝐶 𝑝,𝜙 𝑃ℎ𝑟𝑒𝑒𝑆𝐶𝐴𝐿𝐸 ∑ 𝑚 𝑖 𝑖 ∑ 𝑚 𝑘 𝑘 (37)
The same demonstration can be made on the density. Finally, the more general relationship can be written:

𝑌 𝜙 𝑙𝑖𝑡𝑡 = 𝛿. 𝑌 𝜙 𝑃ℎ𝑟𝑒𝑒𝑆𝐶𝐴𝐿𝐸 (38)
Where:

𝛿 = ∑ 𝑚 𝑖 𝑖 ∑ 𝑚 𝑘 𝑘 (39)
For instance, in the case of Na 2 SO 4 the value computed by PhreeSCALE must be multiplied by 3 to be compared with the value proposed by Phutela and Pitzer (1986a) and [START_REF] Rogers | High-temperature thermodynamic properties of aqueous sodium sulfate solutions[END_REF]:

𝑌 𝜙 𝑙𝑖𝑡𝑡 = ( 𝑚 𝑁𝑎 + + 𝑚 𝑆𝑂 4 2- 𝑚 𝑁𝑎 2 𝑆𝑂 4 ) . 𝑌 𝜙 𝑃ℎ𝑟𝑒𝑒𝑆𝐶𝐴𝐿𝐸 (40)
Since 𝑚 𝑁𝑎 + = 2 * 𝑚 𝑁𝑎 2 𝑆𝑂 4 and 𝑚 𝑆𝑂 4 2-= 𝑚 𝑁𝑎 2 𝑆𝑂 4 , Eq. ( 40) becomes:

𝑌 𝜙 𝑙𝑖𝑡𝑡 = 3. 𝑌 𝜙 𝑃ℎ𝑟𝑒𝑒𝑆𝐶𝐴𝐿𝐸 (41)
At last, the standard deviation of each case is summarized in Table 2. On the binary systems, the largest deviation (2%) is obtained for the apparent molar heat capacities of MgSO 4 -H 2 O system at 20 bar. For more complex system, the properties of solution are well calculated.

Indeed, the largest deviation is equal to 0.415% for the heat capacities of Na-K-Cl-SO 4 -H 2 O system.

In conclusion, the new code PhreeSCALE is able to calculate the thermal and volumetric properties of complex solutions using the Pitzer interaction parameters released by the authors.

4.2.Optimization of interaction parameters

In many studies of brine systems, the standard properties (C 𝑝 0 and V 0 ) are considered as adjustable parameters. The authors use these parameters to optimize their fit of numerical results with experimental data. This approach can lead to inconsistencies of these standard properties with values determined elsewhere in the literature. In order to avoid this, standard properties can be selected from reference studies and thus be excluded from the parameterization procedure. In the present study, we chose to use the HKF formalism [START_REF] Helgeson | Theoretical prediction of the thermodynamic behavior of aqueous electrolytes by high pressures and temperatures; IV, Calculation of activity coefficients, osmotic coefficients, and apparent molal and standard and relative partial molal properties to 600°C[END_REF] to calculate the standard properties. The consequence of using this approach is a loss of coherency between the standard properties and the set of interaction parameters proposed by the authors. A new set of interaction parameters compatible with the calculated C 𝑝 0 and V 0 values must therefore be determined. This is what we did for the two binary systems: Na 2 SO 4 -H 2 O and MgSO 4 -H 2 O.

The Na 2 SO 4 -H 2 O binary system

In Fig. 1, the standard partial molar properties (heat capacity and volume) of the Na 2 SO 4 electrolyte are plotted as a function of temperature at saturation pressure. Symbols represent the C 𝑝 0 and V 0 values optimized by [START_REF] Rogers | High-temperature thermodynamic properties of aqueous sodium sulfate solutions[END_REF] and Phutela and Pitzer (1986b), respectively, and lines represent values calculated with the HKF equations [START_REF] Helgeson | Theoretical prediction of the thermodynamic behavior of aqueous electrolytes by high pressures and temperatures; IV, Calculation of activity coefficients, osmotic coefficients, and apparent molal and standard and relative partial molal properties to 600°C[END_REF]; [START_REF] Johnson | SUPCRT92: A software package for calculating the standard molal thermodynamic properties of minerals, gases, aqueous species, and reactions from 1 to 5000 bar[END_REF]). Despite the fact that results are generally consistent, a nonnegligible difference can be observed for the standard heat capacity over the entire temperature range, which may have significant consequences on the description of the heat capacity of Na 2 SO 4 solutions. This is confirmed in Fig. 2 in which experimental data from the literature (diamonds) are plotted together with values calculated (dotted-line) using the standard heat capacity calculated with the HFK parameters provided by [START_REF] Blanc | Thermoddem a database devoted to waste minerals[END_REF], instead of the heat capacity given by [START_REF] Rogers | High-temperature thermodynamic properties of aqueous sodium sulfate solutions[END_REF], while still using their specific interaction parameters. At 10 °C, the fit with experimental data is not good for Na 2 SO 4 concentrations above 0.25 mol.kg -1 (Fig. 2a). This is consistent with the fact that the model of [START_REF] Rogers | High-temperature thermodynamic properties of aqueous sodium sulfate solutions[END_REF] is valid for temperatures greater than 25 °C. At higher temperatures, the match between experimental and calculated data is satisfactory.

In order to extend the set of interaction parameters of Na 2 SO 4 over the entire temperature and salinity ranges, a new set of interaction parameters was determined using experimental data for the heat capacity and the osmotic coefficient taken from studies reported in Table 3. The new set of interaction parameters is reported in Table 4 and the heat capacity calculated with this new set of interaction parameters is plotted in Fig. 2 and the osmotic coefficient is plotted in Fig. 3 at various temperatures ranging between 0 and 170°C. Data are represented up to saturation with respect to Na 2 SO 4 (s). The a i parameters were optimized using simultaneously osmotic coefficient and heat capacity data.

There are no significant differences between the standard volume of Na 2 SO 4 given by Phutela and Pitzer (1986b) and the one calculated with the HKF equations. Therefore, the density values calculated with the interaction parameters of Phutela and Pitzer (1986b) and the standard volume calculated with the HKF parameters describe very satisfactorily the experimental values (Fig. 4) without further parameter optimization. The errors obtained on the experimental data are 0.068%, 0.049% and 0.153% at 15, 55 and 95 °C respectively.

To describe all of the Na 2 SO 4 -H 2 O binary system's thermodynamic properties, new solubility products of mirabilite (Na 2 SO 4 .10H 2 O) and thenardite (Na 2 SO 4 ) had to be revised in consistency with our new set of interaction parameters (Fig. 5). Experimental solubility data have been reported by several authors and compiled by [START_REF] Linke | Solubilities inorganic and metal-organic compounds[END_REF]. The constants used for the calculation of solubility products according to temperature are given in Table 5.

The MgSO 4 -H 2 O binary system

Fig. 6 presents both the standard heat capacity and the volume of MgSO 4 optimized by Phutela andPitzer (1986a, 1986b) and calculated with the HKF equations [START_REF] Helgeson | Theoretical prediction of the thermodynamic behavior of aqueous electrolytes by high pressures and temperatures; IV, Calculation of activity coefficients, osmotic coefficients, and apparent molal and standard and relative partial molal properties to 600°C[END_REF]; [START_REF] Johnson | SUPCRT92: A software package for calculating the standard molal thermodynamic properties of minerals, gases, aqueous species, and reactions from 1 to 5000 bar[END_REF]). The values of the standard heat capacity obtained by Phutela and Pitzer (1986a) are slightly greater than the values calculated with the HKF equations. For the standard volume, there is no significant difference.

The heat capacity and the density of MgSO 4 solutions were then calculated using the interaction parameters of Phutela andPitzer (1986a, 1986b) and the calculation of standard properties implemented in PhreeSCALE. The interaction parameters (Y) described by Phutela and Pitzer (1986a) were therefore converted to the PhreeSCALE format for comparison (Eq. ( 7)). For the heat capacity, there was a discrepancy between the calculated and experimental values above 50 °C (Fig. 7). Therefore, a new set of interaction parameters was determined using heat capacity and osmotic coefficient data. The new interaction parameters are given in Table 6. The results are plotted in Fig. 7 for the heat capacity (solid lines) and in Fig. 8 for the osmotic coefficient. Data are represented up to saturation with respect to MgSO 4 (s).

For the density, Fig. 9 shows that the parameters of Phutela and Pitzer (1986b) coupled with the standard properties calculated with the HKF equations are in agreement with experimental data. The error obtained on the experimental data is 0.327%, 0.438% and 0.907% at 5, 55 and 95 °C respectively. The solubility of MgSO 4 in H 2 O is plotted in Fig. 10 and can be compared with the experimental data reported by [START_REF] Linke | Solubilities inorganic and metal-organic compounds[END_REF]. The solubility curves were calculated using the solubility products of MgSO 4 .nH 2 O [n = 7, 6 or 1] given in Table 5.

The Na 2 SO 4 -MgSO 4 -H 2 O ternary system

The ternary interaction parameters (𝜃 𝑁𝑎 + /𝑀𝑔 2+ and 𝜓 𝑁𝑎 + /𝑀𝑔 2+ /𝑆𝑂 4 2-) were determined only from the solubility data. A value of 0.07 was chosen for the parameter 𝜃 𝑁𝑎 + /𝑀𝑔 2+ [START_REF] Harvie | The prediction of mineral solubilities in natural waters: The Na-K-Mg-Ca-H-Cl-SO 4 -OH-HCO 3 -CO 3 -CO 2 -H 2 O system to high ionic strengths at 25°C[END_REF]. [START_REF] Pabalan | Thermodynamics of concentrated electrolyte mixtures and the prediction of mineral solubilities to high temperatures for mixtures in the system Na-K-Mg-Cl-SO4-OH-H2O[END_REF] used this value to determine solubility in the NaCl-MgCl 2 -H 2 O ternary system. In order to ensure the consistency of the database, we used their ) were also optimized (Table 5). The solubility diagram in the Na 2 SO 4 -MgSO 4 -H 2 O system was plotted at various temperatures (Fig. 11).

𝜃 𝑁𝑎 + /𝑀𝑔 2+
Solubility data were accurately described up to 75 °C. At 100 °C, however, the invariant points were not correctly described. This discrepancy can be explained by our use of

Conclusion

New equations were derived in the framework of the Pitzer formalism to describe some apparent properties of aqueous solutions, namely enthalpy, heat capacity and density. These equations can theoretically apply to any complex aqueous system since they explicitly include all of the thermal and volumetric specific interactions parameters related to binary and ternary interactions between cations, anions and neutral species. These equations are obtained by temperature or pressure derivatives of the excess Gibbs free energy and are included in the geochemical code PHREEQC V3 [START_REF] Parkhurst | Description of Input and Examples for PHREEQC Version 3-A Computer Program for Speciation, Batch-Reaction, One-Dimensional Transport, and Inverse Geochemical Calculations[END_REF]. A modified version of this new code, PhreeSCALE, was successfully tested on the Na 2 SO 4 -H 2 O, MgSO 4 -H 2 O binary systems and on more complex systems (Na-K-Cl-SO 4 -H 2 O and Na-K-Ca-Mg-Cl-SO 4 -HCO 3 -H 2 O) for which interaction parameters (Y, Y L , Y J and Y V ) and standard properties (𝐶 𝑝,𝑖 0 and 𝑉 𝑖 ° of any compound i) are available in the literature. In this study, the standard properties were calculated directly using PhreeSCALE, which uses the HKF equations [START_REF] Helgeson | Theoretical prediction of the thermodynamic behavior of aqueous electrolytes by high pressures and temperatures; IV, Calculation of activity coefficients, osmotic coefficients, and apparent molal and standard and relative partial molal properties to 600°C[END_REF] instead of the standard properties fitted. Consequently, it can be necessary to reoptimize interaction parameters, which is the case of Na 2 SO 4 -H 2 O and MgSO 4 -H 2 O binary systems. The methodology described for these systems can be applied to many other electrolyte systems if sufficient experimental data are available to enable full parameterization of the systems.

Finally, the main innovation implemented in PhreeSCALE is the calculation of the heat capacity of aqueous solutions using both the Pitzer and the HKF equations. Calculation of the density is also possible with the same formalisms. If the interaction parameters (Y and Y V ) are available PhreeSCALE can be used to calculate directly osmotic coefficient, heat capacity and density; if not they can be optimized by coupling PhreeSCALE and the parameter estimation software PEST [START_REF] Doherty | PEST: Model-Independent Parameter Estimation[END_REF].

Figure captions

ENTHALAPP

The apparent relative molar enthalpy of solution (𝐿 𝜙 ) Eq. ( 6)

CPSTAN("i")

The standard partial molar heat capacity of the aqueous species i at infinite dilution (𝐶 𝑝,𝑖 0 ) HKF model

CPAPP

The apparent molar heat capacity of the solution (𝐶 𝑝,𝜙 ) Eq. ( 15)

CPSOL

The mass heat capacity of the aqueous solution (𝑐 𝑝 ) Eq. ( 18)

VSTAN("i")

The standard partial molar volume of the aqueous species i at infinite dilution (𝑉 𝑖 °)

HKF model

VAPP

The apparent molar volume of solution (𝑉 𝜙 ) Eq. ( 24)

DENSOL

The density of the aqueous solution, ρ Eq. ( 30)

3 𝑉 2 °̅̅̅ , 𝛽 (0)V , 𝛽 (1)V and 𝐶 𝜙𝑉 𝑉 2 °̅̅̅ , 𝛽 (1)V , 𝛽 (2)V and 𝐶 𝜙𝑉

Figure 1 -

 1 Figure 1 -Comparison of the standard properties of Na 2 SO 4 computed by Rogers and Pitzer

Figure 2 -

 2 Figure 2 -Calculated (dotted line for HKF + interaction parameters of Rogers and Pitzer

Figure 4 -

 4 Figure 4 -Calculated (lines) and experimental (open symbol) density for the Na 2 SO 4 -H 2 O

Figure 6 -

 6 Figure 6 -Comparison of the standard properties of MgSO 4 computed by Phutela and Pitzer

Figure 7 -

 7 Figure 7 -Calculated (dotted line for HKF + interaction parameters of Phutela and Pitzer

Figure 9 -

 9 Figure 9 -Calculated (lines) and experimental (open symbol) density for the MgSO 4 -H 2 O

Figure 10 -

 10 Figure 10 -Calculated (lines) and experimental (open symbols) solubility of epsomite

Figure 11 -

 11 Figure 11 -Calculated (lines) and experimental (open symbols) solubility in the MgSO 4 -

  in order to optimize new Pitzer interaction parameters from heat capacity and osmotic coefficient data, or density data. This is on the case of the Na 2 SO 4 -H 2 O and MgSO 4 -H 2 O binary systems, for which we propose a new set of interaction parameters.

  values. Therefore, only 𝜓 𝑁𝑎 + /𝑀𝑔 2+ /𝑆𝑂 4 2-values were optimized and the best value was -0.0181193. The solubility constant of the double salts (Na 2 SO 4 .MgSO 4 .4H 2 O, Na 2 SO 4 .MgSO 4 .2.5H 2 O or 3Na 2 SO 4 .MgSO 4

  doi:10.1039/jr9620000937 Tanger, J.C.,Helgeson, H.C., 1988. Calculation of the thermodynamic and transport properties of aqueous species at high pressures and temperatures; revised equations of state for the standard partial molal properties of ions and electrolytes. American Journal ofScience 288, 19-98. doi:10.2475/ajs.288.1.19

	Thomsen, K., 1997. Aqueous Electrolytes Model Parameters and Process Simulation.
	Technical University of Denmark.
	Wang, P., Anderko, A., Young, R.D., 2002. A speciation-based model for mixed-solvent
	electrolyte systems. Fluid Phase Equilibria 203, 141-176. doi:10.1016/S0378-
	3812(02)00178-4
	Zaytsev, I.D., Aseyev, G.G., 1992. Properties of Aqueous Solutions of Electrolytes. CRC
	Press.

Table 1 .

 1 The absolute average relative deviation for the different systems tested 3

	Keywords

Table 2 .

 2 The absolute average relative deviation for the different systems tested 5

	Systems	Ref.	Parameters used	T	Pressure	Data	∆ 𝑨𝑨𝑫 (%)
				°C
					1 bar or	𝐶 𝑝,𝜙	0.362
					saturated
	Na 2 SO 4 -H 2 O	Rogers and Pitzer, 1981	𝐶 𝑝2 °̅̅̅̅̅ , 𝛽 (0) , 𝛽 (1) and 𝐶 𝜙	25-200	vapor	𝐿 𝜙	0.495
					pressure	𝜙	0.0392
						𝐶 𝑝,𝜙	1.955
	MgSO 4 -H 2 O	Phutela and Pitzer, 1986a	𝐶 𝑝2 °̅̅̅̅̅ , 𝛽 (0) , 𝛽 (1) , 𝛽 (2) and 𝐶 𝜙	25-200	20 bar	𝐿 𝜙	0.0970
						𝜙	0.118
	Na 2 SO 4 -H 2 O
		Phutela and Pitzer, 1986b
	MgSO 4 -H 2 O

Table 3 .

 3 Experimental data used for the optimization of Na 2 SO 4 -H 2 O binary system at saturation pressure 7

	Type of data	T (°C)	m (mol.kg -1 )	Experimental error	Ref.
		25	0.1-4		Robinson et al., 1941
		45	0.76-3.45		Hellams et al., 1965
		60	0.98-3.31		Humphries et al., 1968
		80	0.7-3		Moore et al., 1972
		0	0.1-1.5		Platford, 1973
		25	0.2-1.2		Downes and Pitzer, 1976
	Osmotic coefficient	25	0.46-3.7		Rard and Miller, 1981
		50-250	0.29-3.2		Bhatnagar and Campbell, 1982, 1981
		109.99-225.04	0.5-3		Holmes and Mesmer, 1986
		40	0.2-3.26		Baabor et al., 1998
		25	0.1-3		El Guendouzi et al., 2003
		25; 50	0.1-3.7		Rard et al., 2000
		25-80	0.1-3.4		El Guendouzi and Aboufaris, 2015
		80-180	0.14-0.772	±0.003 cal.g -1 .K -l	Likke and Bromley, 1973
	Heat capacity	25	0.01-1.89		Magalhães, 2002
		50.85-140.75	0-2.63	±0.002 J.g -1 .K -l	Rogers and Pitzer, 1981
		0-90	0.14-1.76	-	Zaytsev and Aseyev, 1992
	8

Table 4 .

 4 Values of the fitting constants for the interaction parameters of Na 2 SO 4 -H 2 O binary system (after optimization) 9

		a 1 (-)	a 2 (T)	a 3 (T²)	a 5 (1/T)	a 6 (lnT)	a 7 (1/(T-263)) a 8 (1/(680-T))
	β 0 (Na + /SO 4 2-)	107.35561	0.04794672	-0.00002034	-2736.492	-19.41903	-0.1243008	0
	β 1 (Na + /SO 4 2-)	1017.4354	0.45031511	-0.00019494	-26137.553	-183.90464	0.00011838	851.963305
	C φ (Na + /SO 4 2-)	-0.13574678	0.00013002	0	30.0944336	0	0.03864754	0
	10
	11

Table 5 .

 5 Values of fitting constants for solubility products 12

		a 1 (-)	a 2 (T)	a 3 (1/T)	a 4 (log T)	a 5 (1/T²)
	Na 2 SO 4 .10H 2 O
	(mirabilite)	-4.86996660E+03 -1.50030705E+00 1.19561391E+05 1.98633627E+03	0
		-3.50862995E+03 -5.18430862E-01 1.99075358E+05 1.26574808E+03	-1.21641701E+07
	Na 2 SO 4 (thenardite)
	MgSO 4 .7H 2 O
		4.08480390E+04	7.49950739E+00 -1.934676E+06 -1.521724E+04	9.39747409E+07
	(epsomite)
	MgSO 4 .6H 2 O
	(hexahydrite)	5.43604369E+03	5.61493175E-01 -3.5461667E+05 -1.8930412E+03	2.38767356E+07
	MgSO 4 .H 2 O (kieserite) 5.39329145E+03	5.18715131E-01 -3.7474065E+05 -1.8602601E+03	2.78132757E+07

Na 2 SO 4 .MgSO 4 .4H 2 O

	(bloedite)	3.56483519E+01 -6.58258168E-02 -5.4674363E+03	0	0
	3Na 2 SO 4 .MgSO 4	1.09286124E+01 -4.26072730E-02	0	0	0
	Na 2 SO 4 .MgSO 4 .2.5H 2 O 6.56596356E+00 -2.73680649E-02	0	0	0
	13
		7

Table 6 .

 6 Experimental data used for the optimization of MgSO 4 -H 2 O binary system at saturation pressure 14

	Type of data	T (°C)	m (mol.kg -1 )	Experimental error	Ref.
		99.6	1.9-4.75		Patterson et al., 1960
		121.1	1.56-5.1		Soldano and Patterson, 1962
		0	0.1-3		Platford, 1973
	Osmotic coefficient	25	0.925-3.61		Rard and Miller, 1981
		109.99, 140.07	1.38-5		Holmes and Mesmer, 1983
		25	0.5-1.3		Archer and Rard, 1998
		25	0.2-3		El Guendouzi et al., 2003
	Heat capacity	80-200	0.18-0.92	±0.003 cal.g -1 .K -1	Likke and Bromley, 1973
		0-90	0.17-1.35	-	Zaytsev and Aseyev, 1992
	15

Table 7 .

 7 Values of the fitting constants for the interaction parameters of MgSO 4 -H 2 O binary system (after optimization) 16

		a 1 (-)	a 2 (T)	a 3 (T²)	a 5 (1/T)	a 6 (lnT)	a 7 (1/(T-263)) a 8 (1/(680-T))
	β 0 (Mg 2+ /SO 4 2-)	1.0245267E+02	4.743914E-02	-2.115E-05	-2.432276E+03 -1.873061E+01 -8.441E-07	1.444294E+02
	β 1 (Mg 2+ /SO 4 2-)	2.2435909E+02	-9.887453E-02	1.43701E-04	-1.049481E+04 -3.016613E+01 -3.5495203	1.1050714E+03
	β 2 (Mg 2+ /SO 4 2-)	7.58048471E+04	6.53134854E+01	-4.388932E-02 -1.05904E+06 -1.544154E+04 -1.36617E-01 4.6838757E+04
	C φ (Mg 2+ /SO 4 2-)	-9.8373402E+01	-7.007396E-02	4.58738E-05	1.7915821E+03 1.93122662E+01 9.40941E-06 -3.159734E+02

Acknowledgements

This work was funded by BRGM (the French Geological Survey). The authors thank the LABEX Voltaire (ANR-10-LABX-100-01).

temperature-independent ternary interaction parameters in order to be coherent with the authors cited above.


A Bourgeat
email: alain.bourgeat@cdcsp.univ-lyon1.fr

M Kern

S Schumacher
email: stephan.schumacher@andra.fr

The COUPLEX Test Cases: Nuclear Waste Disposal Simulation

Keywords: Nuclear waste, numerical simulation, porous media

The models appearing in the COUPLEX benchmark are a set of simplified albeit realistic test cases aimed at simulating the transport of radionuclides around a nuclear waste repository. Three different models were used:

The first test case is related to simulations based on a simplified 2D far-field model close to those used for safety assessments in nuclear waste management. It leads to a classical convection diffusion type problem, but with highly variable parameters in space, highly concentrated sources in space and time, very different time scales and accurate results expected even after millions of years.

The second test case is a simplification of a typical 3D near-field computation, taking into account the glass dissolution of vitrified waste, and the congruent release of several radionuclides (including daughter products), with their migration through the geological barrier.

The aim of the third test case is to use the results of the near-field computation (COUPLEX 2) to drive the behavior of the nuclide source term in the Far Field computation (COUPLEX 1). The modeling of this last case was purposely left rather open, unlike the previous two, leaving the choice to participants of the way the coupling should be made.

Introduction

Deep disposal of nuclear waste raises a number of challenges for numerical simulations : widely differing lengths and time-scales, highly variable physical coefficients and stringent accuracy requirements. It is no wonder therefore if the topic is catching the attention of computational scientists and numerical analysts.

The COUPLEX models are a set of three different test cases, of increasing difficulty.

-The first test case, COUPLEX 1, is a simplified far-field computation. We have chosen a 2D model, with four geological layers, in order to expose the main numerical difficulties: a highly heterogeneous medium (hydraulic conductivities vary over six orders of magnitude), and a highly elongated medium (width-to-height ratio is over 30). In this model, the repository is a simple rectangle acting as a distributed source with known time dependency.

-The second test case is a near-field computation, the aim of which is to provide a better description of the source term in the far-field model. We have incorporated more physics than in the previous model, and a 3D geometry has been used, thus ending up with a more difficult case.

-The last test case was deliberately left more open than the previous two.

Its aim is to couple the near field and the far field, that is to use the results of COUPLEX 2 to drive a model similar to COUPLEX 1, but this time in 3D.

We emphasize once more that all three test cases are simplified with respect to a real safety simulation, in order to concentrate on the numerical difficulties. In all three cases, they derive their approach from the models used by Andra.

The rest of this paper deals with a definition of some general notation, followed by a detailed description of each each test case in turn. We show the geometry, describe in detail the model-defining equations, and provide the physical parameters to be used.

GENERAL TRANSPORT OPERATOR

Each of the physical models we describe is a variant of a convection-diffusionreaction equation. In order to shorten the equations, we found it useful to define a "general" transport operator as follows, in some spatial domain R * :

L(C) = -∇. D∇C -uC in R * .
(1)

The diffusion-dispersion tensor D is a general anisotropic tensor (see Bear [START_REF] Bear | Introduction to modeling of transport phenomena in porous media[END_REF][START_REF] Bear | 2 of Theory and Applications of Transport in Porous Media[END_REF] ), that is given as a function of the Darcy velocity u by:

D = d e I + |u| [α l E(u) + α t (I -E(u))] (2) with E i j (u) = u i u j |u| 2 . ( 3
)
For any species or isotope, d e (m 2 /s) is the effective diffusion coefficient, α l and α t , the longitudinal and transverse dispersion coefficients (both in metres).

couplex.tex; 3/03/2004; 16:37; p.2

COUPLEX 1: Far-field Computation

INTRODUCTION

The repository lies at a depth of 450 m inside a clay layer located between a layer of limestone and a layer of marl above and a layer of Dogger limestone underneath. Water flows slowly (creeping flow) through those porous media and convects the radioactive materials once the containers start to leak. There is also a dilution effect which, in mathematical terms, is similar to diffusion. The problem raises three main difficulties:

1. The geometry of the domain is highly elongated with a width of 25 km and a depth of 700 m.

2. The convection and diffusion constants are very different from one layer to another. In the clay layer, for instance, there is almost no convection while, in the other layers, diffusion and convection are both important.

3. Radioelements leak from containers into the clay within a short period compared with the millions of years over which convection and diffusion remain active.

GEOMETRY

In this first test case, the computation is restricted to a 2D section of the disposal site. Thus, the computational domain lies within a rectangle O = (0, 25000) × (0, 695) in metres. The Dogger, clay, limestone, and marl layers are located as follows (starting from the bottom left corner of the rectangle):

-The Dogger constitutes the bottom layer: 0 < z < 200;

-The clay lies between the horizontal line z = 200 and the line from (0, 295) to (25000, 350);

-The limestone lies between the line from (0, 295) to (25000, 350) and the horizontal line z = 595;

-The marl is near the surface : 595 < z < 695.

The repository, R , is modeled by a uniform rectangular source in the clay layer:

R = {(x, z) ∈ (18440, 21680) × (244, 250)}
The geometry is summarized on Figure 1. For that domain the computation should be carried for t ∈ (0, T ) with T = 10 7 years.

FLOW MODEL

It is assumed that all rock layers are saturated with water and that boundary loads are stationary so that the flow is independent of time. Darcy's Law gives the velocity u in terms of the hydro-dynamic load H = P/ρg + z:

u = -K∇H (4)
where the permeability tensor, K, assumed constant in each layer is given in Table I, P is the pressure and g is Newton's constant. Conservation of mass (∇ • (ρu) = 0, with density, ρ, assumed constant) implies that

∇ • (K∇H) = 0 in O (5)

RADIONUCLIDE MIGRATION

We consider two species of particular interest, iodine 129 and plutonium 242.

Both escape from the repository cave into the water and their concentrations C i , i = 1, 2 are given by two independent convection-diffusion equations:

R i ω ∂C i ∂t + L(C i ) + R i ωλ i C i = f i in O × (0, T ) i = 1, 2. ( 6
)
where L is the general transport operator, as defined in section 1.1, and -R i , the latency retardation factor, with value 1 for 129 I, 10 5 for 242 Pu in the clay and 1 elsewhere for both iodine and plutonium;

-The effective porosity, ω, is equal to 0.001 for 129 I, 0.2 for 242 Pu in the clay layer and 0.1 elsewhere for both;

λ i = log 2/T i where T i is the half-life of the element : T 1 = 1.57 10 7 for 129 I, T 2 = 3.76 10 5 for 242 Pu (in years);

-The effective diffusion/dispersion tensors, D i , have been defined in ( 2) and (3) above. The coefficients are assumed constant in each layer and are given in Table II.

In this test case, the values of the source terms, f i , (i = 1, 2), in the repository R are given in tabulated form in separate data files. The source terms are assumed to be spatially uniformly spread out in all the repository R . It is assumed that there is no source outside the repository (

f i , (i = 1, 2) in O\ R ).
The dependence in time is shown in Figure 2. The structure of the data file is described in Appendix A.

Table II. Dispersion coefficients for the radioelements in the four layers

Iodine Plutonium d e (m 2 /year) α L (m) α T (m) d e (m 2 /year) α L (m) α T (
∂C i ∂n = 0 on {0} × (295, 595) left limestone ∂C i ∂n = 0 on {0} × (0, 200) left Dogger D i ∇C i • n -C i u • n = 0 on (0, 25000) × {0} bottom C i = 0 elsewhere on the boundary.
where n is the outward normal to O.

OUTPUT REQUIREMENTS

The following output quantities are expected from the simulations(both tables and graphical representations): couplex.tex; 3/03/2004; 16:37; p.6 -Contour levels of C i at times 200, 10110, 50110, 10 6 , 10 7 years (the following levels should be used: 10 -12 , 10 -10 , 10 -8 , 10 -6 , 10 -4 );

-Pressure field (10 values uniformly distributed between 180 and 340);

-Darcy velocity field, along the three vertical lines given by x = 50, x = 12500, x = 20000, using 100 points along each line;

-Places where the Darcy velocity is zero;

-Cumulative total fluxes through the top and the bottom clay layer boundaries, as a function of time;

-Cumulative total fluxes through the left boundaries of the Dogger and limestone layers;

-The discretization grid of the domains and the time stepping used in the simulations should also be provided.

COUPLEX 2: Near-field Computation

INTRODUCTION

A real repository is made up of a very large number of disposal modules, each of which includes a large number of cells, linked by a backfilled drift. Since modeling all the details of such a repository still presents difficulties, we follow a usual approach in safety computations: we only consider an "elementary cell", and we add periodicity conditions on the boundaries. That last assumption is valid provided there is a large enough number of similar cells in a disposal module, and that the module boundary layer effects are negligible. The computation in this test case is based on the modeling of three cells linked by a backfilled drift embedded into a homogeneous medium. The dissolution of glass, the ensuing congruent emission of radionuclides, and precipitation of elements with several isotopes are included in the model.

GEOMETRY AND FLOW MODEL

We choose the geometry of the elementary cell to be simple and the host geological medium to be simply made up of the clay layer into which the repository is embedded.

We describe below the elementary cell by successive zooms. The domain ("elementary cell") is a parallelepiped R * , containing three cells and a backfilled drift. Each cell is made up of the container (in this case, a glass parallelepiped, containing the radioactive waste), a buffer and a seal. We assume once again that the flow follows Darcy's Law, so that the piezometric head and the Darcy velocity satisfy Equations ( 4) and ( 5) above.

The boundary conditions are -H is given on both horizontal sides of R * . The difference in piezometric head is taken equal to 50 m;

-Since the container is impervious, the boundary condition on the container boundary is no-flow: ∂H ∂n = 0 there;

-H is supposed to be periodic on the four vertical sides of R * (in order to simulate an entire disposal module). Silica migration is described by a convection-diffusion-reaction type equation:

φR s ∂C s ∂t + L s (C s ) = ρ p ν p λ p 1 - C s S p in R * \ C , ( 7
)
where φ is the porosity, and the retardation factor R s is given by the following equation, (ρ s being the density of the solid phase and the distribution coefficient, K d s , taking adsorption into account in a simplified way):

R s = 1 + ρ s 1 -φ φ K d s .
We have denoted by couplex.tex; 3/03/2004; 16:37; p.9

ρ p : the precipitate concentration (4.2 10 4 mole/m 3 ); -ν p : the precipitation speed (6 10 -9 m/year);

λ p : the inverse of specific surface (between 8 10 -10 m and 10 -3 m);

-S p : the solubility of the precipitated phase (5.4 10 -1 mol/m 3 ).

Silica precipitation follows a first order kinetic law.

We propose two different models for glass dissolution. A detailed model of glass dissolution follows a first order kinetic law. The complete free boundary model associated with the dissolution of the container C is simplified into a Fourier boundary condition on ∂C , the glass-bentonite interface assuming that this interface remains fixed:

-D s ∂C s ∂n = ρ m ν m 1 - C s S m , on ∂C ; ( 8
)
where we have set:

ν m to be the initial dissolution velocity (1.5 10 -5 m/year);

ρ m , the available concentration that may be dissolved in the matrix (4.2 10 3 mol/m 3 );

-S m , the maximum silica concentration (or saturated silica concentration), (8.2 10 -1 mole/m 3 ).

Alternatively, a simplified model of glass dissolution could be used, reduced to a known constant flux on the glass-bentonite interface:

-D s ∂C s ∂n = τ, on ∂C ; (9)
with the release factor τ = 10 -2 mol/m 2 /year. It should be noted that this simplification will in turn make the congruent emission of radionuclides simpler in Equation ( 19) in Section 3.4 below.

The initial silica concentration is considered to be equal to its reference value in the geologic medium, that is 9.9 10 -5 mole/l. The boundary and interface conditions are :

-On the container boundaries, the condition is that given in Equation ( 8) above;

-C s = 0 on the top and bottom boundaries of R * ;

-Periodic boundary conditions on the four outer vertical sides of R * \ C . To make writing the equations easier, we associate to each radionuclide an index k, which will take the values k = s, 0, 1, 2, -1, -2, as shown in Table IV. For k = 1 or k = 2, element k is by definition the parent of element -k. Isotope -1 and -2 correspond to the same chemical element, as isotopes 1 and 2 do.

Each radionuclide satisfies a convection-diffusion-reaction equation analogous to Equation ( 7), obviously with an added radioactive decay term. It may possibly precipitate.

Caesium

As far as caesium (index k = 0) is concerned, we shall assume that adsorption follows a Langmuir isotherm, and that the solubility limit is so high that caesium does not precipitate. Following Equation (2), the corresponding transport equation becomes:

φ ∂ ∂t R 0 C 0 + L 0 (C 0 ) + φR 0 λ 0 C 0 = 0 in R * \ C ( 10
)
where the retardation factor now depends on the concentration: tex; 3/03/2004; 16:37; p.11 through the (concentration dependent) distribution coefficient:

R 0 = 1 + ρ s 1 -φ φ K d 0 (11) couplex.
K d 0 = κ CEC 1 + κC 0 . ( 12
)
In the previous equation, κ = 150 (mol/l) -1 and CEC is the cationic exchange capacity (CEC = 80 mmol/100g).

Plutonium and uranium

On the contrary, since Plutonium and Uranium may precipitate, we introduce the concentration of the solid phase, denoted by F k , k = -2, -1, 1, 2. The transport of each isotope is described by the following equations:

-For k = 1 and 2 (plutonium):

ρφR k ∂C k ∂t + ρL k (C k ) + ρφR k λ k C k = ρφS k in R * \ C (13) ρ s (1 -φ) ∂F k ∂t + ρ s (1 -φ)λ k F k = -φS k in R * \ C ( 14
)
-For k = -1 and -2 (uranium):

ρφR k ∂C k ∂t + ρL k (C k ) + ρφR k λ k C k -ρφR -k λ -k C -k = = ρφS k in R * \ C (15) ρ s (1 -φ) ∂F k ∂t + ρ s (1 -φ)λ k F k -ρ s (1 -φ)λ -k F -k = -φS k in R * \ C (16)
In Equations ( 13) to ( 16), the retardation factor, R k , is still given by an isotherm similar to Equation (11):

R k = 1 + ρ s 1 -φ φ K d k but with a constant coefficient, K d k .
In Equations ( 13) to ( 16), the dissolution-precipitation term, S k , takes into account the precipitation per species, by assuming that each isotope precipitates proportionally to its mass concentration. The dissolution-precipitation term is given by:

S k = σ(C k -C sat k ) if C sat k ≤ C k 0 if C k ≤ C sat k and F k = 0 (17)
couplex.tex; 3/03/2004; 16:37; p.12 with:

C sat k = C sat e(k) φR k C k + ρ s (1 -φ)F k ∑ |k |=1,2 k has same sign as k (φR k C k + ρ s (1 -φ)F k ) , (18)
where C sat e is the solubility limit of the considered species, e ∈ {Pu, U}, with e(1) = e(2) = Pu, e(-1) = e(-2) = U, and σ = 10 4 year -1 .

The dissolution of glass generates a congruent emission of the radionuclides still present in the waste. Taking into account the radioactive decay of the radionuclides leads to a boundary condition with time-dependent coefficients:

D k ∂C k ∂n = N k (t) N 0 s D s ∂C s ∂n , k = -1, 1, 0, 1, 2 (19)
on the glass-bentonite interface ∂ C (supposed fixed as in Equation ( 8)). The constant N 0 s is the initial number of moles of silica in the container, and is equal to N 0 s = 4.56 10 4 (the silica density in the container is equal to the number ρ m which appears in Equation ( 8)).

It should be noted that in case the simplified model of glass dissolution given in Equation ( 9) is used, Equation ( 19) is simplified to:

D k ∂C k ∂n = τ N k (t) N 0 s , k = -1, 1, 0, 1, 2 (20)
with τ = 10 -2 mol/m 2 /year.

The functions N k (t) are the number of moles of the corresponding element at time t. They are solutions of a coupled set of differential equations expressing the decay of the element under consideration. A simple computation gives the expression for N k (t):

N k (t) = N 0 k e -λ k t , k = 0, 1, 2 (21) N k (t) = N 0 -k λ -k λ k -λ -k e -λ -k t + N 0 k + N 0 -k λ -k λ -k -λ k e -λ k t , k = -1, -2 (22)
The constants N 0 k (initial values of the number of mole for each radionuclide) are given in Table V below. They were computed from the (known) activities for each radionuclide, as detailed in Appendix B. Those numbers are valid for one cell.

Initial and boundary conditions

At the beginning of the simulation, radionuclide concentration is zero everywhere. Boundary conditions for the radionuclides are as follows:

-As given in Equation ( 19) above on the container boundary ∂C ;

-C k = 0, k = -2, -1, 1, 2 on the top and bottom sides of the elementary cell R * ;

-Periodic boundary conditions on the four outer vertical sides of R * .

The simulation should cover a period of 10 6 years.

OUTPUT REQUIREMENTS

The following output quantities are expected from the simulations (both tables and graphical representations):

-Total fluxes for all radionuclides on the top and bottom horizontal boundaries of R * , as a function of time;

-Total fluxes through the seal -backfilled drift interfaces, as a function of time;

-Total fluxes through the buffer -seal interfaces, as a function of time;

-Total fluxes through the backfilled drift -geological medium interfaces, as a function of time;

-Concentration profiles for all radionuclides along both symmetry planes, at times 10 2 , 10 3 , 5 × 10 2 , 10 4 , 10 5 , 10 6 years;

-Total cumulative concentration of dissolved silica as a function of time;

-Exact time at which the total initial amount of silica has been dissolved;

-Total cumulative precipitation for each element and isotope, as a function of time. It should be noted that the buffer and the seal are actually the same material.

COUPLEX3: Coupling Near and Far-fields

MOTIVATIONS

The aim of that final test case is to use the results of the near-field computation (COUPLEX 2) to describe the behavior of the radionuclide source term, f i , in the far-field computation (COUPLEX1, Equation ( 6)). In this way, we could include as much as possible of the repository behavior, coming from a near-field computation, inside the radionuclide source term f i necessary for a far-field simulation.

The coupling between the two simulations should then allow possible near-field / far-field feedbacks.

This test case was more open than the previous two, and the way the coupling should be carried out was not specified beforehand, and was left up to the participants.

The results of the COUPLEX2 Test Case on a cell R * (dimension 27 m× 100 m× 49.6 m) have to be used to define the source behavior of the 2D repository R (dimension 3240 m× 6 m), as in figure 1 in COUPLEX1 Section 2.2, and Equation ( 6). The simulation period is T = 10 7 years, as in COUPLEX1.

The cell R * in COUPLEX2 was actually part of a larger repository and we assume for simplicity that the whole repository is made up of 10 parallel lines of containers. Each line of containers is made by repeating 120 times the backfilled drift / cells system described in Figure 4 in COUPLEX 2. The length of a line is then 3240 m. All the lines belong to the same horizontal plane. The distance between two lines is 50 m. The geometry is shown in Figure 6.

The radio-elements to be considered in the repository for the far-field computation are now those used in COUPLEX2, that is an isolated 135 Cs and two radionuclide chains of 238 Pu and 242 Pu, as described in Section 3.4 (COUPLEX 2), following transport Equations ( 7) and ( 13)-( 16) (Sections 3.3 and 3.4 in COUPLEX 2).

We are interested in the far-field simulation of a 2D vertical cross-section of the disposal site, similar to the one in the simulation of COUPLEX1. This cross-section is parallel to the lines of containers.

Data in Table II for COUPLEX 1 have to be completed by the ones in Tables VI-X for COUPLEX 2.

The host geological medium then becomes the clay layer defined for COU-PLEX 1, but with the same porosity of 15% for all radionuclides. Other char-

OUTPUT REQUIREMENTS

-At least one simulation using the same kind of sequential coupling as for COUPLEX 1 (the repository behavior is a simple input to the far-field simulation) should be carried out, and the results should be presented following the requirements mentioned in Section 2.6.

-If another coupling method has been chosen, the results should still be presented according to the requirements metioned in Section 2.6.

Figure 1 .

 1 Figure 1. Geometry of computational domain

Figure 2 .

 2 Figure 2. Release of iodine and plutonium as a function of time

Figure 3 .

 3 Figure 3. Details of the simulation domain R * , (left, front and top views) including the host geological medium, the backfilled drift and the cell

Figure 4 .Figure 5 .

 45 Figure 4. Backfilled drift and three cells included in the host geological medium

Figure 6 .

 6 Figure 6. Geometry of the disposal site

Table I .

 I Permeability tensor in the four rock layers

	Boundary conditions are as follows:
	H = 289	on {25000} × (0, 200),	(right Dogger)
	H = 310	on {25000} × (350, 595), (right limestone)
	H = 180 + 160x/25000 on (0, 25000) × {695},	(top)
	H = 200	on {0} × (295, 595),	(left limestone)
	H = 286	on {0} × (0, 200)	(left Dogger),
	∂H ∂n	= 0	elsewhere.
			Marl	Limestone Clay	Dogger
			K (m/year) 3.1536 10 -5 6.3072	3.1536 10 -6 25.2288
				couplex.tex; 3/03/2004; 16:37; p.4

Table III .

 III We take into account one isolated radionuclide, caesium 135, and two radionuclide chains: 238 Pu → 234 U and 242 Pu → 238 U. Those choices allows us to study the precipitation of each element. The half-lives involved are: Radionuclide half-lives

	3.4. RADIONUCLIDE MIGRATION
	Isotope	135 Cs	238 Pu	234 U	242 Pu	238 U
	Period (years) 2.3 10 6 87.74 2.45 10 5 3.76 10 5 4.47 10 9
			couplex.tex; 3/03/2004; 16:37; p.10

Table IV .

 IV Radionuclide-Index association isotope 135 Cs 238 Pu 234 U 242 Pu 238 U

	index k	0	1	-1	2	-2

Table V .

 V Initial number of moles for the five radionuclides in one cell

	Isotope	135 Cs	238 Pu	234 U	242 Pu	238 U
	Number of moles	41.6	2.01 10 -2 1.40 10 -2 2.03 10 -2	21.1

Table VIII .

 VIII Physical properties for caesiumd e (m 2 /s) α l (m) α t (m) K d (m 3 /kg)

	Table VI. General physical properties 3.6. PHYSICAL COEFFICIENTS
					K (m/s) φ (%)	ρ s
			buffer		10 -13	25	1.8
			backfilled drift 10 -9	40	2.1
			clay		10 -13	15	2.1
	Table VII. Physical properties for silica
			d e (m 2 /s) α l (m) α t (m) K d (m 3 /kg)
	buffer	5 10 -11	0	0	10 -4
	backfilled drift	2 10 -10	0	0	10 -4
	clay		5 10 -11	0	0	10 -4
	buffer	5 10 -11	0	0	cf. (12)
	backfilled drift	2 10 -10	5	1	10 -3
	clay		5 10 -11	5	1	10 -3
	Table IX. Physical properties for plutonium
	d e (m 2 /s) α l (m) α t (m) K d (m 3 /kg) C sat k (mol/l)
	buffer	5 10 -11	0	0	3	4 10 -9
	backfilled drift	2 10 -10	5	1	8.5 10 -3	4 10 -9
	clay	5 10 -11	5	1	8.5 10 -3	4 10 -9
	Table X. Physical properties for uranium
	d e (m 2 /s) α l (m) α t (m) K d (m 3 /kg) C sat k (mol/l)
	buffer	5 10 -11	0	0	3	2 10 -9
	backfilled drift	2 10 -10	5	1	1.7 10 -2	2 10 -9
	clay	5 10 -11	5	1	1.7 10 -2	2 10 -9
			couplex.tex; 3/03/2004; 16:37; p.14 couplex.tex; 3/03/2004; 16:37; p.15

couplex.tex; 3/03/2004; 16:37; p.3

couplex.tex; 3/03/2004; 16:37; p.5

couplex.tex; 3/03/2004; 16:37; p.7

couplex.tex; 3/03/2004; 16:37; p.13

couplex.tex; 3/03/2004; 16:37; p.17

Acknowledgements

We wish to thank all of the teams who took part in the test cases for their very active participations. P. Bastian (IWR, Heidelberg), Ph. Ackerer, R. Mosé and H. Hotteit (IMFS, Strasbourg) made several remarks that helped to clarify the models as presented here.

The original version of this document was an "exercise" which was submitted to the scientific community by Andra (Agence Nationale pour la Gestion des Déchets Radioactifs -French National Radioactive Waste Management Agency) in January 2001, and lasted throughout that year. Results were

Appendix

A. Description of the data file for COUPLEX 1

The file source.dat contains data needed to compute the source term f i in Equation ( 6). These data come from a near-field computation. The file has 212 lines, and each line contains three numbers t p , f p 1 , f p 2 , p = 1, . . . 212, where t p is the time, and the source term f i (t p ) is related to f p i by: f i (t p ) = f p i /S, where S is the surface of the repository.

The times t p are in years, and the numbers f p i are in units of moles /year.

B. Conversion between activities and number of moles

The activity (in becquerels) of a radionuclide is the number of disintegrations per second. The rule for converting the activity to a number of moles is simply:

T being the half-life of the element (in seconds), and N , Avogrado's number (N = 6.02 10 23 mol -1 ).

The initial activities for the various radionuclides are shown in Table XI. Numbers in the table apply to one package, and there are eight packages per cell.


S M Kostritskii
email: skostritskii@optolink.ru

Yu N Korkishko

V A Fedorov

V P Mitrokhin

O G Sevostyanov

I M Chirkova

O Stepanenko

M De Micheli

Subsurface Disorder and Electro-Optical Properties of Proton-Exchanged LiNbO 3 Waveguides Produced by Different Techniques

Keywords: optical waveguides, electro-optic effect, lithium niobate, proton exchange

It has been established, that proton-exchanged LiNbO 3 waveguides have a marked subsurface layer with structural disorder inducing degradation of electro-optical properties of these waveguides. At the same time, such a subsurface disorder is found to be less pronounced in soft proton-exchanged (SPE) waveguides in comparison with annealed proton-exchanged (APE) ones. The experimental samples of phase modulators fabricated by SPE technique exhibit a better electro-optical efficiency compared to the LiNbO 3 modulators produced by the standard and improved APE techniques.

INTRODUCTION

Annealed proton exchange (APE) is now an established technique for fabricating integrated optical devices in LiNbO 3 , as it allows obtaining low-loss waveguides with good electro-optical (EO) performances [START_REF] Rottschalk | Electrooptic behaviour of proton-exchanged LiNbO 3 optical waveguides[END_REF]- [START_REF] Narayan | Electrooptic coefficient variation in proton-exchanged and annealed lithium niobate samples[END_REF]. Note that a significant degradation (more than one order of magnitude) of the EO coefficients occurs immediately after the proton exchange process [START_REF] Rottschalk | Electrooptic behaviour of proton-exchanged LiNbO 3 optical waveguides[END_REF][START_REF] Loni | Very low-loss proton-exchanged waveguides with a substantially restored electro-optic effect[END_REF]. However, after an appropriate annealing in order to obtain LiNbO 3 waveguides in the so-called α-phase [START_REF] Loni | Very low-loss proton-exchanged waveguides with a substantially restored electro-optic effect[END_REF]- [START_REF] Narayan | Electrooptic coefficient variation in proton-exchanged and annealed lithium niobate samples[END_REF], the EO coefficients have almost recovered their initial values. Nevertheless, EO coefficients cannot be completely restored even after a long annealing, when the initial proton exchange has been performed using conditions that are in common use today [START_REF] Rottschalk | Electrooptic behaviour of proton-exchanged LiNbO 3 optical waveguides[END_REF]- [START_REF] Narayan | Electrooptic coefficient variation in proton-exchanged and annealed lithium niobate samples[END_REF]. It was shown [START_REF] Muller | Reduction of lattice defects in proton-exchanged lithium niobate waveguides[END_REF] that an extended annealing in dry air at T ≥ 320 o C, leads to the precipitation of disordered phases in a subsurface layer (depth ≤ 0.5 µm), causing undesirable degradation of the EO properties of the waveguide. Such degradation could be interpreted as the fieldscreening effect in the highly structurally disordered parts of the proton exchanged regions below the surface. Besides, it was recently reported [START_REF] Kostritskii | Leakage of a guided mode caused by static and light-induced inhomogeneities in channel HTPE-LiNbO 3 waveguides[END_REF] that this subsurface disorder provides some marked level of propagation loss related mainly to guided mode scattering within inhomogeneous near-surface part of the waveguide, and also, significant loss of optical nonlinearity [START_REF] Laurell | Loss of optical nonlinearity in proton--exchanged LiNbO 3 waveguides[END_REF]. It has been established [START_REF] Muller | Reduction of lattice defects in proton-exchanged lithium niobate waveguides[END_REF] that the loss of water during annealing is a major contributor to the formation of crystalline disorder and, therefore, annealing under a high pressure of water vapour prevents such a disorder to appear. However, using this method only a partial suppression of the subsurface disorder and the related degradation of waveguide properties were obtained [START_REF] Kostritskii | Electro-Optical Properties of Different H x Li 1-x NbO 3 Phases in Proton-Exchanged LiNbO 3 Waveguides[END_REF]. In this paper we show that, using the soft proton exchange (SPE) process [START_REF] Chanvillard | Soft Proton Exchange on PPLN: a simple waveguide fabrication process for highly efficient nonlinear interactions[END_REF][START_REF] Castaldini | Soft proton exchange tapers for low-insertion loss LiNbO 3 devices[END_REF][START_REF] Mp De | Fabrication and Characterization of Proton Exchanged Waveguides in Periodically Poled Congruent Lithium Niobate[END_REF] recognised as allowing the fabrication of low loss and highly efficient nonlinear components, we obtain waveguides in the α-phase [START_REF] Micheli | Crystalline and optical quality of proton exchanged waveguides[END_REF][START_REF] Yu | Relationships between structural and optical properties of proton-exchanged waveguides on Z-cut lithium niobate[END_REF] and free from any water loss [START_REF] Kostritskii | Leakage of a guided mode caused by static and light-induced inhomogeneities in channel HTPE-LiNbO 3 waveguides[END_REF][START_REF] Kostritskii | Electro-Optical Properties of Different H x Li 1-x NbO 3 Phases in Proton-Exchanged LiNbO 3 Waveguides[END_REF][START_REF] Korkishko | Characterization of α-phase soft proton-exchanged LiNbO 3 optical waveguides[END_REF], and therefore, we avoid completely these undesirable effects.

EXPERIMENT

We started by fabricating a set of planar waveguides on Z-cut LiNbO 3 . The study of these test samples allowed us to determine the fabrication parameters, to use in order to get a low-loss single mode channel waveguide, operating at a chosen wavelength. We choose to optimize our procedure to operate between 1530 and 1560 nm, as the fabrication of EO modulators, operating in the third telecommunication window, has already been developed. The fabrication parameters were chosen in order to obtain a waveguide composed quasi exclusively of a homogeneous α-phase layer. To determine the phase composition of the exchanged layer, we used the structural phase diagram of H x Li 1-x NbO 3 solid solution [START_REF] Kostritskii | Electro-Optical Properties of Different H x Li 1-x NbO 3 Phases in Proton-Exchanged LiNbO 3 Waveguides[END_REF][START_REF] Korkishko | Characterization of α-phase soft proton-exchanged LiNbO 3 optical waveguides[END_REF].

The APE waveguides were fabricated with standard APE technique (APE1 sample was annealed in dry air at 360 °C) [START_REF] Loni | Very low-loss proton-exchanged waveguides with a substantially restored electro-optic effect[END_REF] and improved APE technique (APE2 sample was annealed in wet air at 360 °C) [START_REF] Muller | Reduction of lattice defects in proton-exchanged lithium niobate waveguides[END_REF][START_REF] Kostritskii | Electro-Optical Properties of Different H x Li 1-x NbO 3 Phases in Proton-Exchanged LiNbO 3 Waveguides[END_REF].

The initial proton exchange of these samples was performed in pure benzoic acid at 175 °C for 6 hours.

Parameters of such an as-exchanged sample (denoted as PE) were measured in order to have reference data. Sequent annealing was made for 7 hours. SPE waveguide was fabricated in sealed ampoule [START_REF] Chanvillard | Soft Proton Exchange on PPLN: a simple waveguide fabrication process for highly efficient nonlinear interactions[END_REF][START_REF] Castaldini | Soft proton exchange tapers for low-insertion loss LiNbO 3 devices[END_REF] and with a benzoic acid melt containing 2.9% of lithium benzoate (LB). Subsurface disorder was studied with micro-Raman and IR-reflection spectroscopy technique, as broadening of Raman phonon bands and decrease of reflectivity are well-established clear indications of a marked structural disorder in LiNbO 3 [START_REF] Kostritskii | Electro-Optical Properties of Different H x Li 1-x NbO 3 Phases in Proton-Exchanged LiNbO 3 Waveguides[END_REF][START_REF] Barker | Optical studies of the vibrational properties of disordered solids[END_REF]. Raman spectroscopy was performed on planar and channel waveguides with the aid of a Jobin-Yvon LabRam HR800 spectrometer operating at excitation wavelength of 632.8 nm. Spectra were measured in the autofocus operation mode of the spectrometer, using a ×100 objective and a confocal hole of width = 50 µm. The IR-reflection spectra with the incidence angle Θ varying from 20 to 80° (Figure 1) were measured with an IR-spectrophotometer "Bruker VERTEX 80v".

Transmission spectra in the visible and near-UV ranges were taken with a Shimadzu UV-3101PC spectrophotometer. The optical absorption spectroscopy data in the visible and the near UV ranges can be used to evaluate the electro-optical r 13 coefficients of the proton-exchanged waveguides in LiNbO 3 crystal. Indeed, the shift of fundamental absorption edge, i.e. the reduction of the band-gap energy ΔE g , has been related to the decrease of spontaneous polarization P and, therefore, of the electro-optic coefficient r 13 [START_REF] Ahlfeldt | Structural and optical properties of annealed proton-exchanged waveguides in z-cut LiTaO 3[END_REF][START_REF] Wemple | Relationship between linear and quadratic electrooptic coefficients in LiNbO 3[END_REF].

For direct measurement of the EO effect in the different H x Li 1-x NbO 3 phases, we fabricated several phase modulators consisting in a channel APE waveguide in a 0.5-mm-thick z-cut LiNbO 3 substrate surrounded by two coplanar electrodes. The EO efficiencies of the modulators were measured using a super luminescent diode (central wavelength ~1550 nm and a bandwidth of 50 nm) and a Sagnac fibre interferometer. Therefore, the phase modulator was pigtailed by PM fibre at both ends input and output.

RESULTS AND DISCUSSION

Our data give the direct evidence of high structural disorder, as the strong broadening of the Raman lines (see Figure 2) shows clearly an increase of the phonon damping (related to crystalline disorder [START_REF] Korkishko | Characterization of α-phase soft proton-exchanged LiNbO 3 optical waveguides[END_REF]) in the PE and APE waveguides with any H x Li 1-x NbO 3 phase, even α-phase with low x. The proton exchange degree, x, has been determined using its well known relations with the exchange-induced refractive index change Δn (e) , determined from dark m-lines spectroscopy data and from Raman data on phonon frequencies [START_REF] Kostritskii | Electro-Optical Properties of Different H x Li 1-x NbO 3 Phases in Proton-Exchanged LiNbO 3 Waveguides[END_REF][START_REF] Korkishko | Characterization of α-phase soft proton-exchanged LiNbO 3 optical waveguides[END_REF][START_REF] Domenico | Oxygen-octahedra ferroelectrics. I. Theory of electro-optical and nonlinear optical effects[END_REF].

The region of Raman spectra ranging from 550 cm -1 to 900 cm -1 have been studied in details (Figure 2 and Table 1), as important information on the structural disorder can be obtained from Raman data of this range [START_REF] Kostritskii | Electro-Optical Properties of Different H x Li 1-x NbO 3 Phases in Proton-Exchanged LiNbO 3 Waveguides[END_REF]. In particular the A 1 (LO)-phonon mode was established to be a good marker of structural disorder [START_REF] Kostritskii | Electro-Optical Properties of Different H x Li 1-x NbO 3 Phases in Proton-Exchanged LiNbO 3 Waveguides[END_REF][START_REF] Paz-Pujalt | Characterization of Proton-Exchange Lithium-Niobate Wave-Guides[END_REF]. It has been demonstrated that the shape of the related Raman band is markedly affected by the presence of structural disorder. According to the significant difference in FHWM of this Raman band between waveguides studied, it's evident that subsurface disorder in SPE waveguides fabricated at appropriate conditions, e.g. SPE sample, is much smaller than in any APE samples. The clear broadening of the A 1 (LO) mode peak in APE waveguides is accompanied by an increase of the signal in the extra band cantered at about 660 -690 cm -1 . Such an extra band has been attributed to structural disorder in Nb-O ferroelectric chains [START_REF] Paz-Pujalt | Characterization of Proton-Exchange Lithium-Niobate Wave-Guides[END_REF]. In this energy region the fitting procedure is more difficult due to the superposition of different modes. Nevertheless, the comparative study of the structural disorder in the different waveguides is possible and, thus, one can observe that for the SPE waveguides this extra band is much smaller, demonstrating a less disordered structure, if compared with APE waveguides. 1). This samples order corresponds to spectra sequence from bottom to top at 871 cm -1 (the wavenumber of A 1 (LO)-phonon band). Besides, Raman intensity of this band was measured focusing the laser beam at different depth in the waveguides. The surface Raman intensity of this extra band is much higher than that obtained when the laser spot is at depth greater than 1.5 µm for any APE waveguide fabricated. For the SPE waveguide it is evident that the damaged surface layer is less pronounced and thinner than that of the "best" APE waveguide, as the extra band intensity at the surface is always smaller and decreases more rapidly with depth, reaching the bulk value at a depth as small as 1.0 µm.

Fig. 2. Raman spectra of the subsurface part of the different samples: PE, APE1, APE2, SPE and z0 (notation meaning is given in Table

This result is confirmed by the measurement of the glancing IR-reflection. In this experiment, an attenuation of the reflectivity R within wavenumbers ν range 580 to 880 cm -1 (Figure 3) is an indication of structural disorder. The disorder factor F = 1 -R s /R 0 , with R s and R 0 being the IR-reflection coefficients measured at Θ = 80° in a sample with a H x Li 1-x NbO 3 waveguide and a pure LiNbO 3 , respectively can be evaluated from these data [START_REF] Kostritskii | Electro-Optical Properties of Different H x Li 1-x NbO 3 Phases in Proton-Exchanged LiNbO 3 Waveguides[END_REF][START_REF] Barker | Optical studies of the vibrational properties of disordered solids[END_REF]. According to the data on vibrational spectra of disordered solids [START_REF] Barker | Optical studies of the vibrational properties of disordered solids[END_REF], such an attenuation of R could be interpreted as phonon scattering effects in the highly structurally disordered parts of the proton exchanged regions. Note, that the photometered depth of near-surface area at reflectivity maximum (i.e., within range of the so-called residual rays [START_REF] Barker | Optical studies of the vibrational properties of disordered solids[END_REF]) is about 1 µm for glancing incidence of IR-radiation (i.e., when reflection angle Θ ≥ 70°), Figure 1. Indeed, the photometered depth becomes much more significant at reflection angles that are far from glancing reflection and closer to the normal incidence [START_REF] Micheli | Crystalline and optical quality of proton exchanged waveguides[END_REF]. Our data for all the studied APE and SPE waveguides show, that a relative attenuation of R (normalized to reflectivity value measured in virgin LiNbO 3 at a given Θ) tends to 0 at such a decrease of Θ. Thus IR-reflection factor F is a well adapted parameter to study the subsurface part of the waveguides.

Fig. 3. IR-reflection spectra of the different samples measured at glancing incidence (incidence angle Θ=80°): PE, APE1, APE2, SPE and z0. This samples order corresponds to spectra sequence from bottom

to top at 871 cm -1 . Incident IR light was s-polarized (i.e., it was polarized perpendicular to plane of incidence (Fig. 1)).

It is well known that proton exchange waveguides present a drastic decrease of the spontaneous polarization P 0 [START_REF] Kostritskii | Electro-Optical Properties of Different H x Li 1-x NbO 3 Phases in Proton-Exchanged LiNbO 3 Waveguides[END_REF][START_REF] Barker | Optical studies of the vibrational properties of disordered solids[END_REF] and polarizability [START_REF] Laurell | Loss of optical nonlinearity in proton--exchanged LiNbO 3 waveguides[END_REF]. This variation of P 0 is related to the shift of the absorption edge [START_REF] Ahlfeldt | Structural and optical properties of annealed proton-exchanged waveguides in z-cut LiTaO 3[END_REF][START_REF] Wemple | Relationship between linear and quadratic electrooptic coefficients in LiNbO 3[END_REF], i.e. the change of the band gap energy ΔE g :

ΔE g ≈ aΔP 0 P 0 ( 1
)
with a = -0.35 eV m 4 /C 2 [START_REF] Kostritskii | Electro-Optical Properties of Different H x Li 1-x NbO 3 Phases in Proton-Exchanged LiNbO 3 Waveguides[END_REF] Since linear EO effect is basically a quadratic effect biased by the P 0 , we can write: r 13,S /r 13,0 = (η 0 /η S ) 3 ×{1 + (ΔE g,s )/(aP 0 ) 2 } (2

)
where η is the packing density, subscripts s and 0 indicate values for H x Li 1-x NbO 3 waveguide and virgin LiNbO 3 crystals, respectively. Note, that the values of ΔE g,s = Е g,s -Е g,0 have a negative sign in all the waveguides studied, as we observe a red-shift of the absorption edge in H x Li 1-x NbO 3 waveguide compared to virgin LiNbO 3 (Figure 4). Band-gap energy E g,s can be estimated from data on the apparent fundamental absorption edge E s (m) in the studied sample [START_REF] Kostritskii | Electro-Optical Properties of Different H x Li 1-x NbO 3 Phases in Proton-Exchanged LiNbO 3 Waveguides[END_REF] with the following equation :

Е g,s = Е s (m) + {(Е g,0 -Е s (m)) 2 + (γ/2) 2 ln(hD 0 )/LD s )} 1/2 (3)
where, Е n (m) is the measured energy of the absorption edge, h is the thickness of the waveguiding layer, γ

is the damping coefficient of the electronic oscillator, D 0 is the optical density in pure LiNbO 3 at Е s (m), D s = D s (m) -D 0 is the optical density of the PE layer, D s (m) is the optical density at Е s (m), L is the sample thickness. The values of γ are different for ordinary and extraordinary polarization: γ (o) = 0.547 eV and γ (e) = 0.702 eV [START_REF] Kostritskii | Electro-Optical Properties of Different H x Li 1-x NbO 3 Phases in Proton-Exchanged LiNbO 3 Waveguides[END_REF].

By using the data on the band-gap shift, the so-called local (microscopic) values of the EO coefficient r 13 were evaluated by Eq.( 1) for all the fabricated samples, Table 1. The results obtained demonstrate a significant dependence of the microscopic EO susceptibility on the phase composition of the APE and SPE waveguides. To determine phase composition of SPE and APE waveguides, we use Raman and IRreflection spectroscopy, as each H x Li 1-x NbO 3 phase has a specific phonon spectrum [START_REF] Kostritskii | Leakage of a guided mode caused by static and light-induced inhomogeneities in channel HTPE-LiNbO 3 waveguides[END_REF][START_REF] Kostritskii | Electro-Optical Properties of Different H x Li 1-x NbO 3 Phases in Proton-Exchanged LiNbO 3 Waveguides[END_REF][START_REF] Kuneva | Spectroscopy of optical waveguiding layers[END_REF]. For direct measurement of EO effect in the different H x Li 1-x NbO 3 waveguides, we fabricated the EO phase modulators based on channel APE waveguides in a 0.5-mm-thick x-cut LiNbO 3 substrate. The electrode structure consists in two electrodes coplanar with the waveguide. In our case of Z-cut, Y-propagating geometry of single-mode APE and SPE LiNbO 3 waveguides, the effective value of electro-optic coefficient r 13 (eff) is inversely proportional to the half-wave voltage V of a phase modulator :

r 13 (eff) = (λG/n 3 e VπOL), ( 4
)
where L is the length of the electrodes in each arm, G is the gap width between the electrodes, O is the overlap factor between the field created by the electrodes and the optical mode, λ is the wavelength of the guided light (central wavelength of SLD used is 1534 nm), n e is the effective refractive index of the TM 0 mode.

The local values of r 13 evaluated from ΔE g,s are larger than the effective values derived from the direct measurement of the EO efficiency of the phase modulators (Table 2) that can be related to the fieldscreening effect in a disordered subsurface part of the waveguide. This degradation of the effective values of the EO coefficient r 13 (eff) correlates with the IR-reflection data on structural disorder at subsurface area of waveguides:

r 13 (eff) ~ r 13 (loc)[1 -CF] (5)

CONCLUSION

In this paper we first describe new simple methods to estimate the disorder and the reduction of the EO coefficients induced in LiNbO 3 by the different kind of known proton exchange processes. The experimental characterizations of a set of planar waveguides fabricated by different PE techniques show that the SPE process allows preparing devices with the smaller degradation of the crystalline quality of the subsurface layer and the smaller reduction of the effective electro-optic coefficient. This is confirmed by the characterization of a set of phase modulators prepared by the standard and the improved APE techniques as well as by SPE. SPE devices exhibit the higher electro-optical efficiency and the lower propagation losses. We can also expect that the SPE modulators will present a low drift of the DC bias, as subsurface disorder was previously identified as a key source of such a drift.

Fig. 1 .

 1 Fig. 1. Schematic draw of subsurface layer study with specular IR reflection. d is thickness of the subsurface disordered part of a proton-exchanged LiNbO3 waveguiding layer. Photometered depth, dp, at IR reflection measurements is limited by penetration depth of IR radiation in near-surface layer at glancing incidence. Hence dp ≤ 1 µm that is about d.

Fig. 4 .

 4 Fig. 4. Optical absorption spectra of the different samples: z0, SPE, APE1 and PE. This samples order corresponds to spectra sequence from left to right side at D = 2. Absorption spectrum of z-cut wafer of a pure congruent LiNbO 3 crystal (sample z0) is given as reference data for evaluation of band-gap energy changes.

  Broadening of Raman band correlates with the variation of the disorder factor F, and, therefore, Raman data may be used for the evaluation of the disorder-induced reduction of the effective electro-optic coefficient. To reach easy possibility of such evaluation, we obtain experimental calibration curve for F vs Δν 1/2 (it's a FWHM of A 1 (LO)-phonon band), Figure5. It should be noted, that micro-Raman spectroscopy can be used for study of channel waveguides in contrast to IR-reflection spectroscopy that requires the preparation of planar waveguides.Considering that electro-optic susceptibility is directly related to Raman scattering cross-section[START_REF] Kostritskii | Electro-Optical Properties of Different H x Li 1-x NbO 3 Phases in Proton-Exchanged LiNbO 3 Waveguides[END_REF][START_REF] Kaminow | Quantitative Determination of Sources of electro-optic effect in LiNbO 3 and LiTaO 3[END_REF], micro-Raman measurements were performed to detect any change in the Raman active phonon modes. It is also known that the contribution of the Nb-O bond to the linear and nonlinear susceptibility is larger than the contributions of the Li-O bond. Mode at 870 cm -1 , associated with the Nb-O 6 octahedron belongs to the dominant contributors to the EO effect[START_REF] Kuneva | Spectroscopy of optical waveguiding layers[END_REF][START_REF] Paz-Pujalt | Characterization of Proton-Exchange Lithium-Niobate Wave-Guides[END_REF]. Therefore, a simple correlation between the related Raman intensities and the electro-optic coefficients can be found. In fact, recent findings show that attenuation of Raman intensity (I LO ) of the LO-phonon band at 871 cm -1 are related to the decrease of the local electro-optic response[START_REF] Kostritskii | Electro-Optical Properties of Different H x Li 1-x NbO 3 Phases in Proton-Exchanged LiNbO 3 Waveguides[END_REF][START_REF] Kaminow | Quantitative Determination of Sources of electro-optic effect in LiNbO 3 and LiTaO 3[END_REF]. Moreover, measurement of Raman intensities of certain phonon bands allows for rough estimation of local values of electro-optic coefficients in waveguides in different areas of any proton-exchanged LiNbO 3 channel waveguide.

Fig. 5 .

 5 Fig. 5. Correlation between disorder factor F determined from IR-reflection data and broadening of the Raman band at 871 cm -1 . Δν 1/2 is the change of FWHM of this Raman band measured in subsurface area with confocal micro-Raman spectroscopy.

Table 1

 1

	sample		IR	Raman	Raman	Raman	UV-	r S (loc)/r 0
			F	Δν 1/2 , cm -1	I LO /I TO	I extra /I TO	absorption ΔЕ g,S (o) , eV
	z0,	pure	0	29.7	3.5	0	0	1
	LN
	SPE		0.05	35.3	3.2	0.14	-0.010	0.93
	APE2		0.082	40.5	3.1	0.18	-0.013	0.90
	APE1		0.095	44.1	3.0	0.24	-0.016	0.88
	PE	(no	0.298	37.8	1.5	2.6	-0.119	0.19
	annealing)

Table 1 .

 1 Sample notation, F is the disorder factor, Δν1/2 is the FWHM of the Raman A1(LO)-phonon band at 871 cm-1, ILO and Iextra are the Raman intensities of the A1(LO)-phonon band and extra band, ΔЕg,S(o) is the band-gap shift for the ordinary polarization and rS(loc)/r0 is the normalized values of local electro-optic coefficient r 13 .

Table 2

 2

	sample	APE1 [1]	APE1 our	APE2	SPE	PE**
	r 13 '(loc)	~0.85*	0.88	0.90	0.93	0.19
	r 13 '(eff)	0.69	0.80	0.85	0.90	≤ 0.1
	F	~0.12*	0.095	0.082	0.05	0.22
	*-values evaluated by us for test samples fabricated by technique reported in Ref.[1]
	**-values for a sample proton-exchanged in pure benzoic acid at 175 °C for 6 hours.

Table 2

 2

	. The normalized values of local and effective electro-optic coefficient r13ʹ′ = r13,S/r13,0 , where
	r13,S and r13,0 are values for a sample S studied and pure LiNbO3. Ratio between r13'(eff) and
	r13'(loc) correlates with actual value of disorder factor F.